--[[ Ordnance Survey / Lat Long functions in Lua
****** ORDNANCE SURVEY TO LAT LONG FUNCTIONS ******
Ported to Lua from PHP by Wikipedia User Hike395, 18 Aug 2019
found by RWH at http://www.megalithia.com/search/llfuncshighlight.php
With thanks to Andy, G4JNT for inspiration in GEOG, and to the OSGB for their white paper on coordinate transformation
describing the iterative method used
thanks to the Ordnance survey of Ireland for details of the true and false origins of the Irish grid
You may use and redistribute this code under the terms of the GPL see http://www.gnu.org/copyleft/gpl.html
Written by Richard
www.megalithia.com
v0.something 27/2/2000
v 1.01 28 June 2004
v 1.02 6 Aug 2004 line 89 add "0" to chars in ngr=stripcharsnotinbag Thx Andy
v 1.03 9 Mar 2005 Jan (Klingon) added conversion to WGS84 map datum and removed limitation of digits of the grid ref
v 1.04 10 Aug 2005 Richard correct error trapping (only manifest on malformed ngrs
This code is predicated on the assumption that your are ONLY feeding it Irish or UK Grid references.
It uses the single char prefix of Irish grid refs to tell the difference, UK grid refs have a two letter prefix.
We would like an even number of digits for the rest of the grid ref.
Anything in the NGR other than 0-9, A-Z, a-z is eliminated.
WARNING this assumes there are no decimal points in your NGR components.
The transformation from OSGB36/Ireland 1965 to WGS84 is more precise than 5 m.
The function is case insensitive
Lat/Long to Ordnance Survey conversion is at bottom of file, see further authorship there
]]
local oscoord = {}
local getArgs = require('Module:Arguments').getArgs
local yesno = require('Module:Yesno')
local namespace = mw.title.getCurrentTitle().namespace
local pow = math.pow
local sqrt = math.sqrt
local abs = math.abs
local floor = math.floor
local sin = math.sin
local cos = math.cos
local tan = math.tan
local atan = math.atan
local deg2rad = math.rad
local rad2deg = math.deg
local dr = math.deg(1.0)
-- Transverse Mercator parameters
local scaleGB = 0.9996012717
local eastingGB = 400000.0
local northingGB = -100000.0
local lat0GB = 49.0
local lon0GB = -2.0
--[[ DATUM SHIFT USING HELMERT TRANSFORMATION
*
* height above the ellipsoid is ignored
* latitude and longitude must be given in decimal degrees
*
* no transformation if abs(latitude)>89 or abs(longitude)>89
* (lifting this restriction requires some more lines of code)
*
* Written in 2009 by user Telford at de.wikipedia
* Based on math described in "A Guide to Coordinate Systems in Great Britain"
* by the UK Ordnance Survey
* URL: https://www.ordnancesurvey.co.uk/documents/resources/guide-coordinate-systems-great-britain.pdf
]]
-- Datum parameters
local OSGBglobe = {
semimajor = 6377563.396,
semiminor = 6356256.910,
ecc = 0.006670539761597529073698869358812557054558
}
local WGSglobe = {
semimajor = 6378137.0,
semiminor = 6356752.3141,
ecc = 0.00669438003551279089034150031998869922791
}
local WGS2OSGBparam = {
-- Table 4, Section 6.6, Ordnance Survey document
semimajor_src = WGSglobe.semimajor,
semiminor_src = WGSglobe.semiminor,
ecc_src = WGSglobe.ecc,
semimajor_dst = OSGBglobe.semimajor,
semiminor_dst = OSGBglobe.semiminor,
ecc_dst = OSGBglobe.ecc,
tx = -446.448,
ty = 125.157,
tz = -542.060,
s0 = 0.0000204894,
rx = deg2rad ( -.1502/3600. ),
ry = deg2rad ( -.2470/3600. ),
rz = deg2rad ( -.8421/3600. )
}
local OSGB2WGSparam = {
semimajor_src = OSGBglobe.semimajor,
semiminor_src = OSGBglobe.semiminor,
ecc_src = OSGBglobe.ecc,
semimajor_dst = WGSglobe.semimajor,
semiminor_dst = WGSglobe.semiminor,
ecc_dst = WGSglobe.ecc,
tx = -WGS2OSGBparam.tx,
ty = -WGS2OSGBparam.ty,
tz = -WGS2OSGBparam.tz,
s0 = -WGS2OSGBparam.s0,
rx = -WGS2OSGBparam.rx,
ry = -WGS2OSGBparam.ry,
rz = -WGS2OSGBparam.rz,
}
local function HelmertDatumShift ( latitude , longitude, param )
-- latitude and longitude are in degrees
-- return if latitude or longitude out of range
if abs(latitude) > 89. or abs(longitude) > 89. then
return {lat=latitude, lon=longitude}
end
param = param or WGS2OSGBparam
-- parameters for ellipsoids
-- Annex A.1, Ordnance Survey document --
local a1 = param.semimajor_src
local b1 = param.semiminor_src
local e1 = param.ecc_src
local a2 = param.semimajor_dst
local b2 = param.semiminor_dst
local e2 = param.ecc_dst
-- convert latitude and longitude to cartesian coordinates
-- math in Annex B.1, Ordnance Survey document
local phi = deg2rad ( latitude )
local lambda = deg2rad ( longitude )
local cosphi = cos ( phi )
local sinphi = sin ( phi )
local coslambda = cos ( lambda )
local sinlambda = sin ( lambda )
local ny = a1 / sqrt ( 1. - e1*sinphi*sinphi )
mw.log("New v",ny)
local x1 = ny * cosphi * coslambda
local y1 = ny * cosphi * sinlambda
local z1 = ny * ( 1. - e1 ) * sinphi
mw.log("New xyz1",x1,y1,z1)
-- the helmert transformation proper
-- Equation 3, Section 6.2, Ordnance Survey document
mw.log("New params",param.s0,param.rx,param.ry,param.rz,param.tx,param.ty,param.tz)
local x2 = x1 + param.tx + ( param.s0 * x1 - param.rz * y1 + param.ry * z1 )
local y2 = y1 + param.ty + ( param.rz * x1 + param.s0 * y1 - param.rx * z1 )
local z2 = z1 + param.tz + ( -param.ry *x1 + param.rx * y1 + param.s0 * z1 )
mw.log("New xyz2",x2,y2,z2)
-- convert cartesian coordinates to latitude and longitude
-- math in Annex B.2, of Ordnance Survey document
lambda = atan ( y2 / x2 )
local p2 = sqrt ( x2*x2 + y2*y2 )
phi = atan ( z2 / p2*(1.-e2) )
local n0 = 101
local phi_alt
repeat
phi_alt = phi
ny = a2 / sqrt ( 1. - e2 * sin(phi) * sin(phi) )
phi = atan ( (z2 + e2 * ny * sin(phi)) / p2 )
n0 = n0 - 1
until abs ( phi - phi_alt ) <= 0.000000000001 or n0 <= 0
return {lat=rad2deg(phi),lon=rad2deg(lambda)}
end
-- LAT LONG TO OS GRID LIBRARY RESUMES HERE
local function northeast(lett,num,shift)
-- split into northings and eastings
local le=mw.ustring.len(num)
if le%2 == 1 then
return {err="Malformed numerical part of NGR"}
end
local pr=le/2
local n = mw.ustring.sub(num,pr+1)
local e = mw.ustring.sub(num,1,pr)
-- Hack to move to center of square: append a 5 to northings and eastings
shift = yesno(shift)
if shift then
n = n.."5"
e = e.."5"
pr = pr+1
end
-- end hack
n = n == '' and 0 or n
e = e == '' and 0 or e
pr = pow(10.0,(5.0-pr))
local T1 = mw.ustring.byte(mw.ustring.sub(lett,1,1))-65
if T1>8 then
T1 = T1-1
end
local T2 = nil
if mw.ustring.len(lett)>1 then
T2 = mw.ustring.byte(mw.ustring.sub(lett,2))-65
if T2>8 then
T2 = T2-1
end
end
return {n=n,e=e,pr=pr,T1=T1,T2=T2}
end
local function GBEN2LL(e,n)
-- True Origin is 2 deg W
local phi0uk=lon0GB
-- True Origin is 49 deg N
local lambda0uk=lat0GB
-- scale factor @ central meridian
local F0uk=scaleGB
-- True origin in 400 km E of false origin
local E0uk=eastingGB
--True origin is 100 km S of false origin
local N0uk=northingGB
-- semi-major axis (in line to equator) 0.996012717 is yer scale @ central meridian
local auk= OSGBglobe.semimajor*F0uk
--semi-minor axis (in line to poles)
local buk=OSGBglobe.semiminor*F0uk
-- flatness=a1-b1/(a1+b1)
local n1uk=0.00167322025032508731869331280635710896296
-- first eccentricity squared=2*f-f^2where f=(a1-b1)/a1
local e2uk=OSGBglobe.ecc
local k=(n-N0uk)/auk+lambda0uk/dr
local nextcounter=0
local j3, j4, j5, j6, m
repeat
nextcounter=nextcounter+1
local k3=k-lambda0uk/dr
local k4=k+lambda0uk/dr
j3=(1.0+n1uk+1.25*pow(n1uk,2.0)+1.25*pow(n1uk,3.0))*k3
j4=(3.0*n1uk+3.0*pow(n1uk,2.0)+2.625*pow(n1uk,3.0))*sin(k3)*cos(k4)
j5=(1.875*pow(n1uk,2.0)+1.875*pow(n1uk,3.0))*sin(2.0*k3)*cos(2.0*k4)
j6=35.0/24.0*pow(n1uk,3.0)*sin(3.0*k3)*cos(3.0*k4)
m=buk*(j3-j4+j5-j6)
k=k+(n-N0uk-m)/auk
until abs(n-N0uk-m)<=0.000000001 or nextcounter>=100
local v=auk/sqrt(1.0-e2uk*pow(sin(k),2.0))
local r=v*(1.0-e2uk)/(1.0-e2uk*pow(sin(k),2.0))
local h2=v/r-1.0
local y1=e-E0uk
j3=tan(k)/(2.0*r*v)
j4=tan(k)/(24.0*r*pow(v,3.0))*(5.0+3.0*pow(tan(k),2.0)+h2-9.0*pow(tan(k),2.0)*h2)
j5=tan(k)/(720.0*r*pow(v,5.0))*(61.0+90.0*pow(tan(k),2.0)+45.0*pow(tan(k),4.0))
local k9=k-y1*y1*j3+pow(y1,4.0)*j4-pow(y1,6.0)*j5
j6=1.0/(cos(k)*v)
local j7=1.0/(cos(k)*6.0*pow(v,3.0))*(v/r+2.0*pow(tan(k),2.0))
local j8=1.0/(cos(k)*120.0*pow(v,5.0))*(5.0+28.0*pow(tan(k),2.0)+24.0*pow(tan(k),4.0))
local j9=1.0/(cos(k)*5040.0*pow(v,7.0))
local j9=j9*(61.0+662.0*pow(tan(k),2.0)+1320.0*pow(tan(k),4.0)+720.0*pow(tan(k),6.0))
local long=phi0uk+dr*(y1*j6-y1*y1*y1*j7+pow(y1,5.0)*j8-pow(y1,7.0)*j9)
local lat=k9*dr
local helmert = HelmertDatumShift ( lat, long, OSGB2WGSparam)
local v=OSGBglobe.semimajor/sqrt(1.0-e2uk*pow(sin(k),2.0))
mw.log("Old v",v)
local cartxa=v*cos(k9)*cos(long/dr)
local cartya=v*cos(k9)*sin(long/dr)
local cartza=(1.0-e2uk)*v*sin(k9)
mw.log("Old xyz1",cartxa,cartya,cartza)
-- Helmert-Transformation from OSGB36 to WGS84 map date
local rotx=-0.1502/3600.0/dr
local roty=-0.2470/3600.0/dr
local rotz=-0.8421/3600.0/dr
local scale=-20.4894/1000000.0
mw.log("Old params",scale,-rotx,-roty,-rotz,446.448,-125.157,542.06)
local cartxb=446.448+(1.0+scale)*cartxa+rotz*cartya-roty*cartza
local cartyb=-125.157-rotz*cartxa+(1.0+scale)*cartya+rotx*cartza
local cartzb=542.06+roty*cartxa-rotx*cartya+(1.0+scale)*cartza
mw.log("Old xyz2",cartxb,cartyb,cartzb)
-- convert Cartesian to long/lat
local awgs84=WGSglobe.semimajor
local bwgs84=WGSglobe.semiminor
local e2wgs84=WGSglobe.ecc
local lambdaradwgs84=atan(cartyb/cartxb)
long=lambdaradwgs84*dr
local pxy=sqrt(pow(cartxb,2.0)+pow(cartyb,2.0))
local phiradwgs84
local phinewwgs84=atan(cartzb/pxy/(1.0-e2wgs84))
nextcounter=0
repeat
phiradwgs84=phinewwgs84
nextcounter=nextcounter+1
v=awgs84/sqrt(1.0-e2wgs84*pow(sin(phiradwgs84),2.0))
phinewwgs84=atan((cartzb+e2wgs84*v*sin(phiradwgs84))/pxy)
until abs(phinewwgs84-phiradwgs84)<=0.000000000001 or nextcounter>=100
lat=phinewwgs84*dr
return {region="GB",lat=lat,long=long}
end
local function GB2LL(lett,num)
-- British OS to Lat+Long
-- first caclulate e,n
-- computing e and n exactly, to get SW corner of box
local ne = northeast(lett,num)
if ne.err then
return {region="GB",err=ne.err}
end
-- use British definition of e and n
local e=500000.0*(ne.T1%5)+100000.0*(ne.T2%5)-1000000.0+ne.e*ne.pr
local n=1900000.0-500000.0*floor(ne.T1/5)-100000.0*floor(ne.T2/5)+ne.n*ne.pr
local result = GBEN2LL(e,n)
result.prec = 0.8165*ne.pr
return result
end
local function IrishEN2LL(e,n)
-- True Origin is 8 deg W
local phi0ir=-8.0
-- True Origin is 53.5 deg N
local lambda0ir=53.5
-- scale factor @ central meridian
local F0ir=1.000035
-- True origin in 200 km E of false origin
local E0ir=200000.0
--True origin is 250km N of false origin
local N0ir=250000.0
-- semi-major axis (in line to equator) 1.000035 is yer scale @ central meridian
local air=6377340.189*F0ir
--semi-minor axis (in line to poles)
local bir=6356034.447*F0ir
-- flatness=a1-b1/(a1 + b1)
local n1ir=0.001673220384152058651484728058385228837777
-- first eccentricity squared=2*f-f^2 where f=(a1-b1)/a1
local e2ir=0.006670540293336110419293763349975612794125
local k=(n-N0ir)/air+lambda0ir/dr
local nextcounter=0
local j3,j4,j5,j6,m
repeat
nextcounter=nextcounter+1
local k3=k-lambda0ir/dr
local k4=k+lambda0ir/dr
j3=(1.0+n1ir+1.25*pow(n1ir,2.0)+1.25*pow(n1ir,3.0))*k3
j4=(3.0*n1ir+3.0*pow(n1ir,2.0)+2.625*pow(n1ir,3.0))*sin(k3)*cos(k4)
j5=(1.875*pow(n1ir,2.0)+1.875*pow(n1ir,3.0))*sin(2.0*k3)*cos(2.0*k4)
j6=35.0/24.0*pow(n1ir,3.0)*sin(3.0*k3)*cos(3.0*k4)
m=bir*(j3-j4+j5-j6)
k=k+(n-N0ir-m)/air
until abs(n-N0ir-m)<=0.000000000001 or nextcounter>=10000
local v=air/sqrt(1.0-e2ir*pow(sin(k),2.0))
local r=v*(1.0-e2ir)/(1.0-e2ir*pow(sin(k),2.0))
local h2=v/r-1.0
local y1=e-E0ir
j3=tan(k)/(2.0*r*v)
j4=tan(k)/(24.0*r*pow(v,3.0))*(5.0+3.0*pow(tan(k),2.0)+h2-9.0*pow(tan(k),2.0)*h2)
j5=tan(k)/(720.0*r*pow(v,5.0))*(61.0+90.0*pow(tan(k),2.0)+45.0*pow(tan(k),4.0))
local k9=k-y1*y1*j3+pow(y1,4.0)*j4-pow(y1,6.0)*j5
j6=1.0/(cos(k)*v)
local j7=1.0/(cos(k)*6.0*pow(v,3.0))*(v/r+2.0*pow(tan(k),2.0))
local j8=1.0/(cos(k)*120.0*pow(v,5.0))*(5.0+28.0*pow(tan(k),2.0)+24.0*pow(tan(k),4.0))
local j9=1.0/(cos(k)*5040.0*pow(v,7.0))
local j9=j9*(61.0+662.0*pow(tan(k),2.0)+1320.0*pow(tan(k),4.0)+720.0*pow(tan(k),6.0))
local long=phi0ir+dr*(y1*j6-y1*y1*y1*j7+pow(y1,5.0)*j8-pow(y1,7.0)*j9)
local lat=k9*dr
-- convert long/lat to Cartesian coordinates
v=6377340.189/sqrt(1.0-e2ir*pow(sin(k),2.0))
local cartxa=v*cos(k9)*cos(long/dr)
local cartya=v*cos(k9)*sin(long/dr)
local cartza=(1.0-e2ir)*v*sin(k9)
-- Helmert-Transformation from Ireland 1965 to WGS84 map date
local rotx=1.042/3600.0/dr
local roty=0.214/3600.0/dr
local rotz=0.631/3600.0/dr
local scale=8.15/1000000.0
local cartxb=482.53+(1.0+scale)*cartxa+rotz*cartya-roty*cartza
local cartyb=-130.596-rotz*cartxa+(1.0+scale)*cartya+rotx*cartza
local cartzb=564.557+roty*cartxa-rotx*cartya+(1.0+scale)*cartza
-- convert Cartesian to long/lat
local awgs84=WGSglobe.semimajor
local bwgs84=WGSglobe.semiminor
local e2wgs84=WGSglobe.ecc
local lambdaradwgs84=atan(cartyb/cartxb)
long=lambdaradwgs84*dr
local pxy=sqrt(pow(cartxb,2.0)+pow(cartyb,2.0))
local phinewwgs84=atan(cartzb/pxy/(1.0-e2wgs84))
local phiradwgs84
nextcounter=0
repeat
phiradwgs84=phinewwgs84
nextcounter=nextcounter+1
v=awgs84/sqrt(1.0-e2wgs84*pow(sin(phiradwgs84),2.0))
phinewwgs84=atan((cartzb+e2wgs84*v*sin(phiradwgs84))/pxy)
until abs(phinewwgs84-phiradwgs84)<=0.000000000001 or nextcounter>=10000
lat=phinewwgs84*dr
return {region="IE",lat=lat,long=long}
end
local function Irish2LL(lett,num)
-- Irish OS to Lat+Long
-- first caclulate e,n
-- computing e and n exactly, to get SW corner of box
local ne = northeast(lett,num)
if ne.err then
return {region="IE", err=ne.err}
end
-- use Irish definition of northing and easting
local e=100000.0*(ne.T1%5.0)+ne.e*ne.pr
local n=ne.n*ne.pr+100000.0*(4.0-floor(ne.T1/5.0))
local result = IrishEN2LL(e,n)
result.prec = 0.8165*ne.pr -- useful @ Commons
return result
end
local function empty(s)
return not s or s == ''
end
local function NGR2LL(ngr)
local result = {}
ngr, _ = mw.ustring.gsub(mw.ustring.upper(ngr),"[%s%p]","")
local first, last, lett, num = mw.ustring.find(ngr,"^([A-Z]+)(%d+)$")
if not first or empty(lett) or empty(num) or mw.ustring.len(lett) > 2 then
return {err="Malformed NGR"}
end
if mw.ustring.len(lett) == 1 then
return Irish2LL(lett,num)
end
return GB2LL(lett, num)
end
local function split(s,sep)
-- split a string s into chunks, separated by sep
sep = sep or "%s"
local t = {}
for chunk in mw.ustring.gmatch(s,"([^"..sep.."]+)") do
table.insert(t,chunk)
end
return t
end
local function trim(s)
s, _ = mw.ustring.gsub(s,"^%s+","")
s, _ = mw.ustring.gsub(s,"%s+$","")
return s
end
local function alldigits(s)
return not mw.ustring.find(s,"%D")
end
local function warning(errmsg)
local preview = require('Module:If preview')
local msg = errmsg or 'Empty OS grid ref'
local html = preview._warning({ msg })
if namespace == 0 and errmsg then
html = html..'[[Category:Pages with malformed OS coordinates]]'
end
return html
end
function oscoord._main(args)
local input = args[1]
if empty(input) then
return warning(nil)
end
local linktitle = args[2]
local namearg = args["name"]
local rawurl = yesno(args["rawurl"])
local args = split(input,'_')
local LL
local restargs = 1
local current_page = mw.title.getCurrentTitle()
local pagename = mw.uri.encode( current_page.prefixedText, 'WIKI' )
if #args >= 2 and alldigits(args[2]) then
if mw.ustring.sub(args[1],1,1) == 'i' then
local firstArg = mw.ustring.sub(args[1],2)
if alldigits(firstArg) then
LL = IrishEN2LL(firstArg,args[2])
restargs = 3
if empty(linktitle) then
linktitle=args[1]..'_'..args[2]
end
end
elseif alldigits(args[1]) then
LL = GBEN2LL(args[1],args[2])
restargs = 3
if empty(linktitle) then
linktitle=args[1]..'_'..args[2]
end
end
else
LL = NGR2LL(args[1])
restargs = 2
if empty(linktitle) then
linktitle=args[1]
end
end
linktitle = trim(linktitle)
if not empty(LL.err) then
return linktitle ..warning(LL.err)
end
-- https://geohack.toolforge.org/geohack.php?pagename=Mount_Whitney¶ms=36.578580925_N_118.29199495_W_type:mountain_region:US-CA_scale:100000_source:NGS
local url = ''
if not rawurl then
url = url..'['
end
url = url..'https://geohack.toolforge.org/geohack.php?'
if not empty(pagename) then
url = url..'pagename='..pagename..'&'
end
LL.lat = LL.lat or 0
LL.long = LL.long or 0
url = url..'params='..LL.lat..'_N_'
if LL.long < 0 then
url = url..(-LL.long)..'_W'
else
url = url..LL.long..'_E'
end
for i = restargs,#args do
url = url..'_'..args[i]
end
if not mw.ustring.find(input,"region") and LL.region then
url = url..'_region:'..LL.region
end
if not mw.ustring.find(input,"scale") and
not mw.ustring.find(input,"type") and
not mw.ustring.find(input,"dim") and LL.prec then
url = url..'_dim:'..floor(50*LL.prec+0.5)..'m'
end
if not empty(namearg) then
url = url .. "&title=" .. mw.uri.encode(namearg)
end
if not rawurl then
url = url..' '..linktitle..']'
end
return url
end
function oscoord._oscoord(args)
local output = '<span class="plainlinks nourlexpansion" style="white-space: nowrap">' .. oscoord._main(args) .. '</span>'
if namespace == 0 then
output = output .. '[[Category:Articles with OS grid coordinates]]'
end
return output
end
function oscoord.main(frame)
local args = getArgs(frame)
return oscoord._main(args)
end
function oscoord.oscoord(frame)
local args = getArgs(frame)
return oscoord._oscoord(args)
end
--[[
****** LAT/LONG CONVERSION TO OS GRID REF FUNCTIONS *****
* Uses the WGS-84 ellipsoid and only works for GB grid ref
*
* See also:
* http://www.posc.org/Epicentre.2_2/DataModel/ExamplesofUsage/eu_cs34h.html
* http://kanadier.gps-info.de/d-utm-gitter.htm
* http://www.gpsy.com/gpsinfo/geotoutm/
* http://www.colorado.edu/geography/gcraft/notes/gps/gps_f.html
* http://search.cpan.org/src/GRAHAMC/Geo-Coordinates-UTM-0.05/
* UK Ordnance Survey grid (OSBG36): http://www.gps.gov.uk/guidecontents.asp
* Swiss CH1903: http://www.gps.gov.uk/guidecontents.asp
*
* ----------------------------------------------------------------------
*
* Copyright 2005, Egil Kvaleberg <[email protected]>
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
Converted to Lua by User:Hike395 on 2023-12-15
]]
local function find_M ( lat_rad )
local e = OSGBglobe.ecc
return OSGBglobe.semimajor * ( ( 1 - e/4 - 3 * e * e/64
- 5 * e * e * e/256
) * lat_rad
- ( 3 * e/8 + 3 * e * e/32
+ 45 * e * e * e/1024
) * sin(2 * lat_rad)
+ ( 15 * e * e/256 +
45 * e * e * e/1024
) * sin(4 * lat_rad)
- ( 35 * e * e * e/3072
) * sin(6 * lat_rad) )
end
--[[
* Convert latitude, longitude in decimal degrees to
* Transverse Mercator Easting and Northing based on a GB origin
*
* return nil if problems
]]
local function LatLonOrigin2TM( latitude, longitude )
if longitude < -180 or longitude > 180 or latitude < -80 or latitude > 84 then
return nil
end
local longitude2 = longitude - floor((longitude + 180)/360) * 360
local lat_rad = deg2rad( latitude )
local e = OSGBglobe.ecc
local e_prime_sq = e / (1-e)
local v = OSGBglobe.semimajor / sqrt(1 - e * sin(lat_rad)*sin(lat_rad))
local T = pow( tan(lat_rad), 2)
local C = e_prime_sq * pow( cos(lat_rad), 2)
local A = deg2rad( longitude2 -lon0GB ) * cos(lat_rad)
local M = find_M( lat_rad )
local M0 = 0.0
if latitude_origin ~= 0 then
M0 = find_M(deg2rad( lat0GB ))
end
local northing = northingGB + scaleGB *
( (M - M0) + v*tan(lat_rad) *
( A*A/2
+ (5 - T + 9*C + 4*C*C) * pow(A,4)/24
+ (61 - 58*T + T*T
+ 600*C - 330*e_prime_sq) * pow(A,6)/720 ))
local easting = eastingGB + scaleGB * v *
( A
+ (1-T+C)*pow(A,3)/6
+ (5 - 18*T + pow(T,2) + 72*C
- 58 * e_prime_sq)*pow(A,5)/120 )
return {northing=northing,easting=easting}
end
--[[
* Convert latitude, longitude in decimal degrees to
* OSBG36 Easting and Northing
]]
local function LatLon2OSGB36( latlon, prec )
local tm = LatLonOrigin2TM(latlon.lat, latlon.lon)
if not tm then return '' end
if not tonumber(prec) then prec = 5 end
prec = floor(prec+0.5)
if prec > 5 then prec = 5 end
if prec < 1 then prec = 1 end
-- fix by Roger W Haworth
local grid_x = floor( tm.easting / 100000 )
local grid_y = floor( tm.northing / 100000 )
if grid_x < 0 or grid_x > 6 or grid_y < 0 or grid_y > 12 then return '' end
-- 0000000001111111111222222
-- 1234567890123456789012345
local letters = "ABCDEFGHJKLMNOPQRSTUVWXYZ"
local indx1 = 18-floor(grid_y/5)*5+floor(grid_x/5)
local indx2 = 21-(grid_y%5)*5+grid_x%5
local c1 = mw.ustring.sub(letters,indx1,indx1)
local c2 = mw.ustring.sub(letters,indx2,indx2)
local easting = tm.easting%100000
local northing = tm.northing%100000
local grid = pow(10.0,5.0-prec)
easting = floor(easting/grid+0.5)
northing = floor(northing/grid+0.5)
local fmt = '%0'..prec..'d'
local e = mw.ustring.format(fmt, easting)
local n = mw.ustring.format(fmt, northing)
return c1..c2..e..n
end
function oscoord._WGS2OSGB(lat,lon,prec)
return LatLon2OSGB36(HelmertDatumShift(lat,lon,WGS2OSGBparam),prec)
end
function oscoord.WGS2OSGB(frame)
local args = getArgs(frame)
return args[1] and args[2] and oscoord._WGS2OSGB(args[1],args[2],args[3]) or ''
end
function oscoord.LL2OS(frame)
local args = getArgs(frame)
if not args[1] or not args[2] then return '' end
local gridRef = oscoord._WGS2OSGB(args[1],args[2],args.prec)
if not gridRef or #gridRef == 0 then return '' end
if args[3] then
gridRef = gridRef..'_'..args[3]
end
return oscoord._oscoord({gridRef,args.linktitle,name=args.name})
end
return oscoord
