OS Vector Map District with PostGIS - Take 3

To recap:

• The entire UK's VectorMapDistrict dataset has been loaded into my postgresql/postgis database using my vmd2pgsql script, which is based around shp2pgsql.
• I created some stylesheets for mapnik rendering of the data based on the mapnik tutorial, but they came out with a transparent background and black lines and text, which is not what I wanted.

The rendering problem turned out to be quite easy - it is just that the mapnik tutorial is based on Mapnik 0.7, not mapnik2.  To make it work with mapnik2 you need to:

1. Replace the bgcolor parameter of the map element with background-color.
2. Remove all of the CSSParameter business in the style definitions with simple name="xyz" constructs.

It is a shame that mapnik did not complain about these issues, even in debug mode, but never mind!

Adding the different linear features works ok - you see something of an improvement from the black and white version to something resembling a map:

So now all I need to do is add the areas (water, woodland towns etc.) to make it look extra pretty, and sort out what to plot at different zoom levels so it is reasonable.

Here I hit a problem.   When I tried to add the settlement_area or naturalfeature_area items, mapnik crashed with an "invalid geometry" error, after taking quite a long time thinking about it!   It is something to do with polygons not being closed properly.   I suspect it is rounding errors, but am not too sure.  One possibility is that I should have used the shp2pgsql option to just use integers, which would have avoided rounding, or I could try to fix the database.  As it took 9 hours to import, and I don't know if integers will work, I'll try to fix it.

With much internet searching I discovered that you can check the validity of each geometric feature with the postgis st_isvalid() function.  Running a simple select statement on the database like:

select gid from naturalfeature_area where not st_isvalid(the_geom);

takes a long time (about an hour I think) and gives me a list of invalid geometries, after pages and pages of warnings about things crossing over themselves.

It seems that one trick that people use to force gemetries to be valid is to use the st_buffer() command with the buffer radius set to zero.   I checked it worked by doing:

select st_isvalid(st_buffer(the_geom,0.0)) from naturalfeature_area where gid=2342149;

(2342149 was the first entry in my list of invalid geometries) - Success - it returns 't' for true meaning the st_buffer trick fixes it.

Now to update the database.  To correct the error I had to do:

update naturalfeature_area set the_geom=st_multi(st_buffer(the_geom,0.0)) where gid=2342149;

Note that I do not really know what the st_multi thing does, but I got errors about geometry constraint violations without it - I think that st_buffer returns the simplest type of geometry it can, and st_multi forces it to be a multipolygon, but I could be wrong!.

Now all I need to do is remember enough about subqueries to write a bit of SQL that does the correction for every invalid geometry - back to that well known search engine....

Well, I have set it going using the following to try to repair the geometries:

update naturalfeature_area set the_geom=st_multi(st_buffer(the_geom,0.0)) where gid in (select gid from naturalfeature_area where not st_isvalid(the_geom));

I think I'll leave it for a few hours....

OS Vector Map District with PostGIS - Take 3

To recap:

• The entire UK's VectorMapDistrict dataset has been loaded into my postgresql/postgis database using my vmd2pgsql script, which is based around shp2pgsql.
• I created some stylesheets for mapnik rendering of the data based on the mapnik tutorial, but they came out with a transparent background and black lines and text, which is not what I wanted.

The rendering problem turned out to be quite easy - it is just that the mapnik tutorial is based on Mapnik 0.7, not mapnik2.  To make it work with mapnik2 you need to:

1. Replace the bgcolor parameter of the map element with background-color.
2. Remove all of the CSSParameter business in the style definitions with simple name="xyz" constructs.

It is a shame that mapnik did not complain about these issues, even in debug mode, but never mind!

Adding the different linear features works ok - you see something of an improvement from the black and white version to something resembling a map:

So now all I need to do is add the areas (water, woodland towns etc.) to make it look extra pretty, and sort out what to plot at different zoom levels so it is reasonable.

Here I hit a problem.   When I tried to add the settlement_area or naturalfeature_area items, mapnik crashed with an "invalid geometry" error, after taking quite a long time thinking about it!   It is something to do with polygons not being closed properly.   I suspect it is rounding errors, but am not too sure.  One possibility is that I should have used the shp2pgsql option to just use integers, which would have avoided rounding, or I could try to fix the database.  As it took 9 hours to import, and I don't know if integers will work, I'll try to fix it.

With much internet searching I discovered that you can check the validity of each geometric feature with the postgis st_isvalid() function.  Running a simple select statement on the database like:

select gid from naturalfeature_area where not st_isvalid(the_geom);

takes a long time (about an hour I think) and gives me a list of invalid geometries, after pages and pages of warnings about things crossing over themselves.

It seems that one trick that people use to force gemetries to be valid is to use the st_buffer() command with the buffer radius set to zero.   I checked it worked by doing:

select st_isvalid(st_buffer(the_geom,0.0)) from naturalfeature_area where gid=2342149;

(2342149 was the first entry in my list of invalid geometries) - Success - it returns 't' for true meaning the st_buffer trick fixes it.

Now to update the database.  To correct the error I had to do:

update naturalfeature_area set the_geom=st_multi(st_buffer(the_geom,0.0)) where gid=2342149;

Note that I do not really know what the st_multi thing does, but I got errors about geometry constraint violations without it - I think that st_buffer returns the simplest type of geometry it can, and st_multi forces it to be a multipolygon, but I could be wrong!.

Now all I need to do is remember enough about subqueries to write a bit of SQL that does the correction for every invalid geometry - back to that well known search engine....

Well, I have set it going using the following to try to repair the geometries:

update naturalfeature_area set the_geom=st_multi(st_buffer(the_geom,0.0)) where gid in (select gid from naturalfeature_area where not st_isvalid(the_geom));

I think I'll leave it for a few hours....

Using VectorMapDistrict with PostGIS

Well, vmd2pgsql took about 9 hours to import the entire UK into my postgresql database...No too bad.
Now I need to do something with it, because the mapnik stylesheets that I set up to use the shapefiles will not work - need to change it to use the postgresql database instead.
The things that need changing are:

1. Add the layers and datasources in XML (no need to do it in python now I don't have hundreds of separate datasources).
2. Convert the styles to use lowercase letters for the fields ('featcode' rather than 'FEATCODE') - postgresql seems to be case sensitive, and they have gone into the database lower case, but were upper case in the shapefiles.

Sounds easy,  but I decided I don't want it to look like the OSM stylesheet, which is a bit garbled - I want a nice structure to make it easy to maintain.  I could have used XML entities, but that seems a bit crude - define an entity, then 'use' it to actually 'call' it.   Includes seem much more suitable....But to get includes working I need mapnik2.  I have put a basic structure at code.google.com/p/ntmisc/vmdmap, but still have a few issues:

• I have had to define the database parameters manually because entities do not seem to be working in the include files, which is no good!
• The output is transparent background with black lines (see picture below), which is not what I asked for - not sure if this is a problem with the style file, or my build of mpanik2 - I'll have to run some tests on mapnik2 to make sure it works right.
• But it's too late to fix it now - job for tomorrow!

Using VectorMapDistrict with PostGIS

Well, vmd2pgsql took about 9 hours to import the entire UK into my postgresql database...No too bad.
Now I need to do something with it, because the mapnik stylesheets that I set up to use the shapefiles will not work - need to change it to use the postgresql database instead.
The things that need changing are:

1. Add the layers and datasources in XML (no need to do it in python now I don't have hundreds of separate datasources).
2. Convert the styles to use lowercase letters for the fields ('featcode' rather than 'FEATCODE') - postgresql seems to be case sensitive, and they have gone into the database lower case, but were upper case in the shapefiles.

Sounds easy,  but I decided I don't want it to look like the OSM stylesheet, which is a bit garbled - I want a nice structure to make it easy to maintain.  I could have used XML entities, but that seems a bit crude - define an entity, then 'use' it to actually 'call' it.   Includes seem much more suitable....But to get includes working I need mapnik2.  I have put a basic structure at code.google.com/p/ntmisc/vmdmap, but still have a few issues:

• I have had to define the database parameters manually because entities do not seem to be working in the include files, which is no good!
• The output is transparent background with black lines (see picture below), which is not what I asked for - not sure if this is a problem with the style file, or my build of mpanik2 - I'll have to run some tests on mapnik2 to make sure it works right.
• But it's too late to fix it now - job for tomorrow!

More Rendering of VectorMapDistrict Data

The nice people from Ordnance Survey have sent me a complete set of VectorMapDistrict data (which is on 6 DVDs!).   Surprisingly they didn't even charge me for the DVDs or postage.
Now I have all that data I thought I'd better do something with it!

First step was to get it onto a nice fast disk, so I copied it onto the hard disk of my server...which took a while...

Then I tried to use my vmdmap.py program to render it all, but mapnik bombed out with an error as it was adding the various shapefile layers.  I suspect that this is because every single shapefile is in its own layer, as it is a completely separate datasource, and I think I either ran out of memory or hit some mapnik internal limit.   This means I can't use vmdmap.py to render the whole country, which is a bit of a shame.

To get around this I think I need to merge them into a single datasource.  I don't know enough about manipulating shapefiles to do this, so instead am making use of shp2pgsql which allows you to import a shapefile into a postgresql database.   I have written another program, based on vmdmap.py called vmd2pgsql which will scan through a directory tree looking for the various shapefiles in the vectormap district dataset, and importing them into postgresql.   This gives a much lower number of tables - just one per shapefile name, but each one will have a lot of data.

It is importing now, so will see how long the import takes, then how well it renders.  I suppose it should render ok because I have the whole UK OSM dataset in a single database and that works, but we'll see over the weekend if it ever finishes!  The code is at my google code site.

More Rendering of VectorMapDistrict Data

The nice people from Ordnance Survey have sent me a complete set of VectorMapDistrict data (which is on 6 DVDs!).   Surprisingly they didn't even charge me for the DVDs or postage.
Now I have all that data I thought I'd better do something with it!

First step was to get it onto a nice fast disk, so I copied it onto the hard disk of my server...which took a while...

Then I tried to use my vmdmap.py program to render it all, but mapnik bombed out with an error as it was adding the various shapefile layers.  I suspect that this is because every single shapefile is in its own layer, as it is a completely separate datasource, and I think I either ran out of memory or hit some mapnik internal limit.   This means I can't use vmdmap.py to render the whole country, which is a bit of a shame.

To get around this I think I need to merge them into a single datasource.  I don't know enough about manipulating shapefiles to do this, so instead am making use of shp2pgsql which allows you to import a shapefile into a postgresql database.   I have written another program, based on vmdmap.py called vmd2pgsql which will scan through a directory tree looking for the various shapefiles in the vectormap district dataset, and importing them into postgresql.   This gives a much lower number of tables - just one per shapefile name, but each one will have a lot of data.

It is importing now, so will see how long the import takes, then how well it renders.  I suppose it should render ok because I have the whole UK OSM dataset in a single database and that works, but we'll see over the weekend if it ever finishes!  The code is at my google code site.

Mapnik and OS_OpenData

I decided to compare the quality of OpenStreetMap mapping to that which has been released by Ordnance Survey OpenData.  
I started with VectorMapDistrict, which is a dataset provided as ESRI Shapefiles.   It is provided as a large number of shapefiles, so setting up a mapnik stylesheet manually would have been a pain.  Instead I separated the styles into an XML file, but added the layers using a python program which scans down the directory tree to add the various shapefiles to the map.  The python code is vmdmap.py and the stylesheet is styles.xml.

Then I tried Meridian2.  This is much simpler - just one shapefile for each sort of feature (A-Road, river etc.).  I kept the same structure with the styles defined in an XML file, and the layers in python, but it could have all been done in a single XML file this time.  The python code is md2map.py and the stylesheet is styles_md2.xml.

You can see the results below.

The most notable things are that Meridian2 is much simpler geometry, but it includes more road names than Vector Map District.  The OpenStreetMap rendering is fancier because I used the standard OSM style, rather than the home made ones I used for the OS Data.   I do like the rocks that appear in the VectorMap District rendering though - I will have to import it into OSM!.

VectorMapDistrict

Meridian2

OpenStreetMap

Mapnik and OS_OpenData

I decided to compare the quality of OpenStreetMap mapping to that which has been released by Ordnance Survey OpenData.  
I started with VectorMapDistrict, which is a dataset provided as ESRI Shapefiles.   It is provided as a large number of shapefiles, so setting up a mapnik stylesheet manually would have been a pain.  Instead I separated the styles into an XML file, but added the layers using a python program which scans down the directory tree to add the various shapefiles to the map.  The python code is vmdmap.py and the stylesheet is styles.xml.

Then I tried Meridian2.  This is much simpler - just one shapefile for each sort of feature (A-Road, river etc.).  I kept the same structure with the styles defined in an XML file, and the layers in python, but it could have all been done in a single XML file this time.  The python code is md2map.py and the stylesheet is styles_md2.xml.

You can see the results below.

The most notable things are that Meridian2 is much simpler geometry, but it includes more road names than Vector Map District.  The OpenStreetMap rendering is fancier because I used the standard OSM style, rather than the home made ones I used for the OS Data.   I do like the rocks that appear in the VectorMap District rendering though - I will have to import it into OSM!.

VectorMapDistrict

Meridian2

OpenStreetMap

Reprojecting Ordnance Survey Raster Data

I have been looking at the Ordnance Survey StreetView data that was released to the public on 1st April.
This has got me thinking about a previous project that I failed to make work - I have some old maps of my home town and wanted to be able to display them on the web and switch between them so you can see what changes from year to year - I have maps from 1870 and 1914.  I want to use OpenStreetMap as the recent data.  Now StreetView is available I can use that too for comparison.

The problem is projecting them onto the same coordinate system - I failed miserably last time, but my two scanned old maps were on the same scale, so you can see them at (http://maps.webhop.net/oldmaps/openlayers.html).

The maps are in Ordnance Survey projection, which has the code EPSG:27700.   The downloaded OS data has the origin specified in metres northings and eastings.
You can translate this into the same projection as OpenStreetMap (EPSG:900113) using:

gdalwarp -s_srs EPSG:27700 -t_srs EPSG:900913 *.TIF nz.tif

Check the projection information with

gdalinfo nz.tif | more