Content deleted Content added
fix bug |
start to refactor Helmert |
||
Line 52:
local rad2deg = math.deg
local dr = math.deg(1.0)
-- Transverse Mercator parameters
local scaleGB = 0.9996012717
local eastingGB = 400000.0
Line 59:
local lat0GB = 49.0
local lon0GB = -2.0
--[[ 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 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,
semimajor_dst = OSGBglobe.semimajor,
semiminor_dst = OSGBglobe.semiminor,
tx = -446.448,
ty = 125.157,
tz = -542.060,
s0 = 1.0000204894,
rx = deg2rad ( -.1502/3600. ),
ry = deg2rad ( -.2470/3600. ),
rz = deg2rad ( -.8241/3600. )
}
local function HelmertDatumShift ( latitude , longitude, param )
-- 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 a2 = param.semimajor_dst
local b2 = param.semiminor_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 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 = param.tx + ( param.s0 * x1 - param.rz * y1 + param.ry * z1 )
local y2 = param.ty + ( param.rz * x1 + param.s0 * y1 - param.rx * z1 )
local z2 = 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
aq = a2 * a2
local e2 = ( aq - b2 * b2 ) / aq
lambda = atan ( y2 / x2 )
local p2 = sqrt ( x2*x2 + y2*y2 )
phi = atan ( z2 / p2*(1.-e2) )
local limit = 1. / b2
local n0 = 11
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 ) <= limit 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)
Line 469 ⟶ 581:
]]
local function find_M ( lat_rad )
| |||