I am not expecting any responses but am just checking just in case others have solved this. We have latitude and longitude and make imagery in UTM with the datum as WGS84. I have been getting requests from users now wanting imagery in their local coordinate systems. We have solved a few with routines off the web and for making UTM. Recently got a request now for conversion of latitude and longitudes for the British National Grid. After spending several days, i could not get the correct resources. If any one has solved this or has the "ultimate" resource location it would be appreciated. I have lat/long conversion to UTM routines, Irish National Grid and a few others. I can imagine in the future that I will get inundated by requests for State Plane coordinates, and prefectural coordinate systems from other countries.... what a mess..!!
Announcement
Collapse
No announcement yet.
lat long wgs84 to british national grid
Collapse
X

Does this help?
https://digimap.edina.ac.uk/webhelp/..._and_wgs84.htm
https://www.ordnancesurvey.co.uk/bus...et/coordinates
https://www.ordnancesurvey.co.uk/doc...lculations.xls
In particular, that last link is a spreadsheet which does the conversions so maybe you can extract the formula from it.
It says "To see the Visual Basic code of the functions  Menu = Tools, Macro, Visual Basic Editor, then view the code in "Module 1".
and
"2. The constants in the table above are suitable for transforming from WGS84 GPS coordinates to OSGB36 National Grid."
Paul.

Old code  thimk I have extracted all of it
'=========================================================================
CALLBACK FUNCTION LatLonProc()
LOCAL spare AS STRING
LOCAL n AS LONG
LOCAL temp AS SINGLE
SELECT CASE AS LONG CB.MSG
CASE %WM_INITDIALOG
CONTROL SET OPTION CB.HNDL, 2002, 2001,2002
DECILAT=0:decilon=0
' Initialization handler
CASE %WM_NCACTIVATE
STATIC hWndSaveFocus AS DWORD
IF ISFALSE CB.WPARAM THEN
' Save control focus
hWndSaveFocus = GetFocus()
ELSEIF hWndSaveFocus THEN
' Restore control focus
SetFocus(hWndSaveFocus)
hWndSaveFocus = 0
END IF
CASE %WM_COMMAND
' Process control notifications
SELECT CASE AS LONG CB.CTL
CASE 5001
CONTROL GET TEXT CB.HNDL,1007 TO spare
decilat=VAL(spare)
CONTROL GET TEXT CB.HNDL,1008 TO spare
decilon=VAL(spare)
IF decilat<>0 AND decilon <>0 THEN
DIALOG END CB.HNDL,1
EXIT SELECT
END IF
CONTROL GET TEXT CB.HNDL,1001 TO spare
Nlat=spare+"?"
CONTROL GET TEXT CB.HNDL,1002 TO spare
Nlat=Nlat+spare
CONTROL GET TEXT CB.HNDL,1003 TO spare
temp=VAL(spare)/60
spare =TRIM$(STR$(temp))
Nlat=Nlat+spare+CHR$(34)+"N"+CHR$(34)
CONTROL GET TEXT CB.HNDL,1004 TO spare
Nlon=spare+"?"
CONTROL GET TEXT CB.HNDL,1005 TO spare
NLon=NLon+spare
CONTROL GET TEXT CB.HNDL,1006 TO spare
temp=VAL(spare)/60
spare =TRIM$(STR$(temp))
NLon=NLon+spare+CHR$(34)
CONTROL GET CHECK CB.HNDL,2001 TO n
IF n=1 THEN
NLon=NLon+"E"+CHR$(34)
ELSE
NLon=NLon+"W"+CHR$(34)
END IF
DIALOG END CB.HNDL
CASE 5002
DIALOG END CB.HNDL,0
END SELECT
END SELECT
END FUNCTION
'
FUNCTION GetLatitudeandLongitude(BYVAL hParent AS DWORD) AS LONG
LOCAL lRslt AS LONG
LOCAL hDlg AS LONG
#PBFORMS BEGIN DIALOG %IDD_DIALOG2>>
DIALOG NEW PIXELS, hParent, "Latitude/Longitude conversion", , , 302, 294, %WS_POPUP OR %WS_VISIBLE OR %DS_SYSMODAL OR %DS_3DLOOK OR %DS_NOFAILCREATE OR %DS_SETFONT, %WS_EX_CONTROLPARENT OR %WS_EX_LEFT OR %WS_EX_LTRREADING OR _
%WS_EX_RIGHTSCROLLBAR, TO hDlg
CONTROL ADD FRAME, hDlg, 1, "Latitude", 8, 32, 284, 108
CONTROL ADD FRAME, hDlg, 1, "Longitude", 8, 144, 284, 108
CONTROL ADD BUTTON, hDlg, 5001, "Convert", 24, 260, 88, 21
CONTROL ADD BUTTON, hDlg, 5002, "Cancel", 152, 260, 88, 21
CONTROL ADD LABEL, hDlg, 1, "Enter data in the boxes below and click on Convert.", 4, 8, 292, 21
CONTROL ADD LABEL, hDlg, 1, "Degs.", 16, 49, 36, 14, %WS_CHILD OR %WS_VISIBLE OR %SS_CENTER, %WS_EX_LEFT OR %WS_EX_LTRREADING
CONTROL ADD LABEL, hDlg, 1, "Mins.", 72, 49, 36, 14, %WS_CHILD OR %WS_VISIBLE OR %SS_CENTER, %WS_EX_LEFT OR %WS_EX_LTRREADING
CONTROL ADD LABEL, hDlg, 1, "Secs.", 128, 49, 36, 14, %WS_CHILD OR %WS_VISIBLE OR %SS_CENTER, %WS_EX_LEFT OR %WS_EX_LTRREADING
CONTROL ADD LABEL, hDlg, 1, "N", 176, 67, 20, 21
CONTROL ADD LABEL, hDlg, 1, "Secs.", 128, 162, 36, 17, %WS_CHILD OR %WS_VISIBLE OR %SS_CENTER, %WS_EX_LEFT OR %WS_EX_LTRREADING
CONTROL ADD LABEL, hDlg, 1, "Mins.", 72, 162, 36, 17, %WS_CHILD OR %WS_VISIBLE OR %SS_CENTER, %WS_EX_LEFT OR %WS_EX_LTRREADING
CONTROL ADD LABEL, hDlg, 1, "Degs.", 16, 162, 36, 17, %WS_CHILD OR %WS_VISIBLE OR %SS_CENTER, %WS_EX_LEFT OR %WS_EX_LTRREADING
CONTROL ADD LABEL, hDlg, 1, "Decimal", 16, 92, 138, 16
CONTROL ADD LABEL, hDlg, 1, "Decimal", 16, 204, 138, 16
CONTROL ADD TEXTBOX, hDlg, 1001, "", 16, 63, 36, 21
CONTROL ADD TEXTBOX, hDlg, 1002, "", 72, 63, 36, 21
CONTROL ADD TEXTBOX, hDlg, 1003, "", 128, 63, 44, 21
CONTROL ADD TEXTBOX, hDlg, 1004, "", 16, 179, 36, 21
CONTROL ADD TEXTBOX, hDlg, 1005, "", 72, 179, 36, 21
CONTROL ADD TEXTBOX, hDlg, 1006, "", 128, 179, 44, 21
CONTROL ADD TEXTBOX, hDlg, 1007, "", 16, 112, 148, 21 'lat
CONTROL ADD TEXTBOX, hDlg, 1008, "", 16, 224, 148, 21 'long
'%WS_GROUP
CONTROL ADD OPTION, hDlg, 2001, "East", 176, 179, 44, 21, %WS_CHILD OR %WS_VISIBLE OR %WS_GROUP OR %WS_TABSTOP OR %BS_TEXT OR %BS_AUTORADIOBUTTON OR %BS_LEFT OR %BS_VCENTER, %WS_EX_LEFT OR %WS_EX_LTRREADING
CONTROL ADD OPTION, hDlg, 2002, "West", 232, 179, 48, 21
#PBFORMS END DIALOG
DIALOG SHOW MODAL hDlg, CALL LatLonProc TO lRslt
#PBFORMS BEGIN CLEANUP %IDD_DIALOG2
#PBFORMS END CLEANUP
FUNCTION = lRslt
END FUNCTION
'
'================================================
'########################### entry point ################################
FUNCTION Latlon2grid () AS LONG
GLOBAL deciLat, deg2rad, rad2deg, deciLon AS DOUBLE
LOCAL ngr AS STRING
deg2rad = 3.141592653589793 / 180
rad2deg = 180.0 / 3.141592653589793
GetLatitudeandLongitude gmain
ngr=ParseNMEA (Nlat,Nlon,0)
MSGBOX "Nat. Grid ref. = "+ngr,%MB_TASKMODAL,"Masterdrain
natgrid=ngr
END FUNCTION
'=======================================================================
FUNCTION ParseNMEA(Nlat AS STRING, Nlon AS STRING, height AS DOUBLE) AS STRING
' 'Processes WGS84 lat and lon in NMEA form
' '52?9.1461"N 002?3.3717"W
LOCAL q AS LONG
IF decilat=0 AND decilon=0 THEN
' //grab the bit up to the ?
q=INSTR(Nlat,"?")
deciLat = VAL(LEFT$(Nlat,q1))
' remove that bit from the string now we've used it and the ?symbol
Nlat = RIGHT$(Nlat,LEN(Nlat)q)
q=INSTR(Nlon,"?")
deciLon = VAL(LEFT$(Nlon,q1))
' remove that bit from the string now we've used it and the ?symbol
Nlon = RIGHT$(Nlon,LEN(Nlon)q)
'
' //grab the bit up to the "  divide by 60 to convert to degrees and add it to our double value
q=INSTR(Nlat,CHR$(34))
deciLat=deciLat+ (VAL(LEFT$(Nlat,q1))/60)
q=INSTR(Nlon,CHR$(34))
deciLon=deciLon+ (VAL(LEFT$(Nlon,q1))/60)
' 'check for negative directions
IF INSTR(Nlat,"S")<>0 THEN deciLat = 0  deciLat
IF INSTR(Nlon,"W")<>0 THEN deciLon = 0  deciLon
END IF
' //now we can parse them
FUNCTION= Transform(deciLat, deciLon, height)
END FUNCTION
'=============================================================================
' 'Processes WGS84 lat and lon in decimal form with S and W as ve
FUNCTION Transform(WGlat AS DOUBLE, WGlon AS DOUBLE, height AS DOUBLE) AS STRING
' first off convert to radians
LOCAL radWGlat,radWGlon AS DOUBLE
radWGlat = WGlat * deg2rad
radWGlon = WGlon * deg2rad
' these calculations were derived from the work of
' Roger Muggleton (http://www.carabus.co.uk/)
'
' quoting Roger Muggleton :
' There are many ways to convert data from one system to another, the most accurate
' being the most complex! For this example I shall use a 7 parameter Helmert
' transformation.
' The process is in three parts:
' (1) Convert latitude and longitude to Cartesian coordinates (these also include height data, and have three parameters, X, Y and Z).
' (2) Transform to the new system by applying the 7 parameters and using a little maths.
' (3) Convert back to latitude and longitude.
' For the example we shall transform a GRS80 location to Airy, e.g. a GPS reading to the Airy spheroid.
' The following code converts latitude and longitude to Cartesian coordinates. The input parameters are: WGS84 latitude and longitude,
' axis is the GRS80/WGS84 major axis in metres, ecc is the eccentricity, and height is the height above the ellipsoid.
' v = axis / (Math.sqrt (1  ecc * (Math.pow (Math.sin(lat), 2))));
' x = (v + height) * Math.cos(lat) * Math.cos(lon);
' y = (v + height) * Math.cos(lat) * Math.sin(lon);
' z = ((1  ecc) * v + height) * Math.sin(lat);
' The transformation requires the 7 parameters: xp, yp and zp correct the coordinate origin, xr, yr and zr correct the orientation of
' the axes, and sf deals with the changing scale factors.
'
' these are the values for WGS86(GRS80) to OSGB36(Airy)
LOCAL a,e,h,a2,e2,xp,yp,zp,xr,yr,zr,s AS DOUBLE
a = 6378137 ' WGS84_AXIS
e = 0.00669438037928458 ' WGS84_ECCENTRIC
h = height ' height above datum (from GPS GGA sentence)
a2 = 6377563.396 ' OSGB_AXIS
e2 = 0.0066705397616 ' OSGB_ECCENTRIC
xp = 446.448
yp = 125.157
zp = 542.06
xr = 0.1502
yr = 0.247
zr = 0.8421
s = 20.4894
' convert to cartesian; lat, lon are radians
LOCAL sf,v,x,y,z AS DOUBLE
sf = s * 0.000001
v = a / (SQR(1  (e * (SIN(radWGlat) * SIN(radWGlat)))))
x = (v + h) * COS(radWGlat) * COS(radWGlon)
y = (v + h) * COS(radWGlat) * SIN(radWGlon)
z = ((1  e) * v + h) * SIN(radWGlat)
' // transform cartesian
LOCAL xrot,yrot, zrot,hx,hy,hz AS DOUBLE
xrot = (xr / 3600) * deg2rad
yrot = (yr / 3600) * deg2rad
zrot = (zr / 3600) * deg2rad
hx = x + (x * sf)  (y * zrot) + (z * yrot) + xp
hy = (x * zrot) + y + (y * sf)  (z * xrot) + yp
hz = (1 * x * yrot) + (y * xrot) + z + (z * sf) + zp
'
' Convert back to lat, lon
LOCAL newLon,p,newLat,errvalue,lat0 AS DOUBLE
newLon = ATN(hy / hx)
p = SQR((hx * hx) + (hy * hy))
newLat = ATN(hz / (p * (1  e2)))
v = a2 / (SQR(1  (e2 * SIN(newLat * SIN(newLat)))))
errvalue = 1.0
lat0 = 0
DO
lat0 = ATN((hz + (e2 * v * SIN(newLat))) / p)
errvalue = ABS(lat0  newLat)
newLat = lat0
LOOP WHILE (errvalue > 0.001)
'
' convert back to degrees
newLat = newLat * rad2deg
newLon = newLon * rad2deg
' //convert lat and lon (OSGB36) to OS 6 figure northing and easting
FUNCTION=LLtoNE(newLat, newLon)
END FUNCTION
'==========================================================================
'===========================================================================
' converts lat and lon (OSGB36) to OS 6 figure northing and easting
FUNCTION LLtoNE(lat AS DOUBLE, lon AS DOUBLE) AS STRING
LOCAL phi,lam,a,b,e0,n0,f0,e2,lam0,phi0,af0,bf0 AS DOUBLE
phi = lat * deg2rad ' convert latitude to radians
lam = lon * deg2rad ' convert longitude to radians
a = 6377563.396 ' OSGB semimajor axis
b = 6356256.91 ' OSGB semiminor axis
IF country="eire" THEN
e0 = 200000 ' easting of false origin UK
n0 = 250000 ' northing of false origin UK
ELSE
e0 = 400000 ' easting of false origin UK
n0 = 100000 ' northing of false origin UK
END IF
f0 = 0.9996012717 ' OSGB scale factor on central meridian
e2 = 0.0066705397616 ' OSGB eccentricity squared
lam0 = 0.034906585039886591 ' OSGB false east
phi0 = 0.85521133347722145 ' OSGB false north
af0 = a * f0
bf0 = b * f0
'
' easting
LOCAL slat2,nu,rho,eta2,p,IV,clat3,tlat2,V,clat5,tlat4,VI,east AS DOUBLE
slat2 = SIN(phi)^2
nu = af0 / (SQR(1  (e2 * (slat2))))
rho = (nu * (1  e2)) / (1  (e2 * slat2))
eta2 = (nu / rho)  1
p = lam  lam0
IV = nu * COS(phi)
clat3 = COS(phi)^3
tlat2 = TAN(phi)^2
V = (nu / 6) * clat3 * ((nu / rho)  tlat2)
clat5 = COS(phi)^5
tlat4 = TAN(phi)^4
VI = (nu / 120) * clat5 * ((5  (18 * tlat2)) + tlat4 + (14 * eta2)  (58 * tlat2 * eta2))
east = e0 + (p * IV) + ((p^3) * V) + ((p^5) * VI)
'
' northing
LOCAL n,M,I,II,III,IIIA,north AS DOUBLE
n = (af0  bf0) / (af0 + bf0)
M = Marc(bf0, n, phi0, phi)
I = M + (n0)
II = (nu / 2) * SIN(phi) * COS(phi)
III = ((nu / 24) * SIN(phi) * (COS(phi)^3)) * (5  (TAN(phi)^2) + (9 * eta2))
IIIA = ((nu / 720) * SIN(phi) * clat5) * (61  (58 * tlat2) + tlat4)
north = I + ((p * p) * II) + ((p^4) * III) + ((p^6) * IIIA)
' make whole number values
east = ABS(ROUND(east,0)) ' round to whole number
north = ROUND(north,0) ' round to whole number
?"e "+STR$(east)+" n "+STR$(north)
' convert to nat grid ref
FUNCTION=NE2NGR(east, north)
END FUNCTION
'===============================================================================
FUNCTION Marc(bf0 AS DOUBLE, n AS DOUBLE, phi0 AS DOUBLE, phi AS DOUBLE) AS DOUBLE
LOCAL M AS DOUBLE
M= 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)))))
FUNCTION=M
END FUNCTION
'===============================================================================
' //convert 12 (6e & 6n) figure numeric to letter and number grid system
FUNCTION NE2NGR(east AS DOUBLE, north AS DOUBLE) AS STRING
LOCAL eX,nX,tmp AS DOUBLE
LOCAL eing,ning, ngr AS STRING
' LOCAL lnth AS LONG
eX = east / 500000
nX = north / 500000
tmp = (INT(eX)  (5.0 * INT(nX)) + 17.0) 'INT Returns the largest integer less than or equal to the specified number
nX = 5 * (nX  INT(nX))
eX = 20  (5.0 * INT(nX)) + INT(5.0 * (eX  INT(eX)))
IF (eX > 7.5) THEN eX = eX + 1
IF (tmp > 7.5) THEN tmp = tmp + 1
eing = TRIM$(STR$(east))
eing = RIGHT$("000000"+eing,5)
ning = TRIM$(STR$(north))
ning = RIGHT$("000000"+ning,5)
ngr = CHR$(tmp + 65) + CHR$(eX + 65) + " " + eing + " " + ning
FUNCTION=ngr
END FUNCTION
'###############################################################################################################################
'===========================================================================
[/CODE]“None but those who have experienced them can conceive of the enticements of science”  Mary Shelley
Comment

Iain, I was able to integrate your code for our conversion. Thank you so much for helping out. I did a test with the British National Grid website and the easting is spot on and the northing is off by less than a meter. I know that they do some yearly updates to variables in the BNG conversion to account for seafloor spreading, so I will check to see if all the latest variables and transformation constant are up to date.... Thanks again for your help. Paul, also thanks for your help too... I was just getting started to locate the visual basic routines and then an unexpected gift came!! cheers.
Comment

Originally posted by dean goodman View PostIain, I was able to integrate your code for our conversion. Thank you so much for helping out. I did a test with the British National Grid website and the easting is spot on and the northing is off by less than a meter. I know that they do some yearly updates to variables in the BNG conversion to account for seafloor spreading, so I will check to see if all the latest variables and transformation constant are up to date.... Thanks again for your help. Paul, also thanks for your help too... I was just getting started to locate the visual basic routines and then an unexpected gift came!! cheers.
Dan
Comment

That would explain the difference... First the code from the link that Michael referred I tried and that was off by 100 meters on the eastings etc... The older code from Iain does have the Helmert Transform and that would explain the difference. Wow...a look up table.... So this is becoming nasty.....Thanks Dan....
Comment

The VB code is not as obvious to set up and it also contains the Helmert transformation. For now, i am going to go with Iain' code post and if there are any complaints from our UK users being off a meter or so, I will tell them to import their data already in OSGB 36 coordinates..... Thanks for everyone's help..
Comment

Originally posted by dean goodman View PostThat would explain the difference... First the code from the link that Michael referred I tried and that was off by 100 meters on the eastings etc... The older code from Iain does have the Helmert Transform and that would explain the difference. Wow...a look up table.... So this is becoming nasty.....Thanks Dan....
The OS has downloadable software to do the conversions, called Grid Inquest II. The previous version included a COM library that might've suited your purposes, only I could never get it to work with PB. I don't see any such library (COM or otherwise) with the current software.
Unless something's changed recently, the OS spreadsheet uses only the Helmert transform, and not OSTN15, so if you're using Helmert too you should get a good match with all of your output figures. The OS website meanwhile uses OSTN15 only, so there you'd expect differences of several metres.
Dan
Comment

Well, it turns out that being off a meter or so for British National Grid conversion from WGS84 using the Helmert transform is not good enough for our UK users. I have located a DLL which looks like it can do the conversion as recommend by Dan Soper. There is a snippet of Visual Basic code to implement at Grid Inquest II:
' LibGridInQuestII Coordinates record declaration.
Public Type Coordinates X As Double Y As Double Z As Double End Type
' LibGridInQuestII API declarations.
Declare Function ConvertCoordinates
Lib "GIQ.dll" (ByVal SourceSRID As Integer, ByVal TargetSRID As Integer, ByVal SourceRevision As Integer, ByVal TargetRevision As Integer,
ByRef InputCoordinates As Coordinates, ByRef OutputCoordinates As Coordinates, ByRef Datum As Integer) As Boolean
Does anyone know what this might look like converted to PowerBasic
There is an option to run external bat files which seems to work but would hope to have it internal to the software without DOS windows popping up...
thx.
Comment

The download https://bitbucket.org/PaulFMichell/g...in320101c.zip
contains an Access .mdb. The three modules in it would be trivial to convert to a PB application which does the same thing.
Comment

Thanks Stuart. I had the same link with the files needed for the OSTN 15 lookup table. Now I understand why Dan Soper at post #15 could not get it to work. With my test code here, the execution says that the ConvertCoordinates procedure entry point could not be located. So the GIQ.dll that was posted for this conversion has some error...if not the code below.
#COMPILE EXE
#DIM ALL
#INCLUDE "win32api.inc"
TYPE Coordinates
X AS DOUBLE
Y AS DOUBLE
Z AS DOUBLE
END TYPE
DECLARE FUNCTION ConvertCoordinates _
LIB "GIQ.dll" (BYVAL SourceSRID AS INTEGER,_
BYVAL TargetSRID AS INTEGER,_
BYVAL SourceRevision AS INTEGER,_
BYVAL TargetRevision AS INTEGER,_
BYREF InputCoordinates AS Coordinates,_
BYREF OutputCoordinates AS Coordinates,_
BYREF Datum AS INTEGER) AS LONG
FUNCTION PBMAIN () AS LONG
LOCAL sourcesrid AS INTEGER
LOCAL targetsrid AS INTEGER
LOCAL sourcerevision AS INTEGER
LOCAL targetrevision AS INTEGER
LOCAL datum AS INTEGER
DIM inputcoordinates AS coordinates
DIM outputcoordinates AS coordinates
sourcesrid=4937
targetsrid=27700
sourcerevision=2015
targetrevision=0
datum=30
inputcoordinates.x=1.54
inputcoordinates.y=55.0
inputcoordinates.z=0
CALL convertcoordinates (sourcesrid,targetsrid,sourcerevision,targetrevision,inputcoordinates,outputcoordinates,datum)
MSGBOX STR$(outputcoordinates.x)+STR$(outputcoordinates.y)+STR$(outputcoordinates.z)
END FUNCTION
Comment

I should have included it before: Use the Alias clause:
Code:DECLARE FUNCTION ConvertCoordinates _ LIB "GIQ.dll" ALIAS "ConvertCoordinates" (BYVAL SourceSRID AS INTEGER,_ BYVAL TargetSRID AS INTEGER,_ BYVAL SourceRevision AS INTEGER,_ BYVAL TargetRevision AS INTEGER,_ BYREF InputCoordinates AS Coordinates,_ BYREF OutputCoordinates AS Coordinates,_ BYREF Datum AS INTEGER) AS LONG
The ALIAS clause is very important when importing or exporting Subs and Functions from DLLs. Omitting the ALIAS clause or incorrectly capitalizing the alias name are common causes of DLL load failure problems. "
Comment
Comment