Archive

Posts Tagged ‘DEM’

Parsing USGS BIL Digital Elevation Models in Python

A false colour elevation map of the Hawaiian island of Kahoolawe

The USGS have a digital elevation model covering most of the US, Hawaii and Puerto Rico. If you use the National Map Seamless Server it’s possible to select a region and download a digital elevation model in a number of formats. One of these, BIL, is a simple binary raster file, which can be read quite easily in Python – in this, I’ll base this on code from a previous post.

Getting the data

  • Go to National Map Seamless Server
  • The map is fairly intuitive, but if you get stuck check out their tutorial
  • When you use the download rectangle tool, make sure to click on “Modify data request” and choose the BIL data format (BIL_INT16) rather than the default “Arc_GIS”
  • Choose zip or TGZ
  • You’ll then download a compressed archive – may take a while
  • The archive contains lots of files; look at the contents of output_parameters.txt. This tells you the dimensions of the file. Make a note of the sizes.
  • Extract the .bil file – this is what has the raw data.

Extracting the binary elevation data from the BIL file

# USGS Seamless Server BIL Format Parser
#

import struct

# get these from the output_parameters.txt file
# included in the download
height=437
width=663

# where you put the extracted BIL file
fi=open(r"my-extracted.bil","rb")
contents=fi.read()
fi.close()

# unpack binary data into a flat tuple z
s="<%dH" % (int(width*height),)
z=struct.unpack(s, contents)

After running this, z is now a tuple of width*height elements, with data stored in raster order (entries in row 1 from left-to-right, then row 2 from left-to-right, and so on)

The struct module for reading binary data

struct is a fantastic Python library, it really makes reading binary files pretty easy. Say we have a raster-layout binary file of 50 rows by 100 values per row. We’d expect 5000 values. Each value is an unsigned short (H) and little-endian (<)

In this case, we build up a format specification like this
“<5000H”

Plotting the end result

And optionally, if you have matplotlib and numpy libraries, you can get a quick preview of the resulting array.

from pylab import *
import numpy as np
import matplotlib.cm as cm

heights = np.zeros((height,width))
for r in range(0,height):
    for c in range(0,width):
        elevation=z[((width)*r)+c]
        if (elevation==65535 or elevation<0 or elevation>20000):
            # may not be needed depending on format, and the "magic number"
            # value used for 'void' or missing data
            elevation=0.0
        heights[r][c]=float(elevation)
imshow(heights, interpolation='bilinear',cmap=cm.prism,alpha=1.0)
grid(False)
show()

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

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