Archive

Posts Tagged ‘Elevation’

mapping water drainage with SRTM elevation data in Python

September 18, 2009 Leave a comment

srtm enhanced relief map

With this mini-project, I wanted to try to model the effect of rainfall on the landscape, to see if I could map water drainage – finding watersheds, catchment areas and so on. After lots of tries, I stumbled across a potentially interesting mapping technique by accident, which does a good job of providing a relief map that emphasizes water courses.

The algorithm starts with a square grid of SRTM digital elevation data. (You can read more about this in a previous post)…

For each cell in the grid, a virtual raindrop falls. This drop then repeatedly moves to the lowest neighbouring cell, each time keeping a running total of how many meters it has descended by. Eventually, it falls into a local minima; a ‘well’, a cell with no neighbours lower than itself. This can be the coast, or it might be a lake or a low-lying piece of land in a valley. When it can’t fall any further, the map is coloured according to the total vertical distance travelled.

The height of the cell is relative to the nearest patch of flat land, rather than being relative to sea level. This makes it easier to follow the lie of the land.

The image here is a map of the degree square covering Edinburgh and the Scottish Borders (Dumfries is in the south, Edinburgh in the top-right corner).

Advertisements

Beginning Digital Elevation Model work with Python

September 5, 2009 7 comments

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