From 55112033a4131cc9e5a0f6b2a40cec555d87955a Mon Sep 17 00:00:00 2001 From: Kristian Evers Date: Wed, 5 Feb 2020 07:32:32 +0100 Subject: 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. --- src/conversions/cart.cpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) (limited to 'src') 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); -- cgit v1.2.3