Module:Ordnance Survey coordinates

Documentation for this module may be created at Module:Ordnance Survey coordinates/doc

-- Lat Long functions in Lua 

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

local oscoord = {}
local getArgs = require('Module:Arguments').getArgs
local yesno = require('Module:Yesno')
local preview = require('Module:If preview')
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 dr = math.deg(1.0)

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=-2.0
   -- True Origin is 49 deg N
   local lambda0uk=49.0
   -- 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/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 v=6377563.396/sqrt(1.0-e2uk*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-e2uk)*v*sin(k9)
   -- 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
   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
   -- 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 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=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 {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 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(frame)
  local args = getArgs(frame)
  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&params=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(frame)
	local output = '<span class="plainlinks nourlexpansion" style="white-space: nowrap">' .. oscoord.main(frame) .. '</span>'
	if namespace == 0 then
		output = output .. '[[Category:Articles with OS grid coordinates]]'
	end
	return output
end

return oscoord