PHP and Ordnance Survey Mapping

About a month ago, Ordnance Survey opened up some of their data for public consumption under a brand called OpenData. The data is licenced under a "Creative Commons Attribute"-like license. One of the data sets they provide is Code-Point Open and provides a dataset "that contains postcode units, each of which have a precise geographical location". I've always been a bit of a map-geek, and have always geotagged my pictures and some of my tweets. So I was wondering what cool thing I could do with this newly released data. I decided to map all the postcodes onto the UK map where more postcodes for a specific place would create a "lighter" colour. Each postcode has on average about 15 addresses, so in more densely populated areas you have more "postcodes-per-area". Doing this wasn't very difficult and it resulted in the following map:

You can very clearly see the more densely populated areas such as London.

This generated postcode-density map I wanted to overlay on a real map, such as OpenStreetMap provides. When I tried to align the image that I generated with the OSM map, I ended up with this:

Obviously, the maps don't align. In order to fix this, I ended up learning a lot more about map projections.

The Code-Point Open data uses the UK's National Grid to store location data in. The National Grid consists of 100*100 km squares that are then further subdivided into smaller squares creating grid references such as TQ3012780512, or the numerical version 530128 180512. TQ translates into the "hundred thousands" of the Eastings and Northings according to the grid that you can see here. In this case, it specifies a point 530 128 meters East and 180 512 meters North of the origin. (If you work it out, you'll end up in London).

OpenStreetMap uses a Mercator projection to visualize maps. The Mercator projection is a cylindrical map projection, and it distorts the size and shape of large objects, as the scale increases from the Equator to the poles. Therefore it only works from about 85°N to 85°S (why it is 85° only be came clear after doing all the maths for it). Google Maps uses the same projection.

There were a few challenges plotting the Code-Point Open data on an OpenStreetMap map:

• The National Grid coordinates need to be converted to Latitude/Longitude pairs

• Latitude/Longitude pairs need to be mapped to pixels to align with the Mercator projection of OpenStreetMap maps.

For converting National Grid references to Latitude/Longitude pairs I found some JavaScript code. I converted this code to PHP and to verify whether my code was working properly, I used this site.

Converting Latitude/Longitude pairs to the pixel coordinates of OSM required a bit of maths. OSM provides maps in different zoom levels, from 2 to 17. Each zoom level having twice the amount of pixels horizontally and vertically. Zoom level 2 has 4 times 4 tiles, where each tile is 256*256 pixels. At zoom level two, the whole world fits into 1024*1024 pixels. The number of pixels for each axis for a specific zoom-level can be calculated by:

pow(2^zoom) * 256

For a zoom level of 7 that makes 16384 pixels in each direction. To convert a longitude (in the range -180° to 180°) we can simply apply:

x = ((lon + 180) / 360) * (pow(2^zoom) * 256)

For a Longitude of 0.003117° E at zoom-level 13 that turns out to be pixel 1048594.

Longitude conversions are more difficult due to the Mercator projection itself. To convert we apply:

y = ((atanh(sin(deg2rad(-lat))) / π) + 1) * (pow(2^(zoom-1)))

For a Latitude of 51.502817° N at zoom-level 13 that turns out to be pixel 697399.

Creating an image for the whole world at zoom level 13 is impractical as it is 2097152*2097152 pixels and downloading all the 67 million tiles for this zoom level is probably not liked by the OSM people either. So instead we take a cut out for a specific area only. For the UK (61.37°N -9.49° W, 49.76°N +3.63°E) at zoom level 7 we end up with a 1536x2048 map with the North-Western pixel being 15360,9216. On this map, we can draw the latitudes and longitudes as well as the National Grid lines (full image is on flickr):

The only thing left to do now, is to map the postcode density information to the map. I picked zoom level 6 for this, and the result is (after cropping it to 640 pixels width):

As you can see, this is perfect fit to the outlines of the country. But unfortunately, when we look very closely at the plotted map data, for example for the NW10 3 postcodes, we notice that the mapping is slightly off. The blue dots are what we plotted, and the red dots are what the locations should have been:

The reason for this is that when we converted the National Grid locations to Latitude/Longitude pairs to plot on the OSM maps, I forgot to take into account the different Datums that are used in the projections. The Earth is not a perfect sphere, and an approximation of the ellipsoid of the whole Earth is not necessarily the best fitting for a specific area such as the UK. Therefore, the National Grid uses the OSGB36 Datum which fits more closely to the UK, where as OpenStreetMap uses the WGS84 Datum that is also used by GPS. The Ordnance Survey Ireland has a more thorough explanation on their site. As you can see above, using the wrong Datum can mean locations can be off. In our example about 100 meters. Converting between different Datums is possible, albeit processor intensive.

After I figured out all the maths for this, the only problem that remains that implementing those algorithms in PHP is show—calculating all the positions from the 1.6 million postcode locations takes up to 10 minutes. This is why I am not presenting any code yet. I am planning to implement all the necessary calculations in a PHP extension to speed up the calculations