Sometimes we have to deal with data in a coordinate system that is incompatible with other processes. R (along with the rgdal library) is a great tool to convert between pretty much any coordinate system.

Below is an example that converts the Code-Point database from OSGB36 to WGS84 (so that it can be plotted in Google Maps).

The idea is to find the EPSG  codes for your ‘from’ & ‘to’ coordinate system. The code-point database is in OSGB 1936 which has an EPSG code: 27700. For Google Maps we want EPSG 4326.

Note: Google Maps uses the EPSG: 3857 for Map Tiles, however the coordinates are in latitude/longitude and EPSG: 4326. For a further discussion it may be useful to read this post, where it is mentioned:

All of this further confused by that fact that often even though the map is in Web Mercator(EPSG: 3857), the actual coordinates used are in lat-long (EPSG: 4326).

The x, y coordinates of the data are defined to create a spatial object, which is then ‘projected’ in its original coordinate system. Then the spatial object is transformed to the desired coordinate system.

For a more accurate transformation this post may be useful. Where it is mentioned that:

PROJ.4 contains a large number of transformations and parameters. However by default OGR uses a 3 parameter shift because it covers the largest available area. However, this is only accurate to +/- 20 metres over the whole of the UK. It is equivalent to the ArcGIS transformation called OSGB_1936_To_GS_1984_1 and has parameters of dX: 375, dY-111, dZ 431″

However, since the data deals with postcodes which are polygons, and are approximated with points I don’t think it matters too much if the centroid of the polygon shifts 10 metres right, given that we have around 10-15 households per postcode.

R Script:

########################################################
# OS Code OpenPoint,  Version: 08/2015
# Merge data and convert coordinate system to WGS84
########################################################

require(maptools)
require(rgdal)

setwd('K:/compdata/OS_Open_Point/Data/')
out = 'K:/compdata/OS_Open_Point/'

# Import Raw
files = list.files(pattern="*.csv")
myfiles = do.call("rbind", lapply(files, function(x) read.csv(x, header=FALSE, stringsAsFactors = FALSE)))
names(myfiles) <- c(
  "postcode",
  "positional_quality_indic",
  "eastings",
  "northings",
  "country_code",
  "NHS_regional_HA_code",
  "NHS_HA_code",
  "Admin_county_code",
  "Admin_district_code",
  "Admin_word_code"
  
)

# Remove space in postcode
myfiles$postcode_no_space <- gsub(' ','',myfiles$postcode)

# Process only location data
uk_bng <- myfiles[,c('postcode','postcode_no_space','eastings','northings')]

# Key proj4 transformations for exercise
# USE: http://spatialreference.org/ref/epsg/osgb-1936-british-national-grid/
# FROM (BNG) TO (WGS84)
bng <- "+init=epsg:27700"
wgs84 <- "+init=epsg:4326"
  
# Create spatial object where eastings is x-axis and northings is y-axis
coordinates(uk_bng) <- c("eastings", "northings")
# Set projection of data
proj4string(uk_bng) <- CRS(bng)

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

# Convert SpatialPoints back into a data-frame
converted_coordinates <- as.data.frame(uk_wgs84)
names(converted_coordinates) <- c("postcode", "postcode_no_space", "longitude", "latitude")

########################################################
# Export
########################################################
# CSV
write.table(converted_coordinates, paste(out,'postcode_osgb_082015_v2.csv',sep=''), sep=',', row.names=FALSE)
# Stata
require(foreign)
write.dta(converted_coordinates, paste(out,'postcode_osgb_082015_v2.dta',sep=''))