Inverse tmerc spherical: fix wrong sign of latitude when lat_0 is used (fixes #2468)

Corrected formula given by @evanmiller

2 files changed, 128 insertions, 2 deletions

diff --git a/src/projections/tmerc.cpp b/src/projections/tmerc.cpp
index 8cae9968..2d13b064 100644
--- a/src/projections/tmerc.cpp
+++ b/src/projections/tmerc.cpp
@@ -203,11 +203,13 @@ static PJ_LP approx_s_inv (PJ_XY xy, PJ *P) {
         return proj_coord_error().lp;
     }
     g = .5 * (h - 1. / h);
-    h = cos (P->phi0 + xy.y / Q->esp);
+    /* D, as in equation 8-8 of USGS "Map Projections - A Working Manual" */
+    const double D = P->phi0 + xy.y / Q->esp;
+    h = cos (D);
     lp.phi = asin(sqrt((1. - h * h) / (1. + g * g)));
 
     /* Make sure that phi is on the correct hemisphere when false northing is used */
-    if (xy.y < 0. && -lp.phi+P->phi0 < 0.0) lp.phi = -lp.phi;
+    lp.phi = copysign(lp.phi, D);
 
     lp.lam = (g != 0.0 || h != 0.0) ? atan2 (g, h) : 0.;
     return lp;
diff --git a/test/gie/builtins.gie b/test/gie/builtins.gie
index add5d925..0f86bed2 100644
--- a/test/gie/builtins.gie
+++ b/test/gie/builtins.gie
@@ -5782,6 +5782,130 @@ accept  -200 -100
 expect  -0.001790493 -0.000895247
 
 -------------------------------------------------------------------------------
+# Approx tmerc on a almost spherical ellipsoid, lat_0 north hemisphere
+
+operation +proj=tmerc   +a=6400000 +rf=1e12 +k=0.9 +lat_0=40 +approx
+-------------------------------------------------------------------------------
+tolerance 0.1 mm
+
+accept    0             -30
+expect    0             -7037167.5440
+roundtrip 1
+
+accept    1             -30
+expect    87064.5795    -7037547.4590
+roundtrip 1
+
+accept    -1            -30
+expect    -87064.5795   -7037547.4590
+roundtrip 1
+
+accept    0             30
+expect    0             -1005309.6491
+roundtrip 1
+
+accept    0             40
+expect    0             0
+roundtrip 1
+
+accept    1             41
+expect    75872.2182    100965.3718
+roundtrip 1
+
+-------------------------------------------------------------------------------
+# Approx tmerc on a sphere, lat_0 north hemisphere
+
+operation +proj=tmerc   +R=6400000 +k=0.9 +lat_0=40
+-------------------------------------------------------------------------------
+tolerance 0.1 mm
+
+accept    0             -30
+expect    0             -7037167.5440
+roundtrip 1
+
+accept    1             -30
+expect    87064.5795    -7037547.4590
+roundtrip 1
+
+accept    -1            -30
+expect    -87064.5795   -7037547.4590
+roundtrip 1
+
+accept    0             30
+expect    0             -1005309.6491
+roundtrip 1
+
+accept    0             40
+expect    0             0
+roundtrip 1
+
+accept    1             41
+expect    75872.2182    100965.3718
+roundtrip 1
+
+-------------------------------------------------------------------------------
+# Approx tmerc on a almost spherical ellipsoid, lat_0 south hemisphere
+
+operation +proj=tmerc   +a=6400000 +rf=1e12 +k=0.9 +lat_0=-40 +approx
+-------------------------------------------------------------------------------
+tolerance 0.1 mm
+
+accept    0             -30
+expect    0             1005309.6491
+roundtrip 1
+
+accept    1             -30
+expect    87064.5795    1004929.7341
+roundtrip 1
+
+accept    -1            -30
+expect    -87064.5795   1004929.7341
+roundtrip 1
+
+accept    0             30
+expect    0             7037167.5440
+roundtrip 1
+
+accept    0             -40
+expect    0             0
+roundtrip 1
+
+accept    1             -41
+expect    75872.2182    -100965.3718
+roundtrip 1
+
+-------------------------------------------------------------------------------
+# Approx tmerc on a sphere, lat_0 south hemisphere
+
+operation +proj=tmerc   +R=6400000 +k=0.9 +lat_0=-40
+-------------------------------------------------------------------------------
+tolerance 0.1 mm
+
+accept    0             -30
+expect    0             1005309.6491
+roundtrip 1
+
+accept    1             -30
+expect    87064.5795    1004929.7341
+roundtrip 1
+
+accept    -1            -30
+expect    -87064.5795   1004929.7341
+roundtrip 1
+
+accept    0             30
+expect    0             7037167.5440
+roundtrip 1
+
+accept    0             -40
+expect    0             0
+roundtrip 1
+
+accept    1             -41
+expect    75872.2182    -100965.3718
+roundtrip 1
+
+-------------------------------------------------------------------------------
 operation +proj=tmerc   +R=1
 -------------------------------------------------------------------------------
 direction inverse