183 lines
8.4 KiB
C#
183 lines
8.4 KiB
C#
using System;
|
|
|
|
namespace NZTM_NZGD_Converter
|
|
{
|
|
/// <summary>
|
|
/// Conversion of the c project by Land Information New Zealand
|
|
/// </summary>
|
|
public class Converter
|
|
{
|
|
struct TransverseMercatorProjection
|
|
{
|
|
public double Meridian { get; set; }
|
|
public double ScaleFactor { get; set; }
|
|
public double OriginLatitude { get; set; }
|
|
public double FalseEasting { get; set; }
|
|
public double FalseNorthing { get; set; }
|
|
public double UnitsToMeters { get; set; }
|
|
public double a, rf, f, e2, ep2;
|
|
public double om;
|
|
|
|
/// <summary>
|
|
/// Constructs a new TM Projected point
|
|
/// </summary>
|
|
/// <param name="Meridian"></param>
|
|
/// <param name="ScaleFactor"></param>
|
|
/// <param name="OriginLatitude"></param>
|
|
/// <param name="FalseEasting"></param>
|
|
/// <param name="FalseNorthing"></param>
|
|
/// <param name="UnitsToMeters"></param>
|
|
/// <param name="a"></param>
|
|
/// <param name="rf"></param>
|
|
public TransverseMercatorProjection(double Meridian, double ScaleFactor, double OriginLatitude, double FalseEasting, double FalseNorthing, double UnitsToMeters, double a, double rf)
|
|
{
|
|
this.Meridian = Meridian;
|
|
this.ScaleFactor = ScaleFactor;
|
|
this.OriginLatitude = OriginLatitude;
|
|
this.FalseEasting = FalseNorthing;
|
|
this.FalseNorthing = FalseNorthing;
|
|
this.UnitsToMeters = UnitsToMeters;
|
|
|
|
if (rf != 0.0) f = 1.0 / rf;
|
|
else f = 0.0;
|
|
|
|
this.a = a;
|
|
this.rf = rf;
|
|
e2 = 2.0 * f - f * f;
|
|
ep2 = e2 / (1.0 - e2);
|
|
om = GetMeridianArc(e2, a, OriginLatitude);
|
|
}
|
|
|
|
|
|
|
|
}
|
|
struct GeodeticProjection
|
|
{
|
|
public double Latitude { get; set; }
|
|
public double Longitude { get; set; }
|
|
|
|
}
|
|
/// <summary>
|
|
/// Returns the length of the 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
|
|
/// </summary>
|
|
/// <param name="e2">Source e2</param>
|
|
/// <param name="a">Source a</param>
|
|
/// <param name="lt">Source lt</param>
|
|
/// <returns>Arc-length in meters</returns>
|
|
static double GetMeridianArc(double e2, double a, double lt)
|
|
{
|
|
double e4 = e2 * e2;
|
|
double e6 = e4 * e2;
|
|
|
|
double A0 = 1 - (e2 / 4.0) - (3.0 * e4 / 64.0) - (5.0 * e6 / 256.0);
|
|
double A2 = (3.0 / 8.0) * (e2 + e4 / 4.0 + 15.0 * e6 / 128.0);
|
|
double A4 = (15.0 / 256.0) * (e4 + 3.0 * e6 / 4.0);
|
|
double 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 GetFootPointLatitude(double f, double a, double m)
|
|
{
|
|
double n = f / (2.0 - f);
|
|
double n2 = n * n;
|
|
double n3 = n2 * n;
|
|
double n4 = n2 * n2;
|
|
|
|
double g = a * (1.0 - n) * (1.0 - n2) * (1 + 9.0 * n2 / 4.0 + 225.0 * n4 / 64.0);
|
|
double sig = m / g;
|
|
|
|
return 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);
|
|
}
|
|
|
|
/***************************************************************************/
|
|
/* */
|
|
/* 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 GeodeticProjection TransversMercatorToGeodetic(TransverseMercatorProjection projection, double ce, double cn)
|
|
{
|
|
double cn1 = (cn - projection.FalseNorthing) * projection.UnitsToMeters / projection.ScaleFactor + projection.om;
|
|
double fphi = GetFootPointLatitude(projection.f, projection.a, cn1);
|
|
double slt = Math.Sin(fphi);
|
|
double clt = Math.Cos(fphi);
|
|
|
|
double eslt = (1.0 - projection.e2 * slt * slt);
|
|
double eta = projection.a / Math.Sqrt(eslt);
|
|
double rho = eta * (1.0 - projection.e2) / eslt;
|
|
double psi = eta / rho;
|
|
|
|
double E = (ce - projection.FalseEasting) * projection.UnitsToMeters;
|
|
double x = E / (eta * projection.ScaleFactor);
|
|
double x2 = x * x;
|
|
|
|
|
|
double t = slt / clt;
|
|
double t2 = t * t;
|
|
double t4 = t2 * t2;
|
|
|
|
double trm1 = 1.0 / 2.0;
|
|
|
|
double trm2 = ((-4.0 * psi
|
|
+ 9.0 * (1 - t2)) * psi
|
|
+ 12.0 * t2) / 24.0;
|
|
|
|
double 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;
|
|
|
|
double trm4 = (((1575.0 * t2 + 4095.0) * t2 + 3633.0) * t2 + 1385.0) / 40320.0;
|
|
|
|
GeodeticProjection geodeticProjection = new ();
|
|
geodeticProjection.Latitude = fphi + (t * x * E / (projection.ScaleFactor * 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;
|
|
|
|
geodeticProjection.Longitude = projection.Meridian - (x / clt) * (((trm4 * x2 - trm3) * x2 + trm2) * x2 - trm1);
|
|
|
|
return geodeticProjection;
|
|
}
|
|
|
|
static TransverseMercatorProjection GeodeticToTransverseMercator(GeodeticProjection geodeticProjection)
|
|
{
|
|
TransverseMercatorProjection projection = new();
|
|
|
|
double dlon = geodeticProjection.Longitude -
|
|
}
|
|
}
|
|
}
|