From 4c86bb57358d59f306fb43a02a100b2e65410966 Mon Sep 17 00:00:00 2001 From: Brychan Dempsey Date: Thu, 18 Nov 2021 13:47:54 +1300 Subject: [PATCH] Converted some of it --- NZTM-NZGD Converter/Class1.cs | 182 ++++++++++++++++++ .../NZTM-NZGD Converter.csproj | 8 + NZTM-NZGD Converter/NZTM-NZGD Converter.sln | 25 +++ 3 files changed, 215 insertions(+) create mode 100644 NZTM-NZGD Converter/Class1.cs create mode 100644 NZTM-NZGD Converter/NZTM-NZGD Converter.csproj create mode 100644 NZTM-NZGD Converter/NZTM-NZGD Converter.sln diff --git a/NZTM-NZGD Converter/Class1.cs b/NZTM-NZGD Converter/Class1.cs new file mode 100644 index 0000000..9e6d623 --- /dev/null +++ b/NZTM-NZGD Converter/Class1.cs @@ -0,0 +1,182 @@ +using System; + +namespace NZTM_NZGD_Converter +{ + /// + /// Conversion of the c project by Land Information New Zealand + /// + 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; + + /// + /// Constructs a new TM Projected point + /// + /// + /// + /// + /// + /// + /// + /// + /// + 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; } + + } + /// + /// 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 + /// + /// Source e2 + /// Source a + /// Source lt + /// Arc-length in meters + 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 - + } + } +} diff --git a/NZTM-NZGD Converter/NZTM-NZGD Converter.csproj b/NZTM-NZGD Converter/NZTM-NZGD Converter.csproj new file mode 100644 index 0000000..e20a618 --- /dev/null +++ b/NZTM-NZGD Converter/NZTM-NZGD Converter.csproj @@ -0,0 +1,8 @@ + + + + net5.0 + NZTM_NZGD_Converter + + + diff --git a/NZTM-NZGD Converter/NZTM-NZGD Converter.sln b/NZTM-NZGD Converter/NZTM-NZGD Converter.sln new file mode 100644 index 0000000..9bb193e --- /dev/null +++ b/NZTM-NZGD Converter/NZTM-NZGD Converter.sln @@ -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