Monthly Archives: February 2010

Prettier Maps with Mapnik 0.7.0

Posted on by

[Updated] The Mapnik 0.7.0 release last month has added a number of features, including one that makes choropleth maps, and others with numerous adjacent polygons, more pretty. The PolygonSymbolizer has a new parameter, gamma, which, when set to a value between 0 and 1, causes the polygon edges to bleed into each other slightly, removing distinctive hairlines that allow the background colour to “shine through” the joins.

Here’s a before-and-after of a choropleth of England & Wales. We don’t care much for the boundaries, as they are Middle-Level Super Output Areas (MSOAs), which are administrative boundaries that people are not familiar with, unlike say county boundaries. So, using the gamma parameter, we can hide them.


Mapnik 0.6


Mapnik 0.7

The new parameter is also useful for my OpenOrienteeringMap project, particularly to make large estuaries, which are typically made up of multiple polygons, appear contiguous. However, applying too low a value appears to affect the anti-aliasing of line features, even ones away from the affected polygons. Here, the gamma is just applied to the “water” style. Note the thin white line crossing the water, near the black dotted line that also crosses it, disappears, but the paths and roads start to distort too:


Mapnik 0.6 – thin white line running across the river.


Mapnik 0.7+patch (Gamma 0.8) – the thin white line has gone. Other features are unchanged.

The map tile, by the way, is showing the Greenwich Foot Tunnel crossing the River Thames in east London.

[Update - The crossed-out section highlights a bug that has now been fixed by the Mapnik development team. I've patched my Mapnik build and the gamma changes now only appear where expected.]

OpenLayers 101

Posted on by

Here’s a simple one-file way to get an OpenLayers map on a website. This is all the code you need, including the HTML, Javascript and CSS.

This particular example pulls in some custom tiles from a URL and overlays them on the “standard” OpenstreetMap map. You will need to change the URL (and attribution text) to point to the overlay you are interested in.

<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Strict//EN"
    "http://www.w3.org/TR/xhtml1/DTD/xhtml1-strict.dtd">
<html xmlns="http://www.w3.org/1999/xhtml"
    xml:lang="en" lang="en">
<head>
<meta http-equiv="Content-Type"
    content="text/html;charset=utf-8" />
<title>NPE Scotland</title>
<script type="text/javascript"
    src="http://www.openlayers.org/api/OpenLayers.js">
</script>
<script type="text/javascript">
var map;
var EPSG4326 = new OpenLayers.Projection("EPSG:4326");
var EPSG900913 = new OpenLayers.Projection("EPSG:900913");
var b = 20037508.34;

function init()
{
    map = new OpenLayers.Map ("map",
    {
        controls:[
            new OpenLayers.Control.Navigation(),
            new OpenLayers.Control.PanZoomBar(),
            new OpenLayers.Control.Attribution(),
        new OpenLayers.Control.LayerSwitcher(),
        new OpenLayers.Control.MouseDefaults(),
        new OpenLayers.Control.KeyboardDefaults()],
        maxExtent: new OpenLayers.Bounds(-1*b, -1*b, b, b),
        maxResolution: 156543.0399,
        units: 'm',
        projection: EPSG900913,
        displayProjection: EPSG4326
    });

    layerMapnik = new OpenLayers.Layer.OSM(
    	null, null, { numZoomLevels: 15 });
    layerOverlay = new OpenLayers.Layer.OSM(
        "NPE Scotland",
        "tiles/${z}/${x}/${y}.jpg",
        {
            transitionEffect: "none",
            isBaseLayer: false,
            attribution: "Imagery CC-By-NC-SA NPEMap.org.uk"
        }
    );
    map.addLayer(layerMapnik);
    map.addLayer(layerOverlay);
    var start = new OpenLayers.LonLat(-3.5, 56.5);
    map.setCenter(start.transform(EPSG4326, EPSG900913), 7);
}
</script>
</head>
<body onload="init();" style="margin: 0;">
<div id="map"
    style="position: absolute; width: 100%; height: 100%;">
</div>
</body>
</html>
 

Here’s what it looks like:

Converting OS Eastings/Northings to Grid References in Python

Posted on by

[Updated] I needed to programatically convert a series of Ordnance Survey easting and northings (e.g. 325940, 673060) to “six-figure” grid references (e.g. NT259731) for a project I’m currently working on. It’s a pretty straightforward conversion – no reprojection, just a different way of expressing the same position on the British National Grid.

The specific use is that the 1km x 1km tiles from NPEMap are coding according to eastings & northings, but the calibration files required by TimSC‘s warp-gbos require grid references for the corners and calibration points.

Such a procedure is already well catered for in Perl, PHP and Javascript. Here is a straightforward conversion of some Javascript I found here into Python:

# Derived from
# http://www.movable-type.co.uk/scripts/latlong-gridref.html

def getOSGridReference(e, n):
	import math

	# Note no I
	gridChars = "ABCDEFGHJKLMNOPQRSTUVWXYZ"  	

	# get the 100km-grid indices
	e100k = math.floor(e/100000)
	n100k = math.floor(n/100000)

	if e100k6 or n100k12:
		return ''

	# translate those into numeric equivalents
        # of the grid letters
	l1 = (19-n100k)-(19-n100k)%5+math.floor((e100k+10)/5)
	l2 = (19-n100k)*5%25 + e100k%5

	letPair = gridChars[int(l1)] + gridChars[int(l2)]

	# strip 100km-grid indices from easting & northing,
        # round to 100m
	e100m = math.trunc(round(float(e)/100))
	egr = str(e100m).rjust(4, "0")[1:]
	if n >= 1000000:
		n = n - 1000000 # Fix Shetland northings
	n100m = math.trunc(round(float(n)/100))
	ngr = str(n100m).rjust(4, "0")[1:]

	return letPair + egr + ngr

# test
print getOSGridReference(96000,906000) # NA960060
print getOSGridReference(465149, 1214051) # HP651141
 

[Update - fixed a bug that gave the wrong grid references for the Shetland Islands - their northings have seven figures, and similarly places in the far west or south that only have five figure easting or northings respectively :oops: ]