Module:Ordnance Survey coordinates/sandbox: Difference between revisions

Revision as of 22:33, 16 December 2023

Latest revision as of 14:18, 4 November 2024

(30 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 70 ⟶ 71:

-- Datum parameters

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

local OSGBglobe = {

semimajor = 6377563.396,

semiminor = 6356256.910,

ecc = 0.006670539761597529073698869358812557054558,

n1 = 0.00167322025032508731869331280635710896296,

scale = 0.9996012717,

e0 = 400000.0,

n0 = -100000.0,

lat0 = 49.0,

lon0 = -2.0

}

local IEglobeOSGBglobe = {data.OSGBglobe

local IEglobe = data.IEglobe

semimajor = 6377340.189,

local WGSglobe = data.WGSglobe

semiminor = 6356034.447,

local WGS2OSGBparam = data.WGS2OSGBparam

ecc = 0.006670540293336110419293763349975612794125,

local OSGB2WGSparam = data.OSGB2WGSparam

n1 = 0.001673220384152058651484728058385228837777,

local IE2WGSparam = data.IE2WGSparam

scale = 1.000035,

e0 = 200000.0,

n0 = 250000.0,

lat0 = 53.5,

lon0 = -8.0

}

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 IE2WGSparam = {

semimajor_src = IEglobe.semimajor,

semiminor_src = IEglobe.semiminor,

ecc_src = IEglobe.ecc,

semimajor_dst = WGSglobe.semimajor,

semiminor_dst = WGSglobe.semiminor,

ecc_dst = WGSglobe.ecc,

tx = 482.53,

ty = -130.596,

tz = 564.557,

s0 = 8.15/1000000.0,

rx = deg2rad( -1.042/3600.0 ),

ry = deg2rad( -0.214/3600.0 ),

rz = deg2rad( -0.631/3600.0 )

}

local function HelmertDatumShift ( latitude , longitude, param )

Line 252 ⟶ 183:

local b = datum.semiminor*datum.scale

local lat0rad = deg2rad(datum.lat0)

local n1 = datum.n1

local n12 = n1*n1

local n13 = n12*n1

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

local nextcounter=0

Line 259 ⟶ 193:

local k3=k-lat0rad

local k4=k+lat0rad

j3=(1.0+datum.n1+1.25*pow(datum.n1,2.0)n12+1.25*pow(datum.n1,3.0)n13)*k3

j4=(3.0*datum.n1+3.0*pow(datum.n1,2.0)n12+2.625*pow(datum.n1,3.0)n13)*sin(k3)*cos(k4)

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

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

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

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

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

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

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

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

local h2=v/r-1.0

local y1=e-datum.e0

Line 275 ⟶ 210:

local tank6 = tank2*tank4

local cosk = cos(k)

local v3yv = y1/v*v*v -- dimensionless quantity in series expansion

local v5yv2 = v3yv*v*vyv

local y13yv3 = y1yv*y1*y1yv2

local y14yv5 = y13yv3*y1yv2

local y15yv7 = y14yv5*y1yv2

j3=tank/(2.0*r)

local y16 = y13*y13

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

local y17 = y16*y1

j3j5=tank/(2720.0*r)*v(61.0+90.0*tank2+45.0*tank4)

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

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

j6=1.0/(cosk)

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

local k9j7=k-y11.0/(cosk*y16.0)*j3(v/r+y142.0*j4-y16*j5tank2)

j6local j8=1.0/(cosk*v120.0)*(5.0+28.0*tank2+24.0*tank4)

local j7j9=1.0/(cosk*65040.0*v3)*(v/r+2.0*tank2)

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

local j9=1.0/(cos(k)*5040.0*pow(v,7.0))

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

local long=datum.lon0+rad2deg(y1yv*j6-y13yv3*j7+y15yv5*j8-y17yv7*j9)

local lat=rad2deg(k9)

return {lat=lat,lon=long}

Line 347 ⟶ 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 369 ⟶ 303:

local function trim(s)

local _ = 0

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

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

Line 390 ⟶ 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 430 ⟶ 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 523 ⟶ 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 556 ⟶ 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 -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

Line 567 ⟶ 552:

local northing = OSGBglobe.n0 + OSGBglobe.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 = OSGBglobe.e0 + OSGBglobe.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 609 ⟶ 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)