Module:Ordnance Survey coordinates/sandbox: Difference between revisions

Revision as of 18:53, 16 December 2023

Latest revision as of 14:18, 4 November 2024

(39 intermediate revisions by the same user not shown)

Line 45:

local abs = math.abs

local floor = math.floor

local ceil = math.ceil

local sin = math.sin

local cos = math.cos

Line 51 ⟶ 52:

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

Line 77 ⟶ 71:

-- Datum parameters

local data = mw.loadData('Module:Ordnance Survey coordinates/data')

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

Line 222 ⟶ 179:

end

local function GBEN2LLEN2LL(e,n,datum)

local a = datum.semimajor*datum.scale

-- True Origin is 2 deg W

local phi0ukb =lon0GB datum.semiminor*datum.scale

local lat0rad = deg2rad(datum.lat0)

-- True Origin is 49 deg N

local lambda0ukn1 =lat0GB datum.n1

local n12 = n1*n1

-- scale factor @ central meridian

local F0ukn13 =scaleGB n12*n1

local k=(n-datum.n0)/a+lat0rad

-- 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/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 helmert = HelmertDatumShift ( lat, long, OSGB2WGSparam)

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

local k9=k-y1*(yv*j3+yv3*j4-yv5*j5)

j6=1.0/(cosk)

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

local j8=1.0/(cosk*120.0)*(5.0+28.0*tank2+24.0*tank4)

local j9=1.0/(cosk*5040.0)

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

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

local lat=rad2deg(k9)

return {lat=lat,lon=long}

end

local function GBEN2LL(e,n)

local latlong = EN2LL(e,n,OSGBglobe)

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

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

end

Line 292 ⟶ 253:

local function IrishEN2LL(e,n)

local latlong = EN2LL(e,n,IEglobe)

-- True Origin is 8 deg W

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

local phi0ir=-8.0

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

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

Line 394 ⟶ 280:

local function NGR2LL(ngr)

local result = {}

local _ = 0

ngr, _ = mw.ustring.gsub(mw.ustring.upper(ngr),"[%s%p]","")

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

Line 416 ⟶ 303:

local function trim(s)

local _ = 0

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

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

Line 437 ⟶ 325:

end

-- generate URL of geohack map

function oscoord._main(args)

local function geohack(main_args, other_args, LL)

local input = args[1]

-- create geohack link. Example:

-- 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 = '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 rawurldirect = yesno(argsmain_args["rawurldirect"])

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

Line 477 ⟶ 423:

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

urlLL.lat = url..'params='..ceil(LL.lat..'_N_'*1000000)/1000000

if LL.long <= 0 thenceil(LL.long*1000000)/1000000

local other_args = {}

url = url..(-LL.long)..'_W'

for i = restargs, #args do

table.insert(other_args, args[i])

end

local url

if not direct then

url = geohack(main_args, other_args, LL)

else

url = directLink(main_args, other_args, LL)

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

Line 570 ⟶ 499:

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 * e * ee2/64

- 5 * e * e * ee3/256

) * lat_rad

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

+ 45 * e * e * ee3/1024

) * sin(2 * lat_rad)

+ ( 15 * e * ee2/256 +

45 * e * e * ee3/1024

) * sin(4 * lat_rad)

- ( 35 * e * e * ee3/3072

) * sin(6 * lat_rad) )

end

Line 603 ⟶ 534:

local v = OSGBglobe.semimajor / sqrt(1 - e * sin(lat_rad)*sin(lat_rad))

local Ttank = pow( tan(lat_rad), 2)

local T = tank*tank

local T2 = T*T

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

local A = deg2rad( longitude2 -lon0GBOSGBglobe.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( lat0GBOSGBglobe.lat0 ))

end

local northing = northingGBOSGBglobe.n0 + scaleGBOSGBglobe.scale *

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

( A*AA2/2

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

+ (61 - 58*T + T*T T2

+ 600*C - 330*e_prime_sq) * pow(A,6)A6/720 ))

local easting = eastingGBOSGBglobe.e0 + scaleGBOSGBglobe.scale * v *

( A

+ (1-T+C)*pow(A,3)A3/6

+ (5 - 18*T + pow(T,2)T2 + 72*C

- 58 * e_prime_sq)*pow(A,5)A5/120 )

return {northing=northing,easting=easting}

Line 656 ⟶ 594:

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)