GeoTIFF-TIFF – How to Read Elevation from 1/3 Arc-Second NED .tif or .img Files

elevationgeotiff-tiff

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!

Best Answer

The solution for me was to convert the files to BIL files using GDAL:

gdal_translate -a_srs NAD83 -of EHdr file.img file.bil

This created .bil and .hdr files. The .hdr file contains the details you need to determine where to read the file. I was then able to get the elevation using the PHP unpack function. Performance is lightning fast.

Related Question