Google Mapping, Markers, OS Grid Refs, Latitude Longitude and VBscript

Google Mapping, Markers, OS Grid Refs, Latitude Longitude and ASP VBscript

How to get from SM856127 to 51.7713932226388, -5.10791492926723

A brief history of time-wasting :

I'm employed by the Pembrokeshire Coast National Park (and actually, despite the pay I quite like it), a project came my way to put our National Trail walks on the web. This lead me to Google Maps, as I had some prior knowledge as well as the fact that a lot of info could be displayed from our databases. What do walkers need after a hard day of walking ? - a kip and what do you know, we have a database full of accommodation with OS Grid Refs !

I needed to convert these OS grid Refs to Latitude & Longitude, so that Google could use them. In the past for converting the odd OS ref i'd used a JavaScript code from - http://www.carabus.co.uk/ngr_ll.html. Worth a look if your starting out down a similar path. However I needed to convert over 100 on the fly, in order to plot them all in one go. This method is out because mixing Client Side and Server Side scripting is a pain the ar...

I set out and converted the JavaScript functions into ASP VBScript...
Firstly convert the first two characters. for example SM. Have a goosie here
Secondly, all the maths bit. Have a goosie here

Instructions for use :

1) Save the below in an include file e.g. NEtoLL_inc.asp or download zip file.
2) In your main asp file include it 
3) Then call by passing the text or variable with GridRef e.g. GridRefConvert("SM856127")
Note: remember to validate before passing e.g. GridRefConvert( replace(MyGridRef," ","") )

DOWNLOAD ZIPPED INCLUDE FILE HERE

Shortcomings :

1) It only does SN, SM, SS & SR squares (use http://www.carabus.co.uk/jscalculators.html to find your required numbers)
2) 6 digit OS GridRefs aren't all that accurate - easy enough to update the code to change to 10 though (I may do that one day and re-post here) but as all my data is 6 figure....

Addendum / Disclaimer :

I'm just starting out down this road myself, so if you have any info please have a look at the Google Groups discussion I started on this subject. If this code causes any negative effect on your system(s), on your head be it.

Email me (Pete) skibutt@yahoo.co.uk if you have questions or to point out glaring mistakes i've made.

Some handy links :
The source (Where I nabbed the code from)
LatLongFinder (handy for plotting polylines etc)

The Code

function NEtoLL(east, north)

' converts NGR easting and nothing to lat, lon.
' input metres, output radians
nX = north
eX = east
a = 6377563.396 ' OSI semi-major
b = 6356256.91 ' OSI semi-minor
e0 = 400000 ' easting of false origin
n0 = -100000 ' northing of false origin
f0 = 0.9996012717 ' OSI scale factor on central meridian
e2 = 0.0066705397616 ' OSI eccentricity squared
lam0 = -0.034906585039886591 ' OSI false east
phi0 = 0.85521133347722145 ' OSI false north
af0 = a * f0
bf0 = b * f0
n = (af0 - bf0) / (af0 + bf0)
Et = east - e0
phid = InitialLat(north, n0, af0, phi0, n, bf0)
nu = af0 / (sqr(1 - (e2 * (sin(phid) * sin(phid)))))
rho = (nu * (1 - e2)) / (1 - (e2 * (sin(phid)) * (sin(phid))))
eta2 = (nu / rho) - 1
tlat2 = tan(phid) * tan(phid)
tlat4 = (tan(phid)^ 4)
tlat6 = (tan(phid)^ 6)
clatm1 = (cos(phid)^ -1)
VII = tan(phid) / (2 * rho * nu)
VIII = (tan(phid) / (24 * rho * (nu * nu * nu))) * (5 + (3 * tlat2) + eta2 - (9 * eta2 * tlat2))
IX = ((tan(phid)) / (720 * rho * (nu^ 5))) * (61 + (90 * tlat2) + (45 * (tan(phid)^ 4) ))
phip = (phid - ((Et * Et) * VII) + ((Et^ 4) * VIII) - ((Et^ 6) * IX))
X = (cos(phid)^ -1) / nu
XI = (clatm1 / (6 * (nu * nu * nu))) * ((nu / rho) + (2 * (tlat2)))
XII = (clatm1 / (120 * (nu^ 5))) * (5 + (28 * tlat2) + (24 * tlat4))
XIIA = clatm1 / (5040 * (nu^ 7)) * (61 + (662 * tlat2) + (1320 * tlat4) + (720 * tlat6))
lambdap = (lam0 + (Et * X) - ((Et * Et * Et) * XI) + ((Et^ 5) * XII) - ((Et^ 7) * XIIA))
geo = phip &","& lambdap
NEtoLL = geo
end function

function Marc(bf0, n, phi0, phi)
vMarc = bf0 * (((1 + n + ((5 / 4) * (n * n)) + ((5 / 4) * (n * n * n))) * (phi - phi0)) _
- (((3 * n) + (3 * (n * n)) + ((21 / 8) * (n * n * n))) * (sin(phi - phi0)) * (cos(phi + phi0))) _
+ ((((15 / 8) * (n * n)) + ((15 / 8) * (n * n * n))) * (sin(2 * (phi - phi0))) * (cos(2 * (phi + phi0)))) _
- (((35 / 24) * (n * n * n)) * (sin(3 * (phi - phi0))) * (cos(3 * (phi + phi0)))))
Marc = vMarc
end function

function InitialLat(north, n0, af0, phi0, n, bf0)

phi1 = ((north - n0) / af0) + phi0
M = Marc(bf0, n, phi0, phi1)
phi2 = ((north - n0 - M) / af0) + phi1
ind = 0
do while ((abs(north - n0 - M) > 0.00001) and (ind < 20)) ' max 20 iterations in case of error
ind = ind + 1
phi2 = ((north - n0 - M) / af0) + phi1
M = Marc(bf0, n, phi0, phi2)
phi1 = phi2
loop

InitialLat = phi2
end function

Function rad2deg(Rad)
rad2deg = cdbl((Rad * 180) / 3.1415926535897932)
End Function

The character conversion bit

Function GridRefConvert(osref)
char2=Left(osref,2)
e1=""
n1=""
Select case(char2)
Case "SR"
e1="1"
n1="1"
Case "SM"
e1="1"
n1="2"
Case "SS"
e1="2"
n1="1"
Case "SN"
e1="2"
n1="2"
End Select
If e1 = "" OR n1 = "" Then '03
'do nothing or report error if you're not confident in data
Else
erest=mid(osref,3,3)
nrest=mid(osref,6,3)
east=e1 & erest & "00"
north=n1 & nrest & "00"
end if
geo = NEtoLL(east, north)
lat = left(geo,instr(geo,",")-1) '* rad2deg
lon = right(geo,len(geo)-instr(geo,",")) '* rad2deg
lat = rad2deg(lat)
lon = rad2deg(lon)
GridRefConvert = lat & ", " & lon
end function