From 1a575b61a18eb480d4482c843e4a42db66c81ba5 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Sat, 7 Apr 2018 15:35:08 +0200 Subject: [PATCH 1/3] Pipeline: make sure geocentric/cartesian space transform is done with original ellipsoid parameters (before any projection mess with them) --- src/PJ_pipeline.c | 2 ++ src/proj_4D_api.c | 2 +- 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/src/PJ_pipeline.c b/src/PJ_pipeline.c index b60e05a89a..aa7d76f86b 100644 --- a/src/PJ_pipeline.c +++ b/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); diff --git a/src/proj_4D_api.c b/src/proj_4D_api.c index 11a56ac087..da6a7d1b5b 100644 --- a/src/proj_4D_api.c +++ b/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; From 65d294ea144caceacd8b6e08a1199a22659f5bda Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Sat, 7 Apr 2018 12:24:31 +0200 Subject: [PATCH 2/3] Mercator projections: add a aux_sphere_type=0/1/2 parameter (cleaner modeling of WebMercator) Directly based on the auxiliary sphere type parameter of http://desktop.arcgis.com/en/arcmap/10.3/guide-books/map-projections/mercator.htm Implemented: - 0: use semimajor axis or radius of the geographic coordinate system (WebMercator style) - 1: use semiminor axis or radius - 2: calculate and use authalic radius Unimplemented: - 3: use authalic radius and convert geodetic latitudes to authalic latitudes --- docs/source/operations/projections/merc.rst | 17 +++++++ src/PJ_merc.c | 53 ++++++++++++++++++++- test/gie/4D-API_cs2cs-style.gie | 39 +++++++++++++++ 3 files changed, 108 insertions(+), 1 deletion(-) diff --git a/docs/source/operations/projections/merc.rst b/docs/source/operations/projections/merc.rst index 8895abc8e0..9c6b692ff0 100644 --- a/docs/source/operations/projections/merc.rst +++ b/docs/source/operations/projections/merc.rst @@ -24,6 +24,18 @@ 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 Google Web Mercator:: | +| | | +| | proj=merc +datum=WGS84 +aux_sphere_type=0 | +| | * aux_sphere_type=1 will use the semi-minor radius. | +| | * aux_sphere_type=2 will calculate the authalic radius | +| | and use it. | ++---------------------+----------------------------------------------------------+ @@ -52,6 +64,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``. diff --git a/src/PJ_merc.c b/src/PJ_merc.c index 18303c4586..a7efc2f6ce 100644 --- a/src/PJ_merc.c +++ b/src/PJ_merc.c @@ -1,8 +1,9 @@ #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="; #define EPS10 1.e-10 @@ -48,11 +49,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) diff --git a/test/gie/4D-API_cs2cs-style.gie b/test/gie/4D-API_cs2cs-style.gie index bbd9ee3918..bcbc451614 100644 --- a/test/gie/4D-API_cs2cs-style.gie +++ b/test/gie/4D-API_cs2cs-style.gie @@ -215,6 +215,45 @@ 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 +datum parameters are handled correctly in pipelines. See #872 for details. From 646cfa0c928f1d9d987831739e6fd5d648b68497 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Sat, 7 Apr 2018 13:41:12 +0200 Subject: [PATCH 3/3] Add webmerc projection --- docs/source/operations/projections/index.rst | 1 + docs/source/operations/projections/merc.rst | 4 +- .../source/operations/projections/webmerc.rst | 78 +++++++++++++++++++ src/PJ_merc.c | 19 +++++ src/pj_list.h | 1 + test/gie/4D-API_cs2cs-style.gie | 13 ++++ 6 files changed, 115 insertions(+), 1 deletion(-) create mode 100644 docs/source/operations/projections/webmerc.rst diff --git a/docs/source/operations/projections/index.rst b/docs/source/operations/projections/index.rst index 9255d9fdb5..77decd6c37 100644 --- a/docs/source/operations/projections/index.rst +++ b/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 diff --git a/docs/source/operations/projections/merc.rst b/docs/source/operations/projections/merc.rst index 9c6b692ff0..964cd5d5c3 100644 --- a/docs/source/operations/projections/merc.rst +++ b/docs/source/operations/projections/merc.rst @@ -29,9 +29,11 @@ The projection is conformal which makes it suitable for navigational purposes. | | datum is ellipsoidal. | | | | | | * aux_sphere_type=0 will use the semi-major radius. | -| | This is the value to use for Google Web Mercator:: | +| | 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. | diff --git a/docs/source/operations/projections/webmerc.rst b/docs/source/operations/projections/webmerc.rst new file mode 100644 index 0000000000..d475455649 --- /dev/null +++ b/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 `_: + + 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 `_ + + + diff --git a/src/PJ_merc.c b/src/PJ_merc.c index a7efc2f6ce..d20a09e1cc 100644 --- a/src/PJ_merc.c +++ b/src/PJ_merc.c @@ -4,6 +4,7 @@ #include "projects.h" 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 @@ -127,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; +} diff --git a/src/pj_list.h b/src/pj_list.h index 09a66034da..440f797298 100644 --- a/src/pj_list.h +++ b/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") diff --git a/test/gie/4D-API_cs2cs-style.gie b/test/gie/4D-API_cs2cs-style.gie index bcbc451614..7f3e0c4e30 100644 --- a/test/gie/4D-API_cs2cs-style.gie +++ b/test/gie/4D-API_cs2cs-style.gie @@ -254,6 +254,19 @@ 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.