diff options
| author | Kristian Evers <kristianevers@gmail.com> | 2020-02-05 07:32:32 +0100 |
|---|---|---|
| committer | Kristian Evers <kristianevers@gmail.com> | 2020-02-06 17:53:11 +0100 |
| commit | 55112033a4131cc9e5a0f6b2a40cec555d87955a (patch) | |
| tree | b5faf4c6cc02578a00cb640c8a57612a20b58610 | |
| parent | bd6148db5f5b9a0a13f3dc2d2c35a066c15a57ed (diff) | |
| download | PROJ-55112033a4131cc9e5a0f6b2a40cec555d87955a.tar.gz PROJ-55112033a4131cc9e5a0f6b2a40cec555d87955a.zip | |
cart: Avoid discontinuity at poles in the inverse case
This should avoid issues with numerical stability as uncovered in
https://github.com/OSGeo/PROJ/issues/1906.
Practically speaking this change isn't going to affect real life
scenarios since the position of the center of the Earth is rarely
expressed in geodetic coordinates.
| -rw-r--r-- | src/conversions/cart.cpp | 5 | ||||
| -rw-r--r-- | test/gie/more_builtins.gie | 2 |
2 files changed, 4 insertions, 3 deletions
diff --git a/src/conversions/cart.cpp b/src/conversions/cart.cpp index a7817443..537fc29f 100644 --- a/src/conversions/cart.cpp +++ b/src/conversions/cart.cpp @@ -165,8 +165,9 @@ static PJ_LPZ geodetic (PJ_XYZ cart, PJ *P) { if( fabs(lpz.phi) > M_HALFPI ) { // this happen on non-sphere ellipsoid when x,y,z is very close to 0 // there is no single solution to the cart->geodetic conversion in - // that case, so arbitrarily pickup phi = 0. - lpz.phi = 0; + // that case, clamp to -90/90 deg and avoid a discontinuous boundary + // near the poles + lpz.phi = copysign(M_HALFPI, lpz.phi); } lpz.lam = atan2 (cart.y, cart.x); N = normal_radius_of_curvature (P->a, P->es, lpz.phi); diff --git a/test/gie/more_builtins.gie b/test/gie/more_builtins.gie index 10d3ea47..424f941d 100644 --- a/test/gie/more_builtins.gie +++ b/test/gie/more_builtins.gie @@ -851,7 +851,7 @@ expect 180 0 0 # Center of Earth ! accept 0 0 0 -expect 0 0 -6378137 +expect 0 90 -6356752.314140356 </gie> |