British National Grid transformation and re-projection in R

Georeferenced UK statistical data typically comes in one of two spatial coordinate reference systems – either latitude/longitude points based on the World Geodetic System 1984 (‘WGS 84′), or in the eastings/northings system of the British National Grid (based on the 1936 Ordnance Survey Great Britain datum ‘OSGB 36′). The latter system has the benefit of being measured in meters, which enables more accurate spatial analysis to be undertaken. Converting data between these coordinate systems is a tricky mathematically process, but those facing this for the first time need not worry because clever minds have thought about this since the days of Ancient Greece.

This R script is a minimal example of how to take data in British National Grid (eastings/northings) or WGS84 (latitude/longitude) and convert to the other. The heavy lifting is undertaken by rgdal library, R’s implementation of the cross-platform Geospatial Data Abstraction Library. The transformations themselves are described by the proj.4 strings included.

The script also utilises the useful OpenStreetMap package to download map tiles to add context to data. Online map tools mostly base these on the Mercator projection, which many geographers dislike for its areal distortion towards the poles and the global social inequality this engenders. Mercator does however have the benefit of preserving angularity, meaning the shapes of localities are much less distorted. While the map tiles can also be re-projected to, for example British National Grid, (as demonstrate in the script) this reduces image quality. For this reason I’ve also shown how to reproject data into Mercator coordinates so it can be directly plotted onto these map tiles.

Anyone keen to customise map tiles to suit their data should check out cloudmade.

require(maps)
require(maptools)
require(rgdal)
require(OpenStreetMap)

data(world.cities)
uk = world.cities[world.cities$country.etc == 'UK' & world.cities$pop > 100000,]
head(uk)

# coerce to a spatial object
coordinates(uk) = c('long','lat')
plot(uk)

# key proj4 transformations for exercise
wgs84 = '+proj=longlat +datum=WGS84'
bng = '+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'
mrc = '+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs'

# useful table of proj4 transformations
epsg = make_EPSG()
View(epsg[grep("OSGB", epsg$note),])

# specify coordinate system as WGS84 (typical latlons)
uk@proj4string   # slot will be empty
uk@proj4string = CRS(wgs84)

# reproject data to Mercator coordinate system
uk_merc = spTransform(uk, CRS(mrc))

# download Mercator tiles covering the data area
bbox = uk@bbox
map = openmap(c(bbox[2,2],bbox[1,1]), c(bbox[2,1],bbox[1,2]), 7, 'osm')

plot(map)
plot(uk_merc, add=T)

# conversion to British National Grid (OSGB36)
map_bng = projectRaster(raster(openproj(map)), crs = bng)
uk_bng = spTransform(uk, CRS(bng))

plotRGB(map_bng)
plot(uk_bng, add=T)

#  BNG to WGS84 conversion
uk_wgs84 = spTransform(uk_bng, CRS(wgs84))
coordinates(uk_wgs84)

# custom cloudemade map tiles
plot(openmap(c(bbox[2,2],bbox[1,1]), c(bbox[2,1],bbox[1,2]), 6, 'cloudmade-115496'))