I'm trying to read multiple elevations from a USGS 1/3 arc-second IMG / GeoTIFF file.
A sample of the files is available here (5 MB): http://tdds3.cr.usgs.gov/Ortho9/ned/ned_13/img/n43w108.zip
I can query the .img file and get the elevation using GDAL: gdallocationinfo -wgs84 imgn60w145_13.img -144.267361111152354 59.9981944444444082
This returns: Report: Location: (7918P,25L) Band 1: Value: 1.2593731880188
Because a GeoTIFF can be read quite a bit faster than an IMG file (in my benchmarks) I convert the IMG to a WGS84 GeoTIFF: gdal_translate -a_srs WGS84 -of GTiff imgn60w145_13.img imgn60w145_13.tif
gdallocationinfo only supports querying one point at a time. I need to get the elevation every 10 meters between 2 points (along a straight path), so I'm trying to read the file directly, with the following PHP code (but the approach is the same in C or Python):
$STRIPOFFSETS = 230; // 0xe6
$LEN_OFFSET = 4;
$DATA_VOID = 0x8000; // data void ( = signed int -32768)
$LEN_DATA = 4; // the number of bytes containing each item of elevation data // ( = BitsPerSample tag value / 8)
$fp = fopen("imgn60w145_13.tif", 'rb');
if ($fp === false) {
echo "Could not open the file\n";
}
// first data offset
fseek($fp, $STRIPOFFSETS);
// find the location of the required data row in the StripOffsets data
$dataOffset = $STRIPOFFSETS + ($row * $LEN_OFFSET);
fseek($fp, $dataOffset);
$dataBytes = fread($fp, $LEN_OFFSET);
$data = unpack('VdataOffset', $dataBytes);
echo print_r($data, true);
// this is the offset of the 1st column in the required data row
$firstColOffset = $data['dataOffset'];
// now work out the required column offset relative to the 1st column
$requiredColOffset = $col * $LEN_DATA;
// combine the two and read the elevation data at that address
fseek($fp, $firstColOffset + $requiredColOffset);
$dataBytes = fread($fp, $LEN_DATA);
echo "1: " . $firstColOffset . " + " . $requiredColOffset . "\n";
echo "2: " . $LEN_DATA . "\n";
echo "3: " . $dataBytes . "\n";
$data = unpack('velevation', $dataBytes);
$elevation = $data['elevation'];
if ($elevation == $DATA_VOID) {
$elevation = 0;
}
echo $elevation . "\n";
I don't know that I have the right values for $STRIPOFFSETS, $LEN_OFFSET, or $LEN_DATA. The code above is a limited snippet from SRTMGeoTIFFReader (http://www.osola.org.uk/elevations/index.htm), which I'm not able to get to work either.
I'd really prefer to be able to do this using a direct approach similar to the above (no external libraries, etc), but I can live with a Python / GDAL approach as well.
Thanks!
