bonne: fix inverse map projection computations when lat_1 < 0 (fixes OSGeo#3848)

The fixes are directly following hints of Snyder's "Map projections: A
working manual"

Author: rouault

src/projections/bonne.cpp

@@ -55,35 +55,45 @@ static PJ_XY bonne_s_forward(PJ_LP lp, PJ *P) { /* Spheroidal, forward */
 static PJ_LP bonne_s_inverse(PJ_XY xy, PJ *P) { /* Spheroidal, inverse */
     PJ_LP lp = {0.0, 0.0};
     struct pj_opaque *Q = static_cast<struct pj_opaque *>(P->opaque);
-    double rh;
 
     xy.y = Q->cphi1 - xy.y;
-    rh = hypot(xy.x, xy.y);
+    const double rh = copysign(hypot(xy.x, xy.y), Q->phi1);
     lp.phi = Q->cphi1 + Q->phi1 - rh;
-    if (fabs(lp.phi) > M_HALFPI) {
+    const double abs_phi = fabs(lp.phi);
+    if (abs_phi > M_HALFPI) {
         proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN);
         return lp;
     }
-    if (fabs(fabs(lp.phi) - M_HALFPI) <= EPS10)
+    if (M_HALFPI - abs_phi <= EPS10)
         lp.lam = 0.;
-    else
-        lp.lam = rh * atan2(xy.x, xy.y) / cos(lp.phi);
+    else {
+        const double lm = rh / cos(lp.phi);
+        if (Q->phi1 > 0) {
+            lp.lam = lm * atan2(xy.x, xy.y);
+        } else {
+            lp.lam = lm * atan2(-xy.x, -xy.y);
+        }
+    }
     return lp;
 }
 
 static PJ_LP bonne_e_inverse(PJ_XY xy, PJ *P) { /* Ellipsoidal, inverse */
     PJ_LP lp = {0.0, 0.0};
     struct pj_opaque *Q = static_cast<struct pj_opaque *>(P->opaque);
-    double s, rh;
 
     xy.y = Q->am1 - xy.y;
-    rh = hypot(xy.x, xy.y);
+    const double rh = copysign(hypot(xy.x, xy.y), Q->phi1);
     lp.phi = pj_inv_mlfn(Q->am1 + Q->m1 - rh, Q->en);
-    if ((s = fabs(lp.phi)) < M_HALFPI) {
-        s = sin(lp.phi);
-        lp.lam =
-            rh * atan2(xy.x, xy.y) * sqrt(1. - P->es * s * s) / cos(lp.phi);
-    } else if (fabs(s - M_HALFPI) <= EPS10)
+    const double abs_phi = fabs(lp.phi);
+    if (abs_phi < M_HALFPI) {
+        const double sinphi = sin(lp.phi);
+        const double lm = rh * sqrt(1. - P->es * sinphi * sinphi) / cos(lp.phi);
+        if (Q->phi1 > 0) {
+            lp.lam = lm * atan2(xy.x, xy.y);
+        } else {
+            lp.lam = lm * atan2(-xy.x, -xy.y);
+        }
+    } else if (abs_phi - M_HALFPI <= EPS10)
         lp.lam = 0.;
     else {
         proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN);

test/gie/builtins.gie

@@ -686,6 +686,18 @@ expect  -0.001796699 0.500904369
 accept  -200 -100
 expect  -0.001796698 0.499095631
 
+-------------------------------------------------------------------------------
+operation +proj=bonne   +ellps=GRS80  +lat_1=-0.5
+-------------------------------------------------------------------------------
+tolerance 0.1 mm
+accept  2 1
+expect  222605.2961 165827.6478
+roundtrip   1
+
+accept  2 -1
+expect  222605.2961 -55321.1396
+roundtrip   1
+
 -------------------------------------------------------------------------------
 operation +proj=bonne   +ellps=GRS80  +lat_1=90
 -------------------------------------------------------------------------------
@@ -697,6 +709,17 @@ direction inverse
 accept  0 0
 expect  0 90
 
+-------------------------------------------------------------------------------
+operation +proj=bonne   +ellps=GRS80  +lat_1=-90
+-------------------------------------------------------------------------------
+tolerance 0.1 mm
+accept  0 -90
+expect  0 0
+
+direction inverse
+accept  0 0
+expect  0 -90
+
 -------------------------------------------------------------------------------
 operation +proj=bonne   +R=6400000    +lat_1=0.5
 -------------------------------------------------------------------------------
@@ -720,6 +743,18 @@ expect  -0.001790562 0.500895246
 accept  -200 -100
 expect  -0.001790561 0.499104753
 
+-------------------------------------------------------------------------------
+operation +proj=bonne   +R=6400000    +lat_1=-0.5
+-------------------------------------------------------------------------------
+tolerance 0.1 mm
+accept  2 1
+expect  223368.1156 167517.5994
+roundtrip   1
+
+accept  2 -1
+expect  223368.1156 -55884.5552
+roundtrip   1
+
 -------------------------------------------------------------------------------
 operation +proj=bonne   +R=6400000  +lat_1=90
 -------------------------------------------------------------------------------
@@ -731,6 +766,17 @@ direction inverse
 accept  0 0
 expect  0 90
 
+-------------------------------------------------------------------------------
+operation +proj=bonne   +R=6400000  +lat_1=-90
+-------------------------------------------------------------------------------
+tolerance 0.1 mm
+accept  0 -90
+expect  0 0
+
+direction inverse
+accept  0 0
+expect  0 -90
+
 ===============================================================================
 # Cal Coop Ocean Fish Invest Lines/Stations
 # 	Cyl, Sph&Ell