Skip to content

Instantly share code, notes, and snippets.

@sancarn

sancarn/CoordinateConversion.m

Created March 1, 2024 21:18
Show Gist options
  • Save sancarn/17fdd8358ca2540b534cf9cfacc76f01 to your computer and use it in GitHub Desktop.
Save sancarn/17fdd8358ca2540b534cf9cfacc76f01 to your computer and use it in GitHub Desktop.
Download ZIP
PowerQuery M language for converting from Eastings Northings to Latitude Longitude (and back)
Raw
CoordinateConversion.m
let
//Computed Helmert transformed X coordinate.
//@param {number} X - Cartesian X coord in meters
//@param {number} Y - Cartesian Y coord in meters
//@param {number} Z - Cartesian Z coord in meters
//@param {number} DX - Cartesian X translation in meters
//@param {number} Y_Rot - Y rotation in seconds of arc
//@param {number} Z_Rot - Z rotation in seconds of arc
//@param {number} s - Scale in ppm
//@returns {number} Helmert transformed X coordinate
Helmert_X = (x as number, y as number, z as number, DX as number, Y_Rot as number, Z_Rot as number, s as number) as number =>
(
let
//Convert rotations to radians and ppm scale to a factor
sfactor = s * 0.000001,
RadY_Rot = (Y_Rot / 3600) * Number.PI / 180,
RadZ_Rot = (Z_Rot / 3600) * Number.PI / 180,
//Compute transformed X coordinate
result = x + (x * sfactor) - (y * RadZ_Rot) + (z * RadY_Rot) + DX
in
result
),
//Computed Helmert transformed Y coordinate.
//@param {number} X - Cartesian X coord in meters
//@param {number} Y - Cartesian Y coord in meters
//@param {number} Z - Cartesian Z coord in meters
//@param {number} DY - Cartesian Y translation in meters
//@param {number} X_Rot - X rotation in seconds of arc
//@param {number} Z_Rot - Z rotation in seconds of arc
//@param {number} s - Scale in ppm
//@returns {number} Helmert transformed Y coordinate
Helmert_Y = (x as number, y as number, z as number, DY as number, X_Rot as number, Z_Rot as number, s as number) as number =>
(
let
//Convert rotations to radians and ppm scale to a factor
sfactor = s * 0.000001,
RadX_Rot = (X_Rot / 3600) * Number.PI / 180,
RadZ_Rot = (Z_Rot / 3600) * Number.PI / 180,
//Compute transformed Y coordinate
result = (x * RadZ_Rot) + y + (y * sfactor) - (z * RadX_Rot) + DY
in
result
),
//Computed Helmert transformed Z coordinate.
//@param {number} X - Cartesian X coord in meters
//@param {number} Y - Cartesian Y coord in meters
//@param {number} Z - Cartesian Z coord in meters
//@param {number} DZ - Cartesian Z translation in meters
//@param {number} X_Rot - X rotation in seconds of arc
//@param {number} Y_Rot - Y rotation in seconds of arc
//@param {number} s - Scale in ppm
//@returns {number} Helmert transformed Z coordinate
Helmert_Z = (x as number, y as number, z as number, DZ as number, X_Rot as number, Y_Rot as number, s as number) as number =>
(
let
//Convert rotations to radians and ppm scale to a factor
sfactor = s * 0.000001,
RadX_Rot = (X_Rot / 3600) * Number.PI / 180,
RadY_Rot = (Y_Rot / 3600) * Number.PI / 180,
//Compute transformed Z coordinate
result = - (x * RadY_Rot) + (y * RadX_Rot) + z + (z * sfactor) + DZ
in
result
),
//Convert XYZ to Latitude (PHI) in Dec Degrees.
//@param {number} X - Cartesian X coord in meters
//@param {number} Y - Cartesian Y coord in meters
//@param {number} Z - Cartesian Z coord in meters
//@param {number} a - Ellipsoid semi-major axis in meters
//@param {number} b - Ellipsoid semi-minor axis in meters
//@returns {number} Latitude in decimal (TODO:degrees or radians?)
//@dependencies [Iterate_XYZ_to_Lat]
XYZ_to_Lat = (x as number, y as number, z as number, a as number, b as number) as number =>
(
let
RootXYSqr = Number.Sqrt(x * x + y * y),
e2 = (a * a - b * b) / (a * a),
PHI1 = Number.Atan(z / (RootXYSqr * (1 - e2))),
PHI = Iterate_XYZ_to_Lat(a, e2, PHI1, z, RootXYSqr)
in
PHI * 180 / Number.PI
),
//Iteratively computes Latitude (PHI).
//@param {number} a - Elipsoid semi-major axis in meters
//@param {number} e2 - eta squared
//@param {number} PHI1 - estimated value for latitude in radians
//@param {number} Z - Cartesian coordinate in meters
//@param {number} RootXYSqr - Computed from X and Y in meters
//@returns {number} Latitude in decimal (TODO: degrees or radians?)
Iterate_XYZ_to_Lat = (a as number, e2 as number, PHI1 as number, Z as number, RootXYSqr as number) as number =>
(
let
V = a / (Number.Sqrt(1 - e2 * Number.Sin(PHI1) * Number.Sin(PHI1))),
PHI2 = Number.Atan((Z + e2 * V * Number.Sin(PHI1)) / RootXYSqr),
iterate = (PHI1) =>
(
let
V = a / (Number.Sqrt(1 - e2 * Number.Sin(PHI1) * Number.Sin(PHI1))),
PHI2 = Number.Atan((Z + e2 * V * Number.Sin(PHI1)) / RootXYSqr),
result = if PHI1 - PHI2 > 0.000000001 then @iterate(PHI2) else PHI2
in
result
)
in
iterate(PHI2)
),
//Convert XYZ to Longitude (LAM) in Decimal Degrees.
//@param {number} X - Cartesian X coord in meters
//@param {number} Y - Cartesian Y coord in meters
//@returns {number} Longitude in decimal degrees
XYZ_to_Long = (x as number, y as number) as number =>
(
//tailor the output to fit the equatorial quadrant as determined by the signs of X and Y
let
result =
if x >= 0 then
//longitude is in the W90 thru 0 to E90 hemisphere
Number.Atan(y / x) * 180 / Number.PI
else if x < 0 and y >= 0 then
//longitude is in the E90 to E180 quadrant
Number.Atan(y / x) * 180 / Number.PI + 180
else
//longitude is in the E180 to W90 quadrant
Number.Atan(y / x) * 180 / Number.PI - 180
in
result
),
//Convert XYZ to Ellipsoidal Height (H).
//@param {number} X - Cartesian X coord in meters
//@param {number} Y - Cartesian Y coord in meters
//@param {number} Z - Cartesian Z coord in meters
//@param {number} a - Ellipsoid semi-major axis in meters
//@param {number} b - Ellipsoid semi-minor axis in meters
//@returns {number}
//@dependencies [XYZ_to_Lat()]
XYZ_to_H = (x as number, y as number, z as number, a as number, b as number) as number =>
(
let
//Compute PHI (Dec Degrees) first
phi = XYZ_to_Lat(x, y, z, a, b),
//Convert PHI to radians
RadPHI = phi * Number.PI / 180,
//Compute Elipsoidal Height H
RootXYSqr = Number.Sqrt(x * x + y * y),
e2 = (a * a - b * b) / (a * a),
V = a / (Number.Sqrt(1 - e2 * Number.Sin(RadPHI) * Number.Sin(RadPHI))),
H = (RootXYSqr / Number.Cos(RadPHI)) - V
in
H
),
//Convert geodetic coords lat (PHI), long (LAM) and height (H) to cartesian X coordinate.
//@param {number} PHI - Latitude in decimal degrees
//@param {number} LAM - Longitude in decimal degrees
//@param {number} H - Ellipsoid height in meters
//@param {number} a - Ellipsoid semi-major axis in meters
//@param {number} b - Ellipsoid semi-minor axis in meters
//@returns {number} - Cartesian X coord in meters
Lat_Long_H_to_X = (phi as number, lam as number, H as number, a as number, b as number) as number =>
(
let
//Convert angle measures to radians
RadPHI = phi * Number.PI / 180,
RadLAM = lam * Number.PI / 180,
//Compute eccentricity squared and nu
e2 = (a * a - b * b) / (a * a),
V = a / (Number.Sqrt(1 - e2 * Number.Sin(RadPHI) * Number.Sin(RadPHI))),
//Compute X
X = (V + H) * Number.Cos(RadPHI) * Number.Cos(RadLAM)
in
X
),
//Convert geodetic coords lat (PHI), long (LAM) and height (H) to cartesian Y coordinate.
//@param {number} PHI - Latitude in decimal degrees
//@param {number} LAM - Longitude in decimal degrees
//@param {number} H - Ellipsoid height in meters
//@param {number} a - Ellipsoid semi-major axis in meters
//@param {number} b - Ellipsoid semi-minor axis in meters
//@returns {number} - Cartesian Y coord in meters
Lat_Long_H_to_Y = (phi as number, lam as number, H as number, a as number, b as number) as number =>
(
let
//Convert angle measures to radians
RadPHI = phi * Number.PI / 180,
RadLAM = lam * Number.PI / 180,
//Compute eccentricity squared and nu
e2 = (a * a - b * b) / (a * a),
V = a / (Number.Sqrt(1 - e2 * Number.Sin(RadPHI) * Number.Sin(RadPHI))),
//Compute Y
Y = (V + H) * Number.Cos(RadPHI) * Number.Sin(RadLAM)
in
Y
),
//Convert geodetic coords lat (PHI) and height (H) to cartesian Z coordinate.
//@param {number} PHI - Latitude in decimal degrees
//@param {number} H - Ellipsoid height in meters
//@param {number} a - Ellipsoid semi-major axis in meters
//@param {number} b - Ellipsoid semi-minor axis in meters
//@returns {number} - Cartesian Z coord in meters
Lat_H_to_Z = (phi as number, H as number, a as number, b as number) as number =>
(
let
//Convert angle measures to radians
RadPHI = phi * Number.PI / 180,
//Compute eccentricity squared and nu
e2 = (a * a - b * b) / (a * a),
V = a / (Number.Sqrt(1 - e2 * Number.Sin(RadPHI) * Number.Sin(RadPHI))),
//Compute Z
Z = ((V * (1 - e2)) + H) * Number.Sin(RadPHI)
in
Z
),
//Project Latitude and longitude to Transverse Mercator Eastings.
//@param {number} PHI - Latitude in decimal degrees
//@param {number} LAM - Longitude in decimal degrees
//@param {number} a - Ellipsoid semi-major axis in meters
//@param {number} b - Ellipsoid semi-minor axis in meters
//@param {number} e0 - Eastings of false origin in meters
//@param {number} f0 - Central meridian scale factor
//@param {number} PHI0 - Latitude of false origin in decimal degrees
//@param {number} LAM0 - Longitude of false origin in decimal degrees
Lat_Long_to_East = (
phi as number, lam as number, a as number, b as number, E0 as number, F0 as number, phi0 as number, lam0 as number
) as number =>
(
let
//Convert angle measures to radians
RadPhi = phi * Number.PI / 180,
RadLam = lam * Number.PI / 180,
RadPhi0 = phi0 * Number.PI / 180,
RadLam0 = lam0 * Number.PI / 180,
af0 = a * F0,
bf0 = b * F0,
e2 = (af0 * af0 - bf0 * bf0) / (af0 * af0),
n = (af0 - bf0) / (af0 + bf0),
nu = af0 / (Number.Sqrt(1 - e2 * Number.Sin(RadPhi) * Number.Sin(RadPhi))),
rho = (nu * (1 - e2)) / (1 - e2 * Number.Sin(RadPhi) * Number.Sin(RadPhi)),
eta2 = (nu / rho) - 1,
p = RadLam - RadLam0,
IV = nu * Number.Cos(RadPhi),
V = (nu / 6) * Number.Power(Number.Cos(RadPhi), 3) * ((nu / rho) - Number.Power(
Number.Tan(RadPhi), 2
)),
VI = (nu / 120) * Number.Power(Number.Cos(RadPhi), 5) * (
5 - (18 * Number.Power(Number.Tan(RadPhi), 2)) + Number.Power(Number.Tan(RadPhi), 4) + 14 * eta2 - 58 * Number.Power(
Number.Tan(RadPhi), 2
) * eta2
)
in
E0 + p * IV + Number.Power(p, 3) * V + Number.Power(p, 5) * VI
),
//Project Latitude and longitude to Transverse Mercator Northings.
//@param {number} PHI - Latitude in decimal degrees
//@param {number} LAM - Longitude in decimal degrees
//@param {number} a - Ellipsoid semi-major axis in meters
//@param {number} b - Ellipsoid semi-minor axis in meters
//@param {number} e0 - Eastings of false origin in meters
//@param {number} n0 - Northings of false origin in meters
//@param {number} f0 - Central meridian scale factor
//@param {number} PHI0 - Latitude of false origin in decimal degrees
//@param {number} LAM0 - Longitude of false origin in decimal degrees
//@dependencies [Marc()]
Lat_Long_to_North = (
phi as number,
lam as number,
a as number,
b as number,
E0 as number,
N0 as number,
F0 as number,
phi0 as number,
lam0 as number
) as number =>
(
let
//Convert angle measures to radians
RadPhi = phi * Number.PI / 180,
RadLam = lam * Number.PI / 180,
RadPhi0 = phi0 * Number.PI / 180,
RadLam0 = lam0 * Number.PI / 180,
af0 = a * F0,
bf0 = b * F0,
e2 = (af0 * af0 - bf0 * bf0) / (af0 * af0),
n = (af0 - bf0) / (af0 + bf0),
nu = af0 / (Number.Sqrt(1 - e2 * Number.Sin(RadPhi) * Number.Sin(RadPhi))),
rho = (nu * (1 - e2)) / (1 - e2 * Number.Sin(RadPhi) * Number.Sin(RadPhi)),
eta2 = (nu / rho) - 1,
p = RadLam - RadLam0,
M = Marc(bf0, n, RadPhi0, RadPhi),
I = M + N0,
II = (nu / 2) * Number.Sin(RadPhi) * Number.Cos(RadPhi),
III = (nu / 24) * Number.Sin(RadPhi) * Number.Power(Number.Cos(RadPhi), 3) * (
5 - Number.Power(Number.Tan(RadPhi), 2) + 9 * eta2
),
IIIA = (nu / 720) * Number.Sin(RadPhi) * Number.Power(Number.Cos(RadPhi), 5) * (
61 - 58 * Number.Power(Number.Tan(RadPhi), 2) + Number.Power(Number.Tan(RadPhi), 4)
)
in
I + Number.Power(p, 2) * II + Number.Power(p, 4) * III + Number.Power(p, 6) * IIIA
),
//Un-project Transverse Mercator eastings and northings back to latitude.
//@param {number} East - Eastings in meters
//@param {number} North - Northings in meters
//@param {number} a - Ellipsoid semi-major axis in meters
//@param {number} b - Ellipsoid semi-minor axis in meters
//@param {number} e0 - Eastings of false origin in meters
//@param {number} n0 - Northings of false origin in meters
//@param {number} f0 - Central meridian scale factor
//@param {number} PHI0 - Latitude of false origin in decimal degrees
//@param {number} LAM0 - Longitude of false origin in decimal degrees
//@returns {number} - Latitude in decimal degrees
//@dependencies [Marc(),InitialLat()]
E_N_to_Lat = (
East as number,
North as number,
a as number,
b as number,
E0 as number,
N0 as number,
F0 as number,
phi0 as number,
lam0 as number
) as number =>
(
let
//Convert angle measures to radians
RadPhi0 = phi0 * Number.PI / 180,
RadLam0 = lam0 * Number.PI / 180,
//Compute af0, bf0, e squared (e2), n and Et
af0 = a * F0,
bf0 = b * F0,
e2 = (af0 * af0 - bf0 * bf0) / (af0 * af0),
n = (af0 - bf0) / (af0 + bf0),
Et = East - E0,
//Compute initial value for latitude (PHI) in radians
PHId = InitialLat(North, N0, af0, RadPhi0, n, bf0),
//Compute nu, rho and eta2 using value for PHId
nu = af0 / (Number.Sqrt(1 - e2 * Number.Sin(PHId) * Number.Sin(PHId))),
rho = (nu * (1 - e2)) / (1 - e2 * Number.Sin(PHId) * Number.Sin(PHId)),
eta2 = (nu / rho) - 1,
//Compute Latitude
VII = Number.Tan(PHId) / (2 * rho * nu),
VIII = (Number.Tan(PHId) / (24 * rho * Number.Power(nu, 3))) * (
5 + 3 * Number.Power(Number.Tan(PHId), 2) + eta2 - 9 * eta2 * Number.Power(Number.Tan(PHId), 2)
),
IX = (Number.Tan(PHId) / (720 * rho * Number.Power(nu, 5))) * (
61 + 90 * Number.Power(Number.Tan(PHId), 2) + 45 * Number.Power(Number.Tan(PHId), 4)
)
in
(180 / Number.PI) * (
PHId - (Number.Power(Et, 2) * VII) + (Number.Power(Et, 4) * VIII) - (Number.Power(Et, 6) * IX)
)
),
//Un-project Transverse Mercator eastings and northings back to longitude.
//@param {number} East - Eastings in meters
//@param {number} North - Northings in meters
//@param {number} a - Ellipsoid semi-major axis in meters
//@param {number} b - Ellipsoid semi-minor axis in meters
//@param {number} e0 - Eastings of false origin in meters
//@param {number} n0 - Northings of false origin in meters
//@param {number} f0 - Central meridian scale factor
//@param {number} PHI0 - Latitude of false origin in decimal degrees
//@param {number} LAM0 - Longitude of false origin in decimal degrees
//@returns {number} - Longitude in decimal degrees
//@dependencies [Marc(),InitialLat()]
E_N_to_Long = (
East as number,
North as number,
a as number,
b as number,
E0 as number,
N0 as number,
F0 as number,
phi0 as number,
lam0 as number
) as number =>
(
let
//Convert angle measures to radians
RadPhi0 = phi0 * Number.PI / 180,
RadLam0 = lam0 * Number.PI / 180,
//Compute af0, bf0, e squared (e2), n and Et
af0 = a * F0,
bf0 = b * F0,
e2 = (af0 * af0 - bf0 * bf0) / (af0 * af0),
n = (af0 - bf0) / (af0 + bf0),
Et = East - E0,
//Compute initial value for latitude (PHI) in radians
PHId = InitialLat(North, N0, af0, RadPhi0, n, bf0),
//Compute nu, rho and eta2 using value for PHId
nu = af0 / (Number.Sqrt(1 - e2 * Number.Sin(PHId) * Number.Sin(PHId))),
rho = (nu * (1 - e2)) / (1 - e2 * Number.Sin(PHId) * Number.Sin(PHId)),
eta2 = (nu / rho) - 1,
//Compute Longitude
X = Number.Power(Number.Cos(PHId), -1) / nu,
XI = (Number.Power(Number.Cos(PHId), -1) / (6 * Number.Power(nu, 3))) * (
nu / rho + 2 * Number.Power(Number.Tan(PHId), 2)
),
XII = (Number.Power(Number.Cos(PHId), -1) / (120 * Number.Power(nu, 5))) * (
5 + 28 * Number.Power(Number.Tan(PHId), 2) + 24 * Number.Power(Number.Tan(PHId), 4)
),
XIIA = (Number.Power(Number.Cos(PHId), -1) / (5040 * Number.Power(nu, 7))) * (
61 + 662 * Number.Power(Number.Tan(PHId), 2) + 1320 * Number.Power(Number.Tan(PHId), 4) + 720 * Number.Power(
Number.Tan(PHId), 6
)
)
in
(180 / Number.PI) * (
RadLam0 + Et * X - Number.Power(Et, 3) * XI + Number.Power(Et, 5) * XII - Number.Power(Et, 7) * XIIA
)
),
//Compute initial value for Latitude (PHI) IN RADIANS.
//@param {number} North - Northings of point in meters
//@param {number} n0 - Northings of false origin in meters
//@param {number} afo - semi major axis multiplied by central meridian scale factor in meters
//@param {number} PHI0 - Latitude of false origin in radians
//@param {number} n - (computed from a, b and f0)
//@param {number} bfo - ellipsoid semi major axis multiplied by central meridian scale factor in meters
//@returns {number} initial value for Latitude in radians
//@dependencies [Marc()]
InitialLat = (North as number, N0 as number, afo as number, PHI0 as number, n as number, bfo as number) as number =>
(
let
//First PHI value (PHI1)
PHI1 = (North - N0) / afo + PHI0,
//Calculate M
M = Marc(bfo, n, PHI0, PHI1),
//Iterate to get final value for InitialLat
iterate = (PHI1, M) =>
(
let
PHI2 = (North - N0 - M) / afo + PHI1,
Mx = Marc(bfo, n, PHI0, PHI2),
result = if Number.Abs(North - N0 - Mx) > 0.00001 then @iterate(PHI2, Mx) else PHI2
in
result
)
in
iterate(PHI1, M)
),
//Compute meridional arc.
//@param {number} bfo - ellipsoid semi major axis multiplied by central meridian scale factor in meters
//@param {number} n - (computed from a, b and f0)
//@param {number} PHI0 - Latitude of false origin in radians
//@param {number} PHI - Final latitude of point in radians
//@returns {number} meridional arc
Marc = (bf0 as number, n as number, phi0 as number, phi as number) as number =>
(
let
//Compute M
M = bf0 * (
((1 + n + ((5 / 4) * Number.Power(n, 2)) + ((5 / 4) * Number.Power(n, 3))) * (phi - phi0)) - (
((3 * n) + (3 * Number.Power(n, 2)) + ((21 / 8) * Number.Power(n, 3))) * (
Number.Sin(phi - phi0)
) * (Number.Cos(phi + phi0))
) + (
(((15 / 8) * Number.Power(n, 2)) + ((15 / 8) * Number.Power(n, 3))) * (
Number.Sin(2 * (phi - phi0))
) * (Number.Cos(2 * (phi + phi0)))
) - (((35 / 24) * Number.Power(n, 3)) * (Number.Sin(3 * (phi - phi0))) * (
Number.Cos(3 * (phi + phi0))
))
)
in
M
),
//get Latitude and Longitude from Eastings and Northings
//@param {number} East - Eastings in meters
//@param {number} North - Northings in meters
//@returns {Array} - [Latitude in decimal degrees, Longitude in decimal degrees]
BNG_to_LL = (eastings as number, northings as number) as record =>
(
let
//Get base lat and long
height = 0,
lat1 = E_N_to_Lat(
eastings, northings, 6377563.396, 6356256.91, 400000, -100000, 0.999601272, 49, -2
),
lng1 = E_N_to_Long(
eastings, northings, 6377563.396, 6356256.91, 400000, -100000, 0.999601272, 49, -2
),
//Get cartesian coords
X1 = Lat_Long_H_to_X(lat1, lng1, height, 6377563.396, 6356256.91),
Y1 = Lat_Long_H_to_Y(lat1, lng1, height, 6377563.396, 6356256.91),
Z1 = Lat_H_to_Z(lat1, height, 6377563.396, 6356256.91),
//Offset coordinates to real coordinates
X2 = Helmert_X(X1, Y1, Z1, 446.448, 0.247, 0.8421, -20.4894),
Y2 = Helmert_Y(X1, Y1, Z1, -125.157, 0.1502, 0.8421, -20.4894),
Z2 = Helmert_Z(X1, Y1, Z1, 542.06, 0.1502, 0.247, -20.4894),
//Convert back to latitude, longitude
fLatitude = XYZ_to_Lat(X2, Y2, Z2, 6378137, 6356752.313),
fLongitude = XYZ_to_Long(X2, Y2)
in
[Latitude = fLatitude, Longitude = fLongitude]
),
//Get Eastings and Northings from Latitude and Longitude
//@param {number} Latitude - Latitude in decimal degrees
//@param {number} Longitude - Longitude in decimal degrees
//@returns {Array} - [Eastings in meters, Northings in meters]
LL_to_BNG = (latitude as number, longitude as number) as record =>
(
let
//Get cartesian coords
height = 0,
X1 = Lat_Long_H_to_X(latitude, longitude, height, 6378137, 6356752.313),
Y1 = Lat_Long_H_to_Y(latitude, longitude, height, 6378137, 6356752.313),
Z1 = Lat_H_to_Z(latitude, height, 6378137, 6356752.313),
//Offset coordinates to real coordinates
X2 = Helmert_X(X1, Y1, Z1, -446.448, -0.247, -0.8421, 20.4894),
Y2 = Helmert_Y(X1, Y1, Z1, 125.157, -0.1502, -0.8421, 20.4894),
Z2 = Helmert_Z(X1, Y1, Z1, -542.06, -0.1502, -0.247, 20.4894),
//Convert back to latitude, longitude
latitude2 = XYZ_to_Lat(X2, Y2, Z2, 6377563.396, 6356256.91),
longitude2 = XYZ_to_Long(X2, Y2),
//Get Eastings and Northings
fEastings = Lat_Long_to_East(
latitude2, longitude2, 6377563.396, 6356256.91, 400000, 0.999601272, 49, -2
),
fNorthing = Lat_Long_to_North(
latitude2, longitude2, 6377563.396, 6356256.91, 400000, -100000, 0.999601272, 49, -2
)
in
[Eastings = fEastings, Northings = fNorthing]
),
ApproxEq = (a as number, b as number) as logical => Number.Abs(a - b) < 0.01,
//Test function
Test = [
test_BNG_to_LL = BNG_to_LL(651409.903, 313177.270),
test_LL_to_BNG = LL_to_BNG(52.657570301, 1.717921589),
test_Iterate_XYZ_to_Lat = ApproxEq(
Iterate_XYZ_to_Lat(6378137, 0.00669438037928458, 0.889141291630974, 4929664.92344439, 4026899.36734566),
0.889141264840318
),
test_Marc = ApproxEq(
Marc(6353722.49239479, 1.67322025032507E-03, 0.855211333477221, 0.890263142424236), 223275.25990605
),
test_InitialLat = ApproxEq(
InitialLat(123456, -100000, 6375020.48290224, 0.855211333477221, 1.67322025032507E-03, 6353722.49239479),
0.890291511758942
),
test_E_N_to_Lat = ApproxEq(
E_N_to_Lat(123456, 123456, 6377563.396, 6356256.91, 400000, -100000, 0.999601272, 49, -2),
50.9435419443282
),
test_E_N_to_Long = ApproxEq(
E_N_to_Long(123456, 123456, 6377563.396, 6356256.91, 400000, -100000, 0.999601272, 49, -2),
-5.93731960157555
),
test_Lat_Long_H_to_X = ApproxEq(
Lat_Long_H_to_X(50.9435419443282, -5.93731960157555, 0, 6377563.396, 6356256.91), 4004918.97974917
),
test_Lat_Long_H_to_Y = ApproxEq(
Lat_Long_H_to_Y(50.9435419443282, -5.93731960157555, 0, 6377563.396, 6356256.91), -416504.755830707
),
test_Lat_H_to_Z = ApproxEq(Lat_H_to_Z(50.9435419443282, 0, 6377563.396, 6356256.91), 4929228.95953249),
test_Helmert_X = ApproxEq(
Helmert_X(4004918.97974917, -416504.755830707, 4929228.95953249, 446.448, 0.247, 0.8421, -20.4894),
4005290.97249257
),
test_Helmert_Y = ApproxEq(
Helmert_Y(4004918.97974917, -416504.755830707, 4929228.95953249, -125.157, 0.1502, 0.8421, -20.4894),
-416608.617767794
),
test_Helmert_Z = ApproxEq(
Helmert_Z(4004918.97974917, -416504.755830707, 4929228.95953249, 542.06, 0.1502, 0.247, -20.4894),
4929664.92344439
),
test_XYZ_to_Lat = ApproxEq(
XYZ_to_Lat(4005290.97249257, -416608.617767794, 4929664.92344439, 6378137, 6356752.313), 50.944041866274
),
test_XYZ_to_Long = ApproxEq(XYZ_to_Long(4005290.97249257, -416608.617767794), -5.93824195865157)
]
in
[
BNG_to_LL = BNG_to_LL,
LL_to_BNG = LL_to_BNG,
Tests = Test
]
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment