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 +++-- 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 -- cgit v1.2.3