Beginning Digital Elevation Model work with Python
Digital Elevation Models (DEMs) are height maps of the Earth, taken by satellites. One of the most accessible sources is the Shuttle Radar Topography Mission (SRTM).
Here’s an example,
The raw data files are very basic binary files; I’ve used the 3 arc-seconds data, as .hgt (height) files. You can find these files here.
For the 3 arc-second dataset, an hgt file consists of a grid of 1201×1201 cells representing the heights across a 1 degree by 1 degree cell. They’re held in big-endian format, 2 bytes per cell.
Luckily, Python has a very easy way to parse binary files; the struct module.
Here’s a simple parser class…
import struct
class srtmParser(object):
def parseFile(self,filename):
# read 1,442,401 (1201x1201) high-endian
# signed 16-bit words into self.z
fi=open(filename,"rb")
contents=fi.read()
fi.close()
self.z=struct.unpack(">1442401H", contents)
def writeCSV(self,filename):
if self.z :
fo=open(filename,"w")
for row in range(0,1201):
offset=row*1201
thisrow=self.z[offset:offset+1201]
rowdump = ",".join([str(z) for z in thisrow])
fo.write("%s\n" % rowdump)
fo.close()
else:
return None
That will dump a CSV file ready for use in OpenOffice or Excel; in the case of Excel, a Surface Plot should show a contour map.
For those who have installed matplotlib and numpy, it’s possible to get a quick image…
from pylab import *
import numpy as np
import matplotlib.cm as cm
import matplotlib.pyplot as plt
if __name__ == '__main__':
f = srtmParser()
f.parseFile(r"c:\N55W004.hgt")
zzz = np.zeros((1201,1201))
for r in range(0,1201):
for c in range(0,1201):
va=f.z[(1201*r)+c]
if (va==65535 or va<0 or va>2000):
va=0.0
zzz[r][c]=float(va)
# logarithm color scale
zz=np.log1p(zzz)
imshow(zz, interpolation='bilinear',cmap=cm.gray,alpha=1.0)
grid(False)
show()
Hi I try to execute the first code above, it gave me this error
Traceback (most recent call last):
File “p.py”, line 14, in
rowdump = “,”.join([str(z) for z in thisrow])
TypeError: ‘int’ object is not iterable
I don’t know much about python so can you help me out here?
this is the version of python I am running
Python 2.5.1 (r251:54863, Feb 6 2009, 19:02:12)
[GCC 4.0.1 (Apple Inc. build 5465)] on darwin
I change the code to following, basically took the class structure out
import struct
# read 1,442,401 (1201×1201) high-endian
# signed 16-bit words into self.z
fi=open(‘/Users/sleeplessboy/Sites/py/1.hgt’,”rb”)
contents=fi.read()
fi.close()
z=struct.unpack(“>1442401H”, contents)
if z :
fo=open(‘/Users/sleeplessboy/Sites/py/1.csv’,”w”)
for row in range(0,1201):
offset=row*1201
thisrow = z[offset+1201]
rowdump = “,”.join([str(z) for z in thisrow])
fo.write(“%s\n” % rowdump)
fo.close()
i see the problem… it’s this line :-
You need to change this back to the original
Although the array can be thought of as a 1201 x 1201 array, it is actually a 1442401 x 1 array.
The first 1201 elements are the first row, the next 1201 are the second row, and so on. Offset is a marker which points to the start of the current row.
the z[start:end] is called a split – it creates an array from part of
another array. In this case, it’s lifting a whole row’s worth of values
(1201 values, starting at offset).
The line which fell over…
…failed because ‘thisrow’ was a single int value, not an array/list.
Hope that helps!
Sorry I was trying something else and gave you the wrong error message, this is the one that I could not figure out what’s wrong with it.
Traceback (most recent call last):
File “p.py”, line 13, in
thisrow = z[offset:offset+1201]
TypeError: ‘int’ object is unsubscriptable
import struct
# read 1,442,401 (1201×1201) high-endian
# signed 16-bit words into self.z
fi=open(‘/Users/sleeplessboy/Sites/py/1.hgt’,”rb”)
contents=fi.read()
fi.close()
z=struct.unpack(“>1442401H”, contents)
if z :
fo=open(‘/Users/sleeplessboy/Sites/py/1.csv’,”w”)
for row in range(0,1201):
offset=row*1201
thisrow = z[offset:offset+1201]
rowdump = “,”.join([str(z) for z in thisrow])
fo.write(“%s\n” % rowdump)
fo.close()
Okay, managed to reproduce this error!
It’s a combination of two things – removing the use of classes, and a bad choice of variable names by me
The problem is that the variable name z is used twice; once as an array (the list of elevations), and again as a variable to step through the array.
In the original, it’s not a problem. The two z’s are different because of variable scope; one is self.z, the other is z. Python can tell the two apart; the first belongs to the SRTMParser object, the second is a local variable.
Change this line
to
and it should work.
Or you could re-use the original class as it stands, like this…
I found this python script very helpful any chance you can come up with a version for the NED data format?
Thanks Andrew,
The USGS NED (National Elevation Dataset) is available in a number of formats. BIL is a raw binary file in raster format, very similar to .hgt
Unlike hgt, these files can be any size, and not necessarily square, if you download them from The National Map Seamless Server
It only takes a few tweaks to get this code to read .BIL – the first is to change the way the bytes are read, the second is to allow the code to support arbitrary width and height.
I’ll post a follow-up shortly.
Update
This post shows how to parse BIL files. Hope that helps!