Intense queries on Elevation
Bevan Rudge
I'm a software engineer at CivicActions and am doing some research on
behalf of Conservation Strategy Fund for a tool which will be
available to general public to calculate data about an arbitrary
number of proposed hydro-electric power stations.
The idea is that the tool will accept a line representing the dam wall
(which is probably imprecise), and an elevation of the surface of the
reservoir. The tool will then accurately find the line representing
the dam wall (making the assumption that there is only one dam wall
and that it is approximately straight) and use SRTM3 and/or GEOTOPO30
data to "follow a contour line" and find the approximate polygon that
most accurately represents the surface of the reservoir. From this,
the tool can make estimates about the volume of the proposed
reservoir, the bio-mass of the vegetation that will be inundated in
the reservoir and the environmental costs.
Geoname's SRTM web service
(http://ws.geonames.org/srtm3?lat=-1.086809&lng=-80.134964) provides
the data in an appropriate and accessible format, however the
algorithm will need to make thousands of sequential requests and would
take days to process one reservoir if this data is only available as a
web service. We would need to have a service similar to
ws.geonames.org/srtm3 available on the same machine as the algorithm
that finds the polygon.
How might I go about re-creating this service on our own machines?
Are there specialists who would be able to do this for us in a more
efficient manner? Who should I contact to get an estimate?
Currently we don't have funding but are in the process of applying for
a grant from Google and are confident that we have a good chance at
winning it. This aspect of the system is currently the highest-risk
part and makes it difficult to budget the project and know what size
grant we should request.
I have reviewed the SRTM data at
ftp://e0srp01u.ecs.nasa.gov/srtm/version2/, and my first steps, if we
can not find someone more experienced in this domain to do it for us,
would be to write code that parses this data into a Postres or MySQL
database with 3 indexed columns; latitude, longitude, elevation. Then
I could then query the data for pixels in a containing box at the
specified elevation +/- 100m, or find the line representing the edge
of the lake by interpolating the elevations of neighbouring pixels and
finding the direction which gives the targeted elevation.
I'm sure that this, parts of this or similar problems have been solved
before and am looking to avoid re-inventing the wheel. Please let me
know of any services, individuals or resources that may be of value to
us.
Thanks,
Bevan/
--
Changing the world, one node at a time | CivicActions | http://CivicActions.com
Drupal, Usability, Open Source, Tech | http://Drupal.geek.nz |
http://Twitter.com/BevanR
Gtalk be...@civicactions.net | YIM rudgy_m_nz | AIM b.rudge | skype b.rudge
Bevan Rudge
some ballpark estimates. Anyone interested and available? Anyone who
is able to estimate it now will probably be called on for the work
once it gets approved.
Thanks,
Bevan/
Bevan Rudge
some ballpark estimates. Anyone interested and available? Anyone who
is able to estimate it now will probably be called on for the work
once it gets approved.
Thanks,
Bevan/
On Tue, Dec 9, 2008 at 3:14 PM, Bevan Rudge <b.r...@gmail.com> wrote:
Anthony Cartmell
>> How might I go about re-creating this service on our own machines?
>> Are there specialists who would be able to do this for us in a more
>> efficient manner? Who should I contact to get an estimate?
I've re-created the SRTM elevation look-up on my machine, using the raw
SRTM binary files and PHP. The files take up quite a bit of disk space,
but the time taken to find the elevation of a point is pretty small. Since
the data is well-organised in the files, it's very quick to work out which
file to read, and which part of it to read, then to get the elevation
value. I don't think putting it in a DB will help, unless you want to run
queries more complicated than "how high is this point?".
I don't have a contour-generating algorithm written, but might have some
useful pointers.
HTH,
Anthony
--
www.fonant.com - Quality web sites
Bevan Rudge
> I've re-created the SRTM elevation look-up on my machine, using the raw
> SRTM binary files and PHP. The files take up quite a bit of disk space,
> but the time taken to find the elevation of a point is pretty small. Since
> the data is well-organised in the files, it's very quick to work out which
> file to read, and which part of it to read, then to get the elevation
> value.
Would you be willing to share this code on a GPL license or otherwise?
The fact it's PHP is a bonus for us.
> I don't think putting it in a DB will help, unless you want to run
> queries more complicated than "how high is this point?".
It's likely that it will be convenient to query the data in other
ways, such as; "what points in a containing box are at elevation X +/-
Y?" And "what are the coordinates and elevations of the Z points
closest to point W?" where W is an arbitrary point not where an SRTM
data point is.
Thanks!
Anthony Cartmell
> Would you be willing to share this code on a GPL license or otherwise?
> The fact it's PHP is a bonus for us.
Here's my PHP class to extract elevation data from the (uncompressed) SRTM
files:
<?php
class SrtmReader
{
public static function heightForPoint(GeoPoint $point)
{
$latinput = $point->latitude;
$lnginput = $point->longitude;
$lat_degree = floor($latinput);
$lng_degree = floor($lnginput);
$lat_letter = ($latinput>=0)?'N':'S';
$lng_letter = ($lnginput>=0)?'E':'W';
$lng_tilenum = abs($lng_degree);
$lat_tilenum = abs($lat_degree);
$filename =
$lat_letter.sprintf('%02d',$lat_tilenum).$lng_letter.sprintf('%03d',$lng_tilenum).".hgt";
$dirname = substr($filename,0,2);
$filename = '/home/sites/dem/'.$dirname.'/'.$filename;
//echo $filename;
if (!file_exists($filename)) return;
$col = round(($lnginput-$lng_degree)*1200);
$row = 1200 - round(($latinput-$lat_degree)*1200);
$fp = fopen($filename,'r');
fseek($fp,(2402*$row)+($col*2));
$h = fread($fp,2);
$heights = unpack('n',$h); // 'n' returns unsigned short (always
16 bit, big endian byte order)
$height = $heights[1]; // result of unpack() indexes from 1,
not zero.
if ($height>32767) $height -= 65536; // make signed
if ($height==-32768) return;
return $height;
}
}
?>
The function returns NULL if the SRTM file for the location doesn't
exists, or if the height for the point is -32768m (e.g. no data).
Please feel free to make use of it, tidy it, find bugs and fix them, etc :)
...and let me know of any useful changes! And perhaps mention
ajcar...@fonant.com in the credits somewhere?
>> I don't think putting it in a DB will help, unless you want to run
>> queries more complicated than "how high is this point?".
>
> It's likely that it will be convenient to query the data in other
> ways, such as; "what points in a containing box are at elevation X +/-
> Y?" And "what are the coordinates and elevations of the Z points
> closest to point W?" where W is an arbitrary point not where an SRTM
> data point is.
In which case a DB may well be useful. There are an awfully large number
of data points though!
Cheers!
Bevan Rudge
Thanks very much! I'll let you know if/when we use it to see how we
can collaborate.
I did some quick math on the numbers for the database; there are, of
the order of, 10 billion data points. I'll also need to do some
research and/or testing IRT the efficiency of queries on such a large
data set. Luckily every field can be an int (if coordinates are
stored as arc-seconds instead of degrees + decimal, or degrees +
minutes + seconds)!
I'm curious about how many of the points have the -32768 flag – since
these would be excluded from the database.
Cheers,
Bevan/
Anthony Cartmell
> these would be excluded from the database.
Not sure, but they tend to appear in mountainous regions, I think where a
peak has hidden a valley point from the Space Shuttle's radar.
XCRunner
on earth were scanned 2-3 times to produce this dataset. The
limitation however was in geographic areas of high slope there could
be shadows from the synthetic aperture radar. For example, as the
radar scanned back and forth, at areas of great eccentricity, a valley
could be hidden behind a high mountain peak (line of sight to the
shuttle). In most places the repeat scans from different angles
provided the necessary information to correct for this. Synthetic
aperture radar also requires many Fourier transforms which do not
always have a solution, which would also yield no elevation
information. In all these cases the point in question would be marked
null.
Rather than exclude these points, it would be better to use several
utilities available to replace these points with the mean of the
points around them. This is an approximation, but sufficient for your
purposes. Or, you could modify the above class file to seek right/left
or up/down until a non-null point is found and return that value.
Scanning further to the east (right) would be easiest given the class
as is.
For your application you could get the data from many locations such
as the USGS elevation web API, which has a greater granularity and
accuracy than even the SRTM-1 v2 data (or Geonames). The limitation
here is they can only process a single point per connection instance
taking ~600ms per point.
You can insert all datapoints into a mysql db and pull each out for
about 1-2ms overhead time. Using this manner you can rather easily and
quickly produce the results you are seeking in seconds or minutes
rather than days. The limitation here is that you will take up much
more space with the indicies than the data being stored in the DB. The
DB can be compressed, but since it is static data, a db is overkill.
Or, you can decompress all the files (.hgt, a.k.a. "height" files) to
the ~28GB mass and run them using the PHP class, but this will allow
only singular access to the points again and you must load the 25MB
file for each point, again reducing speed. If you are pulling a ton of
points, this will be a limitation. However, using the above class with
fclose() added, it still takes only 0.8ms to pull a point. For your
application then, this is likely optimal.
Joseph
On Dec 12 2008, 4:27 am, "Anthony Cartmell" <ajcartm...@fonant.com>
wrote:
Kirk
currently store our DEM in a ArcInfo grid. I have a new project where
I will to move this off site open source solution, where I would
require batch processing as well as transactional. I am quite new to
this.
On Dec 10 2008, 4:19 pm, "Bevan Rudge" <be...@civicactions.com> wrote:
> Anthony,
>
> Thanks very much! I'll let you know if/when we use it to see how we
> cancollaborate.
>
> I did some quick math on the numbers for the database; there are, of
> the order of, 10 billion data points. I'll also need to do some
> research and/or testing IRT the efficiency of queries on such a large
> data set. Luckily every field can be an int (if coordinates are
>
> I'm curious about how many of the points have the -32768 flag – since
> these would be excluded from the database.
>
> Cheers,
> Bevan/
>
>
>
> On Thu, Dec 11, 2008 at 4:20 AM,Anthony Cartmell<ajcartm...@fonant.com> wrote:
>
> > Bevan,
>
> >> Would you be willing to share this code on aGPL licenseor otherwise?
>
> > Here's myPHP classto extract elevation data from the (uncompressed) SRTM
> > files:
>
> > <?php
>
> > class SrtmReader
> > {
> > publicstatic functionheightForPoint(GeoPoint $point)
> > $latinput = $point->latitude;
> > $lnginput = $point->longitude;
>
> > $lat_degree = floor($latinput);
> > $lng_degree = floor($lnginput);
>
> > $lat_letter = ($latinput>=0)?'N':'S';
> > $lng_letter = ($lnginput>=0)?'E':'W';
>
> > $lng_tilenum = abs($lng_degree);
> > $lat_tilenum = abs($lat_degree);
>
> > $filename =
> > $lat_letter.sprintf('%02d',$lat_tilenum).$lng_letter.sprintf('%03d',$lng_tilenum).".hgt";
> > $dirname = substr($filename,0,2);
> > $filename = '/home/sites/dem/'.$dirname.'/'.$filename;
>
> > //echo $filename;
>
> > if (!file_exists($filename)) return;
>
>
> > $fp = fopen($filename,'r');
>
> > fseek($fp,(2402*$row)+($col*2));
>
> > $h = fread($fp,2);
> > $heights = unpack('n',$h); // 'n' returnsunsigned short(always
> > not zero.
> > if ($height>32767) $height -= 65536; // make signed
> > if ($height==-32768) return;
> > return $height;
> > }
> > }
>
> > ?>
>
> > The function returns NULL if the SRTM file for the location doesn't
> > exists, or if the height for the point is -32768m (e.g. no data).
>
> > Please feel free to make use of it, tidy it, find bugs and fix them, etc :)
>
> > ...and let me know of any useful changes! And perhaps mention
> >>> I don't think putting it in a DB will help, unless you want to run
> >>> queries more complicated than "how high is this point?".
>
> >> Y?" And "what are the coordinates and elevations of the Z points
> >> closest to point W?" where W is an arbitrary point not where an SRTM
> >> data point is.
>
> > In which case a DB may well be useful. There are an awfully large number
> > of data points though!
>
> > Cheers!
>
> > Anthony
> > --
> Drupal, Usability,Open Source, Tech |http://Drupal.geek.nz|http://Twitter.com/BevanR
Steve
On Dec 10 2008, 4:20 pm, "Anthony Cartmell" <ajcartm...@fonant.com>
wrote:
>
> > Would you be willing to share this code on a GPL license or otherwise?
> > The fact it's PHP is a bonus for us.
>
> Here's my PHP class to extract elevation data from the (uncompressed) SRTM
> files:
>
Javascript but without success.
> <?php
>
> class SrtmReader
> {
> public static function heightForPoint(GeoPoint $point)
Regards,
Steve
Steve
On Jan 30, 3:19 pm, Steve <stephen.jo...@googlemail.com> wrote:
> On Dec 10 2008, 4:20 pm, "Anthony Cartmell" <ajcartm...@fonant.com>
> wrote:> Bevan,
>
> Could you give an example of how a call to this function would look?
>
What I was trying to ask is what is the format of the GeoPoint $point?
From looking at the code I assume that it is simply an object with two
properties, latitude and longitude, such as the Google Maps API
GLatLng object. Is this assumption correct?
Regards,
Steve
Anthony Cartmell
> What I was trying to ask is what is the format of the GeoPoint $point?
>> From looking at the code I assume that it is simply an object with two
> properties, latitude and longitude, such as the Google Maps API
> GLatLng object. Is this assumption correct?
Yes :)
Sorry, have been busy over the last few days, but I was planning to answer
your question at some point...
I should probably re-write the function to just have $latitude and
$longitude arguments, which would be usable without having to construct
and object first. This is just how it is in my code, where I like to use
objects to get some encapsulation and argument type-checking in PHP.
Steve
SRTM Data at work.
Regards,
Steve.
Kirk Webb
Anthony Cartmell
There is no interpolation :)
Kirk Webb
I have switched gears and am now trying to use the jasper JP2 library with GDAL to retrieve the elevation for a point. I will still need some kind of interpolation. As well as figuring out the optimal image size for the fastest response.