Heights & OS

Working with Digital Elevation Models (DEM) is an interesting extension to creating maps, which are usually a flat representation of part of the world. I really want to find a way to show elevation in a way that is a bit different from a flat map. Working with the OS DEM data has whet my appetite to try something new, but first I need a map to work with.

The most detailed DEM data I have is based on Ordnance Survey OpenData, so creating a map in the OS projection will be useful. I use TileMill to create maps from OSM data.

Firstly I needed OSM data in the Ordnance Survey projection. That means loading some OSM data into a fresh PostgreSQL database. I created a PostgreSQL database and, as usual, add the extension for PostGIS. this creates a table called spacial_ref_sys that includes the OS projection, amongst many others. I often add the hstore extension too, but this is a simple map so I didn't need it.

createdb -E UTF8 EYOS
echo "CREATE EXTENSION postgis;" | psql -d EYOS

I loaded an extract of OSM data using the usual osm2pgsql utility except the projection was needed too to convert the data to OS projection as it is loaded.

osm2pgsql --slim -d EYOS -C 1024 ey.osm.pbf --proj 27700

I decided to add a coastline, so that needed to be in OS projection too. OSM coastlines are handled differently from all other data. They are extracted from the main DB, checked for consistency and created into a shapefile for the world. This is known as processed_p.shp. I have my own copy with a cut-down version with only the British Isles in it to make rendering a bit quicker. I reprojected that to a copy in OS projection using OGR2OGR, part of the Geospatial Data Abstraction Library. 

ogr2ogr -t_srs 'EPSG:27700' -s_srs 'EPSG:3857' coast_bi_os.shp coast_bi.shp

Armed with all of this I could now start TileMill and add the layers I need for the map. Each of the layers, including the coast shapefile, needed a custom projection. This is:

+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.999601 +x_0=400000 +y_0=-100000 +ellps=airy  +units=m +towgs84=446.448,-125.157,542.060,0.1502,0.2470,0.8421,-20.4894 +units=m +nodefs

 I got this from the PostgreSQL table postgis created above. Once I had designed the map as I wanted it I exported the Mapnik XML and ran it through Mapnik. I discovered that the Mapnik XML was not quite right. It needed to have the third line changed so the srs part matches the custom projection above. There doesn't seem to be a way to set this in TileMill, so a manual edit was needed.

To run the Mapnik XML through Mapnik I used the following python code:

#!/usr/bin/python

# generate a map image in OS projection epsg:27700

import mapnik
import sys, os
def drawMap(filename, west,south,east,north):
    print(filename)
    sz = 5000
    ossrs = "+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +towgs84=446.448,-125.157,542.06,0.1502,0.247,0.8421,-20.4894 +units=m +no_defs"
    m = mapnik.Map(sz,sz,ossrs)
    mapnik.load_map(m,"osproj.xml")
    bbox = mapnik.Envelope(west,south,east,north)
    m.zoom_to_box(bbox)
    im = mapnik.Image(sz,sz)
    mapnik.render(m, im)
    view = im.view(0,0,sz,sz) # x,y,width,height
    view.save(filename,'png')

if __name__ == '__main__':
    drawMap('cott.png',500000,430000,510000,440000) 

 That long-winded projection was needed again. Notice the coordinates in the drawMap function are OS coordinates, not longitude and latitude. Everything must match the chosen projection.

This gives me an image of the map in the OS projection, but the style could be any style you choose, though I'd be wary of copying the OS style too closely. This will now match the DEM data if they are combined. My style is still a bit stark and only renders a few objects, but it is something to work with.

Next I need to use it imaginatively.