Ortho: add visibility condition for ellipsoidal case. Credits to @cffk

2 files changed, 56 insertions, 0 deletions

diff --git a/src/projections/ortho.cpp b/src/projections/ortho.cpp
index ff9d9341..0594205b 100644
--- a/src/projections/ortho.cpp
+++ b/src/projections/ortho.cpp
@@ -50,6 +50,13 @@ static PJ_XY ortho_s_forward (PJ_LP lp, PJ *P) {           /* Spheroidal, forwar
         break;
     case OBLIQ:
         sinphi = sin(lp.phi);
+
+        // Is the point visible from the projection plane ?
+        // From https://lists.osgeo.org/pipermail/proj/2020-September/009831.html
+        // this is the dot product of the normal of the ellipsoid at the center of
+        // the projection and at the point considered for projection.
+        // [cos(phi)*cos(lambda), cos(phi)*sin(lambda), sin(phi)]
+        // Also from Snyder's Map Projection - A working manual, equation (5-3), page 149
         if (Q->sinph0 * sinphi + Q->cosph0 * cosphi * coslam < - EPS10)
             return forward_error(P, lp, xy);
         xy.y = Q->cosph0 * sinphi - Q->sinph0 * cosphi * coslam;
@@ -131,6 +138,14 @@ static PJ_XY ortho_e_forward (PJ_LP lp, PJ *P) {           /* Ellipsoidal, forwa
     const double sinphi = sin(lp.phi);
     const double coslam = cos(lp.lam);
     const double sinlam = sin(lp.lam);
+
+    // Is the point visible from the projection plane ?
+    // Same condition as in spherical case
+    if (Q->sinph0 * sinphi + Q->cosph0 * cosphi * coslam < - EPS10) {
+        xy.x = HUGE_VAL; xy.y = HUGE_VAL;
+        return forward_error(P, lp, xy);
+    }
+
     const double nu = 1.0 / sqrt(1.0 - P->es * sinphi * sinphi);
     xy.x = nu * cosphi * sinlam;
     xy.y = nu * (sinphi * Q->cosph0 - cosphi * Q->sinph0 * coslam) +
diff --git a/test/gie/builtins.gie b/test/gie/builtins.gie
index 48afc00a..9d5b6644 100644
--- a/test/gie/builtins.gie
+++ b/test/gie/builtins.gie
@@ -4434,13 +4434,23 @@ accept      2.12955      53.80939444444444
 expect      -189011.711  -128640.567
 roundtrip   1
 
+-------------------------------------------------------------------------------
 # Equatorial
+-------------------------------------------------------------------------------
 operation   +proj=ortho +ellps=WGS84
 tolerance   0.1 mm
 accept      0               0
 expect      0               0
 roundtrip   1
 
+accept      89.99           0
+expect      6378136.9029    0
+roundtrip   1
+
+accept      -89.99          0
+expect      -6378136.9029   0
+roundtrip   1
+
 accept      0               89.99
 expect      0               6356752.2167
 roundtrip   1
@@ -4449,12 +4459,33 @@ accept      0               -89.99
 expect      0               -6356752.2167
 roundtrip   1
 
+# Point not visible from the projection plane
+accept      90.00001        0
+expect      failure         errno tolerance_condition
+
+# Point not visible from the projection plane
+accept      -90.00001       0
+expect      failure         errno tolerance_condition
+
+# Consistant with WGS84 semi-major axis
+# The inverse transformation doesn't converge due to properties of the projection
+accept      90              0
+expect      6378137         0
+
+accept      -90             0
+expect      -6378137        0
+
 # Consistant with WGS84 semi-minor axis
 # The inverse transformation doesn't converge due to properties of the projection
 accept      0               90
 expect      0               6356752.3142
 
+accept      0               -90
+expect      0               -6356752.3142
+
+-------------------------------------------------------------------------------
 # North pole tests
+-------------------------------------------------------------------------------
 operation   +proj=ortho +ellps=WGS84 +lat_0=90
 tolerance   0.1 mm
 accept      0               90
@@ -4465,12 +4496,18 @@ accept      1               89.999
 expect      1.9493         -111.6770
 roundtrip   1
 
+# Point not visible from the projection plane
+accept      0               -0.0000001
+expect      failure         errno tolerance_condition
+
 # Consistant with WGS84 semi-major axis
 # The inverse transformation doesn't converge due to properties of the projection
 accept      0               0
 expect      0               -6378137
 
+-------------------------------------------------------------------------------
 # South pole tests
+-------------------------------------------------------------------------------
 operation   +proj=ortho +ellps=WGS84 +lat_0=-90
 tolerance   0.1 mm
 accept      0               -90
@@ -4481,6 +4518,10 @@ accept      1               -89.999
 expect      1.9493          111.6770
 roundtrip   1
 
+# Point not visible from the projection plane
+accept      0               0.0000001
+expect      failure         errno tolerance_condition
+
 # Consistant with WGS84 semi-major axis
 # The inverse transformation doesn't converge due to properties of the projection
 accept      0               0