Home > Python > Beginning Digital Elevation Model work with Python

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,

SRTM DEM

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()
Advertisements
  1. sleeplessboy
    September 26, 2009 at 4:32 am

    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()

    • September 26, 2009 at 1:46 pm

      i see the problem… it’s this line :-

      thisrow = z[offset+1201]
      

      You need to change this back to the original

      thisrow = z[offset:offset+1201]
      

      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…

      rowdump = “,”.join([str(z) for z in thisrow])
      

      …failed because ‘thisrow’ was a single int value, not an array/list.

      Hope that helps!

  2. sleeplessboy
    September 26, 2009 at 3:38 pm

    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()

    • September 27, 2009 at 7:18 am

      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

      rowdump = “,”.join([str(z) for z in thisrow])
      

      to

      rowdump = “,”.join([str(foo) for foo in thisrow])
      

      and it should work.

      Or you could re-use the original class as it stands, like this…

      srtmparser = srtmParser()
      srtmparser.parseFile(r"/Users/sleeplessboy/Sites/py/1.hgt") 
      srtmparser.writeCSV(r"/Users/sleeplessboy/Sites/py/1.csv")
      
  3. andrew
    May 26, 2010 at 2:25 am

    I found this python script very helpful any chance you can come up with a version for the NED data format?

    • May 29, 2010 at 5:27 pm

      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!

  1. May 29, 2010 at 6:51 pm

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

%d bloggers like this: