Merge pull request OSGeo#2853 from rouault/fix_inverse_ortho_e

Fix error in implementation of Inverse ellipsoidal orthographic projection that cause convergence to sometimes fail (fixes OSGeo#2844)
a0x8o · Sep 15, 2021 · b506e1a · b506e1a
1 parent f8d0f50
commit b506e1a
Show file tree

Hide file tree

Showing 2 changed files with 83 additions and 3 deletions.
diff --git a/src/projections/ortho.cpp b/src/projections/ortho.cpp
@@ -254,15 +254,23 @@ static PJ_LP ortho_e_inverse (PJ_XY xy, PJ *P) {           /* Ellipsoidal, inver
         const double J11 = -rho * sinphi * sinlam;
         const double J12 = nu * cosphi * coslam;
         const double J21 = rho * (cosphi * Q->cosph0 + sinphi * Q->sinph0 * coslam);
-        const double J22 = nu * Q->sinph0 * Q->cosph0 * sinlam;
+        const double J22 = nu * Q->sinph0 * cosphi * sinlam;
         const double D = J11 * J22 - J12 * J21;
         const double dx = xy.x - xy_new.x;
         const double dy = xy.y - xy_new.y;
         const double dphi = (J22 * dx - J12 * dy) / D;
         const double dlam = (-J21 * dx + J11 * dy) / D;
         lp.phi += dphi;
-        if( lp.phi > M_PI_2) lp.phi = M_PI_2;
-        else if( lp.phi < -M_PI_2) lp.phi = -M_PI_2;
+        if( lp.phi > M_PI_2)
+        {
+            lp.phi = M_PI_2 - (M_PI_2 - lp.phi);
+            lp.lam = adjlon(lp.lam + 180);
+        }
+        else if( lp.phi < -M_PI_2)
+        {
+            lp.phi = -M_PI_2 + (-M_PI_2 - lp.phi);
+            lp.lam = adjlon(lp.lam + 180);
+        }
         lp.lam += dlam;
         if( fabs(dphi) < 1e-12 && fabs(dlam) < 1e-12 )
         {

diff --git a/test/gie/builtins.gie b/test/gie/builtins.gie
@@ -4531,6 +4531,78 @@ tolerance   0.1 mm
 accept      0               -6343600
 expect      0               -59.966377950099655436
 
+# At pole or very close to it
+direction   forward
+tolerance   0.1 mm
+accept      0              90
+expect      0              5523613.1150
+
+direction   inverse
+tolerance   0.1 mm
+accept      0              5523613.1150
+expect      0              90
+
+direction   forward
+tolerance   0.1 mm
+accept      0              89.99999999
+expect      0              5523613.1145
+roundtrip   1
+
+direction   forward
+tolerance   0.1 mm
+accept      180            89.99999999
+expect      0              5523613.1156
+roundtrip   1
+
+direction   forward
+tolerance   0.1 mm
+accept      90             89.99999999
+expect      0.0011         5523613.1150
+# Roundrip doesn't work on 32bit builds
+#roundtrip   1
+
+direction   forward
+tolerance   0.1 mm
+accept      -90            89.99999999
+expect      -0.0011        5523613.1150
+roundtrip   1
+
+# Southern hemisphere
+operation   +proj=ortho +ellps=WGS84 +lat_0=-30
+
+# At pole or very close to it
+direction   forward
+tolerance   0.1 mm
+accept      0              -90
+expect      0              -5523613.1150
+
+direction   inverse
+tolerance   0.1 mm
+accept      0              -5523613.1150
+expect      0              -90
+
+direction   forward
+tolerance   0.1 mm
+accept      0              -89.99999999
+expect      0              -5523613.1145
+roundtrip   1
+
+direction   forward
+tolerance   0.1 mm
+accept      180            -89.99999999
+expect      0              -5523613.1156
+roundtrip   1
+
+-------------------------------------------------------------------------------
+# Oblique
+# Test case from https://github.com/OSGeo/PROJ/issues/2844
+-------------------------------------------------------------------------------
+operation   +proj=ortho +ellps=WGS84 +lon_0=23.0 +lat_0=37.0
+tolerance   0.1 mm
+accept      120.0           84.0
+expect      663929.0678     5118338.2423
+roundtrip   1
+
 
 -------------------------------------------------------------------------------
 # Equatorial