Module:Ordnance Survey coordinates/sandbox: Difference between revisions

Revision as of 23:25, 14 February 2021

Latest revision as of 14:18, 4 November 2024

(90 intermediate revisions by 2 users not shown)

Line 1:

--[[ 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

Ported to Lua from PHP by Wikipedia User Hike395, 18 Aug 2019

-- found by RWH at http://www.megalithia.com/search/llfuncshighlight.php

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

With thanks to Andy, G4JNT for inspiration in GEOG, and to the OSGB for their white paper on coordinate transformation

-- You may use and redistribute this code under the terms of the GPL see http://www.gnu.org/copyleft/gpl.html

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

Written by Richard

-- This code is predicated on the assumption that your are ONLY feeding it Irish or UK Grid references.

www.megalithia.com

-- It uses the single char prefix of Irish grid refs to tell the difference, UK grid refs have a two letter prefix.

v0.something 27/2/2000

-- We would like an even number of digits for the rest of the grid ref.

v 1.01 28 June 2004

-- Anything in the NGR other than 0-9, A-Z, a-z is eliminated.

v 1.02 6 Aug 2004 line 89 add "0" to chars in ngr=stripcharsnotinbag Thx Andy

-- WARNING this assumes there are no decimal points in your NGR components.

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.

-- The transformation from OSGB36/Ireland 1965 to WGS84 is more precise than 5 m.

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

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 namespace = mw.title.getCurrentTitle().namespace

local pow = math.pow

local sqrt = math.sqrt

local abs = math.abs

local floor = math.floor

local ceil = math.ceil

local sin = math.sin

local cos = math.cos

local tan = math.tan

local atan = math.atan

local deg2rad = math.rad

local rad2deg = math.deg

--[[ 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 data = mw.loadData('Module:Ordnance Survey coordinates/data')

local OSGBglobe = data.OSGBglobe

local IEglobe = data.IEglobe

local WGSglobe = data.WGSglobe

local WGS2OSGBparam = data.WGS2OSGBparam

local OSGB2WGSparam = data.OSGB2WGSparam

local IE2WGSparam = data.IE2WGSparam

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 )

local x1 = ny * cosphi * coslambda

local y1 = ny * cosphi * sinlambda

local z1 = ny * ( 1. - e1 ) * sinphi

-- the helmert transformation proper

-- Equation 3, Section 6.2, Ordnance Survey document

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 )

-- 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=stringmw.ustring.len(num)

if le%2 == 1 then

return {err="Malformed numerical part of NGR"}

end

local pr=le/2

local n = stringmw.ustring.sub(num,pr+1)

local e = stringmw.ustring.sub(num,1,pr)

-- Hack to move to center of square: append a 5 to northings and eastings

if shift ~= nil and yesno(shift > 0 then)

if shift then

n = n.."5"

e = e.."5"

Line 50 ⟶ 162:

end

-- end hack

ifn = string.len(n) == '' and 0 thenor n

e = ne == '' and 0 or e

pr = pow(10.0,(5.0-pr))

end

local T1 = mw.ustring.byte(mw.ustring.sub(lett,1,1))-65

if string.len(e) == 0 then

e = 0

end

pr = math.pow(10.0,(5.0-pr))

local T1 = string.byte(string.sub(lett,1,1))-65

if T1>8 then

T1 = T1-1

end

local T2 = nil

if stringmw.ustring.len(lett)>1 then

T2 = stringmw.ustring.byte(stringmw.ustring.sub(lett,2))-65

if T2>8 then

T2 = T2-1

end

end

return nil,{n=n,e=e,pr=pr,T1=T1,T2=T2}

end

local function GBEN2LLEN2LL(e,n,datum)

local pow,sqrt,absa =math datum.pow,mathsemimajor*datum.sqrt,math.absscale

local b = datum.semiminor*datum.scale

local sin,cos,tan,atan=math.sin,math.cos,math.tan,math.atan

local drlat0rad = math.degdeg2rad(1datum.0lat0)

local n1 = datum.n1

-- True Origin is 2 deg W

local phi0ukn12 =-2.0 n1*n1

local n13 = n12*n1

-- True Origin is 49 deg N

local lambda0ukk=49(n-datum.0n0)/a+lat0rad

-- scale factor @ central meridian

local F0uk=0.9996012717

-- True origin in 400 km E of false origin

local E0uk=400000.0

--True origin is 100 km S of false origin

local N0uk=-100000.0

-- semi-major axis (in line to equator) 0.996012717 is yer scale @ central meridian

local auk=6377563.396*F0uk

--semi-minor axis (in line to poles)

local buk=6356256.91*F0uk

-- flatness=a1-b1/(a1+b1)

local n1uk=0.00167322025032508731869331280635710896296

-- first eccentricity squared=2*f-f^2where f=(a1-b1)/a1

local e2uk=0.006670539761597529073698869358812557054558

local k=(n-N0uk)/auk+lambda0uk/dr

local nextcounter=0

local j3, j4, j5, j6, m

repeat

nextcounter=nextcounter+1

local k3=k-lambda0uk/drlat0rad

local k4=k+lambda0uk/drlat0rad

j3=(1.0+n1ukn1+1.25*pow(n1uk,2.0)n12+1.25*pow(n1uk,3.0)n13)*k3

j4=(3.0*n1ukn1+3.0*pow(n1uk,2.0)n12+2.625*pow(n1uk,3.0)n13)*sin(k3)*cos(k4)

j5=(1.875*pow(n1uk,2.0)n12+1.875*pow(n1uk,3.0)n13)*sin(2.0*k3)*cos(2.0*k4)

j6=35.0/24.0*pow(n1uk,3.0)n13*sin(3.0*k3)*cos(3.0*k4)

m=bukb*(j3-j4+j5-j6)

k=k+(n-N0ukdatum.n0-m)/auka

until abs(n-N0ukdatum.n0-m)<=0.000000001 or nextcounter>=100

local vx =auk/sqrt( 1.0-e2ukdatum.ecc*pow(sin(k),2.0))

local rv=v*(1.0-e2uk)a/sqrt(1.0-e2uk*pow(sin(k),2.0)x)

local r=v*(1.0-datum.ecc)/x

local h2=v/r-1.0

local y1=e-E0ukdatum.e0

j3local tank = tan(k)/(2.0*r*v)

local tank2 = tank*tank

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)

local tank4 = tank2*tank2

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 tank6 = tank2*tank4

local k9=k-y1*y1*j3+pow(y1,4.0)*j4-pow(y1,6.0)*j5

j6local cosk =1.0/( cos(k)*v)

local yv = y1/v -- dimensionless quantity in series expansion

local j7=1.0/(cos(k)*6.0*pow(v,3.0))*(v/r+2.0*pow(tan(k),2.0))

local yv2 = yv*yv

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 j9yv3 =1.0/(cos(k) yv*5040.0*pow(v,7.0))yv2

local yv5 = yv3*yv2

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 yv7 = yv5*yv2

local long=phi0uk+dr*(y1*j6-y1*y1*y1*j7+pow(y1,5.0)*j8-pow(y1,7.0)*j9)

j3=tank/(2.0*r)

local lat=k9*dr

j4=tank/(24.0*r)*(5.0+3.0*tank2+h2-9.0*tank2*h2)

local v=6377563.396/sqrt(1.0-e2uk*pow(sin(k),2.0))

j5=tank/(720.0*r)*(61.0+90.0*tank2+45.0*tank4)

local cartxa=v*cos(k9)*cos(long/dr)

local cartyak9=vk-y1*cos(k9)yv*sin(long/drj3+yv3*j4-yv5*j5)

local cartzaj6=(1.0-e2uk)*v*sin/(k9cosk)

local j7=1.0/(cosk*6.0)*(v/r+2.0*tank2)

-- Helmert-Transformation from OSGB36 to WGS84 map date

local rotxj8=-01.15020/3600(cosk*120.0/dr)*(5.0+28.0*tank2+24.0*tank4)

local rotyj9=-1.0.2470/3600(cosk*5040.0/dr)

j9=j9*(61.0+662.0*tank2+1320.0*tank4+720.0*tank6)

local rotz=-0.8421/3600.0/dr

local long=datum.lon0+rad2deg(yv*j6-yv3*j7+yv5*j8-yv7*j9)

local scale=-20.4894/1000000.0

local lat=rad2deg(k9)

local cartxb=446.448+(1.0+scale)*cartxa+rotz*cartya-roty*cartza

return {lat=lat,lon=long}

local cartyb=-125.157-rotz*cartxa+(1.0+scale)*cartya+rotx*cartza

end

local cartzb=542.06+roty*cartxa-rotx*cartya+(1.0+scale)*cartza

-- convert Cartesian to long/lat

local function GBEN2LL(e,n)

local awgs84=6378137.0

local bwgs84latlong =6356752.3141 EN2LL(e,n,OSGBglobe)

local helmert = HelmertDatumShift ( latlong.lat, latlong.lon, OSGB2WGSparam)

local e2wgs84=0.00669438003551279089034150031998869922791

return {region="GB",lat=helmert.lat,long=helmert.lon}

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 "GB",nil,lat,long

end

Line 159 ⟶ 240:

-- first caclulate e,n

-- computing e and n exactly, to get SW corner of box

local err, n, e, pr, T1, T2ne = northeast(lett,num,0)

if ne.err ~= nil then

return {region="GB",err,0=ne.0,0.0err}

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*math.floor(ne.T1/5)-100000.0*math.floor(ne.T2/5)+ne.n*ne.pr

returnlocal result = GBEN2LL(e,n)

result.prec = 0.8165*ne.pr

return result

end

local function IrishEN2LL(e,n)

local pow,sqrt,abslatlong =math.pow EN2LL(e,math.sqrtn,math.absIEglobe)

local helmert = HelmertDatumShift ( latlong.lat, latlong.lon, IE2WGSparam)

local sin,cos,tan,atan=math.sin,math.cos,math.tan,math.atan

return {region="IE",lat=helmert.lat,long=helmert.lon}

local dr=math.deg(1.0)

-- 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=6378137.0

local bwgs84=6356752.3141

local e2wgs84=0.00669438003551279089034150031998869922791

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 "IE",nil,lat,long

end

Line 257 ⟶ 262:

-- first caclulate e,n

-- computing e and n exactly, to get SW corner of box

local err, n, e, pr, T1ne = northeast(lett,num,0)

if ne.err ~= nil then

return {region="IE", err,0=ne.0,0.0err}

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-math.floor(ne.T1/5.0))

returnlocal 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 = {}

-- returns a country,error,lat,long list

local _ = 0

ngr = string.gsub(string.upper(ngr),"[^%d%u]","")

localngr, lett_ = stringmw.ustring.gsub(mw.ustring.upper(ngr),"[^%us%p]","")

local first, last, lett, num = stringmw.gsubustring.find(ngr,"[^([A-Z]+)(%d]","+)$")

if stringnot first or empty(lett) or empty(num) or mw.ustring.len(lett) ==> 12 then

return {err="Malformed NGR"}

return Irish2LL(lett,num)

end

if stringmw.ustring.len(lett) == 21 then

return GB2LLIrish2LL(lett,num)

end

return nilGB2LL(lett,"Malformed NGR",0.0,0.0num)

end

local function split(s,sep)

-- split a string s into chunks, separated by sep

if sep == nilsep or then"%s"

sep = "%s"

end

local t = {}

for chunk in stringmw.ustring.gmatch(s,"([^"..sep.."]+)") do

table.insert(t,chunk)

end

Line 294 ⟶ 303:

local function trim(s)

local _ = 0

s = string.gsub(s,"^%s+","")

s, _ = stringmw.ustring.gsub(s,"^%s+$","")

s, _ = mw.ustring.gsub(s,"%s+$","")

return s

end

local function alldigits(s)

return stringnot mw.ustring.find(s,"%D") == nil

end

local function warning(errmsg)

local framepreview = mw.getCurrentFramerequire()'Module:getParent(If preview')

local htmlmsg = ""errmsg or 'Empty OS grid ref'

local msg = errmsg

local html = preview._warning({ msg })

if errmsg == nil then

msg = "Empty OS grid ref"

if namespace == 0 and errmsg then

end

if frame:preprocess( "{{REVISIONID}}" ) == "" then

html = '<div style="color:red"><strong>Warning:</strong> '..msg

html = html..' <small>(this message is shown only in preview)</small></div>'

end

if namespace == 0 and errmsg ~= nil then

html = html..'[[Category:Pages with malformed OS coordinates]]'

end

Line 320 ⟶ 325:

end

-- generate URL of geohack map

function oscoord.main(frame)

local function geohack(main_args, other_args, LL)

local args = getArgs(frame,{parentFirst=true,parentOnly=false,frameOnly=false})

-- create geohack link. Example:

local input = args[1]

-- https://geohack.toolforge.org/geohack.php?pagename=Mount_Whitney&params=36.578580925_N_118.29199495_W_type:mountain_region:US-CA_scale:100000_source:NGS

if input == nil or string.len(input)==0 then

local url = 'https://geohack.toolforge.org/geohack.php?'

local input = main_args[1]

local namearg = main_args["name"]

local current_page = mw.title.getCurrentTitle()

local pagename = mw.uri.encode( current_page.prefixedText, 'WIKI' )

if not empty(pagename) then

url = url..'pagename='..pagename..'&'

end

url = url..'params='..mw.ustring.format('%.6f',LL.lat)..'_N_'

if LL.long < 0 then

url = url..mw.ustring.format('%.6f',-LL.long)..'_W'

else

url = url..mw.ustring.format('%.6f',LL.long)..'_E'

end

for _, a in ipairs(other_args) do

url = url..'_'..a

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

return url

end

-- function to generate direct link to OS map

local function directLink(main_args, other_args, LL)

-- create link to Bing server for OS maps. Example:

-- https://www.bing.com/maps/?mkt=en-gb&v=2&cp=56.796026%7E-5.01307&lvl=16.0&sp=Point.56.796029_-5.004711_Ben+Nevis&sty=s&style=s

local current_page = mw.title.getCurrentTitle()

local namearg = mw.uri.encode( main_args["name"] or current_page.prefixedText, 'QUERY' )

local args = {}

for _, a in ipairs(other_args) do

local splitOut = mw.text.split(a, ':', true)

args[splitOut[1]] = splitOut[2]

end

if not args.scale and not args.type and not args.dim then

args.dim = LL.prec and tostring(floor(50*LL.prec+0.5))..'m'

end

args.viewport_cm = 10

local zoom = require('Module:Infobox dim')._zoom

local lvl = zoom(args) or 12

local url = mw.ustring.format('https://www.bing.com/maps/?mkt=en-gb&v=2&cp=%.6f~%.6f&lvl=%d&sp=Point.%.6f_%.6f',LL.lat,LL.long,lvl,LL.lat,LL.long)

if not empty(namearg) then

url = url..'_'..namearg

end

url = url..'&sty=s&style=s'

return url

end

function oscoord._main(main_args)

local input = main_args[1]

if empty(input) then

return warning(nil)

end

local linktitle = argsmain_args[2]

local nameargrawurl = argsyesno(main_args["namerawurl"])

local direct = yesno(main_args["direct"])

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 stringmw.ustring.sub(args[1],1,1) == 'i' then

local firstArg = stringmw.ustring.sub(args[1],2)

if alldigits(firstArg) then

LL = {IrishEN2LL(firstArg,args[2])}

restargs = 3

if linktitle == nil or string.lenempty(linktitle) == 0 then

linktitle=args[1]..'_'..args[2]

end

end

elseif alldigits(args[1]) then

LL = {GBEN2LL(args[1],args[2])}

restargs = 3

if linktitle == nil or string.lenempty(linktitle) == 0 then

linktitle=args[1]..'_'..args[2]

end

end

else

LL = {NGR2LL(args[1])}

restargs = 2

if linktitle == nil or string.lenempty(linktitle) == 0 then

linktitle=args[1]

end

end

linktitle = trim(linktitle)

if not empty(LL[2] ~= nil.err) then

return linktitle ..warning(LL.err)

local html = linktitle

html = html..warning(LL[2])

return html

end

LL.lat = LL.lat or 0

-- https://geohack.toolforge.org/geohack.php?pagename=Mount_Whitney&params=36.578580925_N_118.29199495_W_type:mountain_region:US-CA_scale:100000_source:NGS

LL.long = LL.long or 0

local url = ''

LL.lat = ceil(LL.lat*1000000)/1000000

if not args['rawurl'] then

LL.long = ceil(LL.long*1000000)/1000000

url = url..'['

local other_args = {}

for i = restargs, #args do

table.insert(other_args, args[i])

end

local url

url = url..'https://geohack.toolforge.org/geohack.php?'

if pagenamenot ~= nildirect then

url = geohack(main_args, other_args, LL)

url = url..'pagename='..pagename..'&'

end

url = url..'params='..LL[3]..'_N_'

if LL[4] < 0 then

url = url..(-LL[4])..'_W'

else

url = directLink(main_args, other_args, LL)

url = url..LL[4]..'_E'

end

if not rawurl then

for i = restargs,#args do

url = '['..url..'_ '..args[ilinktitle..']'

end

if string.find(input,"region") == nil then

url = url..'_region:'..LL[1]

end

if namearg then

url = url .. "&title=" .. mw.uri.encode(namearg)

end

if not args['rawurl'] then

url = url..' '..linktitle..']'

end

return url

end

function oscoord.oscoord_oscoord(frameargs)

local output = '<span class="plainlinks nourlexpansion" style="white-space: nowrap">' .. oscoord.main_main(frameargs) .. '</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 <[e-mail]>

*

* 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

local e2 = e*e

local e3 = e*e*e

return OSGBglobe.semimajor * ( ( 1 - e/4 - 3 * e2/64

- 5 * e3/256

) * lat_rad

- ( 3 * e/8 + 3 * e2/32

+ 45 * e3/1024

) * sin(2 * lat_rad)

+ ( 15 * e2/256 +

45 * e3/1024

) * sin(4 * lat_rad)

- ( 35 * e3/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 tank = tan(lat_rad)

local T = tank*tank

local T2 = T*T

local C = e_prime_sq * pow( cos(lat_rad), 2)

local A = deg2rad( longitude2 -OSGBglobe.lon0 ) * cos(lat_rad)

local A2 = A*A

local A3 = A2*A

local A4 = A2*A2

local A5 = A3*A2

local A6 = A3*A3

local M = find_M( lat_rad )

local M0 = 0.0

if latitude_origin ~= 0 then

M0 = find_M(deg2rad( OSGBglobe.lat0 ))

end

local northing = OSGBglobe.n0 + OSGBglobe.scale *

( (M - M0) + v*tan(lat_rad) *

( A2/2

+ (5 - T + 9*C + 4*C*C) * A4/24

+ (61 - 58*T + T2

+ 600*C - 330*e_prime_sq) * A6/720 ))

local easting = OSGBglobe.e0 + OSGBglobe.scale * v *

( A

+ (1-T+C)*A3/6

+ (5 - 18*T + T2 + 72*C

- 58 * e_prime_sq)*A5/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)

northing = floor(northing/grid)

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