Merge pull request #4386 from rouault/fix_4385

Inverse +proj=cass: fix non-convergence on inputs where easting=false_easting or northing=false_northing
OSGeo · Jan 18, 2025 · e1733c3 · e1733c3
2 parents c524721 + 90ad86d
commit e1733c3
Show file tree

Hide file tree

Showing 2 changed files with 31 additions and 22 deletions.
diff --git a/src/generic_inverse.cpp b/src/generic_inverse.cpp
@@ -87,29 +87,23 @@ PJ_LP pj_generic_inverse_2d(PJ_XY xy, PJ *P, PJ_LP lpInitial,
             }
         }
 
-        if (xy.x != 0) {
-            // Limit the amplitude of correction to avoid overshoots due to
-            // bad initial guess
-            const double delta_lam = std::max(
-                std::min(deltaX * deriv_lam_X + deltaY * deriv_lam_Y, 0.3),
-                -0.3);
-            lp.lam -= delta_lam;
-            if (lp.lam < -M_PI)
-                lp.lam = -M_PI;
-            else if (lp.lam > M_PI)
-                lp.lam = M_PI;
-        }
+        // Limit the amplitude of correction to avoid overshoots due to
+        // bad initial guess
+        const double delta_lam = std::max(
+            std::min(deltaX * deriv_lam_X + deltaY * deriv_lam_Y, 0.3), -0.3);
+        lp.lam -= delta_lam;
+        if (lp.lam < -M_PI)
+            lp.lam = -M_PI;
+        else if (lp.lam > M_PI)
+            lp.lam = M_PI;
 
-        if (xy.y != 0) {
-            const double delta_phi = std::max(
-                std::min(deltaX * deriv_phi_X + deltaY * deriv_phi_Y, 0.3),
-                -0.3);
-            lp.phi -= delta_phi;
-            if (lp.phi < -M_HALFPI)
-                lp.phi = -M_HALFPI;
-            else if (lp.phi > M_HALFPI)
-                lp.phi = M_HALFPI;
-        }
+        const double delta_phi = std::max(
+            std::min(deltaX * deriv_phi_X + deltaY * deriv_phi_Y, 0.3), -0.3);
+        lp.phi -= delta_phi;
+        if (lp.phi < -M_HALFPI)
+            lp.phi = -M_HALFPI;
+        else if (lp.phi > M_HALFPI)
+            lp.phi = M_HALFPI;
     }
     proj_context_errno_set(P->ctx,
                            PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN);

diff --git a/test/gie/builtins.gie b/test/gie/builtins.gie
@@ -920,6 +920,21 @@ accept  179.99433652777776 -16.841456527777776
 expect  16015.28901692     13369.66005367
 roundtrip 1
 
+
+-------------------------------------------------------------------------------
+# Scenario of https://github.com/OSGeo/PROJ/issues/4385
+-------------------------------------------------------------------------------
+operation +proj=cass +lat_0=50.6177 +lon_0=-1.19725 +x_0=500000 +y_0=100000 +ellps=airy +units=m
+
+tolerance 0.1 mm
+direction inverse
+
+accept 300000 100000
+expect -4.022094267169   50.583438725252
+
+accept 500000 100000
+expect -1.19725 50.6177
+
 ===============================================================================
 # Central Conic
 # 	Sph