Python Archive

0

Converting Latitude and Longitude to British National grid

This code reads in a .csv file called LatLon, expecting two columns with headers – Latitude and Longitude (in WGS84, decimal form). If the script is run in the same directory as LatLon.csv, it will spit out a second file called LatLonandBNG.csv, with two additional columns: OSGB36 Eastings and Northings respectively.

For the inverse transform - OSGB36 to WGS84 – please refer to this post, where a you can find the relevant python script, and more details on the algorithms involved.

If you don’t have python, you can find instructions for the bits you need on another of my previous posts.

#This code converts lat lon (WGS84) to british national grid (OSBG36)

from scipy import *

import csv

def WGS84toOSGB36(lat, lon):

    #First convert to radians
    #These are on the wrong ellipsoid currently: GRS80. (Denoted by _1)
    lat_1 = lat*pi/180
    lon_1 = lon*pi/180

    #Want to convert to the Airy 1830 ellipsoid, which has the following:
    a_1, b_1 =6378137.000, 6356752.3141 #The GSR80 semi-major and semi-minor axes used for WGS84(m)
    e2_1 = 1- (b_1*b_1)/(a_1*a_1)   #The eccentricity of the GRS80 ellipsoid
    nu_1 = a_1/sqrt(1-e2_1*sin(lat_1)**2)

    #First convert to cartesian from spherical polar coordinates
    H = 0 #Third spherical coord.
    x_1 = (nu_1 + H)*cos(lat_1)*cos(lon_1)
    y_1 = (nu_1+ H)*cos(lat_1)*sin(lon_1)
    z_1 = ((1-e2_1)*nu_1 +H)*sin(lat_1)

    #Perform Helmut transform (to go between GRS80 (_1) and Airy 1830 (_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, b = 6377563.396, 6356256.909 #The GSR80 semi-major and semi-minor axes used for WGS84(m)
    e2 = 1- (b*b)/(a*a)   #The eccentricity of the Airy 1830 ellipsoid
    p = sqrt(x_2**2 + y_2**2)

    #Lat is obtained by an iterative proceedure:
    lat = arctan2(z_2,(p*(1-e2))) #Initial value
    latold = 2*pi
    while abs(lat - latold)>10**-16:
        lat, latold = latold, lat
        nu = a/sqrt(1-e2*sin(latold)**2)
        lat = arctan2(z_2+e2*nu*sin(latold), p)

    #Lon and height are then pretty easy
    lon = arctan2(y_2,x_2)
    H = p/cos(lat) - nu

    #E, N are the British national grid coordinates - eastings and northings
    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)
    n = (a-b)/(a+b)

    #meridional radius of curvature
    rho = a*F0*(1-e2)*(1-e2*sin(lat)**2)**(-1.5)
    eta2 = nu*F0/rho-1

    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)          

    I = M + N0
    II = nu*F0*sin(lat)*cos(lat)/2
    III = nu*F0*sin(lat)*cos(lat)**3*(5- tan(lat)**2 + 9*eta2)/24
    IIIA = nu*F0*sin(lat)*cos(lat)**5*(61- 58*tan(lat)**2 + tan(lat)**4)/720
    IV = nu*F0*cos(lat)
    V = nu*F0*cos(lat)**3*(nu/rho - tan(lat)**2)/6
    VI = nu*F0*cos(lat)**5*(5 - 18* tan(lat)**2 + tan(lat)**4 + 14*eta2 - 58*eta2*tan(lat)**2)/120

    N = I + II*(lon-lon0)**2 + III*(lon- lon0)**4 + IIIA*(lon-lon0)**6
    E = E0 + IV*(lon-lon0) + V*(lon- lon0)**3 + VI*(lon- lon0)**5 

    #Job's a good'n.
    return E,N

#Read in from a file
LatLon = csv.reader(open('LatLon.csv', 'rU'), delimiter = ',')
LatLon.next()

#Get the output file ready
outputFile = open('LatLonandBNG.csv', 'wb')
output=csv.writer(outputFile,delimiter=',')
output.writerow(['Lat', 'Lon', 'E', 'N'])

#Loop through the data
for line in LatLon:
    lat = line[0]
    lon = line[1]
    E, N = WGS84toOSGB36(float(lat), float(lon))
    output.writerow([str(lat), str(lon), str(E), str(N)])

#Close the output file
outputFile.close()
Tags:
0

Converting British National Grid to Latitude and Longitude II

A few months ago, I wrote a python script to convert British National grid coordinates (OSGB36) to latitude and longitude (WGS84).  A fellow blogger Andrzej Bieniek very kindly pointed out that the algorithm was only accurate to around 100m because of an additional subtlety which I hadn’t taken into account.

I’ll have a quick bash at explaining the reasoning behind this difference,  but if you want to skip to the new version of the code (now accurate to 5m) I’ve posted it at the bottom of the page in all its delicious glory. Alternatively, if you are looking for a similarly accurate script applying the inverse transform – ie. WGS84 lat, lon to OSGB36 Eastings, Northings – you can find it here.

Because the Earth’s surface is not smooth, any attempt to represent a geographical area in two dimensions  - on paper, on screen etc – will naturally result in some distortions and discontinuities. One way to get round this, is to approximate the Earth using a ’3D ellipse’ (or ‘ellipsoid’, if we’re being proper) try and fit it as best you can to the Earth’s squashed sphere shape, cross your fingers & not get too worried about ‘land’, ‘oceans’ and other such trivial inconveniences.

If you’re looking to consider large distances, or perhaps the whole globe at once – as is the case with WGS84 –  you’ll need to fit the ellipsoid to the entire Earth’s surface at once.  If, however, you’re only interested in mapping a small area of the globe –  as is the case with British National Grid  - it may be that a much smaller ellipse would work much better:

Of course, the picture above – taken from the Ordinance Survey website – is a greatly exaggerated illustration of this idea, but I think it demonstrates the concept pretty well.  This is exactly what I hadn’t realised when doing my previous post: the ellipsoids that are used for the two systems of British National grid and WGS84 latitude and longitude are ever so slightly different – hence why the results could be up to 100m out. British National grid is based on the Airy* 1830 ellipsoid, while latitude and longitude use the GSR80 (or WGS 1984).

The origins of these two, A1830 and GSR80 are in slightly different positions, their axes are a little rotated, and one is a tiny bit larger than the other – as in the picture below.

It’s not sufficient to just find  ‘lat/lon style’ coordinates from Eastings and Northings. You’ve also got to take the differences between the two ellipsoids before you can reproduce the latitude and longitude coordinates.

So, on top of the original post, to convert between the two ellipsoids you must scale, translate and rotate every point on the surface.

Actually, once you get down to it, the maths isn’t that bad – although the transform is a great deal simpler if you first convert to Cartesian coordinates, and that’s the only annoying bit. The full process is actually all contained on the Ordinance Survey guide  (just annoyingly not all in one place) and can be summarised as follows:

1) Use your Eastings and Northings to find lat/lon style coordinates on the Airy 1830 ellipsoid. (As in the previous post, or on pages 41&42 of the guide).

2) Convert those lat/lon coords to 3D Cartesian coordinates. (As on page 38 of the guide). Take the third spherical polar coordinate – height –   as zero in your calculations. (As though the point sits on the surface of the ellipsoid).

3) Use Helmert’s datum transform to convert between the two ellipsoids, and anchor your coordinates on the GSR80 system. (As on page 29 of the guide).

4) Your answer will be in Cartesian coordinates, so the final step is to convert back to spherical polars to get the final lat/lon. (As on page 39 of the guide).

All the constants you will need are also there in the guide – although you might have to hunt for them a little. Or just nick ‘em from my script below.

As before, the code reads in a .csv file called BNG, expecting two columns with headers – East and North. If the script is run in the same directory as BNG.csv, it will spit out a second file called BNGandLatLon.csv, with both sets of coordinate (lat and lon will be in decimal form). If you don’t have python, you can find instructions for the bits you need on another of my previous posts. Aren’t I nice to you?

*The same chap as Airy’s function, which I virtually have tattooed on the insides of my eyelids from my fluid dynamics days. Not that this is remotely interesting, I just think he’s a top fella is all.

#This code converts british national grid to lat lon

from math import sqrt, pi, sin, cos, tan, atan2 as arctan2
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 >= 0.00001: #Accurate to 0.01mm
        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**-16: 
        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

    #Uncomment this line if you want to print the results
    #print [(lat-lat_1)*180/pi, (lon - lon_1)*180/pi]

    #Convert to degrees
    lat = lat*180/pi
    lon = lon*180/pi

    #Job's a good'n. 
    return lat, lon

#Read in from a file
BNG = csv.reader(open('BNG.csv', 'rU'), delimiter = ',')
BNG.next()

#Get the output file ready
outputFile = open('BNGandLatLon.csv', 'wb')
output=csv.writer(outputFile,delimiter=',')
output.writerow(['Lat', 'Lon', 'E', 'N'])

#Loop through the data
for E,N in BNG:
    lat, lon = OSGB36toWGS84(float(E), float(N))
    output.writerow([str(lat), str(lon), str(E), str(N)])
#Close the output file
outputFile.close()
Tags:
0

Bezier Curves

I wonder how long I can carry on working at the Centre for Advanced Spatial Analysis and still manage to avoid learning any GIS.

This week: Bezier curves, how to draw them in Python and particularly, how to decide where to put the control points.

A Bezier curve, is a special type of parametric curve used frequently in computer graphics. You may have seen them in Powerpoint, or using the pen tool in illustrator, but for my purpose, they are a lovely way to visualise flow paths on a map, like this:



This particular map shows the home addresses and riot locations of a subset of suspects involved in the London riots. The data was kindly given to us by Alex Singleton.

The curves themselves are anchored at both ends (origin P0, destination P3) and get their shape from two control points (blue P1 and red P2), as in the following picture:

The curve itself is created from the fairly simple expression (everything in bold is a vector):

 \textbf{B}(t)=(1-t)^3\textbf{P}_0+3t(1-t)^2\textbf{P}_1+3t^2(1-t)\textbf{P}_2+t^3\textbf{P}_3,\;\;\;t\in(0,1)

The parameterisation t starts from zero, so that the curve begins at P0, and as it increases bends the line towards the two control points, first P1 and then P2. Finally, when t  is 1, the line finishes at P3.

It’s all very clever stuff, and there is lots and lots of literature available on the web about them. But – as far as I can see – no one ever tells you how to decide where to put your control points, and certainly doesn’t offer an algorithm to automatically do so.

But fear not people! I have written one. And ugly and fudgy as it is, it does the job.

The code below takes 6 inputs:

-origin points,

-destination points,

-a length for the blue line,

-a length for the red line (both as a fraction of the length of the green line)

-the angle between the blue and green lines

-the angle between the red and green lines (both in radians. I am a mathematician after all).

And it will spit out a totally delicious plot of your curves.

It’s set up to bend the curves clockwise (as in the image at the top of the page), which generally looks nicer when plotting more than one bezier curve.

Have a play with the last four inputs: the longer the length, the more the curve will deviate from the green line. The higher the angle, the steeper the curve will roll in to it’s destination. Indeed, choosing the angles to be on the other side of the green line will result in an s-shaped curve. Anyway – the idea is for you to hack out the bits of the code that are useful.

import scipy as sp
import pylab as plt

#### Inputs

#A list of P0's and P3's. Must be the same length
origins = [[1,0],[-0.5,sp.sqrt(3)/2], [-0.5,-sp.sqrt(3)/2]]
destinations = [[0,0],[0,0],[0,0]]

#The angle the control point will make with the green line
blue_angle = sp.pi/6
red_angle = sp.pi/4

#And the lengths of the lines (as a fraction of the length of the green one)
blue_len = 1./5
red_len = 1./3

### Workings

#Generate the figure
fig = plt.figure()
ax = fig.add_subplot(111)
ax.hold(True)

#Setup the parameterisation
t = sp.linspace(0,1,100)

for i in xrange(len(origins)):
    #Read in the origin & destination points
    POx,POy = origins[i][0], origins[i][1]
    P3x,P3y = destinations[i][0], destinations[i][1]

    #Add those to the axes
    ax.plot(POx,POy, 'ob')
    ax.plot(P3x,P3y, 'or')
    ax.plot((POx,P3x),(POy,P3y), 'g')

    #Work out r and theta (as if based at P3)
    r = ((POx-P3x)**2 + (POy-P3y)**2)**0.5
    theta = sp.arctan2((POy-P3y),(POx-P3x))

    #Find the relevant angles for the control points
    aO =  theta + blue_angle+ sp.pi
    aD = theta - red_angle

    #Work out the control points
    P1x, P1y = POx+ blue_len*r*sp.cos(aO), POy + blue_len*r*sp.sin(aO)
    P2x, P2y = P3x+ red_len*r*sp.cos(aD), P3y + red_len*r*sp.sin(aD)

    #Plot the control points and their vectors
    ax.plot((P3x,P2x),(P3y,P2y), 'r')
    ax.plot((POx,P1x),(POy,P1y), 'b')
    ax.plot(P1x, P1y, 'ob')
    ax.plot(P2x, P2y, 'or')

    #Use the Bezier formula
    Bx = (1-t)**3*POx + 3*(1-t)**2*t*P1x + 3*(1-t)*t**2*P2x + t**3*P3x
    By = (1-t)**3*POy + 3*(1-t)**2*t*P1y + 3*(1-t)*t**2*P2y + t**3*P3y

    #Plot the Bezier curve
    ax.plot(Bx, By, 'k')

#Save it
plt.savefig('totally-awesome-bezier.png')
plt.show()
#Bosch.
Tags:
0

Converting British National Grid to Latitude and Longitude I

EDIT: Andrzej Bieniek brought to my attention that this version is correct to 100m. For a more accurate script (accurate to 5m) see my new post.

I have recently started to deal with a lot of geographical data in my research and and begun to realise it is difficult to go to long without stumbling across the sticky world of map projections – something I knew almost nothing about a year ago.

There is a nice blog post explaining the background by James Cheshire which I shan’t attempt to reproduce, but explanations aside, I found myself today trying to convert a long list of British National grid coordinates into Latitude and Longitude.

Of course, if you’ve got a week spare, you can do these conversions using QGis or ArcGis, but all I was looking for was just a quick .csv to .csv converter  of one coordinate system to the other. After a fruitless half hour of googling, I decided it would just be quicker to write a python script myself to get the job done. And because I know I’ll use it again, here it is.

The code reads in a .csv file called BNG, expecting two columns with headers – East and North. If the script is run in the same directory as BNG.csv, it will spit out a second file called BNGandLatLon.csv, with both sets of coordinate. All the maths behind the conversion can be found on page 41 of a pdf written by the chaps at ordnance survey. Sorry about the lack of highlighting – I’ll sort it out when I have a spare moment.

#This code converts british national grid to lat lon

from scipy import *
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 >= 0.00001: #Accurate to 0.01mm
        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

    lat = lat - VII*dE**2 + VIII*dE**4 - IX*dE**6
    lon = lon0 + X*dE - XI*dE**3 + XII*dE**5 - XIIA*dE**7

    #Convert to degrees
    lat = lat*180/pi
    lon = lon*180/pi

    return lat, lon

#Read in from a file
BNG = csv.reader(open('BNG.csv', 'rU'), delimiter = ',')
BNG.next()

#Get the output file ready
outputFile = open('BNGandLatLon.csv', 'wb')
output=csv.writer(outputFile,delimiter=',')
output.writerow(['Lat', 'Lon', 'E', 'N'])

#Loop through the data
for E,N in BNG:
    lat, lon = OSGB36toWGS84(float(E), float(N))
    output.writerow([str(lat), str(lon), str(E), str(N)])
#Close the output file
outputFile.close()
Tags: