From d209e4ab62e9ab6f1c3825acd518bec421aebf12 Mon Sep 17 00:00:00 2001
From: Even Rouault <even.rouault@spatialys.com>
Date: Wed, 15 Sep 2021 06:41:23 +0200
Subject: Inverse ellipsoidal orthographic projection (oblique case): fix
 convergence at pole

---
 src/projections/ortho.cpp | 12 +++++++--
 test/gie/builtins.gie     | 62 +++++++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 72 insertions(+), 2 deletions(-)

diff --git a/src/projections/ortho.cpp b/src/projections/ortho.cpp
index d2d113ef..c334f398 100644
--- a/src/projections/ortho.cpp
+++ b/src/projections/ortho.cpp
@@ -261,8 +261,16 @@ static PJ_LP ortho_e_inverse (PJ_XY xy, PJ *P) {           /* Ellipsoidal, inver
         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
index 9cebec84..569ca812 100644
--- a/test/gie/builtins.gie
+++ b/test/gie/builtins.gie
@@ -4531,6 +4531,68 @@ 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
-- 
cgit v1.2.3