Content deleted Content added
lazy load module |
start on inverse code |
||
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
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
Line 43 ⟶ 49:
local tan = math.tan
local atan = math.atan
local deg2rad = math.rad
local rad2deg = math.deg
local dr = math.deg(1.0)
local wgs84radius = 6378137.0
local wgs84ecc = 0.00669438003551279089034150031998869922791
local scaleUK = 0.9996012717
local function northeast(lett,num,shift)
Line 85 ⟶ 96:
local lambda0uk=49.0
-- scale factor @ central meridian
local F0uk=
-- True origin in 400 km E of false origin
local E0uk=400000.0
Line 140 ⟶ 151:
local cartzb=542.06+roty*cartxa-rotx*cartya+(1.0+scale)*cartza
-- convert Cartesian to long/lat
local awgs84=
local bwgs84=6356752.3141
local e2wgs84=
local lambdaradwgs84=atan(cartyb/cartxb)
long=lambdaradwgs84*dr
Line 335 ⟶ 346:
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 408 ⟶ 419:
return output
end
--[[
* Convert latitude longitude geographical coordinates to
* Transverse Mercator coordinates.
*
* Uses the WGS-84 ellipsoid by default
*
* 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 <egil@kvaleberg.no>
*
* 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
]]
--[[
* datum shift using Helmert transformation
* parameters for transformation WGS84 --> OSGB36 hard coded
* possibility to implement other parameters should be evident
*
* 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
]]
local function HelmertDatumShift ( latitude , longitude )
-- return if latitude or longitude out of range
if abs(latitude) > 89. or abs(longitude) > 89. then
return {lat=latitude, lon=longitude}
end
-- parameters for ellipsoids
-- Annex A.1, Ordnance Survey document --
local a1 = 6378137.0
local b1 = 6356752.3141
local a2 = 6377563.396
local b2 = 6356256.910
-- parameters for transformation WGS84 --> OSGB36
-- Table 4, Section 6.6, Ordnance Survey document
local tx = -446.448
local ty = 125.157
local tz = -542.060
local s0 = 1.0000204894
local rx = deg2rad ( -.1502/3600. )
local ry = deg2rad ( -.2470/3600. )
local rz = deg2rad ( -.8241/3600. )
-- 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 aq = a1 * a1
local e1 = ( aq - b1 * b1 ) / aq
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 = tx + ( s0 * x1 - rz * y1 + ry * z1 )
local y2 = ty + ( rz * x1 + s0 * y1 - rx * z1 )
local z2 = tz + ( -ry *x1 + rx * y1 + s0 * z1 )
-- convert cartesian coordinates to latitude and longitude
-- math in Annex B.2, of Ordnance Survey document
local aq = a2 * a2
local e2 = ( aq - b2 * b2 ) / aq
lambda = atan2 ( y2 , x2 )
local p2 = sqrt ( x2*x2 + y2*y2 )
phi = atan2 ( z2 , p2*(1.-e2) )
limit = 1. / b2
n0 = 11
repeat
local phi_alt = phi
ny = a2 / sqrt ( 1. - e2 * sin(phi) * sin(phi) )
phi = atan2 ( z2 + e2 * ny * sin(phi) , p2 )
n0 = no - 1
until abs ( phi - phi_alt ) <= limit or n0 <= 0
return {lat=rad2deg(phi),lon=rad2deg(lambda)}
end
return oscoord
| |||