[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
]
Do you have or know of a python script to do conversion the other way, OS Grid -> eastings/northings
Tried to reverse engineer the linked javascript but failing miserably
Not on me, I’m afraid. However adapting parts of phpcoord/JScoord for python may be fruitful – see http://www.jstott.me.uk/phpcoord/
This may be too late for Tony but here’s a python script to convert the other way.
#!/bin/env python
'''
Convert OS refs with letters to cartesian coords in meters.
'''
import math
def os_cart(inref):
''' Return a tuple of two 7 digit OS numeric grid refs from an XXNNNN type reference.
Also returns the accuracy of the original OS ref since this is lost when
the data are padded to two 6 digit grid refs.
Example: os_cart("SN109112") returns (210900, 211200, 100)
'''
#print inref
'''get numeric values of letter references, mapping A->0, B->1, C->2, etc:'''
inref = inref.replace(" ","") #Strip all spaces out
# Deal with the letters
l1 = ord(inref[0].upper())-ord("A")
l2 = ord(inref[1].upper())-ord("A")
# shuffle down letters after 'I' since 'I' is not used in grid:
if l1 > 7: l1 -= 1
if l2 > 7: l2 -= 1
e = str(((l1-2)%5)*5 + (l2%5)) #easting
n = str(int((19-math.floor(l1/5)*5) - math.floor(l2/5))) #northing
# Now the numbers
gridref = inref[2:]
e += gridref[:int(len(gridref)/2)]
n += gridref[int(len(gridref)/2):]
# Pad short refs with correct number of zeros for postGIS
# This does imply a greater accuracy then the original data suggest
# however and should be used with caution, see next line.
a = 1 * int(math.pow(10,6-len(n))) # Calculate the accuracy of original coordinates in metres.
n = int(math.pow(10,6-len(n))) * int(n)
e = int(math.pow(10,6-len(e))) * int(e)
return int(e), int(n), int(a)
if __name__ == "__main__":
''' Run some test code id called directly from command line'''
OSREF = "SN109112"
point = os_cart(OSREF)
print OSREF, point