cass: fix forward computation of easting (fixes OSGeo#3432), and use …

…iterative generic inversion to make GIGS 5108 roundtripping tests succeed The error on the sign after "nu * A * (1. - A2 * T * (C1 ...") dates back to PROJ origins. All other independent implementations (Apache SIS cf /~https://github.com/apache/sis/blob/a04bd4b29816c91fc4dfc66135b8de5e962933a0/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/CassiniSoldner.java, CSMAP cf /~https://github.com/normanb/CS-Map/blob/1d8f2db7f13194a99dc218d4cd91754c64a2a774/Source/CS_csini.c) use the same formulas as EPSG guidance note which is consistent with Snyder's Map projections - A working manual). So align on them, so that everyone is consistent (and hopefully right). While doing it, I've modified the inverse method to use the generic inversion method to improve its accuracy, which make it possible to enable the 5108.gie tests that were previously failing on the 1000 iterations roundtripping.
rouault · Nov 4, 2022 · c34f791 · c34f791
1 parent 4586998
commit c34f791
Show file tree

Hide file tree

Showing 8 changed files with 48 additions and 29 deletions.
diff --git a/src/generic_inverse.cpp b/src/generic_inverse.cpp
@@ -41,7 +41,7 @@
  *
  * Starts with initial guess provided by user in lpInitial
  */
-PJ_LP pj_generic_inverse_2d(PJ_XY xy, PJ *P, PJ_LP lpInitial) {
+PJ_LP pj_generic_inverse_2d(PJ_XY xy, PJ *P, PJ_LP lpInitial, double deltaXYTolerance) {
     PJ_LP lp = lpInitial;
     double deriv_lam_X = 0;
     double deriv_lam_Y = 0;
@@ -51,7 +51,7 @@ PJ_LP pj_generic_inverse_2d(PJ_XY xy, PJ *P, PJ_LP lpInitial) {
         PJ_XY xyApprox = P->fwd(lp, P);
         const double deltaX = xyApprox.x - xy.x;
         const double deltaY = xyApprox.y - xy.y;
-        if (fabs(deltaX) < 1e-10 && fabs(deltaY) < 1e-10) {
+        if (fabs(deltaX) < deltaXYTolerance && fabs(deltaY) < deltaXYTolerance) {
             return lp;
         }
 

diff --git a/src/proj_internal.h b/src/proj_internal.h
@@ -875,7 +875,7 @@ void pj_clear_vgridshift_knowngrids_cache();
 
 void pj_clear_sqlite_cache();
 
-PJ_LP pj_generic_inverse_2d(PJ_XY xy, PJ *P, PJ_LP lpInitial);
+PJ_LP pj_generic_inverse_2d(PJ_XY xy, PJ *P, PJ_LP lpInitial, double deltaXYTolerance);
 
 
 

diff --git a/src/projections/adams.cpp b/src/projections/adams.cpp
@@ -296,7 +296,9 @@ static PJ_LP adams_inverse(PJ_XY xy, PJ *P)
     lp.lam = fabs(lp.phi) >= M_HALFPI ? 0 :
         std::max(std::min(xy.x / 2.62205760 / cos(lp.phi), 1.0), -1.0) * M_PI;
 
-    return pj_generic_inverse_2d(xy, P, lp);
+
+    constexpr double deltaXYTolerance = 1e-10;
+    return pj_generic_inverse_2d(xy, P, lp, deltaXYTolerance);
 }
 
 static PJ_LP peirce_q_square_inverse(PJ_XY xy, PJ *P)
@@ -339,7 +341,9 @@ static PJ_LP peirce_q_square_inverse(PJ_XY xy, PJ *P)
     }
     else /* if( xy.x <= 0 && xy.y <= 0 ) */
         lp.lam = -M_PI / 2;
-    return pj_generic_inverse_2d(xy, P, lp);
+
+    constexpr double deltaXYTolerance = 1e-10;
+    return pj_generic_inverse_2d(xy, P, lp, deltaXYTolerance);
 }
 
 static PJ_LP peirce_q_diamond_inverse(PJ_XY xy, PJ *P)
@@ -382,7 +386,8 @@ static PJ_LP peirce_q_diamond_inverse(PJ_XY xy, PJ *P)
         lp.phi = -M_PI / 4;
     }
 
-    return pj_generic_inverse_2d(xy, P, lp);
+    constexpr double deltaXYTolerance = 1e-10;
+    return pj_generic_inverse_2d(xy, P, lp, deltaXYTolerance);
 }
 
 static PJ *setup(PJ *P, projection_type mode) {

diff --git a/src/projections/cass.cpp b/src/projections/cass.cpp
@@ -43,7 +43,7 @@ static PJ_XY cass_e_forward (PJ_LP lp, PJ *P) {          /* Ellipsoidal, forward
     const double A2 = A * A;
 
     xy.x = nu * A * (1. - A2 * T *
-        (C1 - (8. - T + 8. * C) * A2 * C2));
+        (C1 + (8. - T + 8. * C) * A2 * C2));
     xy.y = M - Q->m0 + nu * tanphi * A2 *
         (.5 + (5. - T + 6. * C) * A2 * C3);
     if( Q->hyperbolic )
@@ -82,13 +82,13 @@ static PJ_LP cass_e_inverse (PJ_XY xy, PJ *P) {          /* Ellipsoidal, inverse
     lp.lam = D * (1. + T1 * D2 *
         (-C4 + (1. + 3. * T1) * D2 * C5)) / cos (phi1);
 
-    if( Q->hyperbolic )
-    {
-        // EPSG guidance note 7-2 suggests a custom approximation for the
-        // 'Vanua Levu 1915 / Vanua Levu Grid' case, but better use the
-        // generic inversion method
-        lp = pj_generic_inverse_2d(xy, P, lp);
-    }
+    // EPSG guidance note 7-2 suggests a custom approximation for the
+    // 'Vanua Levu 1915 / Vanua Levu Grid' case, but better use the
+    // generic inversion method
+    // Actually use it in the non-hyperbolic case. It enables to make the
+    // 5108.gie roundtripping tests to success, with at most 2 iterations.
+    constexpr double deltaXYTolerance = 1e-12;
+    lp = pj_generic_inverse_2d(xy, P, lp, deltaXYTolerance);
 
     return lp;
 }

diff --git a/src/projections/wink2.cpp b/src/projections/wink2.cpp
@@ -48,7 +48,8 @@ static PJ_LP wink2_s_inverse(PJ_XY xy, PJ *P)
     lpInit.phi = xy.y;
     lpInit.lam = xy.x;
 
-    return pj_generic_inverse_2d(xy, P, lpInit);
+    constexpr double deltaXYTolerance = 1e-10;
+    return pj_generic_inverse_2d(xy, P, lpInit, deltaXYTolerance);
 }
 
 

diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
@@ -80,7 +80,7 @@ proj_add_gie_test("GIGS-5104" "gigs/5104.gie")
 proj_add_gie_test("GIGS-5105.2" "gigs/5105.2.gie")
 proj_add_gie_test("GIGS-5106" "gigs/5106.gie")
 proj_add_gie_test("GIGS-5107" "gigs/5107.gie")
-#proj_add_gie_test("GIGS-5108" "gigs/5108.gie")
+proj_add_gie_test("GIGS-5108" "gigs/5108.gie")
 proj_add_gie_test("GIGS-5109" "gigs/5109.gie")
 #proj_add_gie_test("GIGS-5110" "gigs/5110.gie")
 proj_add_gie_test("GIGS-5111.1" "gigs/5111.1.gie")

diff --git a/test/gie/builtins.gie b/test/gie/builtins.gie
@@ -850,6 +850,18 @@ expect  -0.001790493 0.000895247
 accept  -200 -100
 expect  -0.001790493 -0.000895247
 
+-------------------------------------------------------------------------------
+# test point from EPSG Guidance Note 7.2
+operation +proj=cass +lat_0=10.4416666666667 +lon_0=-61.3333333333333 \
+          +x_0=86501.46392052 +y_0=65379.0134283 \
+          +a=6378293.64520876 +b=6356617.98767984 +to_meter=0.201166195164
+-------------------------------------------------------------------------------
+
+tolerance 0.1 mm
+accept  -62                10
+expect  66644.94040882     82536.21873655
+roundtrip 1
+
 -------------------------------------------------------------------------------
 # Hyperbolic variant: test point from EPSG Guidance Note 7.2
 operation +proj=cass +hyperbolic +a=6378306.376305601 +rf=293.466307 \
@@ -859,7 +871,7 @@ operation +proj=cass +hyperbolic +a=6378306.376305601 +rf=293.466307 \
 
 tolerance 0.1 mm
 accept  179.99433652777776 -16.841456527777776
-expect  16015.289017555102 13369.660053668682
+expect  16015.28901692     13369.66005367
 roundtrip 1
 
 ===============================================================================

diff --git a/test/gigs/5108.gie.failing → test/gigs/5108.gie b/test/gigs/5108.gie.failing → test/gigs/5108.gie
@@ -4,18 +4,19 @@ Test 5108, Cassini-Soldner, v2-0_2011-06-28.
 
 --------------------------------------------------------------------------------
 
-# GDM2000
-<4742> +proj=longlat +ellps=GRS80  <>
+<gie-strict>
 
-# GDM2000 / Johor Grid
-<3377> +proj=cass +lat_0=2.121679744444445 +lon_0=103.4279362361111 +x_0=-14810.562 +y_0=8758.32 +ellps=GRS80 +units=m  <>
+use_proj4_init_rules true
 
+# GDM2000
+# <4742> +proj=longlat +ellps=GRS80  <>
 
-<gie>
+# GDM2000 / Johor Grid
+# <3377> +proj=cass +lat_0=2.121679744444445 +lon_0=103.4279362361111 +x_0=-14810.562 +y_0=8758.32 +ellps=GRS80 +units=m  <>
 
 --------------------------------------------------------------------------------
-operation  +proj=pipeline
-           +step +init=epsg:4742 +inv
+operation  +proj=pipeline \
+           +step +init=epsg:4742 +inv \
            +step +init=epsg:3377
 --------------------------------------------------------------------------------
 tolerance  0.05 m
@@ -87,8 +88,8 @@ accept    104 5
 expect    48630.563 327067.097
 
 --------------------------------------------------------------------------------
-operation  +proj=pipeline
-           +step +init=epsg:3377 +inv
+operation  +proj=pipeline \
+           +step +init=epsg:3377 +inv \
            +step +init=epsg:4742
 --------------------------------------------------------------------------------
 tolerance  0.05 m
@@ -160,8 +161,8 @@ accept    48630.563 327067.097
 expect    104 5
 
 --------------------------------------------------------------------------------
-operation  +proj=pipeline
-           +step +init=epsg:4742 +inv
+operation  +proj=pipeline \
+           +step +init=epsg:4742 +inv \
            +step +init=epsg:3377
 --------------------------------------------------------------------------------
 tolerance  0.006 m
@@ -232,4 +233,4 @@ tolerance  0.006 m
 accept    104 5
 roundtrip 1000
 
-</gie>
+</gie-strict>
Original file line number	Diff line number	Diff line change
Expand Up		@@ -875,7 +875,7 @@ void pj_clear_vgridshift_knowngrids_cache();

		void pj_clear_sqlite_cache();

		PJ_LP pj_generic_inverse_2d(PJ_XY xy, PJ *P, PJ_LP lpInitial);
		PJ_LP pj_generic_inverse_2d(PJ_XY xy, PJ *P, PJ_LP lpInitial, double deltaXYTolerance);



Expand Down