1

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!

  • GDAL would be good for this, you can open the raster and read the block (extent) and then index the rows from the returned array. I don't know PHP so can't help any further.. sorry. Another option is BIL format which can be read directly as a binary file. The data is stored as just data and can be indexed as such if it's in the right coordinate system. – Michael Stimson May 12 '14 at 21:57
  • Thanks, BIL was just the nudge I needed. Much easier than GeoTiff (no offsets to deal with), and lightning fast. – user1517922 May 13 '14 at 14:19
  • You're welcome. I've used BIL format for just that reason myself, also the other way around - it's easy to write to the binary at file offsets and the header is just text... that was before GDAL though. – Michael Stimson May 14 '14 at 1:06
2

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.

  • Are you able to give a code sample of your solution using unpack()? I am trying to get the elevation for a lat/lon from a NED img file. – fooquency Sep 4 '15 at 15:40
  • There is a lot more to this (determining which row, which column, getting NCOLS from the .hdr, but this should work: $bil = fopen('file.bil', 'rb'); $NCOLS = 3601; $col = 1; $row = 1; fseek($bil, ($col * 4) + ($row * $NCOLS * 4)); $dataBytes = fread(bil, 4); $data = unpack('f', $dataBytes); echo $data['1'] . " meters\n"; – user1517922 Sep 4 '15 at 19:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.