Extracting elevation from .HGT file?

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.


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.


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

Tags:

Elevation

Srtm

Recent Posts