This link is pretty close to what I am aftering (Converting NZMG (or NZTM) to latitude / longitude for use with R map library).
Basically I am looking for a formula that can convert NZTM coordinates to lat/long (to be run on C# platform). I need this implemented in a C# code that can take in a pair of NZTM coordinates and return lat/long result. In other words:
public double[] (double nztmLat, double nztmLong){
//Implementation details
//return lat/long in a double array.
}
Converting from NZTM to Lat/Long is not a simple task, and requires either a projection library, such as Proj4Net, or a few well-tuned functions.
LINZ has a free nztm.zip program written in ANSI C code for converting coordinates between the New Zealand Transverse Merctator and latitude and longitude on the New Zealand Geodetic Datum 2000. You can port parts of this into C#, I suggest doing this manually.
Looking at nztm.h there are two functions that do the forward and inverse conversions:
void nztm_geod( double n, double e, double *lt, double *ln );
void geod_nztm( double lt, double ln, double *n, double *e );
A quick demo of compiling and running the program:
$ gcc -lm -o nztm nztm.c
$ echo '1576041.15 6188574.24' | ./nztm
Enter NZTM easting, northing:
Input NZTM e,n: 1576041.150 6188574.240
Output Lat/Long: -34.444066 172.739194
Output NZTM e,n: 1576041.150 6188574.240
Difference: 0.000 -0.000
In case anyone is interested, I've implemented this in Python to run over a list of sites with NZTM coordinates. Here's my code:
def nztm_to_lat_long(input_filename, output_filename):
import math
infile = open(input_filename)
lines = infile.readlines()
infile.close()
outfile = open(output_filename, 'w')
outfile.write("site_id, latitude, longitude\n")
for line in lines[1:]: #skip header row
site_id, nztm_e, nztm_n = line.split(',')
#Common variables for NZTM2000
a = 6378137;
f = 1 / 298.257222101;
phizero = 0;
lambdazero = 173;
Nzero = 10000000;
Ezero = 1600000;
kzero = 0.9996;
#input Northing(Y); Easting(X) variables
N = int(nztm_n);
E = int(nztm_e);
#Calculation: From NZTM to lat/Long
b = a * (1 - f);
esq = 2 * f - f ** 2;
Z0 = 1 - esq / 4 - 3 * (esq ** 2) / 64 - 5 * (esq ** 3) / 256;
A2 = 0.375 * (esq + esq ** 2 / 4 + 15 * (esq ** 3) / 128);
A4 = 15 * ((esq ** 2) + 3 * (esq ** 3) / 4) / 256;
A6 = 35 * (esq ** 3) / 3072;
Nprime = N - Nzero;
mprime = Nprime / kzero;
smn = (a - b) / (a + b);
G = a * (1 - smn) * (1 - (smn ** 2)) * (1 + 9 * (smn ** 2) / 4 + 225 * (smn ** 4) / 64) * math.pi/ 180.0;
sigma = mprime * math.pi / (180 * G);
phiprime = sigma + (3 * smn / 2 - 27 * (smn ** 3) / 32) * math.sin(2 * sigma) + (21 * (smn ** 2) / 16 - 55 * (smn ** 4) / 32) * math.sin(4 * sigma) + (151 * (smn ** 3) / 96) * math.sin(6 * sigma) + (1097 * (smn ** 4) / 512) *math.sin(8 * sigma);
rhoprime = a * (1 - esq) / ((1 - esq * ((math.sin(phiprime)) ** 2)) ** 1.5);
upsilonprime = a / math.sqrt(1 - esq * ((math.sin(phiprime)) ** 2));
psiprime = upsilonprime / rhoprime;
tprime = math.tan(phiprime);
Eprime = E - Ezero;
chi = Eprime / (kzero * upsilonprime);
term_1 = tprime * Eprime * chi / (kzero * rhoprime * 2);
term_2 = term_1 * (chi ** 2) / 12 * (-4 * (psiprime ** 2) + 9 * psiprime * (1 - (tprime ** 2)) + 12 * (tprime ** 2));
term_3 = tprime * Eprime * (chi ** 5) / (kzero * rhoprime * 720) * (8 * (psiprime ** 4) * (11 - 24 * (tprime ** 2)) - 12 * (psiprime ** 3) * (21 - 71 * (tprime ** 2)) + 15 * (psiprime ** 2) * (15 - 98 * (tprime ** 2) + 15 * (tprime ** 4)) + 180 * psiprime * (5 * (tprime ** 2) - 3 * (tprime ** 4)) + 360 * (tprime ** 4));
term_4 = tprime * Eprime * (chi ** 7) / (kzero * rhoprime * 40320) * (1385 + 3633 * (tprime ** 2) + 4095 * (tprime ** 4) + 1575 * (tprime ** 6));
term1 = chi * (1 / math.cos(phiprime));
term2 = (chi ** 3) * (1 / math.cos(phiprime)) / 6 * (psiprime + 2 * (tprime ** 2));
term3 = (chi ** 5) * (1 / math.cos(phiprime)) / 120 * (-4 * (psiprime ** 3) * (1 - 6 * (tprime ** 2)) + (psiprime ** 2) * (9 - 68 * (tprime ** 2)) + 72 * psiprime * (tprime ** 2) + 24 * (tprime ** 4));
term4 = (chi ** 7) * (1 / math.cos(phiprime)) / 5040 * (61 + 662 * (tprime ** 2) + 1320 * (tprime ** 4) + 720 * (tprime ** 6));
latitude = (phiprime - term_1 + term_2 - term_3 + term_4) * 180 / math.pi;
longitude = lambdazero + 180 / math.pi * (term1 - term2 + term3 - term4);
outfile.write("{}, {}, {}\n".format(site_id, latitude, longitude))
outfile.close()
And my call to the command line: nztm_to_lat_long('NZTM.csv', 'lat_long.csv'). To use this command just make sure to create the blank output CSV before running the code and that both CSVs are closed and saved in the same location as your .py file.
Here is a C# port of LINZ free nztm.zip download for anyone that's interested. I have tried to stick as close as possible to the original code by using the unsafe flag to handle memory pointers in original C code. You will need to enable -unsafe flag on c# compiler.
Program.cs
using System;
using System.Threading;
namespace NZGD_to_WGS84
{
unsafe class Program
{
public const double PI = 3.1415926535898;
public const double TWOPI = 2.0 * PI;
public const double rad2deg = 180 / PI;
static void Main(string[] args)
{
double e, n, lt, ln, e1, n1;
NZTM nztm = new NZTM();
Console.WriteLine("Enter NZTM easting, northing: ");
string line = Console.ReadLine();
var parts = line.Split(' ');
e = Convert.ToDouble(parts[0]);
n = Convert.ToDouble(parts[1]);
nztm.nztm_geod(n, e, <, &ln);
nztm.geod_nztm(lt, ln, &n1, &e1);
Console.WriteLine("Output Lat/Long: " + lt * rad2deg + " " + ln * rad2deg);
Console.WriteLine("Output NZTM e,n: " + e1 + " " + n1);
Console.WriteLine("Difference: " + (e1 - e) + " " + (n1 - n));
Thread.Sleep(20000);
}
}
}
NZTM.cs
using System;
namespace NZGD_to_WGS84
{
/* Structure used to define a TM projection */
public unsafe struct tmprojection
{
public double meridian; /* Central meridian */
public double scalef; /* Scale factor */
public double orglat; /* Origin latitude */
public double falsee; /* False easting */
public double falsen; /* False northing */
public double utom; /* Unit to metre conversion */
public double a, rf, f, e2, ep2; /* Ellipsoid parameters */
public double om; /* Intermediate calculation */
}
public unsafe class NZTM
{
public const double PI = 3.1415926535898;
public const double TWOPI = 2.0 * PI;
public const double rad2deg = 180 / PI;
public const double NZTM_A = 6378137;
public const double NZTM_RF = 298.257222101;
public const double NZTM_CM = 173.0;
public const double NZTM_OLAT = 0.0;
public const double NZTM_SF = 0.9996;
public const double NZTM_FE = 1600000.0;
public const double NZTM_FN = 10000000.0;
public static tmprojection nztm_projection;
public static Boolean initiallized = false;
//constructor
public NZTM()
{
}
/***************************************************************************/
/* */
/* meridian_arc */
/* */
/* Returns the length of meridional arc (Helmert formula) */
/* Method based on Redfearn's formulation as expressed in GDA technical */
/* manual at http://www.anzlic.org.au/icsm/gdatm/index.html */
/* */
/* Parameters are */
/* projection */
/* latitude (radians) */
/* */
/* Return value is the arc length in metres */
/* */
/***************************************************************************/
static double meridian_arc(tmprojection* tm, double lt)
{
double e2 = tm->e2;
double a = tm->a;
double e4;
double e6;
double A0;
double A2;
double A4;
double A6;
e4 = e2 * e2;
e6 = e4 * e2;
A0 = 1 - (e2 / 4.0) - (3.0 * e4 / 64.0) - (5.0 * e6 / 256.0);
A2 = (3.0 / 8.0) * (e2 + e4 / 4.0 + 15.0 * e6 / 128.0);
A4 = (15.0 / 256.0) * (e4 + 3.0 * e6 / 4.0);
A6 = 35.0 * e6 / 3072.0;
return a * (A0 * lt - A2 * Math.Sin(2 * lt) + A4 * Math.Sin(4 * lt) - A6 * Math.Sin(6 * lt));
}
/*************************************************************************/
/* */
/* foot_point_lat */
/* */
/* Calculates the foot point latitude from the meridional arc */
/* Method based on Redfearn's formulation as expressed in GDA technical*/
/* manual at http://www.anzlic.org.au/icsm/gdatm/index.html */
/* */
/* Takes parameters */
/* tm definition (for scale factor) */
/* meridional arc (metres) */
/* */
/* Returns the foot point latitude (radians) */
/* */
/*************************************************************************/
static double foot_point_lat(tmprojection* tm, double m)
{
double f = tm->f;
double a = tm->a;
double n;
double n2;
double n3;
double n4;
double g;
double sig;
double phio;
n = f / (2.0 - f);
n2 = n * n;
n3 = n2 * n;
n4 = n2 * n2;
g = a * (1.0 - n) * (1.0 - n2) * (1 + 9.0 * n2 / 4.0 + 225.0 * n4 / 64.0);
sig = m / g;
phio = sig + (3.0 * n / 2.0 - 27.0 * n3 / 32.0) * Math.Sin(2.0 * sig)
+ (21.0 * n2 / 16.0 - 55.0 * n4 / 32.0) * Math.Sin(4.0 * sig)
+ (151.0 * n3 / 96.0) * Math.Sin(6.0 * sig)
+ (1097.0 * n4 / 512.0) * Math.Sin(8.0 * sig);
return phio;
}
/***************************************************************************/
/* */
/* tmgeod */
/* */
/* Routine to convert from Tranverse Mercator to latitude and longitude. */
/* Method based on Redfearn's formulation as expressed in GDA technical */
/* manual at http://www.anzlic.org.au/icsm/gdatm/index.html */
/* */
/* Takes parameters */
/* input easting (metres) */
/* input northing (metres) */
/* output latitude (radians) */
/* output longitude (radians) */
/* */
/***************************************************************************/
static void tm_geod(tmprojection* tm,
double ce, double cn, double* ln, double* lt)
{
double fn = tm->falsen;
double fe = tm->falsee;
double sf = tm->scalef;
double e2 = tm->e2;
double a = tm->a;
double cm = tm->meridian;
double om = tm->om;
double utom = tm->utom;
double cn1;
double fphi;
double slt;
double clt;
double eslt;
double eta;
double rho;
double psi;
double E;
double x;
double x2;
double t;
double t2;
double t4;
double trm1;
double trm2;
double trm3;
double trm4;
cn1 = (cn - fn) * utom / sf + om;
fphi = foot_point_lat(tm, cn1);
slt = Math.Sin(fphi);
clt = Math.Cos(fphi);
eslt = (1.0 - e2 * slt * slt);
eta = a / Math.Sqrt(eslt);
rho = eta * (1.0 - e2) / eslt;
psi = eta / rho;
E = (ce - fe) * utom;
x = E / (eta * sf);
x2 = x * x;
t = slt / clt;
t2 = t * t;
t4 = t2 * t2;
trm1 = 1.0 / 2.0;
trm2 = ((-4.0 * psi
+ 9.0 * (1 - t2)) * psi
+ 12.0 * t2) / 24.0;
trm3 = ((((8.0 * (11.0 - 24.0 * t2) * psi
- 12.0 * (21.0 - 71.0 * t2)) * psi
+ 15.0 * ((15.0 * t2 - 98.0) * t2 + 15)) * psi
+ 180.0 * ((-3.0 * t2 + 5.0) * t2)) * psi + 360.0 * t4) / 720.0;
trm4 = (((1575.0 * t2 + 4095.0) * t2 + 3633.0) * t2 + 1385.0) / 40320.0;
*lt = fphi + (t * x * E / (sf * rho)) * (((trm4 * x2 - trm3) * x2 + trm2) * x2 - trm1);
trm1 = 1.0;
trm2 = (psi + 2.0 * t2) / 6.0;
trm3 = (((-4.0 * (1.0 - 6.0 * t2) * psi
+ (9.0 - 68.0 * t2)) * psi
+ 72.0 * t2) * psi
+ 24.0 * t4) / 120.0;
trm4 = (((720.0 * t2 + 1320.0) * t2 + 662.0) * t2 + 61.0) / 5040.0;
*ln = cm - (x / clt) * (((trm4 * x2 - trm3) * x2 + trm2) * x2 - trm1);
}
/***************************************************************************/
/* */
/* geodtm */
/* */
/* Routine to convert from latitude and longitude to Transverse Mercator.*/
/* Method based on Redfearn's formulation as expressed in GDA technical */
/* manual at http://www.anzlic.org.au/icsm/gdatm/index.html */
/* Loosely based on FORTRAN source code by J.Hannah and A.Broadhurst. */
/* */
/* Takes parameters */
/* input latitude (radians) */
/* input longitude (radians) */
/* output easting (metres) */
/* output northing (metres) */
/* */
/***************************************************************************/
static void geod_tm(tmprojection* tm,
double ln, double lt, double* ce, double* cn)
{
double fn = tm->falsen;
double fe = tm->falsee;
double sf = tm->scalef;
double e2 = tm->e2;
double a = tm->a;
double cm = tm->meridian;
double om = tm->om;
double utom = tm->utom;
double dlon;
double m;
double slt;
double eslt;
double eta;
double rho;
double psi;
double clt;
double w;
double wc;
double wc2;
double t;
double t2;
double t4;
double t6;
double trm1;
double trm2;
double trm3;
double gce;
double trm4;
double gcn;
dlon = ln - cm;
while (dlon > PI) dlon -= TWOPI;
while (dlon < -PI) dlon += TWOPI;
m = meridian_arc(tm, lt);
slt = Math.Sin(lt);
eslt = (1.0 - e2 * slt * slt);
eta = a / Math.Sqrt(eslt);
rho = eta * (1.0 - e2) / eslt;
psi = eta / rho;
clt = Math.Cos(lt);
w = dlon;
wc = clt * w;
wc2 = wc * wc;
t = slt / clt;
t2 = t * t;
t4 = t2 * t2;
t6 = t2 * t4;
trm1 = (psi - t2) / 6.0;
trm2 = (((4.0 * (1.0 - 6.0 * t2) * psi
+ (1.0 + 8.0 * t2)) * psi
- 2.0 * t2) * psi + t4) / 120.0;
trm3 = (61 - 479.0 * t2 + 179.0 * t4 - t6) / 5040.0;
gce = (sf * eta * dlon * clt) * (((trm3 * wc2 + trm2) * wc2 + trm1) * wc2 + 1.0);
*ce = gce / utom + fe;
trm1 = 1.0 / 2.0;
trm2 = ((4.0 * psi + 1) * psi - t2) / 24.0;
trm3 = ((((8.0 * (11.0 - 24.0 * t2) * psi
- 28.0 * (1.0 - 6.0 * t2)) * psi
+ (1.0 - 32.0 * t2)) * psi
- 2.0 * t2) * psi
+ t4) / 720.0;
trm4 = (1385.0 - 3111.0 * t2 + 543.0 * t4 - t6) / 40320.0;
gcn = (eta * t) * ((((trm4 * wc2 + trm3) * wc2 + trm2) * wc2 + trm1) * wc2);
*cn = (gcn + m - om) * sf / utom + fn;
return;
}
static void define_tmprojection(tmprojection* tm, double a, double rf, double cm, double sf, double lto, double fe, double fn, double utom)
{
double f;
tm->meridian = cm;
tm->scalef = sf;
tm->orglat = lto;
tm->falsee = fe;
tm->falsen = fn;
tm->utom = utom;
if (rf != 0.0) f = 1.0 / rf; else f = 0.0;
tm->a = a;
tm->rf = rf;
tm->f = f;
tm->e2 = 2.0 * f - f * f;
tm->ep2 = tm->e2 / (1.0 - tm->e2);
tm->om = meridian_arc(tm, tm->orglat);
}
/* Define a static implementation of tmprojection */
/* Note: for some implementations it may be better to create this
dynamically and develop modified versions of the transformation
functions to take this as a parameter */
static tmprojection* get_nztm_projection()
{
unsafe
{
fixed (tmprojection* p = &nztm_projection) //fix struct memory address
{
if (!initiallized)
{
define_tmprojection(p, NZTM_A, NZTM_RF,
NZTM_CM / rad2deg, NZTM_SF, NZTM_OLAT / rad2deg, NZTM_FE, NZTM_FN,
1.0);
}
initiallized = true;
return p;
}
}
}
/* Functions implementation the TM projection specifically for the
NZTM coordinate system
*/
public void nztm_geod(double n, double e, double *lt, double *ln)
{
tmprojection* nztm = get_nztm_projection();
tm_geod(nztm, e, n, ln, lt);
}
public void geod_nztm(double lt, double ln, double* n, double* e)
{
tmprojection* nztm = get_nztm_projection();
geod_tm(nztm, ln, lt, e, n);
}
} //end class
} //end namespace
I know this wasn't requested, but for those of us working on the web, I've transformed the equation to JS too.
// Convert NZTM2000 to Latitude Longitude
function convert(X, Y) {
var a = 6378137;
var f = 1 / 298.257222101;
var phizero = 0;
var lambdazero = 173;
var Nzero = 10000000;
var Ezero = 1600000;
var kzero = 0.9996;
var N = Y;
var E = X;
var b = a * (1 - f);
var esq = 2 * f - Math.pow(f,2);
var Z0 = 1 - esq / 4 - 3 * Math.pow(esq, 2) / 64 - 5 * Math.pow(esq, 3) / 256;
var A2 = 0.375 * (esq + Math.pow(esq, 2) / 4 + 15 * Math.pow(esq, 3) / 128);
var A4 = 15 * (Math.pow(esq, 2) + 3 * Math.pow(esq, 2) / 4) / 256;
var A6 = 35 * Math.pow(esq, 3) / 3072;
var Nprime = N - Nzero;
var mprime = Nprime / kzero;
var smn = (a - b) / (a + b);
var G = a * (1 - smn) * (1 - Math.pow(smn, 2)) * (1 + 9 * Math.pow(smn, 2) / 4 + 225 * Math.pow(smn, 4) / 64) * Math.PI / 180.0;
var sigma = mprime * Math.PI / (180 * G);
var phiprime = sigma + (3 * smn / 2 - 27 * Math.pow(smn, 3) / 32) * Math.sin(2 * sigma) + (21 * Math.pow(smn, 2) / 16 - 55 * Math.pow(smn, 4) / 32) * Math.sin(4 * sigma) + (151 * Math.pow(smn, 3) / 96) * Math.sin(6 * sigma) + (1097 * Math.pow(smn, 4) / 512) * Math.sin(8 * sigma);
var rhoprime = a * (1 - esq) / Math.pow(Math.pow(1 - esq * Math.sin(phiprime),2),1.5);
var upsilonprime = a / Math.sqrt(1 - esq * Math.pow(Math.sin(phiprime),2));
var psiprime = upsilonprime / rhoprime;
var tprime = Math.tan(phiprime);
var Eprime = E - Ezero;
var chi = Eprime / (kzero * upsilonprime);
var term_1 = tprime * Eprime * chi / (kzero * rhoprime * 2);
var term_2 = term_1 * Math.pow(chi,2) / 12 * (-4 * Math.pow(psiprime,2) + 9 * psiprime * (1 - Math.pow(tprime,2)) + 12 * Math.pow(tprime,2));
var term_3 = tprime * Eprime * Math.pow(chi,5) / (kzero * rhoprime * 720) * (8 * Math.pow(psiprime,4) * (11 - 24 * Math.pow(tprime,2)) - 12 * Math.pow(psiprime,3) * (21 - 71 * Math.pow(tprime,2)) + 15 * Math.pow(psiprime,2) * (15 - 98 * Math.pow(tprime,2) + 15 * Math.pow(tprime,4)) + 180 * psiprime * (5 * Math.pow(tprime,2) - 3 * Math.pow(tprime,4)) + 360 * Math.pow(tprime,4));
var term_4 = tprime * Eprime * Math.pow(chi,7) / (kzero * rhoprime * 40320) * (1385 + 3633 * Math.pow(tprime,2) + 4095 * Math.pow(tprime,4) + 1575 * Math.pow(tprime,6));
var term1 = chi * (1 / Math.cos(phiprime));
var term2 = Math.pow(chi,3) * (1 / Math.cos(phiprime)) / 6 * (psiprime + 2 * Math.pow(tprime,2));
var term3 = Math.pow(chi,5) * (1 / Math.cos(phiprime)) / 120 * (-4 * Math.pow(psiprime,3) * (1 - 6 * Math.pow(tprime,2)) + Math.pow(psiprime,2) * (9 - 68 * Math.pow(tprime,2)) + 72 * psiprime * Math.pow(tprime,2) + 24 * Math.pow(tprime,4));
var term4 = Math.pow(chi,7) * (1 / Math.cos(phiprime)) / 5040 * (61 + 662 * Math.pow(tprime,2) + 1320 * Math.pow(tprime,4) + 720 * Math.pow(tprime,6));
var latitude = (phiprime - term_1 + term_2 - term_3 + term_4) * 180 / Math.PI;
var longitude = lambdazero + 180 / Math.PI * (term1 - term2 + term3 - term4);
console.log(longitude,latitude);
}
I've just done something similar in SQL. The Formula for the conversion is provided here by LINZ.
We implemented this first in Excel, and then SQL. We have shared my SQL code in this post.
Edit: I got bored and had a play. Here you go.This works.
static void Main(string[] args)
{
//Common variables for NZTM2000
double a = 6378137;
double f = 1 / 298.257222101;
double phizero = 0;
double lambdazero = 173;
double Nzero = 10000000;
double Ezero = 1600000;
double kzero = 0.9996;
//input Northing(Y); Easting(X) variables
double N = 5427502.0;
double E = 1749165.0;
//Calculation: From NZTM to lat/Long
double b = a * (1 - f);
double esq = 2 * f - Math.Pow(f, 2);
double Z0 = 1 - esq / 4 - 3 * Math.Pow(esq, 2) / 64 - 5 * Math.Pow(esq, 3) / 256;
double A2 = 0.375 * (esq + Math.Pow(esq, 2) / 4 + 15 * Math.Pow(esq, 3) / 128);
double A4 = 15 * (Math.Pow(esq, 2) + 3 * Math.Pow(esq, 3) / 4) / 256;
double A6 = 35 * Math.Pow(esq, 3) / 3072;
double Nprime = N - Nzero;
double mprime = Nprime / kzero;
double smn = (a - b) / (a + b);
double G = a * (1 - smn) * (1 - Math.Pow(smn, 2)) * (1 + 9 * Math.Pow(smn, 2) / 4 + 225 * Math.Pow(smn, 4) / 64) * Math.PI/ 180.0;
double sigma = mprime * Math.PI / (180 * G);
double phiprime = sigma + (3 * smn / 2 - 27 * Math.Pow(smn, 3) / 32) * Math.Sin(2 * sigma) + (21 * Math.Pow(smn, 2) / 16 - 55 * Math.Pow(smn, 4) / 32) * Math.Sin(4 * sigma) + (151 * Math.Pow(smn, 3) / 96) * Math.Sin(6 * sigma) + (1097 * Math.Pow(smn, 4) / 512) *Math.Sin(8 * sigma);
double rhoprime = a * (1 - esq) / Math.Pow((1 - esq * Math.Pow((Math.Sin(phiprime)), 2)), 1.5);
double upsilonprime = a / Math.Sqrt(1 - esq * Math.Pow((Math.Sin(phiprime)), 2));
double psiprime = upsilonprime / rhoprime;
double tprime = Math.Tan(phiprime);
double Eprime = E - Ezero;
double chi = Eprime / (kzero * upsilonprime);
double term_1 = tprime * Eprime * chi / (kzero * rhoprime * 2);
double term_2 = term_1 * Math.Pow(chi, 2) / 12 * (-4 * Math.Pow(psiprime, 2) + 9 * psiprime * (1 - Math.Pow(tprime, 2)) + 12 * Math.Pow(tprime, 2));
double term_3 = tprime * Eprime * Math.Pow(chi, 5) / (kzero * rhoprime * 720) * (8 * Math.Pow(psiprime, 4) * (11 - 24 * Math.Pow(tprime, 2)) - 12 * Math.Pow(psiprime, 3) * (21 - 71 * Math.Pow(tprime, 2)) + 15 * Math.Pow(psiprime, 2) * (15 - 98 * Math.Pow(tprime, 2) + 15 * Math.Pow(tprime, 4)) + 180 * psiprime * (5 * Math.Pow(tprime, 2) - 3 * Math.Pow(tprime, 4)) + 360 * Math.Pow(tprime, 4));
double term_4 = tprime * Eprime * Math.Pow(chi ,7) / (kzero * rhoprime * 40320) * (1385 + 3633 * Math.Pow(tprime, 2) + 4095 * Math.Pow(tprime, 4) + 1575 * Math.Pow(tprime, 6));
double term1 = chi * (1 / Math.Cos(phiprime));
double term2 = Math.Pow(chi, 3) * (1 / Math.Cos(phiprime)) / 6 * (psiprime + 2 * Math.Pow(tprime, 2));
double term3 = Math.Pow(chi, 5) * (1 / Math.Cos(phiprime)) / 120 * (-4 * Math.Pow(psiprime, 3) * (1 - 6 * Math.Pow(tprime, 2)) + Math.Pow(psiprime, 2) * (9 - 68 * Math.Pow(tprime, 2)) + 72 * psiprime * Math.Pow(tprime, 2) + 24 * Math.Pow(tprime, 4));
double term4 = Math.Pow(chi, 7) * (1 / Math.Cos(phiprime)) / 5040 * (61 + 662 * Math.Pow(tprime, 2) + 1320 * Math.Pow(tprime, 4) + 720 * Math.Pow(tprime, 6));
double latitude = (phiprime - term_1 + term_2 - term_3 + term_4) * 180 / Math.PI;
double longitude = lambdazero + 180 / Math.PI * (term1 - term2 + term3 - term4);
Console.WriteLine("Lat/Long: " + latitude.ToString() + ',' + longitude.ToString());
Console.ReadLine();
}
Here is a C# class which can convert between WGS and NZTM
It can be called like so:
var coords = NZTM.NztmToWgs(northing, easting);
or:
var coords = NZTM.NztmToWgs(lat, lon);
It is largely based on @ozzyzig's answer, but I've refactored it to bring it into the 21st century.
public class NZTM
{
/* Structure used to define a TM projection */
protected struct tmprojection
{
internal double meridian; /* Central meridian */
internal double scalef; /* Scale factor */
internal double orglat; /* Origin latitude */
internal double falsee; /* False easting */
internal double falsen; /* False northing */
internal double utom; /* Unit to metre conversion */
internal double a, rf, f, e2, ep2; /* Ellipsoid parameters */
internal double om; /* Intermediate calculation */
}
protected const double PI = 3.1415926535898;
protected const double TWOPI = 2.0 * PI;
protected const double rad2deg = 180 / PI;
protected const double NZTM_A = 6378137;
protected const double NZTM_RF = 298.257222101;
protected const double NZTM_CM = 173.0;
protected const double NZTM_OLAT = 0.0;
protected const double NZTM_SF = 0.9996;
protected const double NZTM_FE = 1600000.0;
protected const double NZTM_FN = 10000000.0;
protected static tmprojection nztm_projection;
/***************************************************************************/
/* */
/* meridian_arc */
/* */
/* Returns the length of meridional arc (Helmert formula) */
/* Method based on Redfearn's formulation as expressed in GDA technical */
/* manual at http://www.anzlic.org.au/icsm/gdatm/index.html */
/* */
/* Parameters are */
/* projection */
/* latitude (radians) */
/* */
/* Return value is the arc length in metres */
/* */
/***************************************************************************/
protected static double meridian_arc(tmprojection tm, double lt)
{
double e2 = tm.e2;
double a = tm.a;
double e4;
double e6;
double A0;
double A2;
double A4;
double A6;
e4 = e2 * e2;
e6 = e4 * e2;
A0 = 1 - (e2 / 4.0) - (3.0 * e4 / 64.0) - (5.0 * e6 / 256.0);
A2 = (3.0 / 8.0) * (e2 + e4 / 4.0 + 15.0 * e6 / 128.0);
A4 = (15.0 / 256.0) * (e4 + 3.0 * e6 / 4.0);
A6 = 35.0 * e6 / 3072.0;
return a * (A0 * lt - A2 * Math.Sin(2 * lt) + A4 * Math.Sin(4 * lt) - A6 * Math.Sin(6 * lt));
}
/*************************************************************************/
/* */
/* foot_point_lat */
/* */
/* Calculates the foot point latitude from the meridional arc */
/* Method based on Redfearn's formulation as expressed in GDA technical*/
/* manual at http://www.anzlic.org.au/icsm/gdatm/index.html */
/* */
/* Takes parameters */
/* tm definition (for scale factor) */
/* meridional arc (metres) */
/* */
/* Returns the foot point latitude (radians) */
/* */
/*************************************************************************/
protected static double foot_point_lat(tmprojection tm, double m)
{
double f = tm.f;
double a = tm.a;
double n;
double n2;
double n3;
double n4;
double g;
double sig;
double phio;
n = f / (2.0 - f);
n2 = n * n;
n3 = n2 * n;
n4 = n2 * n2;
g = a * (1.0 - n) * (1.0 - n2) * (1 + 9.0 * n2 / 4.0 + 225.0 * n4 / 64.0);
sig = m / g;
phio = sig + (3.0 * n / 2.0 - 27.0 * n3 / 32.0) * Math.Sin(2.0 * sig)
+ (21.0 * n2 / 16.0 - 55.0 * n4 / 32.0) * Math.Sin(4.0 * sig)
+ (151.0 * n3 / 96.0) * Math.Sin(6.0 * sig)
+ (1097.0 * n4 / 512.0) * Math.Sin(8.0 * sig);
return phio;
}
/***************************************************************************/
/* */
/* tmgeod */
/* */
/* Routine to convert from Tranverse Mercator to latitude and longitude. */
/* Method based on Redfearn's formulation as expressed in GDA technical */
/* manual at http://www.anzlic.org.au/icsm/gdatm/index.html */
/* */
/* Takes parameters */
/* input easting (metres) */
/* input northing (metres) */
/* output latitude (radians) */
/* output longitude (radians) */
/* */
/***************************************************************************/
protected static(double lat, double lon) tm_geod(tmprojection tm,
double ce, double cn)
{
double lat;
double lon;
double fn = tm.falsen;
double fe = tm.falsee;
double sf = tm.scalef;
double e2 = tm.e2;
double a = tm.a;
double cm = tm.meridian;
double om = tm.om;
double utom = tm.utom;
double cn1;
double fphi;
double slt;
double clt;
double eslt;
double eta;
double rho;
double psi;
double E;
double x;
double x2;
double t;
double t2;
double t4;
double trm1;
double trm2;
double trm3;
double trm4;
cn1 = (cn - fn) * utom / sf + om;
fphi = foot_point_lat(tm, cn1);
slt = Math.Sin(fphi);
clt = Math.Cos(fphi);
eslt = (1.0 - e2 * slt * slt);
eta = a / Math.Sqrt(eslt);
rho = eta * (1.0 - e2) / eslt;
psi = eta / rho;
E = (ce - fe) * utom;
x = E / (eta * sf);
x2 = x * x;
t = slt / clt;
t2 = t * t;
t4 = t2 * t2;
trm1 = 1.0 / 2.0;
trm2 = ((-4.0 * psi
+ 9.0 * (1 - t2)) * psi
+ 12.0 * t2) / 24.0;
trm3 = ((((8.0 * (11.0 - 24.0 * t2) * psi
- 12.0 * (21.0 - 71.0 * t2)) * psi
+ 15.0 * ((15.0 * t2 - 98.0) * t2 + 15)) * psi
+ 180.0 * ((-3.0 * t2 + 5.0) * t2)) * psi + 360.0 * t4) / 720.0;
trm4 = (((1575.0 * t2 + 4095.0) * t2 + 3633.0) * t2 + 1385.0) / 40320.0;
lat = fphi + (t * x * E / (sf * rho)) * (((trm4 * x2 - trm3) * x2 + trm2) * x2 - trm1);
trm1 = 1.0;
trm2 = (psi + 2.0 * t2) / 6.0;
trm3 = (((-4.0 * (1.0 - 6.0 * t2) * psi
+ (9.0 - 68.0 * t2)) * psi
+ 72.0 * t2) * psi
+ 24.0 * t4) / 120.0;
trm4 = (((720.0 * t2 + 1320.0) * t2 + 662.0) * t2 + 61.0) / 5040.0;
lon = cm - (x / clt) * (((trm4 * x2 - trm3) * x2 + trm2) * x2 - trm1);
return (lat * rad2deg, lon * rad2deg);
}
/***************************************************************************/
/* */
/* geodtm */
/* */
/* Routine to convert from latitude and longitude to Transverse Mercator.*/
/* Method based on Redfearn's formulation as expressed in GDA technical */
/* manual at http://www.anzlic.org.au/icsm/gdatm/index.html */
/* Loosely based on FORTRAN source code by J.Hannah and A.Broadhurst. */
/* */
/* Takes parameters */
/* input latitude (radians) */
/* input longitude (radians) */
/* output easting (metres) */
/* output northing (metres) */
/* */
/***************************************************************************/
protected static(double n, double e) geod_tm(tmprojection tm,
double ln, double lt)
{
double fn = tm.falsen;
double fe = tm.falsee;
double sf = tm.scalef;
double e2 = tm.e2;
double a = tm.a;
double cm = tm.meridian;
double om = tm.om;
double utom = tm.utom;
double dlon;
double m;
double slt;
double eslt;
double eta;
double rho;
double psi;
double clt;
double w;
double wc;
double wc2;
double t;
double t2;
double t4;
double t6;
double trm1;
double trm2;
double trm3;
double gce;
double trm4;
double gcn;
dlon = ln - cm;
while (dlon > PI) dlon -= TWOPI;
while (dlon < -PI) dlon += TWOPI;
m = meridian_arc(tm, lt);
slt = Math.Sin(lt);
eslt = (1.0 - e2 * slt * slt);
eta = a / Math.Sqrt(eslt);
rho = eta * (1.0 - e2) / eslt;
psi = eta / rho;
clt = Math.Cos(lt);
w = dlon;
wc = clt * w;
wc2 = wc * wc;
t = slt / clt;
t2 = t * t;
t4 = t2 * t2;
t6 = t2 * t4;
trm1 = (psi - t2) / 6.0;
trm2 = (((4.0 * (1.0 - 6.0 * t2) * psi
+ (1.0 + 8.0 * t2)) * psi
- 2.0 * t2) * psi + t4) / 120.0;
trm3 = (61 - 479.0 * t2 + 179.0 * t4 - t6) / 5040.0;
gce = (sf * eta * dlon * clt) * (((trm3 * wc2 + trm2) * wc2 + trm1) * wc2 + 1.0);
double e = gce / utom + fe;
trm1 = 1.0 / 2.0;
trm2 = ((4.0 * psi + 1) * psi - t2) / 24.0;
trm3 = ((((8.0 * (11.0 - 24.0 * t2) * psi
- 28.0 * (1.0 - 6.0 * t2)) * psi
+ (1.0 - 32.0 * t2)) * psi
- 2.0 * t2) * psi
+ t4) / 720.0;
trm4 = (1385.0 - 3111.0 * t2 + 543.0 * t4 - t6) / 40320.0;
gcn = (eta * t) * ((((trm4 * wc2 + trm3) * wc2 + trm2) * wc2 + trm1) * wc2);
double n = (gcn + m - om) * sf / utom + fn;
return (n, e);
}
protected static tmprojection define_tmprojection(double a, double rf, double cm, double sf, double lto, double fe, double fn, double utom)
{
var tm = new tmprojection();
double f;
tm.meridian = cm;
tm.scalef = sf;
tm.orglat = lto;
tm.falsee = fe;
tm.falsen = fn;
tm.utom = utom;
if (rf != 0.0) f = 1.0 / rf; else f = 0.0;
tm.a = a;
tm.rf = rf;
tm.f = f;
tm.e2 = 2.0 * f - f * f;
tm.ep2 = tm.e2 / (1.0 - tm.e2);
tm.om = meridian_arc(tm, tm.orglat);
return tm;
}
protected static tmprojection nztmProjection = define_tmprojection(NZTM_A, NZTM_RF,
NZTM_CM / rad2deg, NZTM_SF, NZTM_OLAT / rad2deg, NZTM_FE, NZTM_FN,
1.0);
/* Functions implementation the TM projection specifically for the
NZTM coordinate system
*/
public static (double lat, double lon) NztmToWgs(double n, double e)
{
return tm_geod(nztmProjection, e, n);
}
public static (double n, double e) WgsToNztm(double lat, double lon)
{
return geod_tm(nztmProjection, lon / rad2deg, lat / rad2deg);
}
}
Created a javascript module "nztm2000-latlng" and released it to NPM repository. The precision is comparable to that of the C implementation.
Javascript is not precise enough in floating point calculations, so the module employs the decimal.js library.
Here is the NPM module. Source code can be found from this Github project.
I have taken the code from Sir Swears-a-lot's answer and modified the syntax into ArcGIS Arcade code, ready to be adapted for expressions in ArcGIS environments:
function NZTM_to_latlong(NZTM_EASTING, NZTM_NORTHING)
{
//input Northing(Y); Easting(X) variables
var N = NZTM_NORTHING;
var E = NZTM_EASTING;
//Common variables for NZTM2000
var a = 6378137;
var f = 1/298.257222101;
var phizero = 0;
var lambdazero = 173;
var Nzero = 10000000;
var Ezero = 1600000;
var kzero = 0.9996;
var pi = 3.14159265359;
//Calculation: From NZTM to lat/Long
var b = a * (1 - f);
var esq = 2 * f - Pow(f, 2);
var Z0 = 1 - esq / 4 - 3 * Pow(esq, 2) / 64 - 5 * Pow(esq, 3) / 256;
var A2 = 0.375 * (esq + Pow(esq, 2) / 4 + 15 * Pow(esq, 3) / 128);
var A4 = 15 * (Pow(esq, 2) + 3 * Pow(esq, 3) / 4) / 256;
var A6 = 35 * Pow(esq, 3) / 3072;
var Nprime = N - Nzero;
var mprime = Nprime / kzero;
var smn = (a - b) / (a + b);
var G = a * (1 - smn) * (1 - Pow(smn, 2)) * (1 + 9 * Pow(smn, 2) / 4 + 225 * Pow(smn, 4) / 64) * pi/ 180.0;
var sigma = mprime * pi / (180 * G);
var phiprime = sigma + (3 * smn / 2 - 27 * Pow(smn, 3) / 32) * Sin(2 * sigma) + (21 * Pow(smn, 2) / 16 - 55 * Pow(smn, 4) / 32) * Sin(4 * sigma) + (151 * Pow(smn, 3) / 96) * Sin(6 * sigma) + (1097 * Pow(smn, 4) / 512) *Sin(8 * sigma);
var rhoprime = a * (1 - esq) / Pow((1 - esq * Pow((Sin(phiprime)), 2)), 1.5);
var upsilonprime = a / Sqrt(1 - esq * Pow((Sin(phiprime)), 2));
var psiprime = upsilonprime / rhoprime;
var tprime = Tan(phiprime);
var Eprime = E - Ezero;
var chi = Eprime / (kzero * upsilonprime);
var term_1 = tprime * Eprime * chi / (kzero * rhoprime * 2);
var term_2 = term_1 * Pow(chi, 2) / 12 * (-4 * Pow(psiprime, 2) + 9 * psiprime * (1 - Pow(tprime, 2)) + 12 * Pow(tprime, 2));
var term_3 = tprime * Eprime * Pow(chi, 5) / (kzero * rhoprime * 720) * (8 * Pow(psiprime, 4) * (11 - 24 * Pow(tprime, 2)) - 12 * Pow(psiprime, 3) * (21 - 71 * Pow(tprime, 2)) + 15 * Pow(psiprime, 2) * (15 - 98 * Pow(tprime, 2) + 15 * Pow(tprime, 4)) + 180 * psiprime * (5 * Pow(tprime, 2) - 3 * Pow(tprime, 4)) + 360 * Pow(tprime, 4));
var term_4 = tprime * Eprime * Pow(chi ,7) / (kzero * rhoprime * 40320) * (1385 + 3633 * Pow(tprime, 2) + 4095 * Pow(tprime, 4) + 1575 * Pow(tprime, 6));
var term1 = chi * (1 / Cos(phiprime));
var term2 = Pow(chi, 3) * (1 / Cos(phiprime)) / 6 * (psiprime + 2 * Pow(tprime, 2));
var term3 = Pow(chi, 5) * (1 / Cos(phiprime)) / 120 * (-4 * Pow(psiprime, 3) * (1 - 6 * Pow(tprime, 2)) + Pow(psiprime, 2) * (9 - 68 * Pow(tprime, 2)) + 72 * psiprime * Pow(tprime, 2) + 24 * Pow(tprime, 4));
var term4 = Pow(chi, 7) * (1 / Cos(phiprime)) / 5040 * (61 + 662 * Pow(tprime, 2) + 1320 * Pow(tprime, 4) + 720 * Pow(tprime, 6));
var latitude = (phiprime - term_1 + term_2 - term_3 + term_4) * 180 / pi;
var longitude = lambdazero + 180 / pi * (term1 - term2 + term3 - term4);
var string = "Lat/Long: " + Text(latitude) + ',' + Text(longitude);
return string;
}
Example of calling the above function on a point feature
var x = Geometry($feature).x;
var y = Geometry($feature).y;
NZTM_to_latlong(x,y);
In python, geopandas, takes care of translating between projections. The trick is that geopandas calls NZTM:2000 by the name EPSG:2193.
import geopandas
series = geopandas.GeoSeries.from_wkt(['POINT(1820605.400590 5530627.89)'], crs="EPSG:2193")
print(series)
0 POINT (1820605.401 5530627.890)
WGS84 = 4326
wgs_series = series.to_crs(epsg=WGS84)
print(wgs_series)
0 POINT (175.59742 -40.34575)
had now you have WGS84 co-ordinates which are recognizable decimal longitude and latitude.
LINZ has an online converter on their website.
