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:
- SRTM1 ( 1-arc second sampling)
- SRTM3 ( 3-arc second sampling)
1 second in angular units is
3-arc second is equivalent to
If we consider N26E092.hgt
file to be a SRTM3 data, then it contains
data for latitudes .
If above file was in 1-arc second resolution, then it would contain data for latitudes .
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:
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:
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
Links
- Official SRTM Download Link
- Coverage Map of latest SRTM data - Jan 2015
- Linked coverage map of elevation data - Easy way to download collection of hgt files - Void Filled
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.