diff options
| -rw-r--r-- | src/projections/ortho.cpp | 15 | ||||
| -rw-r--r-- | test/gie/builtins.gie | 41 |
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 |