Reading SRTM data with Python


Reading SRTM data with Python

2015/10/29

Tags: SRTM GIS HGT

Digital Elevation Model (DEM) provides elevation data for use in programs and analysis. These models have varied practical uses like line-of-sight and viewshed analysis. Shuttle Radar Topography Mission (SRTM), as the name suggests, is a research mission to produce a free DEM for general use. If you are like me, you hadn’t heard about DEM before you actually had to use it. And, you will be amazed where they are being used.

Data Description

File Name

SRTM data file has a .hgt extension (a shortform for height ?) and its name defines its coverage. For example, a filename named N26E092.hgt contains elevation data that stretches from 26°N 92°E to 27°N 93°E

Data Resolution

Two kinds of SRTM data are available:

1 second in angular units is 1 60 minutes = 1 3600 degrees

3-arc second is equivalent to 3 3600 degrees = 1 1200 degrees

If we consider N26E092.hgt file to be a SRTM3 data, then it contains data for latitudes 26, 26 + 11200, 26 + 21200, .... , 26 +  12001200 = 27.

If above file was in 1-arc second resolution, then it would contain data for latitudes 26, 26 + 13600, 26 + 23600, .... , 26 +  36003600 = 27.

Same is the case for Longitude.

Hence, SRTM3 data contains 1201 * 1201 individual elevations and SRTM1 data contains 3601 * 3601 elevations.

I assume you are now aware of the differences between these two resolutions. Now, let’s move on to the physical representation of data in the file.

Data Structure

From now on, I am going to assume a 3-arc second data(SRTM3) named N27E086.hgt. Following figure depicts its general structure:

Physical representation of data in hgt file

But, how do hgt files store 2d array inside a binary file? If you know a bit about array indexing, array data is stored sequentially in memory and a special function maps it to an individual cell. We also need a mapping function to convert a given latitude-longitude combo to data position in file.

SRTM defines data as 16bit signed integer in Big-Endian. So, if you can read binary files, you could read it easily after you account for these details. But, numpy makes our job a lot more easier. Provide numpy the specification of data and ask it to convert that data into a 2D array.

Enough talk, following Python program does this:

import numpy as np
SAMPLES = 1201 # Change this to 3601 for SRTM1
def read_elevation_from_file(hgt_file, lon, lat):
    with open(hgt_file, 'rb') as hgt_data:
        # Each data is 16bit signed integer(i2) - big endian(>)
        elevations = np.fromfile(hgt_data, np.dtype('>i2'), SAMPLES*SAMPLES)\
                                .reshape((SAMPLES, SAMPLES))

        lat_row = int(round((lat - int(lat)) * (SAMPLES - 1), 0))
        lon_row = int(round((lon - int(lon)) * (SAMPLES - 1), 0))

        return elevations[SAMPLES - 1 - lat_row, lon_row].astype(int)

Numpy converts the data to 2d array. To find the index given latitude and longitude, we need to solve for x to get the index:

27+ x 1200 = latitude

86+ x 1200 = longitude

Say, Latitude = 27.006666666666666667. It implies that x = 8. Solving the equation is actually easy as we know for sure that the integer part of latitude is always going to be 27 for a given file. (lat - int(lat)) * 1200 represents this logic. Notice that, latitude starts from lower left and moves upwards. So, we need to account for this. Hence, Latitude index = 1200 - 8. That’s what 1200-lat_row does. Same is the case for longitude - only index correction is not required. Rounding is done as values outside exact multiples of 1/1200 may be present in fractional part of given latitude-longitude combo.

Reading from multiple files

We got the general idea of how to read data from a single data file. Following snippet combines this logic to print out data for any latitude-longitude given the file that contains that latitude-longitude is available for reading:

import os
import json
import numpy as np

SAMPLES = 1201  # Change this to 3601 for SRTM1
HGTDIR = 'hgt'  # All 'hgt' files will be kept here uncompressed


def get_elevation(lon, lat):
    hgt_file = get_file_name(lon, lat)
    if hgt_file:
        return read_elevation_from_file(hgt_file, lon, lat)
    # Treat it as data void as in SRTM documentation
    # if file is absent
    return -32768


def read_elevation_from_file(hgt_file, lon, lat):
    with open(hgt_file, 'rb') as hgt_data:
        # HGT is 16bit signed integer(i2) - big endian(>)
        elevations = np.fromfile(hgt_data, np.dtype('>i2'), SAMPLES*SAMPLES)\
                                .reshape((SAMPLES, SAMPLES))

        lat_row = int(round((lat - int(lat)) * (SAMPLES - 1), 0))
        lon_row = int(round((lon - int(lon)) * (SAMPLES - 1), 0))

        return elevations[SAMPLES - 1 - lat_row, lon_row].astype(int)

def get_file_name(lon, lat):
    """
    Returns filename such as N27E086.hgt, concatenated
    with HGTDIR where these 'hgt' files are kept
    """

    if lat >= 0:
        ns = 'N'
    elif lat < 0:
        ns = 'S'

    if lon >= 0:
        ew = 'E'
    elif lon < 0:
        ew = 'W'

    hgt_file = "%(ns)s%(lat)02d%(ew)s%(lon)03d.hgt" % {'lat': abs(lat), 'lon': abs(lon), 'ns': ns, 'ew': ew}
    hgt_file_path = os.path.join(HGTDIR, hgt_file)
    if os.path.isfile(hgt_file_path):
        return hgt_file_path
    else:
        return None

# Mt. Everest
print get_elevation(86.925278, 27.988056)
# Kanchanjunga
print get_elevation(88.146667, 27.7025)

You can clone the latest source from here.

For non-programmers

If you just like to get the elevation without programming or using this code, use gdallocationinfo. In Linux, you can install it using:

$ sudo apt-get install gdal-bin
$ gdallocationinfo N27E088.hgt -wgs84 88.2123 27.1115

Some thoughts

Things like this is what I love about programming. You spend more time learning than you do programming. Also, I was trying to learn Inkscape. So, I drew the above illustration in Inkscape.

Please feel free to ask or comment.