Merge pull request #2397 from cffk/merc-update

Update Mercator projection, more accurate, faster

docs/source/operations/projections/merc.rst

@@ -158,7 +158,35 @@ Inverse projection
     \lambda = \frac{x}{k_0 a}; \quad \psi = \frac{y}{k_0 a}
 
 The latitude :math:`\phi` is found by inverting the equation for
-:math:`\psi` iteratively.
+:math:`\psi`.  This follows the method given by :cite:`Karney2011tm`.
+Start by introducing the conformal latitude
+
+.. math::
+   \chi = \tan^{-1}\sinh\psi
+
+The tangents of the latitudes :math:`\tau = \tan\phi` and :math:`\tau' =
+\tan\chi = \sinh\psi` are related by
+
+.. math::
+   \tau' = \tau \sqrt{1 + \sigma^2} - \sigma \sqrt{1 + \tau^2}
+
+where
+
+.. math::
+   \sigma = \sinh\bigl(e \tanh^{-1}(e \tau/\sqrt{1 + \tau^2}) \bigr)
+
+This is obtained by taking the :math:`\sinh` of the equation for
+:math:`\psi` and using the multiple argument formula.  The equation for
+:math:`\tau'` can be solved to give :math:`\tau` using Newton's method
+using :math:`\tau = \tau'/(1 - e^2)` as an initial guess and with the
+needed derivative given by
+
+..math::
+    \frac{d\tau'}{d\tau} = \frac{1 - e^2}{1 + (1 - e^2)\tau^2}
+    \sqrt{1 + \tau'^2} \sqrt{1 + \tau^2}
+
+This converges after no more than 2 iterations.  Finally set
+:math:`\phi=\tan^{-1}\tau`.
 
 Further reading
 ###############

src/apps/gie.cpp

@@ -1124,7 +1124,7 @@ static const struct errno_vs_err_const lookup[] = {
     {"pjd_err_invalid_x_or_y"           ,  -15},
     {"pjd_err_wrong_format_dms_value"   ,  -16},
     {"pjd_err_non_conv_inv_meri_dist"   ,  -17},
-    {"pjd_err_non_con_inv_phi2"         ,  -18},
+    {"pjd_err_non_conv_sinhpsi2tanphi"  ,  -18},
     {"pjd_err_acos_asin_arg_too_large"  ,  -19},
     {"pjd_err_tolerance_condition"      ,  -20},
     {"pjd_err_conic_lat_equal"          ,  -21},

src/phi2.cpp

@@ -1,68 +1,134 @@
 /* Determine latitude angle phi-2. */
 
 #include <math.h>
+#include <limits>
+#include <algorithm>
 
 #include "proj.h"
 #include "proj_internal.h"
 
-static const double TOL = 1.0e-10;
-static const int N_ITER = 15;
+double pj_sinhpsi2tanphi(projCtx ctx, const double taup, const double e) {
+  /****************************************************************************
+  * Convert tau' = sinh(psi) = tan(chi) to tau = tan(phi).  The code is taken
+  * from GeographicLib::Math::tauf(taup, e).
+  *
+  * Here
+  *   phi = geographic latitude (radians)
+  * psi is the isometric latitude
+  *   psi = asinh(tan(phi)) - e * atanh(e * sin(phi))
+  *       = asinh(tan(chi))
+  * chi is the conformal latitude
+  *
+  * The representation of latitudes via their tangents, tan(phi) and tan(chi),
+  * maintains full *relative* accuracy close to latitude = 0 and +/- pi/2.
+  * This is sometimes important, e.g., to compute the scale of the transverse
+  * Mercator projection which involves cos(phi)/cos(chi) tan(phi)
+  *
+  * From Karney (2011), Eq. 7,
+  *
+  *   tau' = sinh(psi) = sinh(asinh(tan(phi)) - e * atanh(e * sin(phi)))
+  *        = tan(phi) * cosh(e * atanh(e * sin(phi))) -
+  *          sec(phi) * sinh(e * atanh(e * sin(phi)))
+  *        = tau * sqrt(1 + sigma^2) - sqrt(1 + tau^2) * sigma
+  * where
+  *   sigma = sinh(e * atanh( e * tau / sqrt(1 + tau^2) ))
+  *
+  * For e small,
+  *
+  *    tau' = (1 - e^2) * tau
+  *
+  * The relation tau'(tau) can therefore by reliably inverted by Newton's
+  * method with
+  *
+  *    tau = tau' / (1 - e^2)
+  *
+  * as an initial guess.  Newton's method requires dtau'/dtau.  Noting that
+  *
+  *   dsigma/dtau = e^2 * sqrt(1 + sigma^2) /
+  *                 (sqrt(1 + tau^2) * (1 + (1 - e^2) * tau^2))
+  *   d(sqrt(1 + tau^2))/dtau = tau / sqrt(1 + tau^2)
+  *
+  * we have
+  *
+  *   dtau'/dtau = (1 - e^2) * sqrt(1 + tau'^2) * sqrt(1 + tau^2) /
+  *                (1 + (1 - e^2) * tau^2)
+  *
+  * This works fine unless tau^2 and tau'^2 overflows.  This may be partially
+  * cured by writing, e.g., sqrt(1 + tau^2) as hypot(1, tau).  However, nan
+  * will still be generated with tau' = inf, since (inf - inf) will appear in
+  * the Newton iteration.
+  *
+  * If we note that for sufficiently large |tau|, i.e., |tau| >= 2/sqrt(eps),
+  * sqrt(1 + tau^2) = |tau| and
+  *
+  *   tau' = exp(- e * atanh(e)) * tau
+  *
+  * So
+  *
+  *   tau = exp(e * atanh(e)) * tau'
+  *
+  * can be returned unless |tau| >= 2/sqrt(eps); this then avoids overflow
+  * problems for large tau' and returns the correct result for tau' = +/-inf
+  * and nan.
+  *
+  * Newton's method usually take 2 iterations to converge to double precision
+  * accuracy (for WGS84 flattening).  However only 1 iteration is needed for
+  * |chi| < 3.35 deg.  In addition, only 1 iteration is needed for |chi| >
+  * 89.18 deg (tau' > 70), if tau = exp(e * atanh(e)) * tau' is used as the
+  * starting guess.
+  ****************************************************************************/
+
+  constexpr int numit = 5;
+  // min iterations = 1, max iterations = 2; mean = 1.954
+  static const double rooteps = sqrt(std::numeric_limits<double>::epsilon());
+  static const double tol = rooteps / 10; // the criterion for Newton's method
+  static const double tmax = 2 / rooteps; // threshold for large arg limit exact
+  const double e2m = 1 - e * e;
+  const double stol = tol * std::max(1.0, fabs(taup));
+  // The initial guess.  70 corresponds to chi = 89.18 deg (see above)
+  double tau = fabs(taup) > 70 ? taup * exp(e * atanh(e)) : taup / e2m;
+  if (!(fabs(tau) < tmax))      // handles +/-inf and nan and e = 1
+    return tau;
+  // If we need to deal with e > 1, then we could include:
+  // if (e2m < 0) return std::numeric_limits<double>::quiet_NaN();
+  int i = numit;
+  for (; i; --i) {
+    double tau1 = sqrt(1 + tau * tau);
+    double sig = sinh( e * atanh(e * tau / tau1) );
+    double taupa = sqrt(1 + sig * sig) * tau - sig * tau1;
+    double dtau = ( (taup - taupa) * (1 + e2m * (tau * tau)) /
+                    (e2m * tau1 * sqrt(1 + taupa * taupa)) );
+    tau += dtau;
+    if (!(fabs(dtau) >= stol))  // backwards test to allow nans to succeed.
+      break;
+  }
+  if (i == 0)
+    pj_ctx_set_errno(ctx, PJD_ERR_NON_CONV_SINHPSI2TANPHI);
+  return tau;
+}
 
 /*****************************************************************************/
 double pj_phi2(projCtx ctx, const double ts0, const double e) {
-/******************************************************************************
-Determine latitude angle phi-2.
-Inputs:
-  ts = exp(-psi) where psi is the isometric latitude (dimensionless)
-  e = eccentricity of the ellipsoid (dimensionless)
-Output:
-  phi = geographic latitude (radians)
-Here isometric latitude is defined by
-  psi = log( tan(pi/4 + phi/2) *
-             ( (1 - e*sin(phi)) / (1 + e*sin(phi)) )^(e/2) )
-      = asinh(tan(phi)) - e * atanh(e * sin(phi))
-This routine inverts this relation using the iterative scheme given
-by Snyder (1987), Eqs. (7-9) - (7-11)
-*******************************************************************************/
-    const double eccnth = .5 * e;
-    double ts = ts0;
-#ifdef no_longer_used_original_convergence_on_exact_dphi
-    double Phi = M_HALFPI - 2. * atan(ts);
-#endif
-    int i = N_ITER;
-
-    for(;;) {
-        /*
-         * sin(Phi) = sin(PI/2 - 2* atan(ts))
-         *          = cos(2*atan(ts))
-         *          = 2*cos^2(atan(ts)) - 1
-         *          = 2 / (1 + ts^2) - 1
-         *          = (1 - ts^2) / (1 + ts^2)
-         */
-        const double sinPhi = (1 - ts * ts) / (1 + ts * ts);
-        const double con = e * sinPhi;
-        double old_ts = ts;
-        ts = ts0 * pow((1. - con) / (1. + con), eccnth);
-#ifdef no_longer_used_original_convergence_on_exact_dphi
-        /* The convergence criterion is nominally on exact dphi */
-        const double newPhi = M_HALFPI - 2. * atan(ts);
-        const double dphi = newPhi - Phi;
-        Phi = newPhi;
-#else
-        /* If we don't immediately update Phi from this, we can
-         * change the conversion criterion to save us computing atan() at each step.
-         * Particularly we can observe that:
-         * |atan(ts) - atan(old_ts)| <= |ts - old_ts|
-         * So if |ts - old_ts| matches our convergence criterion, we're good.
-         */
-        const double dphi = 2 * (ts - old_ts);
-#endif
-        if (fabs(dphi) > TOL && --i) {
-            continue;
-        }
-        break;
-    }
-    if (i <= 0)
-        pj_ctx_set_errno(ctx, PJD_ERR_NON_CON_INV_PHI2);
-    return M_HALFPI - 2. * atan(ts);
+  /****************************************************************************
+   * Determine latitude angle phi-2.
+   * Inputs:
+   *   ts = exp(-psi) where psi is the isometric latitude (dimensionless)
+   *        this variable is defined in Snyder (1987), Eq. (7-10)
+   *   e = eccentricity of the ellipsoid (dimensionless)
+   * Output:
+   *   phi = geographic latitude (radians)
+   * Here isometric latitude is defined by
+   *   psi = log( tan(pi/4 + phi/2) *
+   *              ( (1 - e*sin(phi)) / (1 + e*sin(phi)) )^(e/2) )
+   *       = asinh(tan(phi)) - e * atanh(e * sin(phi))
+   *       = asinh(tan(chi))
+   *   chi = conformal latitude
+   *
+   * This routine converts t = exp(-psi) to
+   *
+   *   tau' = tan(chi) = sinh(psi) = (1/t - t)/2
+   *
+   * returns atan(sinpsi2tanphi(tau'))
+   ***************************************************************************/
+  return atan(pj_sinhpsi2tanphi(ctx, (1/ts0 - ts0) / 2, e));
 }

src/proj_internal.h

@@ -626,7 +626,7 @@ struct FACTORS {
 #define PJD_ERR_INVALID_X_OR_Y          -15
 #define PJD_ERR_WRONG_FORMAT_DMS_VALUE  -16
 #define PJD_ERR_NON_CONV_INV_MERI_DIST  -17
-#define PJD_ERR_NON_CON_INV_PHI2        -18
+#define PJD_ERR_NON_CONV_SINHPSI2TANPHI -18
 #define PJD_ERR_ACOS_ASIN_ARG_TOO_LARGE -19
 #define PJD_ERR_TOLERANCE_CONDITION     -20
 #define PJD_ERR_CONIC_LAT_EQUAL         -21
@@ -846,6 +846,7 @@ double  pj_qsfn(double, double, double);
 double  pj_tsfn(double, double, double);
 double  pj_msfn(double, double, double);
 double  PROJ_DLL pj_phi2(projCtx_t *, const double, const double);
+double  pj_sinhpsi2tanphi(projCtx_t *, const double, const double);
 double  pj_qsfn_(double, PJ *);
 double *pj_authset(double);
 double  pj_authlat(double, double *);

src/projections/gstmerc.cpp

@@ -28,9 +28,9 @@ static PJ_XY gstmerc_s_forward (PJ_LP lp, PJ *P) {           /* Spheroidal, forw
     double L, Ls, sinLs1, Ls1;
 
     L = Q->n1*lp.lam;
-    Ls = Q->c + Q->n1 * log(pj_tsfn(-1.0 * lp.phi, -1.0 * sin(lp.phi), P->e));
+    Ls = Q->c + Q->n1 * log(pj_tsfn(-lp.phi, -sin(lp.phi), P->e));
     sinLs1 = sin(L) / cosh(Ls);
-    Ls1 = log(pj_tsfn(-1.0 * asin(sinLs1), 0.0, 0.0));
+    Ls1 = log(pj_tsfn(-asin(sinLs1), -sinLs1, 0.0));
     xy.x = (Q->XS + Q->n2*Ls1) * P->ra;
     xy.y = (Q->YS + Q->n2*atan(sinh(Ls) / cos(L))) * P->ra;
 
@@ -45,9 +45,9 @@ static PJ_LP gstmerc_s_inverse (PJ_XY xy, PJ *P) {           /* Spheroidal, inve
 
     L = atan(sinh((xy.x * P->a - Q->XS) / Q->n2) / cos((xy.y * P->a - Q->YS) / Q->n2));
     sinC = sin((xy.y * P->a - Q->YS) / Q->n2) / cosh((xy.x * P->a - Q->XS) / Q->n2);
-    LC = log(pj_tsfn(-1.0 * asin(sinC), 0.0, 0.0));
+    LC = log(pj_tsfn(-asin(sinC), -sinC, 0.0));
     lp.lam = L / Q->n1;
-    lp.phi = -1.0 * pj_phi2(P->ctx, exp((LC - Q->c) / Q->n1), P->e);
+    lp.phi = -pj_phi2(P->ctx, exp((LC - Q->c) / Q->n1), P->e);
 
     return lp;
 }
@@ -60,13 +60,13 @@ PJ *PROJECTION(gstmerc) {
     P->opaque = Q;
 
     Q->lamc = P->lam0;
-    Q->n1 = sqrt(1.0 + P->es * pow(cos(P->phi0), 4.0) / (1.0 - P->es));
+    Q->n1 = sqrt(1 + P->es * pow(cos(P->phi0), 4.0) / (1 - P->es));
     Q->phic = asin(sin(P->phi0) / Q->n1);
-    Q->c = log(pj_tsfn(-1.0 * Q->phic, 0.0, 0.0))
-         - Q->n1 * log(pj_tsfn(-1.0 * P->phi0, -1.0 * sin(P->phi0), P->e));
-    Q->n2 = P->k0 * P->a * sqrt(1.0 - P->es) / (1.0 - P->es * sin(P->phi0) * sin(P->phi0));
+    Q->c = log(pj_tsfn(-Q->phic, -sin(P->phi0) / Q->n1, 0.0))
+         - Q->n1 * log(pj_tsfn(-P->phi0, -sin(P->phi0), P->e));
+    Q->n2 = P->k0 * P->a * sqrt(1 - P->es) / (1 - P->es * sin(P->phi0) * sin(P->phi0));
     Q->XS = 0;
-    Q->YS = -1.0 * Q->n2 * Q->phic;
+    Q->YS = -Q->n2 * Q->phic;
 
     P->inv = gstmerc_s_inverse;
     P->fwd = gstmerc_s_forward;

src/projections/lcc.cpp

@@ -106,21 +106,21 @@ PJ *PROJECTION(lcc) {
         double ml1, m1;
 
         m1 = pj_msfn(sinphi, cosphi, P->es);
-        ml1 = pj_tsfn(Q->phi1, sinphi, P->e);
-        if( ml1 == 0 ) {
+        if( fabs(Q->phi1) == M_HALFPI ) {
             return pj_default_destructor(P, PJD_ERR_LAT_1_OR_2_ZERO_OR_90);
         }
+        ml1 = pj_tsfn(Q->phi1, sinphi, P->e);
         if (secant) { /* secant cone */
             sinphi = sin(Q->phi2);
             Q->n = log(m1 / pj_msfn(sinphi, cos(Q->phi2), P->es));
             if (Q->n == 0) {
                 // Not quite, but es is very close to 1...
                 return pj_default_destructor(P, PJD_ERR_INVALID_ECCENTRICITY);
             }
-            const double ml2 = pj_tsfn(Q->phi2, sinphi, P->e);
-            if( ml2 == 0 ) {
-                return pj_default_destructor(P, PJD_ERR_LAT_1_OR_2_ZERO_OR_90);
+            if( fabs(Q->phi2) == M_HALFPI ) {
+              return pj_default_destructor(P, PJD_ERR_LAT_1_OR_2_ZERO_OR_90);
             }
+            const double ml2 = pj_tsfn(Q->phi2, sinphi, P->e);
             const double denom = log(ml1 / ml2);
             if( denom == 0 ) {
                 // Not quite, but es is very close to 1...

src/projections/merc.cpp

@@ -10,45 +10,29 @@
 PROJ_HEAD(merc, "Mercator") "\n\tCyl, Sph&Ell\n\tlat_ts=";
 PROJ_HEAD(webmerc, "Web Mercator / Pseudo Mercator") "\n\tCyl, Ell\n\t";
 
-#define EPS10 1.e-10
-static double logtanpfpim1(double x) {       /* log(tan(x/2 + M_FORTPI)) */
-    if (fabs(x) <= DBL_EPSILON) {
-        /* tan(M_FORTPI + .5 * x) can be approximated by  1.0 + x */
-        return log1p(x);
-    }
-    return log(tan(M_FORTPI + .5 * x));
-}
-
 static PJ_XY merc_e_forward (PJ_LP lp, PJ *P) {          /* Ellipsoidal, forward */
     PJ_XY xy = {0.0,0.0};
-    if (fabs(fabs(lp.phi) - M_HALFPI) <= EPS10) {
-        proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION);
-        return xy;
-    }
     xy.x = P->k0 * lp.lam;
-    xy.y = - P->k0 * log(pj_tsfn(lp.phi, sin(lp.phi), P->e));
+    // Instead of calling tan and sin, call sin and cos which the compiler
+    // optimizes to a single call to sincos.
+    double sphi = sin(lp.phi);
+    double cphi = cos(lp.phi);
+    xy.y = P->k0 * (asinh(sphi/cphi) - P->e * atanh(P->e * sphi));
     return xy;
 }
 
 
 static PJ_XY merc_s_forward (PJ_LP lp, PJ *P) {           /* Spheroidal, forward */
     PJ_XY xy = {0.0,0.0};
-    if (fabs(fabs(lp.phi) - M_HALFPI) <= EPS10) {
-        proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION);
-        return xy;
-}
     xy.x = P->k0 * lp.lam;
-    xy.y = P->k0 * logtanpfpim1(lp.phi);
+    xy.y = P->k0 * asinh(tan(lp.phi));
     return xy;
 }
 
 
 static PJ_LP merc_e_inverse (PJ_XY xy, PJ *P) {          /* Ellipsoidal, inverse */
     PJ_LP lp = {0.0,0.0};
-    if ((lp.phi = pj_phi2(P->ctx, exp(- xy.y / P->k0), P->e)) == HUGE_VAL) {
-        proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION);
-        return lp;
-}
+    lp.phi = atan(pj_sinhpsi2tanphi(P->ctx, sinh(xy.y / P->k0), P->e));
     lp.lam = xy.x / P->k0;
     return lp;
 }