diff options
| author | Even Rouault <even.rouault@spatialys.com> | 2021-09-15 06:41:23 +0200 |
|---|---|---|
| committer | Even Rouault <even.rouault@spatialys.com> | 2021-09-15 06:58:26 +0200 |
| commit | d209e4ab62e9ab6f1c3825acd518bec421aebf12 (patch) | |
| tree | 70ee380c12a202a82b4883783cb40f120ebb533a | |
| parent | 3bd64e8532521d61119df0a6a365a0a32b5c9b15 (diff) | |
| download | PROJ-d209e4ab62e9ab6f1c3825acd518bec421aebf12.tar.gz PROJ-d209e4ab62e9ab6f1c3825acd518bec421aebf12.zip | |
Inverse ellipsoidal orthographic projection (oblique case): fix convergence at pole
| -rw-r--r-- | src/projections/ortho.cpp | 12 | ||||
| -rw-r--r-- | 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 |