Created
September 26, 2016 10:23
Star
You must be signed in to star a gist
Python script to convert a Ordnance Survey Code Point Open files to Redis 'GEOADD' commands
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
#!/usr/bin/python | |
# based on code from http://www.hannahfry.co.uk/blog/2012/02/01/converting-british-national-grid-to-latitude-and-longitude-ii | |
#This code converts lat lon (WGS84) to british national grid (OSBG36) | |
from scipy import * | |
from itertools import imap | |
import sys | |
import csv | |
def OSGB36toWGS84(E,N): | |
#E, N are the British national grid coordinates - eastings and northings | |
a, b = 6377563.396, 6356256.909 #The Airy 180 semi-major and semi-minor axes used for OSGB36 (m) | |
F0 = 0.9996012717 #scale factor on the central meridian | |
lat0 = 49*pi/180 #Latitude of true origin (radians) | |
lon0 = -2*pi/180 #Longtitude of true origin and central meridian (radians) | |
N0, E0 = -100000, 400000 #Northing & easting of true origin (m) | |
e2 = 1 - (b*b)/(a*a) #eccentricity squared | |
n = (a-b)/(a+b) | |
#Initialise the iterative variables | |
lat,M = lat0, 0 | |
while N-N0-M >= 1: #Accurate to 1m | |
lat = (N-N0-M)/(a*F0) + lat; | |
M1 = (1 + n + (5./4)*n**2 + (5./4)*n**3) * (lat-lat0) | |
M2 = (3*n + 3*n**2 + (21./8)*n**3) * sin(lat-lat0) * cos(lat+lat0) | |
M3 = ((15./8)*n**2 + (15./8)*n**3) * sin(2*(lat-lat0)) * cos(2*(lat+lat0)) | |
M4 = (35./24)*n**3 * sin(3*(lat-lat0)) * cos(3*(lat+lat0)) | |
#meridional arc | |
M = b * F0 * (M1 - M2 + M3 - M4) | |
#transverse radius of curvature | |
nu = a*F0/sqrt(1-e2*sin(lat)**2) | |
#meridional radius of curvature | |
rho = a*F0*(1-e2)*(1-e2*sin(lat)**2)**(-1.5) | |
eta2 = nu/rho-1 | |
secLat = 1./cos(lat) | |
VII = tan(lat)/(2*rho*nu) | |
VIII = tan(lat)/(24*rho*nu**3)*(5+3*tan(lat)**2+eta2-9*tan(lat)**2*eta2) | |
IX = tan(lat)/(720*rho*nu**5)*(61+90*tan(lat)**2+45*tan(lat)**4) | |
X = secLat/nu | |
XI = secLat/(6*nu**3)*(nu/rho+2*tan(lat)**2) | |
XII = secLat/(120*nu**5)*(5+28*tan(lat)**2+24*tan(lat)**4) | |
XIIA = secLat/(5040*nu**7)*(61+662*tan(lat)**2+1320*tan(lat)**4+720*tan(lat)**6) | |
dE = E-E0 | |
#These are on the wrong ellipsoid currently: Airy1830. (Denoted by _1) | |
lat_1 = lat - VII*dE**2 + VIII*dE**4 - IX*dE**6 | |
lon_1 = lon0 + X*dE - XI*dE**3 + XII*dE**5 - XIIA*dE**7 | |
#Want to convert to the GRS80 ellipsoid. | |
#First convert to cartesian from spherical polar coordinates | |
H = 0 #Third spherical coord. | |
x_1 = (nu/F0 + H)*cos(lat_1)*cos(lon_1) | |
y_1 = (nu/F0+ H)*cos(lat_1)*sin(lon_1) | |
z_1 = ((1-e2)*nu/F0 +H)*sin(lat_1) | |
#Perform Helmut transform (to go between Airy 1830 (_1) and GRS80 (_2)) | |
s = -20.4894*10**-6 #The scale factor -1 | |
tx, ty, tz = 446.448, -125.157, + 542.060 #The translations along x,y,z axes respectively | |
rxs,rys,rzs = 0.1502, 0.2470, 0.8421 #The rotations along x,y,z respectively, in seconds | |
rx, ry, rz = rxs*pi/(180*3600.), rys*pi/(180*3600.), rzs*pi/(180*3600.) #In radians | |
x_2 = tx + (1+s)*x_1 + (-rz)*y_1 + (ry)*z_1 | |
y_2 = ty + (rz)*x_1 + (1+s)*y_1 + (-rx)*z_1 | |
z_2 = tz + (-ry)*x_1 + (rx)*y_1 + (1+s)*z_1 | |
#Back to spherical polar coordinates from cartesian | |
#Need some of the characteristics of the new ellipsoid | |
a_2, b_2 =6378137.000, 6356752.3141 #The GSR80 semi-major and semi-minor axes used for WGS84(m) | |
e2_2 = 1- (b_2*b_2)/(a_2*a_2) #The eccentricity of the GRS80 ellipsoid | |
p = sqrt(x_2**2 + y_2**2) | |
#Lat is obtained by an iterative proceedure: | |
lat = arctan2(z_2,(p*(1-e2_2))) #Initial value | |
latold = 2*pi | |
while abs(lat - latold)>10**-6: | |
lat, latold = latold, lat | |
nu_2 = a_2/sqrt(1-e2_2*sin(latold)**2) | |
lat = arctan2(z_2+e2_2*nu_2*sin(latold), p) | |
#Lon and height are then pretty easy | |
lon = arctan2(y_2,x_2) | |
H = p/cos(lat) - nu_2 | |
#Convert to degrees | |
lat = lat*180/pi | |
lon = lon*180/pi | |
#Job's a good'n. | |
return lat, lon | |
#Read in from a file | |
LatLon = csv.reader(sys.stdin, delimiter = ',') | |
#Loop through the data | |
for line in LatLon: | |
pc = line[0] # first column is the postcode | |
E = line[2] # third column is the Easting | |
N = line[3] # fourth column is the Northing | |
lat, lon = OSGB36toWGS84(float(E), float(N)) | |
# output as a Redis GEOADD command | |
print "GEOADD postcodes " + str(lon) + " " + str(lat) + " " + pc.replace(" ", "") | |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment