Merge 646cfa0 into b4b15d0

docs/source/operations/projections/index.rst

@@ -141,6 +141,7 @@ Projections map the spherical 3D space to a flat 2D space.
    wag5
    wag6
    wag7
+   webmerc
    weren
    wink1
    wink2

docs/source/operations/projections/merc.rst

@@ -24,6 +24,20 @@ The projection is conformal which makes it suitable for navigational purposes.
 +---------------------+----------------------------------------------------------+
 | `+k_0`              | Scaling factor. Defaults to 1.0                          |
 +---------------------+----------------------------------------------------------+
+| `+aux_sphere_type`  | Auxiliary sphere type. Added in proj 5.1.0.              |
+|                     | When specified, spherical formulas are used even if the  |
+|                     | datum is ellipsoidal.                                    |
+|                     |                                                          |
+|                     | * aux_sphere_type=0 will use the semi-major radius.      |
+|                     |   This is the value to use for Web Mercator::            |
+|                     |                                                          |
+|                     |     proj=merc +datum=WGS84 +aux_sphere_type=0            |
+|                     |                                                          |
+|                     |   You can also use the :ref:`webmerc` projection.        |
+|                     | * aux_sphere_type=1 will use the semi-minor radius.      |
+|                     | * aux_sphere_type=2 will calculate the authalic radius   |
+|                     |   and use it.                                            |
++---------------------+----------------------------------------------------------+
 
 
 
@@ -52,6 +66,11 @@ Example using scaling factor::
     echo 56.35 12.32 | proj +proj=merc +k_0=2
     12545706.61     2746073.80
 
+Example using aux_sphere_type (proj >= 5.1.0)::
+
+    $ echo 2 49 | proj +proj=merc +datum=WGS84 +aux_sphere_type=0
+    222638.98       6274861.39
+
 
 Note that ``+lat_ts`` and ``+k_0`` are mutually exclusive.
 If used together, ``+lat_ts`` takes precedence over ``+k_0``.

docs/source/operations/projections/webmerc.rst

@@ -0,0 +1,78 @@
+.. _webmerc:
+
+********************************************************************************
+Web Mercator / Pseudo Mercator
+********************************************************************************
+
+The Web Mercator / Pseudo Mercator projection is a cylindrical map projection.
+This is a variant of the regular :ref:`merc` projection, except that the computation
+is done on a sphere, using the semi-major axis of the ellipsoid.
+
+From `Wikipedia <https://en.wikipedia.org/wiki/Web_Mercator>`_:
+
+    This projection is widely used by the Web Mercator, Google Web Mercator,
+    Spherical Mercator, WGS 84 Web Mercator[1] or WGS 84/Pseudo-Mercator is a
+    variant of the Mercator projection and is the de facto standard for Web
+    mapping applications. [...]
+    It is used by virtually all major online map providers [...]
+    Its official EPSG identifier is EPSG:3857, although others have been used
+    historically.
+
+
++---------------------+----------------------------------------------------------+
+| **Classification**  | Cylindrical (non conformant if used with ellipsoid)      |
++---------------------+----------------------------------------------------------+
+| **Version**         | Added in proj 5.1.0                                      |
++---------------------+----------------------------------------------------------+
+| **Available forms** | Forward and inverse, spherical projection                |
++---------------------+----------------------------------------------------------+
+| **Defined area**    | Global, but best used near the equator                   |
++---------------------+----------------------------------------------------------+
+| **Options**         | Neither lat_0, lat_ts or k_0 should be used.             |
++---------------------+----------------------------------------------------------+
+
+Usage
+########
+
+Example::
+
+    $ echo 2 49 | proj +proj=webmerc +datum=WGS84
+    222638.98       6274861.39
+
+Mathematical definition
+#######################
+
+The formulas describing the Mercator projection are all taken from G. Evenden's libproj manuals [Evenden2005]_.
+
+Forward projection
+==================
+
+.. math::
+
+    x = \lambda
+
+.. math::
+
+    y = \ln \left[ \tan \left(\frac{\pi}{4} + \frac{\phi}{2} \right) \right]
+
+
+Inverse projection
+==================
+
+.. math::
+
+    \lambda = {x}
+
+.. math::
+
+    \phi = \frac{\pi}{2} - 2 \arctan \left[ e^{-y} \right]
+
+
+
+Further reading
+###############
+
+#. `Wikipedia <https://en.wikipedia.org/wiki/Web_Mercator>`_
+
+
+

src/PJ_merc.c

@@ -1,8 +1,10 @@
 #define PJ_LIB__
+#include "proj_internal.h"
 #include "proj.h"
 #include "projects.h"
 
-PROJ_HEAD(merc, "Mercator") "\n\tCyl, Sph&Ell\n\tlat_ts=";
+PROJ_HEAD(merc, "Mercator") "\n\tCyl, Sph&Ell\n\tlat_ts= aux_sphere_type=";
+PROJ_HEAD(webmerc, "Web Mercator / Pseudo Mercator") "\n\tCyl, Sph\n\t";
 
 #define EPS10 1.e-10
 
@@ -48,11 +50,61 @@ static LP s_inverse (XY xy, PJ *P) {           /* Spheroidal, inverse */
     return lp;
 }
 
+/* Cf http://desktop.arcgis.com/en/arcmap/10.3/guide-books/map-projections/mercator.htm */
+enum AuxSphereType
+{
+    AUX_SPHERE_TYPE_SEMI_MAJOR=0,
+    AUX_SPHERE_TYPE_SEMI_MINOR=1,
+    AUX_SPHERE_TYPE_AUTHALIC_RADIUS=2
+    /* AUX_SPHERE_TYPE_AUTHALIC_LATITUDE=3 */ /* unimplemented for now */
+};
 
 PJ *PROJECTION(merc) {
     double phits=0.0;
     int is_phits;
 
+    if( pj_param(P->ctx, P->params, "taux_sphere_type").i ) {
+        int aux_sphere_type = pj_param(P->ctx, P->params, "iaux_sphere_type").i;
+        switch( aux_sphere_type )
+        {
+            case AUX_SPHERE_TYPE_SEMI_MAJOR:
+            {
+                P->b = P->a;
+                /* Clean up the ellipsoidal parameters to reflect the sphere */
+                P->es = P->e = P->f = 0;
+                pj_calc_ellipsoid_params (P, P->a, 0);
+                break;
+            }
+
+            case AUX_SPHERE_TYPE_SEMI_MINOR:
+            {
+                P->a = P->b;
+                /* Clean up the ellipsoidal parameters to reflect the sphere */
+                P->es = P->e = P->f = 0;
+                pj_calc_ellipsoid_params (P, P->a, 0);
+                break;
+            }
+
+            case AUX_SPHERE_TYPE_AUTHALIC_RADIUS:
+            {
+                double qp = pj_qsfn(1.0, P->e, P->one_es);
+                P->a = P->a*sqrt(0.5*qp); /* Set P->a to authalic radius. */
+                P->b = P->a;
+                /* Clean up the ellipsoidal parameters to reflect the sphere */
+                P->es = P->e = P->f = 0;
+                pj_calc_ellipsoid_params (P, P->a, 0);
+                break;
+            }
+
+            default:
+            {
+                proj_log_debug (P, "merc: aux_sphere_type=%d unsupported",
+                                aux_sphere_type);
+                return pj_default_destructor(P, PJD_ERR_INVALID_ARG);
+            }
+        }
+    }
+
     if( (is_phits = pj_param(P->ctx, P->params, "tlat_ts").i) ) {
         phits = fabs(pj_param(P->ctx, P->params, "rlat_ts").f);
         if (phits >= M_HALFPI)
@@ -76,3 +128,21 @@ PJ *PROJECTION(merc) {
     return P;
 }
 
+PJ *PROJECTION(webmerc) {
+
+    if( pj_param(P->ctx, P->params, "tk_0").i ) {
+        return pj_default_destructor(P, PJD_ERR_INVALID_ARG);
+    }
+
+    if( pj_param(P->ctx, P->params, "tlat_0").i ) {
+        return pj_default_destructor(P, PJD_ERR_INVALID_ARG);
+    }
+
+    P->b = P->a;
+    /* Clean up the ellipsoidal parameters to reflect the sphere */
+    P->es = P->e = P->f = 0;
+    pj_calc_ellipsoid_params (P, P->a, 0);
+    P->inv = s_inverse;
+    P->fwd = s_forward;
+    return P;
+}

src/PJ_pipeline.c

@@ -306,6 +306,8 @@ static void set_ellipsoid(PJ *P) {
         /* the PJ you provided".                                         */
         proj_errno_reset (P);
     }
+    P->a_orig = P->a;
+    P->es_orig = P->es;
 
     pj_calc_ellipsoid_params (P, P->a, P->es);
 

src/pj_list.h

@@ -156,6 +156,7 @@ PROJ_HEAD(wag4, "Wagner IV")
 PROJ_HEAD(wag5, "Wagner V")
 PROJ_HEAD(wag6, "Wagner VI")
 PROJ_HEAD(wag7, "Wagner VII")
+PROJ_HEAD(webmerc, "Web Mercator / Pseudo Mercator")
 PROJ_HEAD(weren, "Werenskiold I")
 PROJ_HEAD(wink1, "Winkel I")
 PROJ_HEAD(wink2, "Winkel II")

src/proj_4D_api.c

@@ -512,7 +512,7 @@ Returns 1 on success, 0 on failure
     /* geocentric/cartesian space or we need to do a Helmert transform.         */
     if (P->is_geocent || P->helmert || do_cart) {
         char def[150];
-        sprintf (def, "break_cs2cs_recursion     proj=cart   a=%40.20g  f=%40.20g", P->a, P->f);
+        sprintf (def, "break_cs2cs_recursion     proj=cart   a=%40.20g  es=%40.20g", P->a_orig, P->es_orig);
         Q = proj_create (P->ctx, def);
         if (0==Q)
             return 0;

test/gie/4D-API_cs2cs-style.gie

@@ -215,6 +215,58 @@ accept      487147.594520173    4934316.46263998    0
 expect      -10370728.80        5552839.74          0
 -------------------------------------------------------------------------------
 
+-------------------------------------------------------------------------------
+Test that Google's Web Mercator with +proj=merc +aux_sphere_type=0
+-------------------------------------------------------------------------------
+operation   proj=pipeline step init=epsg:26915 inv step proj=merc datum=WGS84 aux_sphere_type=0
+-------------------------------------------------------------------------------
+tolerance   20 cm
+accept      487147.594520173    4934316.46263998
+expect      -10370728.80        5552839.74
+
+accept      487147.594520173    4934316.46263998    0
+expect      -10370728.80        5552839.74          0
+-------------------------------------------------------------------------------
+
+-------------------------------------------------------------------------------
+Test +proj=merc +aux_sphere_type=1
+-------------------------------------------------------------------------------
+operation   proj=pipeline step init=epsg:26915 inv step proj=merc datum=WGS84 aux_sphere_type=1
+-------------------------------------------------------------------------------
+tolerance   20 cm
+accept      487147.594520173    4934316.46263998
+expect      -10335957.71        5534222.12
+
+accept      487147.594520173    4934316.46263998    0
+expect      -10335957.71        5534222.12          0
+-------------------------------------------------------------------------------
+
+-------------------------------------------------------------------------------
+Test +proj=merc +aux_sphere_type=2
+-------------------------------------------------------------------------------
+operation   proj=pipeline step init=epsg:26915 inv step proj=merc datum=WGS84 aux_sphere_type=2
+-------------------------------------------------------------------------------
+tolerance   20 cm
+accept      487147.594520173    4934316.46263998
+expect      -10359135.85        5546632.48
+
+accept      487147.594520173    4934316.46263998    0
+expect      -10359135.85        5546632.48          0
+-------------------------------------------------------------------------------
+
+-------------------------------------------------------------------------------
+Test that Google's Web Mercator with +proj=webmerc
+-------------------------------------------------------------------------------
+operation   proj=pipeline step init=epsg:26915 inv step proj=webmerc datum=WGS84
+-------------------------------------------------------------------------------
+tolerance   20 cm
+accept      487147.594520173    4934316.46263998
+expect      -10370728.80        5552839.74
+
+accept      487147.594520173    4934316.46263998    0
+expect      -10370728.80        5552839.74          0
+-------------------------------------------------------------------------------
+
 -------------------------------------------------------------------------------
 Test that +datum parameters are handled correctly in pipelines.
 See #872 for details.