Extracting elevation from .HGT file?

Question

I want to assign a specific long/lat position on a map to the elevation from SRTM3 data files, but have no idea how to find the specific value. So I want some example of how I can find in N50E14.hgt elevation to 50°24'58.888"N, 14°55'11.377"E.

Answer 1

Data format

I'll take it as a little exercise in how to program a data reader. Have a look at the documentation:

SRTM data are distributed in two levels: SRTM1 (for the U.S. and its territories and possessions) with data sampled at one arc-second intervals in latitude and longitude, and SRTM3 (for the world) sampled at three arc-seconds.

Data are divided into one by one degree latitude and longitude tiles in "geographic" projection, which is to say a raster presentation with equal intervals of latitude and longitude in no projection at all but easy to manipulate and mosaic.

File names refer to the latitude and longitude of the lower left corner of the tile - e.g. N37W105 has its lower left corner at 37 degrees north latitude and 105 degrees west longitude. To be more exact, these coordinates refer to the geometric center of the lower left pixel, which in the case of SRTM3 data will be about 90 meters in extent.

Height files have the extension .HGT and are signed two byte integers. The bytes are in Motorola "big-endian" order with the most significant byte first, directly readable by systems such as Sun SPARC, Silicon Graphics and Macintosh computers using Power PC processors. DEC Alpha, most PCs and Macintosh computers built after 2006 use Intel ("little-endian") order so some byte-swapping may be necessary. Heights are in meters referenced to the WGS84/EGM96 geoid. Data voids are assigned the value -32768.

How to proceed

For your position, 50°24'58.888"N 14°55'11.377"E, you already found the correct tile, N50E14.hgt. Let's find out which pixel you are interested in. First latitude, 50°24'58.888"N:

24'58.888" = (24 * 60)" + 58.888" = 1498.888"

arc seconds. Divided by three and rounded to the closest integer gives a grid row of 500. The same calculation for longitude results in grid column 1104.

The quickstart documentation lacks information about how rows and columns are organized in the file, but in the full documentation it is stated that

The data are stored in row major order (all the data for row 1, followed by all the data for row 2, etc.)

The first row in the file is very likely the northernmost one, i.e. if we are interested in row 500 from the lower edge, we actually have to look at row

1201 - 500 = 701

from the start if the file. Our grid cell is number

(1201 * 700) + 1104 = 841804

from the start of the file (i.e. skip 700 rows, and in the 701st one take sample 1104). Two bytes per sample means we have to skip the first 1683606 bytes in the file and then read two bytes in order to get our grid cell. The data is big-endian, which means that you need to swap the two bytes on e.g. Intel platforms.

Sample program

A simplistic Python program to retrieve the right data would look like this (see the docs for use of the struct module):

import struct

def get_sample(filename, n, e):
    i = 1201 - int(round(n / 3, 0))
    j = int(round(e / 3, 0))
    with open(filename, "rb") as f:
        f.seek(((i - 1) * 1201 + (j - 1)) * 2)  # go to the right spot,
        buf = f.read(2)  # read two bytes and convert them:
        val = struct.unpack('>h', buf)  # ">h" is a signed two byte integer
        if not val == -32768:  # the not-a-valid-sample value
            return val
        else:
            return None

if __name__ == "__main__":
    n = 24 * 60 + 58.888
    e = 55 * 60 + 11.377
    tile = "N50E14.hgt"  # Or some magic to figure it out from position
    print get_sample(tile, n, e)

Note that efficient data retrieval would have to look a little more sophisticated (e.g. not opening the file for each and every sample).

Alternatives

You could also use a program which can read the .hgt files out of the box. But that is boring.

Answer 2

GDAL can read/write these raster formats with the SRTMHGT driver. This means you can view the raster with QGIS, ArcGIS, or use GDAL utilities like gdallocationinfo to get values from a point, e.g.:

Convert DMS to DD:

• Lat: 50°24'58.888"N = 50 + (24 / 60) + (58.888 / 3600) = 50.4163577778
• Long: 14°55'11.377"E = 14 + (55 / 60) + (11.377 / 3600) = 14.9198269444

Then from a shell, use gdallocationinfo file.hgt -wgs84 long lat:

$ gdallocationinfo N50E14.hgt -wgs84 14.9198269444 50.4163577778
Report:
  Location: (1104P,700L)
  Band 1:
    Value: 216

The elevation is 216 m.

Answer 3

If you do use QGIS, check if the python plugin "Point Sampling Tool" is installed. You'll find it at -> Enhancements (Python) -> Analyse.

Select your point layer of the required positions, then start the PST, choose the hgt (or whatever raster/polygone file) and choose a new point shape for output.

Thats all :-)

  Chris

Answer 4

Chris' answer indicates it is straightforward to sample points from a layer in QGIS.

However, since your reply to my comment clarifies you are writing your own program to read elevation values from the .hgt files, have another look at the Quickstart PDF in the SRTM docs. It explains how the elevation data is stored. To summarize:

• SRTM3 files contain a sequence of big-endian integer values.
• Each integer value represents an elevation "in meters referenced to the WGS84/EGM96 geoid", except for values of -32768, which indicate no-data pixels.
• There are 1201 lines of 1201 samples, so there should be 1442401 integer values altogether.

You say you can convert between lon/lat coordinates and pixels, so getting the elevation is a matter of reading the integer value from the appropriate offset in the file. Given pixel coordinates x and y relative to the upper left corner of the scene, that's basically offset = (y * 1201) + x. Pixel 0,0 is the first integer in the file and pixel 1200,1200 is the last integer in the file.