Mapnik GB Projection

Although there is loads of documentation in how to generate your own customised maps, most of it assumes you want everything in the 900913 (Google) projection as used by Google Maps and the like.

These instructions are for those of us who need to create some tiles on a different projection - in the example the British National Grid.

'Please Note' Some steps may not be done in the best way possible but it all works. These instructions also assume you already have a working Mapnik install and have successfully created a map in the standard Google projection, and were written based on my experience with 32bit Ubuntu 8.10 on a machine with less than 500MB of RAM.

In theory it should be a simple matter of enabling support for the UK projection ([EPSG:27700]) in the database, amending the generate_image.py file to use the UK geocodes (instead of Lat and Long), and one change in the OSM.xml file - unfortunately I was unable to get this to work - However you may as well give it a try - you'll need to make all these changes anyway and the bugs may have been machine specific.

As a quick background, the UK national grid is based on Geocodes - that is the number of meters North and East of a point south west of the Scilly Isles. This means all of Great Britain has positive coordinates - much of Éire is to the west of the zero line and as such has negative geocodes, the west coast of Ireland also starts to look rather odd as well as the grid is a best fit to Great Britain (the Ordnance Survey of Northern Ireland use a different grid which is a better fit for the whole of the island).

Update the database to support British National Grid

Firstly we need to ensure that the database supports the British national grid, to do this create a sql file

gedit ~/osm2pgsql/27700.sql &

Which contains the following (on one line)

INSERT into spatial_ref_sys (srid, auth_name, auth_srid, proj4text, srtext) values ( 927700, 'epsg', 27700, '+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +datum=OSGB36 +units=m +no_defs ', 'PROJCS["OSGB 1936 / British National Grid",GEOGCS["OSGB 1936",DATUM["OSGB_1936",SPHEROID["Airy 1830",6377563.396,299.3249646,AUTHORITY["EPSG","7001"]],AUTHORITY["EPSG","6277"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4277"]],UNIT["metre",1,AUTHORITY["EPSG","9001"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",49],PARAMETER["central_meridian",-2],PARAMETER["scale_factor",0.9996012717],PARAMETER["false_easting",400000],PARAMETER["false_northing",-100000],AUTHORITY["EPSG","27700"],AXIS["Easting",EAST],AXIS["Northing",NORTH]]');

Now add it to the database

sudo -u postgres -i
psql -d gis -f ~/osm2pgsql/27700.sql

Amend the XML stylesheet

Firstly create a copy of the osm.xml stylesheet file

cp ~/mapnik/osm.xml ~/mapnik/osmgb.xml

Now change the srs attribute on the third line of the xml file to

+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

Create your image generation script

gedit ~/mapnik/generate_british.py &

#!/usr/bin/python
from mapnik import *
import sys, os
def drawMap(filename,west,south,east,north):
    print(filename)
    sz = 250
    m = Map(sz,sz,"+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")
    load_map(m,"osmgb.xml")
    bbox = Envelope(west,south,east,north)
    m.zoom_to_box(bbox)
    im = Image(sz,sz)
    render(m, im)
    render(m, im)
    view = im.view(0,0,sz,sz) # x,y,width,height
    view.save(filename,'png')
drawMap('preston.png',340000,420000,360000,440000)

execute the script

cd ~/mapnik/
python generate_british.py

Assuming everything has worked okay you will have an image looking similar to the following.

It doesn't work!

However when I run this I only saw roads for around the first 20 pixels and a torn coastline - this is either a bug in Mapnik or a lack of memory on your PC. However it's not the end of the world, and we can work round this.

Importing data into the PostGIS database

Firstly we need to import the data into the database, but instead of the standard Spherical projection we want to import it in the UK geocode format.

./osm2pgsql -d gis --proj 27700 ~/osm/great_britain.osm.bz2
./osm2pgsql -d gis -a --proj 27700 ~/osm/ireland.osm.bz2

(the -a on the Ireland import means append this data to what's already in the database).

Possible problem : Buffer Overflow

Unfortunately I was getting some buffer overflow errors, I was able to fix this by recompiling the osm2pgsql script with the debug flags turned on.

Firstly we want to download the source code

$ cd ~/
$ svn co http://svn.openstreetmap.org/applications/utils/export/osm2pgsql/
$ cd osm2pgsql

Now amend the Makefile

gedit Makefile &

Do this by changing both occurrences of

-g -O2

to

-g -ggdb

Now compile the program

make clean
make

Now running the following commands should import all the data in the British national grid projection

~/osm2pgsql/osm2pgsql -d gis --proj 27700 ~/osm/great_britain.osm.bz2
~/osm2pgsql/osm2pgsql -d gis -a --proj 27700 ~/osm/ireland.osm.bz2

Stylesheet Changes

If you try and generate the images now you'll find that you have a coastline but no roads. This is a simply because we need to amend the stylesheet file to tell it our database is now in the British national grid projection.

In the xml file look for all occurrences of <Parameter name="table> then change the associates srs attribute to

srs="+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"

But my coastline is missing sometimes

This appears to be down to a lack of memory as only those areas with a complex coastline are effected - the easiest way is to re-project the shape files you stores in the "world_boundaries" folder.

• Download and install Map Window GIS on a Windows PC
• Copy the world_boundaries_m files and builtup_area files to a windows box.
• Start Map Window GIS
• Click the Green Plus, select Add Layer and select world_boundaries_m
• Click the "Select" icon and draw a box round the British Isles (I included the Faroe Islands, Norway, northern France and Benelux to get a complete map when you are showing the whole of the UK)
• Click GIS Tools, Click Vector, Click Export Selected Shapes to New Shape File
• Give it a sensible name (e.g. british_isles_landmass)
• Repeat for the builtup_area shape file.
• Click the new button
• Click "GIS Tools", "Vector", "Reproject a shape file"
• Select british_isles_landmass.shp
• In the dialogue that appears select
  • Category : Projected Coordinate Systems
  • Group : National Grids
  • Name : British National Grid
• Click OK
• This should take a few seconds, at which point a "Finished, Add to map?" dialogue will be shown - Click "NO" otherwise it screws up.
• Repeat for british_isles_urban.shp
• Now copy the files with _reprojected in the file names back to the world_boundaries folder on your Ubuntu box.

As with the reprojected database we need to amend the stylesheet - search for all occurrences of <Parameter name="shape> and then change the attributes as follows

* Layer -> srs = +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.999601 +x_0=400000 +y_0=-100000 +ellps=airy  +units=m +nodefs
* Parameter  -> file = ** Your Reprojected Shapefile Filenames **

You will also need to amend the Style rules for world and builtup to have the following Scale Denominators

• MaxScaleDenominator = 60000000
• MinScaleDenominator = 0

Now if you try and generate an image everything should work fine - you can now start customising your style sheet to your requirements.

Thanks

With thanks to User:Matt for all the assistance he has given me in getting everything working, and everyone else on the mail lists who helped point me in the right direction