BrychanD/NZTM-NZGD-Converter
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

Converted some of it

Browse Source
This commit is contained in:
Brychan Dempsey 2021-11-18 13:47:54 +13:00
parent b8aa8eefdd
commit 4c86bb5735
3 changed files with 215 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

182
NZTM-NZGD Converter/Class1.cs Normal file
View File

@ -0,0 +1,182 @@
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 -
}
}
}

8
NZTM-NZGD Converter/NZTM-NZGD Converter.csproj Normal file
View File

@ -0,0 +1,8 @@
<Project Sdk="Microsoft.NET.Sdk">
<PropertyGroup>
<TargetFramework>net5.0</TargetFramework>
<RootNamespace>NZTM_NZGD_Converter</RootNamespace>
</PropertyGroup>
</Project>

25
NZTM-NZGD Converter/NZTM-NZGD Converter.sln Normal file
View File

@ -0,0 +1,25 @@

Microsoft Visual Studio Solution File, Format Version 12.00
# Visual Studio Version 17
VisualStudioVersion = 17.0.31808.319
MinimumVisualStudioVersion = 10.0.40219.1
Project("{FAE04EC0-301F-11D3-BF4B-00C04F79EFBC}") = "NZTM-NZGD Converter", "NZTM-NZGD Converter.csproj", "{70329027-2E94-4AC6-B3E2-B8B622C133D5}"
EndProject
Global
GlobalSection(SolutionConfigurationPlatforms) = preSolution
Debug|Any CPU = Debug|Any CPU
Release|Any CPU = Release|Any CPU
EndGlobalSection
GlobalSection(ProjectConfigurationPlatforms) = postSolution
{70329027-2E94-4AC6-B3E2-B8B622C133D5}.Debug|Any CPU.ActiveCfg = Debug|Any CPU
{70329027-2E94-4AC6-B3E2-B8B622C133D5}.Debug|Any CPU.Build.0 = Debug|Any CPU
{70329027-2E94-4AC6-B3E2-B8B622C133D5}.Release|Any CPU.ActiveCfg = Release|Any CPU
{70329027-2E94-4AC6-B3E2-B8B622C133D5}.Release|Any CPU.Build.0 = Release|Any CPU
EndGlobalSection
GlobalSection(SolutionProperties) = preSolution
HideSolutionNode = FALSE
EndGlobalSection
GlobalSection(ExtensibilityGlobals) = postSolution
SolutionGuid = {C2E9734C-BBE4-4E14-A92B-9C829D9361FC}
EndGlobalSection
EndGlobal