Improve Transverse Mercator accuracy

dvdoug · Apr 18, 2022 · aec3885 · aec3885
1 parent 86162fd
commit aec3885
Show file tree

Hide file tree

Showing 5 changed files with 64 additions and 23 deletions.
diff --git a/resources/papers/1002.1417.pdf b/resources/papers/1002.1417.pdf
diff --git a/src/GeographicPoint.php b/src/GeographicPoint.php
@@ -1797,12 +1797,16 @@ public function transverseMercator(
         $f = $ellipsoid->getFlattening();
 
         $n = $f / (2 - $f);
-        $B = ($a / (1 + $n)) * (1 + $n ** 2 / 4 + $n ** 4 / 64);
+        $B = ($a / (1 + $n)) * (1 + $n ** 2 / 4 + $n ** 4 / 64 + $n ** 6 / 256 + (25 / 16384) * $n ** 8);
 
-        $h1 = $n / 2 - (2 / 3) * $n ** 2 + (5 / 16) * $n ** 3 + (41 / 180) * $n ** 4;
-        $h2 = (13 / 48) * $n ** 2 - (3 / 5) * $n ** 3 + (557 / 1440) * $n ** 4;
-        $h3 = (61 / 240) * $n ** 3 - (103 / 140) * $n ** 4;
-        $h4 = (49561 / 161280) * $n ** 4;
+        $h1 = $n / 2 - (2 / 3) * $n ** 2 + (5 / 16) * $n ** 3 + (41 / 180) * $n ** 4 - (127 / 288) * $n ** 5 + (7891 / 37800) * $n ** 6 + (72161 / 387072) * $n ** 7 - (18975107 / 50803200) * $n ** 8;
+        $h2 = (13 / 48) * $n ** 2 - (3 / 5) * $n ** 3 + (557 / 1440) * $n ** 4 + (281 / 630) * $n ** 5 - (1983433 / 1935360) * $n ** 6 + (13769 / 28800) * $n ** 7 + (148003883 / 174182400) * $n ** 8;
+        $h3 = (61 / 240) * $n ** 3 - (103 / 140) * $n ** 4 + (15061 / 26880) * $n ** 5 + (167603 / 181440) * $n ** 6 - (67102379 / 29030400) * $n ** 7 + (79682431 / 79833600) * $n ** 8;
+        $h4 = (49561 / 161280) * $n ** 4 - (179 / 168) * $n ** 5 + (6601661 / 7257600) * $n ** 6 + (97445 / 49896) * $n ** 7 - (40176129013 / 7664025600) * $n ** 8;
+        $h5 = (34729 / 80640) * $n ** 5 - (3418889 / 1995840) * $n ** 6 + (14644087 / 9123840) * $n ** 7 + (2605413599 / 622702080) * $n ** 8;
+        $h6 = (212378941 / 319334400) * $n ** 6 - (30705481 / 10378368) * $n ** 7 + (175214326799 / 58118860800) * $n ** 8;
+        $h7 = (1522256789 / 1383782400) * $n ** 7 - (16759934899 / 3113510400) * $n ** 8;
+        $h8 = (1424729850961 / 743921418240) * $n ** 8;
 
         if ($latitudeOrigin === 0.0) {
             $mO = 0;
@@ -1818,7 +1822,11 @@ public function transverseMercator(
             $xiO2 = $h2 * sin(4 * $xiO0);
             $xiO3 = $h3 * sin(6 * $xiO0);
             $xiO4 = $h4 * sin(8 * $xiO0);
-            $xiO = $xiO0 + $xiO1 + $xiO2 + $xiO3 + $xiO4;
+            $xiO5 = $h5 * sin(10 * $xiO0);
+            $xiO6 = $h6 * sin(12 * $xiO0);
+            $xiO7 = $h7 * sin(14 * $xiO0);
+            $xiO8 = $h8 * sin(16 * $xiO0);
+            $xiO = $xiO0 + $xiO1 + $xiO2 + $xiO3 + $xiO4 + $xiO5 + $xiO6 + $xiO7 + $xiO8;
             $mO = $B * $xiO;
         }
 
@@ -1834,8 +1842,16 @@ public function transverseMercator(
         $eta3 = $h3 * cos(6 * $xi0) * sinh(6 * $eta0);
         $xi4 = $h4 * sin(8 * $xi0) * cosh(8 * $eta0);
         $eta4 = $h4 * cos(8 * $xi0) * sinh(8 * $eta0);
-        $xi = $xi0 + $xi1 + $xi2 + $xi3 + $xi4;
-        $eta = $eta0 + $eta1 + $eta2 + $eta3 + $eta4;
+        $xi5 = $h5 * sin(10 * $xi0) * cosh(10 * $eta0);
+        $eta5 = $h5 * cos(10 * $xi0) * sinh(10 * $eta0);
+        $xi6 = $h6 * sin(12 * $xi0) * cosh(12 * $eta0);
+        $eta6 = $h6 * cos(12 * $xi0) * sinh(12 * $eta0);
+        $xi7 = $h7 * sin(14 * $xi0) * cosh(14 * $eta0);
+        $eta7 = $h7 * cos(14 * $xi0) * sinh(14 * $eta0);
+        $xi8 = $h8 * sin(16 * $xi0) * cosh(16 * $eta0);
+        $eta8 = $h8 * cos(16 * $xi0) * sinh(16 * $eta0);
+        $xi = $xi0 + $xi1 + $xi2 + $xi3 + $xi4 + $xi5 + $xi6 + $xi7+ $xi8;
+        $eta = $eta0 + $eta1 + $eta2 + $eta3 + $eta4 + $eta5 + $eta6 + $eta7+ $eta8;
 
         $easting = $falseEasting->asMetres()->getValue() + $kO * $B * $eta;
         $northing = $falseNorthing->asMetres()->getValue() + $kO * ($B * $xi - $mO);

diff --git a/src/ProjectedPoint.php b/src/ProjectedPoint.php
@@ -1893,12 +1893,16 @@ public function transverseMercator(
         $f = $ellipsoid->getFlattening();
 
         $n = $f / (2 - $f);
-        $B = ($a / (1 + $n)) * (1 + $n ** 2 / 4 + $n ** 4 / 64);
+        $B = ($a / (1 + $n)) * (1 + $n ** 2 / 4 + $n ** 4 / 64 + $n ** 6 / 256 + (25 / 16384) * $n ** 8);
 
-        $h1 = $n / 2 - (2 / 3) * $n ** 2 + (37 / 96) * $n ** 3 - (1 / 360) * $n ** 4;
-        $h2 = (1 / 48) * $n ** 2 + (1 / 15) * $n ** 3 - (437 / 1440) * $n ** 4;
-        $h3 = (17 / 480) * $n ** 3 - (37 / 840) * $n ** 4;
-        $h4 = (4397 / 161280) * $n ** 4;
+        $h1 = $n / 2 - (2 / 3) * $n ** 2 + (5 / 16) * $n ** 3 + (41 / 180) * $n ** 4 - (127 / 288) * $n ** 5 + (7891 / 37800) * $n ** 6 + (72161 / 387072) * $n ** 7 - (18975107 / 50803200) * $n ** 8;
+        $h2 = (13 / 48) * $n ** 2 - (3 / 5) * $n ** 3 + (557 / 1440) * $n ** 4 + (281 / 630) * $n ** 5 - (1983433 / 1935360) * $n ** 6 + (13769 / 28800) * $n ** 7 + (148003883 / 174182400) * $n ** 8;
+        $h3 = (61 / 240) * $n ** 3 - (103 / 140) * $n ** 4 + (15061 / 26880) * $n ** 5 + (167603 / 181440) * $n ** 6 - (67102379 / 29030400) * $n ** 7 + (79682431 / 79833600) * $n ** 8;
+        $h4 = (49561 / 161280) * $n ** 4 - (179 / 168) * $n ** 5 + (6601661 / 7257600) * $n ** 6 + (97445 / 49896) * $n ** 7 - (40176129013 / 7664025600) * $n ** 8;
+        $h5 = (34729 / 80640) * $n ** 5 - (3418889 / 1995840) * $n ** 6 + (14644087 / 9123840) * $n ** 7 + (2605413599 / 622702080) * $n ** 8;
+        $h6 = (212378941 / 319334400) * $n ** 6 - (30705481 / 10378368) * $n ** 7 + (175214326799 / 58118860800) * $n ** 8;
+        $h7 = (1522256789 / 1383782400) * $n ** 7 - (16759934899 / 3113510400) * $n ** 8;
+        $h8 = (1424729850961 / 743921418240) * $n ** 8;
 
         if ($latitudeOrigin === 0.0) {
             $mO = 0;
@@ -1914,10 +1918,23 @@ public function transverseMercator(
             $xiO2 = $h2 * sin(4 * $xiO0);
             $xiO3 = $h3 * sin(6 * $xiO0);
             $xiO4 = $h4 * sin(8 * $xiO0);
-            $xiO = $xiO0 + $xiO1 + $xiO2 + $xiO3 + $xiO4;
+            $xiO5 = $h5 * sin(10 * $xiO0);
+            $xiO6 = $h6 * sin(12 * $xiO0);
+            $xiO7 = $h7 * sin(14 * $xiO0);
+            $xiO8 = $h8 * sin(16 * $xiO0);
+            $xiO = $xiO0 + $xiO1 + $xiO2 + $xiO3 + $xiO4 + $xiO5 + $xiO6 + $xiO7 + $xiO8;
             $mO = $B * $xiO;
         }
 
+        $h1 = $n / 2 - (2 / 3) * $n ** 2 + (37 / 96) * $n ** 3 - (1 / 360) * $n ** 4 - (81 / 512) * $n ** 5 + (96199 / 604800) * $n ** 6 - (5406467 / 38707200) * $n ** 7 + (7944359 / 67737600) * $n ** 8;
+        $h2 = (1 / 48) * $n ** 2 + (1 / 15) * $n ** 3 - (437 / 1440) * $n ** 4 + (46 / 105) * $n ** 5 - (1118711 / 3870720) * $n ** 6 + (51841 / 1209600) * $n ** 7 + (24749483 / 348364800) * $n ** 8;
+        $h3 = (17 / 480) * $n ** 3 - (37 / 840) * $n ** 4 - (209 / 4480) * $n ** 5 + (5569 / 90720) * $n ** 6 + (9261899 / 58060800) * $n ** 7 - (6457463 / 17740800) * $n ** 8;
+        $h4 = (4397 / 161280) * $n ** 4 - (11 / 504) * $n ** 5 - (830251 / 7257600) * $n ** 6 + (466511 / 2494800) * $n ** 7 + (324154477 / 7664025600) * $n ** 8;
+        $h5 = (4583 / 161280) * $n ** 5 - (108847 / 3991680) * $n ** 6 - (8005831 / 63866880) * $n ** 7 + (22894433 / 124540416) * $n ** 8;
+        $h6 = (20648693 / 638668800) * $n ** 6 - (16363163 / 518918400) * $n ** 7 - (2204645983 / 12915302400) * $n ** 8;
+        $h7 = (219941297 / 5535129600) * $n ** 7 - (497323811 / 12454041600) * $n ** 8;
+        $h8 = (191773887257 / 3719607091200) * $n ** 8;
+
         $eta = $easting / ($B * $kO);
         $xi = ($northing + $kO * $mO) / ($B * $kO);
         $xi1 = $h1 * sin(2 * $xi) * cosh(2 * $eta);
@@ -1928,8 +1945,16 @@ public function transverseMercator(
         $eta3 = $h3 * cos(6 * $xi) * sinh(6 * $eta);
         $xi4 = $h4 * sin(8 * $xi) * cosh(8 * $eta);
         $eta4 = $h4 * cos(8 * $xi) * sinh(8 * $eta);
-        $xi0 = $xi - ($xi1 + $xi2 + $xi3 + $xi4);
-        $eta0 = $eta - ($eta1 + $eta2 + $eta3 + $eta4);
+        $xi5 = $h5 * sin(10 * $xi) * cosh(10 * $eta);
+        $eta5 = $h5 * cos(10 * $xi) * sinh(10 * $eta);
+        $xi6 = $h6 * sin(12 * $xi) * cosh(12 * $eta);
+        $eta6 = $h6 * cos(12 * $xi) * sinh(12 * $eta);
+        $xi7 = $h7 * sin(14 * $xi) * cosh(14 * $eta);
+        $eta7 = $h7 * cos(14 * $xi) * sinh(14 * $eta);
+        $xi8 = $h8 * sin(16 * $xi) * cosh(16 * $eta);
+        $eta8 = $h8 * cos(16 * $xi) * sinh(16 * $eta);
+        $xi0 = $xi - $xi1 - $xi2 - $xi3 - $xi4 - $xi5 - $xi6 - $xi7 - $xi8;
+        $eta0 = $eta - $eta1 - $eta2 - $eta3 - $eta4 - $eta5 - $eta6 - $eta7 - $eta8;
 
         $beta = self::asin(sin($xi0) / cosh($eta0));
 

diff --git a/tests/CoordinateOperation/AutoConversionTest.php b/tests/CoordinateOperation/AutoConversionTest.php
@@ -208,11 +208,11 @@ public function testBritishNationalGridToUTM(): void
         $to = $from->convert($toCRS);
 
         if (class_exists(OSTN15OSGM15Provider::class)) {
-            self::assertEqualsWithDelta(31326368.174145, $to->getEasting()->getValue(), 0.0001);
-            self::assertEqualsWithDelta(5708455.8991285, $to->getNorthing()->getValue(), 0.0001);
+            self::assertEqualsWithDelta(31326368.093447, $to->getEasting()->getValue(), 0.0001);
+            self::assertEqualsWithDelta(5708454.7196684, $to->getNorthing()->getValue(), 0.0001);
         } else {
-            self::assertEqualsWithDelta(31326366.159092, $to->getEasting()->getValue(), 0.0001);
-            self::assertEqualsWithDelta(5708455.7715515, $to->getNorthing()->getValue(), 0.0001);
+            self::assertEqualsWithDelta(31326366.078970, $to->getEasting()->getValue(), 0.0001);
+            self::assertEqualsWithDelta(5708454.600444621, $to->getNorthing()->getValue(), 0.0001);
         }
     }
 

diff --git a/tests/GeographicPointTest.php b/tests/GeographicPointTest.php
@@ -289,7 +289,7 @@ public function testDistanceDifferentCRSs(): void
         $to = GeographicPoint::create(Geographic2D::fromSRID(Geographic2D::EPSG_PZ_90), new Degree(51.507977), new Degree(-0.124588), null);
 
         if (class_exists(OSTN15OSGM15Provider::class)) {
-            self::assertEqualsWithDelta(3734.135, $from->calculateDistance($to)->getValue(), 0.001);
+            self::assertEqualsWithDelta(3735.308, $from->calculateDistance($to)->getValue(), 0.001);
         } else {
             self::assertEqualsWithDelta(3735.156, $from->calculateDistance($to)->getValue(), 0.001);
         }
@@ -1117,8 +1117,8 @@ public function testOSTN15(): void
         $toCRS = Projected::fromSRID(Projected::EPSG_OSGB36_BRITISH_NATIONAL_GRID);
         $to = $from->OSTN15($toCRS, (new OSTN15OSGM15Provider())->provideGrid());
 
-        self::assertEqualsWithDelta(651409.80373330, $to->getEasting()->asMetres()->getValue(), 0.00000001);
-        self::assertEqualsWithDelta(313177.44988696, $to->getNorthing()->asMetres()->getValue(), 0.00000001);
+        self::assertEqualsWithDelta(651409.804, $to->getEasting()->asMetres()->getValue(), 0.001);
+        self::assertEqualsWithDelta(313177.450, $to->getNorthing()->asMetres()->getValue(), 0.001);
     }
 
     public function testNADCON5Forward3DTransform3DPoint(): void