Archive

Posts Tagged ‘Mapping’

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

Plotting points on an OpenStreetMap Export

February 24, 2010 Leave a comment

This map shows the location of pictures in a Flickr group superimposed on an OpenStreetMap (OSM) export.

If you try plotting points directly on an OSM map, you’ll find that points are all over the shop. The reason is that OSM exports use the Mercator projection; you need to change the latitude and longitude coordinates into the Mercator projection.

Basemap to the rescue!

If you use the code in the previous post, you can change the plotting code.

from basemap import Basemap 
import matplotlib.pyplot as plt
import matplotlib.pylab as pylab
import numpy as np
import string
import matplotlib.cm as cm

x=[] #longitudes
y=[] #latitudes

fi=open(r'C:\infoviz\scotland_photos.csv','r')

linenum=0
for line in fi:
    if linenum>0:
        line=string.replace(line, "\n","")
        try:
            fields=string.split(line,",")
            lon,lat=fields[0:2]
            x.append(float(lon))
            y.append(float(lat))
        except:
            print "Error!"
    linenum+=1
fi.close()

m = Basemap(llcrnrlon=-8.0,llcrnrlat=54.5,urcrnrlon=1.5,urcrnrlat=59.5,lat_ts=20,
            resolution='h',projection='merc',lon_0=-4.36,lat_0=54.5)
x1,y1=m(x,y)
m.drawmapboundary(fill_color='white') # fill to edge
m.scatter(x1,y1,s=5,c='r',marker="o",cmap=cm.jet,alpha=1.0)

Rather than using basemap to draw the outline of Scotland, this script simply creates a scatter plot on a white background, like so:-

Now, you need to export the map from OpenStreetMap.

The corners have been set as follows..

… llcrnrlon=-8.0,llcrnrlat=54.5,urcrnrlon=1.5,urcrnrlat=59.5 …

llcrn stands for the lower-left coordinate, and and urcrn for the upper-right coordinate.

So if you if you export from OpenStreetMap using these coordinates…

.. you’ll have a map of Scotland with the Mercator projection.

The two images can then be composited in Photoshop or another graphics app like the Gimp.

Creating a pinboard map of geotagged photos in a flickr pool

January 14, 2010 1 comment

In this post I’ll show how to produce a simple pinboard map of geotagged photos in a flickr group pool, using Python and Basemap/Matplotlib. You’ll need:-

There are two short scripts here:-

  • A script to find the longitude and latitude of geotagged photos in the group pool
  • A script to generate the plot

The first script produces a CSV file; the second uses this CSV file to produce the plot.

Here’s the script to produce the CSV file with photo locations:-

# -*- coding: UTF8 -*-
'''
Created on 12 May 2009
Based on beejs flickr API
Produce a list of photo locations for a given
group's pool on flickr
@author: Steven Kay
'''

import flickrapi
import string
import datetime
import string
import time

# Enter your API key below
# You can apply for an API key at 
# http://www.flickr.com/services/apps/create/apply
api_key = '' 

# paste group NSID below
group = '1124494@N22'

# sample... fetch your latest images with 
# a count of the views, faves and comments

if __name__ == '__main__':
    flickr = flickrapi.FlickrAPI(api_key)
    response_photos = flickr.groups_pools_getPhotos(group_id=group,per_page=500,extras='geo')
    root=response_photos.findall('.//photos')
    pages=int(root[0].get('pages'))
    if pages>8:
        # stop after 8 pages of 500 images
        # not sure if groups.pools.getPhotos has the same
        # 4000 image limit as photos.search..?
        pages=8
    
    fo=open(r"C:\infoviz\scotland_photos.csv","w")
    print "Longitude,Latitude"
    fo.write("Longitude,Latitude\n")
    for page in range(0,pages):
        response_photos = flickr.groups_pools_getPhotos(group_id=group,per_page=500,page=str(page),extras='geo') 
        for photo in response_photos.findall(".//photos/photo"):
            try:
                lat=photo.get('latitude')
                lon=photo.get('longitude')
                st="%s,%s" %(lon,lat)
                if not st=="0,0":
                    # ignore the odd buggy 0,0 coords
                    print "%s,%s" %(lon,lat)
                    fo.write("%s,%s\n" %(lon,lat))
            except:
                pass
        time.sleep(1)
    fo.close()

You’ll need to find the NSID of the group as an input; you can find this with the flickr API call flickr.group.search.

Now, you have a simple CSV file with the latitude and longitude of each geotagged image in the pool.

Longitude,Latitude
-5.167792,58.352519
-4.024359,57.675544
-4.230251,57.497356
-4.2348,57.501045
-4.84703,56.646034
-4.306168,55.873986
-3.586263,56.564732
...

This demo uses the Photography Guide to Scotland pool.

The next step is to plot the map.

'''
Simple Matplotlib/Basemap pinboard map for
Flickr Groups.

Need to provide a CSV file in following format

Longitude,Latitude
20.1,-3.25
20.225,-3.125
.. etc..

Created on 10 Oct 2009

@author: Steven Kay
'''

from basemap import Basemap 
import matplotlib.pyplot as plt
import matplotlib.pylab as pylab
import numpy as np
import string
import matplotlib.cm as cm

x=[] #longitudes
y=[] #latitudes

fi=open(r'C:\infoviz\scotland_photos.csv','r')

linenum=0
for line in fi:
    if linenum>0:
        line=string.replace(line, "\n","")
        try:
            fields=string.split(line,",")
            lon,lat=fields[0:2]
            x.append(float(lon))
            y.append(float(lat))
        except:
            pass
    linenum+=1
fi.close()

# cass projection centred on scotland
# will need to replace with a projection more suited
# to the group you're plotting

m = Basemap(llcrnrlon=-8.0,llcrnrlat=54.5,urcrnrlon=1.5,urcrnrlat=59.5,
            resolution='h',projection='cass',lon_0=-4.36,lat_0=54.5)
x1,y1=m(x,y)
m.drawmapboundary(fill_color='cyan') # fill to edge
m.drawcountries()
m.drawrivers() # you may want to turn this off for larger areas like continents
m.fillcontinents(color='white',lake_color='cyan',zorder=0)
m.scatter(x1,y1,s=5,c='r',marker="o",cmap=cm.jet,alpha=1.0)

plt.title("Photography Guide to Scotland in FlickR") # might want to change this!
plt.show()

This script uses a projection centred around scotland; you’ll need to change the following line…

m = Basemap(llcrnrlon=-8.0,llcrnrlat=54.5,urcrnrlon=1.5,urcrnrlat=59.5,
            resolution='h',projection='cass',lon_0=-4.36,lat_0=54.5)

…to something more suitable for your needs. Basemap provides an intimidating list of projections which should meet your needs.

Scatter plots with Basemap and Matplotlib

October 12, 2009 Leave a comment

flickr-geotagging-with-base
A while back I used the flickr api to map 24 hours worth of geotagged photos.

My previous attempts needed some manual Photoshop work to superimpose the plots on a map. The next logical step is to do the whole process – from start to finish – in code, and remove the manual steps.

To do this, I tried the awesome Basemap toolkit. This library allows all sorts of cartographic projections…

Installing Basemap

Basemap is an extention available with Matplotlib. You can download it here (under matplotlib-toolkits)

I installed the version for Python 2.5 on Windows; this missed out a dependency to httplib2 which I needed to install separately from here.

Getting started

Let’s assume you have 3 arrays – x, y and z. These contain the longitudes, latitudes, and data values at each point. In this case, I binned the geotagged photos into a grid of degree points (360×180), so that each degree square contained the number of photos tagged in that degree square.

Setting up

from basemap import Basemap 
import matplotlib.pyplot as plt
import numpy as np
import string
import matplotlib.cm as cm

x=[]
y=[]
z=[]

Now, you need to populate the x,y and z arrays with values. I’ll leave that an exercise to you 🙂 All three arrays need to be the same length.

Now, you need to decide which projection to use. Here, I’ve used the Orthographic projection.

m = Basemap(projection='ortho',lon_0=-50,lat_0=60,resolution='l')

Here is the secret sauce I took a while to work out. That’ll teach me not to R the FM. This line transforms all the lat, lon coordinates into the appropriate projection.

x1,y1=m(x,y)

The next bit, you can decide which bits you want to plot – land masses, country boundaries etc.

m.drawmapboundary(fill_color='black') # fill to edge
m.drawcountries()
m.fillcontinents(color='white',lake_color='black',zorder=0)

Finally, the scatter plot.

m.scatter(x1,y1,s=sizes,c=cols,marker="o",cmap=cm.cool,alpha=0.7)
plt.title("Flickr Geotagging Counts with Basemap")
plt.show()

visualising the nationality of Nobel Peace Prize Winners

October 9, 2009 Leave a comment

Visualizing the nationality of Nobel Peace Prize winners over time

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