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!

user1517922
  • 141
  • 4
  • 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 01:06

1 Answers1

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.

user1517922
  • 141
  • 4
  • 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 04 '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 04 '15 at 19:59