/
ll2nztm.js
64 lines (53 loc) · 3.7 KB
/
ll2nztm.js
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
// Convert NZGD2000 Latitude Longitude to NZTM2000
// Formulae from "https://www.linz.govt.nz/data/geodetic-services/coordinate-conversion/projection-conversions/transverse-mercator-transformation-formulae"
// Written to run as a function in Google Sheets, but could be used elsewhere
// Pass in Lat, Long in the format -41.000000, 174.000000
function LL2NZTM(LL) {
// Split LL variable into an array
var LLarray = LL.split(',');
Logger.log(LLarray[0]);
Logger.log(LLarray[1]);
// Enter constants
var a = 6378137.0; //Semi-major axis of reference ellipsoid
var f = 1 / 298.257222101; //Ellipsoidal flattening, this could be one of two values, but makes vary little difference
var phizero = 0.0001 * Math.PI / 180; //Origin latitude
var lambdazero = 173.0001 * Math.PI / 180; //Origin longitude
var Nzero = 10000000.0; //False Northing
var Ezero = 1600000.0; //False Easting
var kzero = 0.9996; //Central meridian scale factor
// Pass Lat and Long variables in and convert to radians
var phi = LLarray[0] * Math.PI / 180; //latitude; //Latitude of computation point
var lambda = LLarray[1] * Math.PI / 180; //longitude; //Longitude of computation point
// Work out projection constants
var esq = 2 * f - Math.pow(f,2); //Projection constant
var A0 = 1 - esq / 4 - 3 * Math.pow(esq, 2) / 64 - 5 * Math.pow(esq, 3) / 256; //Projection constant
var A2 = 0.375 * (esq + Math.pow(esq, 2) / 4 + 15 * Math.pow(esq, 3) / 128); //Projection constant
var A4 = 15 / 256 * (Math.pow(esq, 2) + 3 * Math.pow(esq, 3) / 4); //Projection constant
var A6 = 35 * Math.pow(esq, 3) / 3072; //Projection constant
var m = a * (A0 * phi - A2 * Math.sin(2 * phi) + A4 * Math.sin(4 * phi) - A6 * Math.sin(6 * phi)); //Projection constant
var mzero = a * (A0 * phizero - A2 * Math.sin(2 * phizero) + A4 * Math.sin(4 * phizero) - A6 * Math.sin(6 * phizero)); //Projection constant
// Determine variables at coordinates supplied
var rho = a * (1 - esq) / Math.pow(1 - esq * Math.pow(Math.sin(phi),2),1.5);
var upsilon = a / Math.sqrt(1 - esq * Math.pow(Math.sin(phi),2));
var psi = upsilon / rho;
var t = Math.tan(phi);
var omega = lambda - lambdazero;
// Calculate Northing
var Nterm1 = Math.pow(omega,2) / 2 * upsilon * Math.sin(phi) * Math.cos(phi);
var Nterm2 = Math.pow(omega,4) / 24 * upsilon * Math.sin(phi) * Math.pow(Math.cos(phi),3) * (4 * Math.pow(psi,2) + psi - Math.pow(t,2));
var Nterm3 = Math.pow(omega,6) / 720 * upsilon * Math.sin(phi) * Math.pow(Math.cos(phi),5) * (8 * Math.pow(psi,4) * (11 - 24 * Math.pow(t,2)) - 28 * Math.pow(psi,3) * (1 - 6 * Math.pow(t,2)) + Math.pow(psi,2) * (1 - 32 * Math.pow(t,2)) - psi * (2 * Math.pow(t,2)) + Math.pow(t,4));
var Nterm4 = Math.pow(omega,8) / 40320 * upsilon * Math.sin(phi) * Math.pow(Math.cos(phi),7) * (1385 - 3111 * Math.pow(t,2) + 543 * Math.pow(t,4) - Math.pow(t,6));
var N = Nzero + kzero * (m - mzero + Nterm1 + Nterm2 + Nterm3 + Nterm4); //projection northing
// Calculate Easting
var Eterm1 = Math.pow(omega,2) / 6 * Math.pow(Math.cos(phi),2) * (psi - Math.pow(t,2));
var Eterm2 = Math.pow(omega,4) / 120 * Math.pow(Math.cos(phi),4) * (4 * Math.pow(psi,3) * (1 - 6 * Math.pow(t,2)) + Math.pow(psi,2) * (1 + 8 * Math.pow(t,2)) - psi * 2 * Math.pow(t,2) + Math.pow(t,4));
var Eterm3 = Math.pow(omega,6) / 5040 * Math.pow(Math.cos(phi),6) * (61 - 479 * Math.pow(t,2) + 179 * Math.pow(t,4) - Math.pow(t,6));
var E = Ezero + kzero * upsilon * omega * Math.cos(phi) * (1 + Eterm1 + Eterm2 + Eterm3);
Logger.log(N);
Logger.log(E);
// Round to 2dp and concatenate into a single variable
var N = Math.round((N + Number.EPSILON) * 100) / 100;
var E = Math.round((E + Number.EPSILON) * 100) / 100;
var NE = N + ", " + E;
return (NE);
}