Posts Tagged ‘Python’

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

# where you put the extracted BIL file

# 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

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 as cm

heights = np.zeros((height,width))
for r in range(0,height):
    for c in range(0,width):
        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
imshow(heights, interpolation='bilinear',cmap=cm.prism,alpha=1.0)

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 as cm

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


for line in fi:
    if linenum>0:
        line=string.replace(line, "\n","")
            print "Error!"

m = Basemap(llcrnrlon=-8.0,llcrnrlat=54.5,urcrnrlon=1.5,urcrnrlat=59.5,lat_ts=20,
m.drawmapboundary(fill_color='white') # fill to edge

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 
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')
    if pages>8:
        # stop after 8 pages of 500 images
        # not sure if groups.pools.getPhotos has the same
        # 4000 image limit as
    print "Longitude,Latitude"
    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"):
                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))

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

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


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

.. 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 as cm

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


for line in fi:
    if linenum>0:
        line=string.replace(line, "\n","")

# 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,
m.drawmapboundary(fill_color='cyan') # fill to edge
m.drawrivers() # you may want to turn this off for larger areas like continents

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

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,

…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

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 as cm


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.


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

Finally, the scatter plot.

plt.title("Flickr Geotagging Counts with Basemap")

Image Steganography with PIL

October 7, 2009 Leave a comment

Steganography is greek for ‘hidden writing‘; the act of hiding a message inside another message.

In this case, hiding an image inside another image, without it being obvious to the viewer. The example I’ll give here is only a ‘toy’ implementation, for two reasons:-

  • easily cracked: it wouldn’t take the authorities long to spot the hidden message, not least because the algorithm is described on wikipedia 😉
  • fragile: the hidden image-within-an-image can easily be broken, if the image has its colours changed afterwards.

But it does illustrate how to do bitwise-manipulation of images in PIL using the ImageMath module, which is the purpose of the post.

How it works

The watermark – the image we wish to hide – is a bitonal image, with black and white pixels only. It’s then resized to be the same size as the original image.

We ‘smuggle’ the watermark inside the original by replacing the LSB (least significant bit) of each colour channel (R,G and B) in the original with the corresponding pixel in the watermark – either 1 for white, or 0 for black.

This image shows the binary arithmetic…least significant bit on the right.


Hiding the watermark image inside our image

For this, we’ll need these imports…

from PIL import Image, ImageMath

and open the two files. The watermark is scaled to match the size of the original image."c:\watermark.png")"c:\original.jpg")

ImageMath only works with single channel (greyscale) images, so we need to split the two images into their three channels (Red, Green and Blue) using the split() method.

red, green, blue = original.split()
wred, wgreen, wblue = watermark.split()

Now, using ImageMath. ImageMath lets you write simple expressions using values from one or more images. Here, ‘a’ and ‘b’ are bound to the values in the original and watermarked images, respectively. The convert() call is needed to prevent problems later; we need to cast the results back to a greyscale image (mode ‘L’).

red2 = ImageMath.eval("convert(a&0xFE|b&0x1,'L')", a=red, b=wred)
green2 = ImageMath.eval("convert(a&0xFE|b&0x1,'L')", a=green, b=wgreen)
blue2 = ImageMath.eval("convert(a&0xFE|b&0x1,'L')", a=blue, b=wblue)

Okay, so now we have three channels whose LSBs have been replaced with the LSB of the watermark.

But we need to combine the 3 channels back to get an RGB image ready for saving.

out = Image.merge("RGB", (red2, green2, blue2))"c:\merged.png")

Open the original and the processed images; can you see any difference?

Extracting the hidden image

All this is for nought if you can’t extract the hidden image afterwards.

This is simpler, as we only need to produce a black/white image from the LSB of the image. Here, I’ve only bothered with the Red channel."c:\merged.png")
red, green, blue = stegged.split()
watermark=ImageMath.eval("(a&0x1)*255",a=red) # convert to 0 or 255

Using Flickr API to get the views, faves and comments of your most popular images

September 25, 2009 Leave a comment

One of the first things I wanted to do with the Flickr API was to get some stats on my most popular images.

You can get this info through the web front end, but there’s no option to download the stats in delimited format (such as CSV) so it can be analysed in a spreadsheet.

I wanted to work out if there was a pattern emerging in the key stats for my 200 most popular images…

  1. Number of Views
  2. Number of Favourites
  3. Number of Groups posted to
  4. Number of sets an image is in

Using a Python script (v2.5) and Beej’s FlickR API, this is fairly straightforward. It doesn’t require authentication.

The script runs slowly as it ‘plays nice’, leaving a seconds pause between calls, courtesy of the time.sleep() function. I don’t want to thrash the server.

# -*- coding: UTF8 -*-

import flickrapi
import datetime
import time
import string

# enter your api key below
api_key = 'PUT_YOUR_API_KEY_HERE' 

# enter the user id below (you can use flickr.people.findByUsername to get this for any user)
# it'll look something like 99999999@N99

# delimiter. Use comma if you want, I tend to use ~

# dump number of views in delimited format

if __name__ == '__main__':
    #output format : "photoid,title,views,faves,groups,sets"
    flickr = flickrapi.FlickrAPI(api_key)
    photos = flickr.photos_Search(user_id= userid,extras='views', per_page='200', page='1', sort='interestingness-desc')
    for photo in photos.find('photos'):

        title = string.replace(photo.get('title'),",","") #in case you want to use comment as a delimiter ;0)
        # number of views
        id = photo.get('id')
        # fave count (up to 50)
        faves = flickr.photos_getFavorites(photo_id=id,per_page=50)
        # pools and sets posted to
        # comments
        # output as delimited text
        print converted

This script dumps to the console, rather than a file; but it’s easily modified to write to a file. It should work with comma as a separator (for CSV use) as the title tag is stripped of commas…

Here’s some sample output from my photostream…

The format is :
photo id~title~views~favourites~groups~sets

3598511429~paris photo heatmap~736~6~12~5~3
3609118442~heart texture~861~10~26~0~2
2304836447~persistence de-motivator~4964~1~4~1~1
3688253826~St Anthony's Chapel Edinburgh~124~15~15~9~1
2717978614~st marys~88~9~7~3~2

Once you have the output saved to a text file, you can import it into a spreadsheet (like OpenOffice or Excel) and play around with the figures 🙂

mapping flickr group activity

September 20, 2009 Leave a comment

flickr group activity visualization

Mapping the activity levels of approx 1,500 Flickr groups against the number of members of each group, using the Flickr API.

The x axis is the group size, the y axis is the number of seconds an image can expect to stay on the ‘front page’. This was measured as the timestamp difference between the 1st and 13th images in the pool (the landing page for a group shows 12 images; the user has to follow a link to show more).

Note that the image is a log-log plot, as the groups follow a power-law distribution.

Groups in the bottom right are the busiest – this includes the B&W and FlickrCentral groups. Groups in the top right have lots of members but are more vigorously moderated, so images are added to the pool more slowly (or new submissions are deleted). The group touching the bottom of the graph is the “10 million images” group, where users are encouraged to “dump and run”.

This is a hexbin plot – the colour represents the number of groups falling within a certain range of values. Red=Lots, Green=Fewer, Blue=Few.

As you’d expect, larger groups tend to have a higher turnaround of images, but there’s a lot of variation.

The most common group size seems to be around 2000-3000 members; a group this size, you can expect an image to stay on the front page for around 2-3 hours. With the largest groups, this drops to around 5 minutes.

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

Creating a Unicode character map in Python

September 17, 2009 Leave a comment

When I were a lad, you had ASCII, and that was it. Or EBCDIC.. but the less said about that the better.

Unfortunately ASCII is useless when you want to start adding characters from non-latin scripts, such as… well, pretty much every language which isn’t English. That’s where Unicode comes in.


The following Python script creates a web page showing a table of Unicode characters. Invalid Unicode characters appear as grey squares; hover over a character to get more info (a description of the character and its hex code).

# -*- coding: UTF8 -*-
Created on 15 Sep 2009

@author: Steven Kay
import unicodedata

# generate table of Unicode characters in Python
# using the unicodedata module

fo.write("<table border=1 style=\"text-align:center\">")
fo.write("".join("<td>%s</td>"%x for x in " .".split("."))) #header
for page in range(0,1024):
    s="<td>%s</td>" % hex(page)
    for cell in range(0,16):
            # valid Unicode char
            s=s+"<td>"+"<span title=\"0x%0X\n%s\">" % (x,name) +unichr(x)+"</span></td>"
            # not a valid Unicode char
            s=s+"<td style=\"background:#A99\">&nbsp;</td>"

You won’t be able to see all characters; many fonts only support the Latin characters (A-Z,a-z,0-9 and punctuation characters).

The script also stops after 16k, so it won’t pick up all the characters in some eastern languages. Easy enough to change though!

If you’re using Windows XP, you’ll need to install Eastern fonts (using Regional and Language Options/Languages on the Control Panel) to see characters in languages such as Japanese, Thai, Chinese etc. I normally have this turned off as it takes up a fair amount of memory…

Categories: Python Tags: , , , ,

generating ASCII art from photographs in Python

September 8, 2009 Leave a comment

I’m old school – I have vague recollections of seeing an ASCII art printout of Snoopy printed on mainframe feed paper as a child back in the 70s. I think I drew on top of this with Crayola crayons.

I thought it would be fun to try to recreate this with photos.

ASCII Art works by using bog-standard ASCII characters to build up an image. It’s the creative use of positive and negative space in typography; some characters use more ink (like #) and some use less (like .), and by using the ‘weight’ of characters, you can build up an image.


What we need to do is to shrink the image to thumbnail size, then convert the image to monochrome.

With Python, the standard image processing library is PIL. Creating a monochrome thumbnail is trivial:-"c:\test.jpg")
im=im.resize((75, 75), Image.ANTIALIAS)
im=im.convert("L") # convert to mono

Then, it’s a matter of splitting the luminosity values into 7 bands, and assigning each pixel a random character from a group of ascii characters of similar “optical weight”.

Here’s the script…

ASCII Art maker
Creates an ascii art image from an arbitrary image
Created on 7 Sep 2009

@author: Steven Kay

from PIL import Image
import random
from bisect import bisect

# greyscale.. the following strings represent
# 7 tonal ranges, from lighter to darker.
# for a given pixel tonal level, choose a character
# at random from that range.

greyscale = [
            " ",
            " ",

# using the bisect class to put luminosity values
# in various ranges.
# these are the luminosity cut-off points for each
# of the 7 tonal levels. At the moment, these are 7 bands
# of even width, but they could be changed to boost
# contrast or change gamma, for example.


# open image and resize
# experiment with aspect ratios according to font"c:\test.jpg")
im=im.resize((160, 75),Image.BILINEAR)
im=im.convert("L") # convert to mono

# now, work our way over the pixels
# build up str

for y in range(0,im.size[1]):
    for x in range(0,im.size[0]):

print str

The text can then be pasted into a text editor or web page, but you’ll need to use a monospaced font, where each character takes the same amount of space horizontally. And depending on your editor, you may need to turn off word wrap 🙂

It’s amazing how much difference there is between monospaced fonts; some work better than others. Some fonts are ‘wider’ than others, so it’s worth experimenting with aspect ratios.

Categories: Python Tags: , ,