diff --git a/docs/source/operations/projections/eqc.rst b/docs/source/operations/projections/eqc.rst index 47184752fe..8d2c2965a6 100644 --- a/docs/source/operations/projections/eqc.rst +++ b/docs/source/operations/projections/eqc.rst @@ -6,10 +6,13 @@ Equidistant Cylindrical (Plate Carrée) The simplest of all projections. Standard parallels (0° when omitted) may be specified that determine latitude of true scale (k=h=1). +This projection is available in both spherical (EPSG:1029) and ellipsoidal (EPSG:1028) forms. +When using an ellipsoid, PROJ uses the ellipsoidal formulas from IOGP Guidance Note 7-2. + +---------------------+----------------------------------------------------------+ | **Classification** | Conformal cylindrical | +---------------------+----------------------------------------------------------+ -| **Available forms** | Forward and inverse | +| **Available forms** | Forward and inverse, spherical and ellipsoidal | +---------------------+----------------------------------------------------------+ | **Defined area** | Global, but best used near the equator | +---------------------+----------------------------------------------------------+ @@ -92,16 +95,19 @@ Parameters Mathematical definition ####################### -The formulas describing the Equidistant Cylindrical projection are all taken from :cite:`Snyder1987`. - :math:`\phi_{ts}` is the latitude of true scale, i.e., the standard parallel where the scale of the projection is true. It can be set with ``+lat_ts``. :math:`\phi_0` is the latitude of origin that match the center of the map. It can be set with ``+lat_0``. -Forward projection +Spherical formulas ================== +The spherical formulas are taken from :cite:`Snyder1987` and correspond to EPSG method 1029. + +Forward projection +------------------ + .. math:: x = \lambda \cos \phi_{ts} @@ -111,21 +117,58 @@ Forward projection y = \phi - \phi_0 Inverse projection -================== +------------------ .. math:: - \lambda = x / cos \phi_{ts} + \lambda = x / \cos \phi_{ts} .. math:: \phi = y + \phi_0 +Ellipsoidal formulas +==================== + +The ellipsoidal formulas are from IOGP Guidance Note 7-2, Section 3.2.5, and correspond to EPSG method 1028. + +Forward projection +------------------ + +.. math:: + + x = \nu_1 \cos \phi_{ts} \cdot \lambda + +.. math:: + + y = M(\phi) - M(\phi_0) + +where :math:`\nu_1` is the radius of curvature in the prime vertical at the standard parallel: + +.. math:: + + \nu_1 = \frac{a}{\sqrt{1 - e^2 \sin^2 \phi_{ts}}} + +and :math:`M(\phi)` is the meridional arc distance from the equator to latitude :math:`\phi`. + +Inverse projection +------------------ + +.. math:: + + \lambda = x / (\nu_1 \cos \phi_{ts}) + +.. math:: + + \phi = M^{-1}(y + M(\phi_0)) + + Further reading ############### #. `Wikipedia `_ #. `Wolfram Mathworld `_ +#. `IOGP Guidance Note 7-2 `_ - Section 3.2.5 diff --git a/src/projections/eqc.cpp b/src/projections/eqc.cpp index c5564521dc..2d7adccf7d 100644 --- a/src/projections/eqc.cpp +++ b/src/projections/eqc.cpp @@ -1,19 +1,42 @@ +/****************************************************************************** + * Project: PROJ + * Purpose: Implementation of the eqc (Equidistant Cylindrical) projection. + * Also known as Plate Carree when lat_ts=0. + * + * This file implements both the spherical (EPSG:1029) and ellipsoidal + * (EPSG:1028) forms of the Equidistant Cylindrical projection. + * + * Spherical formulas (EPSG:1029): + * E = FE + R cos(lat_ts) (lon - lon_0) + * N = FN + R lat + * + * Ellipsoidal formulas (EPSG:1028) from IOGP Guidance Note 7-2: + * E = FE + nu1 cos(lat_ts) (lon - lon_0) + * N = FN + M + * + * + * Reference: IOGP Publication 373-7-2, Section 3.2.5 + * + *****************************************************************************/ #include #include +#include #include "proj.h" #include "proj_internal.h" namespace { // anonymous namespace struct pj_eqc_data { - double rc; + double rc; // Spherical: cos(lat_ts); Ellipsoidal: nu1 * cos(lat_ts) + double M0; // Meridional arc at latitude of origin (lat_0) + double *en; // Coefficients for meridional arc computation }; } // anonymous namespace PROJ_HEAD(eqc, "Equidistant Cylindrical (Plate Carree)") -"\n\tCyl, Sph\n\tlat_ts=[, lat_0=0]"; +"\n\tCyl, Sph&Ell\n\tlat_ts=[, lat_0=0]"; static PJ_XY eqc_s_forward(PJ_LP lp, PJ *P) { /* Spheroidal, forward */ PJ_XY xy = {0.0, 0.0}; @@ -35,21 +58,88 @@ static PJ_LP eqc_s_inverse(PJ_XY xy, PJ *P) { /* Spheroidal, inverse */ return lp; } +// Ellipsoidal forward (EPSG method 1028) +// E = FE + nu1 cos(lat_ts) (lon - lon_0) +// N = FN + M - M0 +// where M is the meridional arc from equator to latitude +static PJ_XY eqc_e_forward(PJ_LP lp, PJ *P) { /* Ellipsoidal, forward */ + PJ_XY xy = {0.0, 0.0}; + struct pj_eqc_data *Q = static_cast(P->opaque); + + const double sinphi = sin(lp.phi); + const double cosphi = cos(lp.phi); + + xy.x = Q->rc * lp.lam; + xy.y = pj_mlfn(lp.phi, sinphi, cosphi, Q->en) - Q->M0; + + return xy; +} + +// Ellipsoidal inverse (EPSG method 1028) +// lon = lon_0 + (E - FE) / (nu1 cos lat_ts) +// lat is found by inverting M = N - FN + M0 +static PJ_LP eqc_e_inverse(PJ_XY xy, PJ *P) { /* Ellipsoidal, inverse */ + PJ_LP lp = {0.0, 0.0}; + struct pj_eqc_data *Q = static_cast(P->opaque); + + lp.lam = xy.x / Q->rc; + lp.phi = pj_inv_mlfn(xy.y + Q->M0, Q->en); + + return lp; +} + +static PJ *pj_eqc_destructor(PJ *P, int errlev) { /* Destructor */ + if (nullptr == P) + return nullptr; + + if (nullptr == P->opaque) + return pj_default_destructor(P, errlev); + + free(static_cast(P->opaque)->en); + return pj_default_destructor(P, errlev); +} + PJ *PJ_PROJECTION(eqc) { struct pj_eqc_data *Q = static_cast( calloc(1, sizeof(struct pj_eqc_data))); if (nullptr == Q) return pj_default_destructor(P, PROJ_ERR_OTHER /*ENOMEM*/); P->opaque = Q; + P->destructor = pj_eqc_destructor; - if ((Q->rc = cos(pj_param(P->ctx, P->params, "rlat_ts").f)) <= 0.) { + const double phi1 = pj_param(P->ctx, P->params, "rlat_ts").f; + const double cos_phi1 = cos(phi1); + + if (cos_phi1 <= 0.) { proj_log_error( P, _("Invalid value for lat_ts: |lat_ts| should be <= 90°")); return pj_default_destructor(P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE); } - P->inv = eqc_s_inverse; - P->fwd = eqc_s_forward; - P->es = 0.; + + if (P->es != 0.0) { + // Ellipsoidal case (EPSG:1028) + const double sin_phi1 = sin(phi1); + // nu1 = a / sqrt(1 - e^2 sin^2(lat_ts)), normalized by a + const double nu1 = 1.0 / sqrt(1.0 - P->es * sin_phi1 * sin_phi1); + Q->rc = nu1 * cos_phi1; + + Q->en = pj_enfn(P->n); + if (nullptr == Q->en) + return pj_default_destructor(P, PROJ_ERR_OTHER /*ENOMEM*/); + + Q->M0 = pj_mlfn(P->phi0, sin(P->phi0), cos(P->phi0), Q->en); + + P->inv = eqc_e_inverse; + P->fwd = eqc_e_forward; + } else { + // Spheroidal case (EPSG:1029) + Q->rc = cos_phi1; + Q->en = nullptr; + Q->M0 = 0.0; + + P->inv = eqc_s_inverse; + P->fwd = eqc_s_forward; + } return P; } diff --git a/test/cli/test_cs2cs_ignf.yaml b/test/cli/test_cs2cs_ignf.yaml index 0b36df8793..51f52fe2d0 100644 --- a/test/cli/test_cs2cs_ignf.yaml +++ b/test/cli/test_cs2cs_ignf.yaml @@ -35,52 +35,52 @@ tests: out: "1238837.253\t5057451.037 0.000" - args: +init=IGNF:LAMBE +to +init=IGNF:GEOPORTALFXX -f %.3f in: 600000.0000 2600545.4523 0.0000 - out: "179040.148\t5610495.275 0.000" + out: "179356.306\t5585333.153 0.000" - args: +init=IGNF:LAMBE +to +init=IGNF:GEOPORTALFXX -f %.3f in: 135638.3592 2418760.4094 0.0000 - out: "-303729.363\t5410118.356 0.000" + out: "-304265.704\t5385136.401 0.000" - args: +init=IGNF:LAMBE +to +init=IGNF:GEOPORTALFXX -f %.3f in: 998137.3947 2413822.2844 0.0000 - out: "592842.792\t5410120.554 0.000" + out: "593889.664\t5385138.596 0.000" - args: +init=IGNF:LAMBE +to +init=IGNF:GEOPORTALFXX -f %.3f in: 600000.0000 2200000.0000 0.0000 - out: "179041.670\t5209746.080 0.000" + out: "179357.831\t5185007.149 0.000" - args: +init=IGNF:LAMBE +to +init=IGNF:GEOPORTALFXX -f %.3f in: 311552.5340 1906457.4840 0.0000 - out: "-96825.465\t4909184.136 0.000" + out: "-96996.444\t4884928.293 0.000" - args: +init=IGNF:LAMBE +to +init=IGNF:GEOPORTALFXX -f %.3f in: 960488.4138 1910172.8812 0.0000 - out: "523880.019\t4909191.141 0.000" + out: "524805.113\t4884935.286 0.000" - args: +init=IGNF:LAMBE +to +init=IGNF:GEOPORTALFXX -f %.3f in: 600000.0000 1699510.8340 0.0000 - out: "179047.633\t4708817.007 0.000" + out: "179363.804\t4684962.279 0.000" - args: +init=IGNF:LAMBE +to +init=IGNF:GEOPORTALFXX -f %.3f in: 1203792.5981 626873.17210 0.0000 - out: "658259.467\t3623786.764 0.000" + out: "659421.855\t3603178.891 0.000" - args: +init=IGNF:RGF93G +to +init=IGNF:GEOPORTALFXX -f %.3f in: 2d20'11.4239243" 50d23'59.7718445" 0.0 - out: "179040.151\t5610495.281 0.000" + out: "179356.309\t5585333.158 0.000" - args: +init=IGNF:RGF93G +to +init=IGNF:GEOPORTALFXX -f %.3f in: -3d57'49.4051448" 48d35'59.7121716" 0.0 - out: "-303729.365\t5410118.352 0.000" + out: "-304265.706\t5385136.397 0.000" - args: +init=IGNF:RGF93G +to +init=IGNF:GEOPORTALFXX -f %.3f in: 7d44'12.1439796" 48d35'59.7832558" 0.0 - out: "592842.794\t5410120.550 0.000" + out: "593889.666\t5385138.593 0.000" - args: +init=IGNF:RGF93G +to +init=IGNF:GEOPORTALFXX -f %.3f in: 2d20'11.4951975" 46d47'59.8029841" 0.0 - out: "179041.668\t5209746.077 0.000" + out: "179357.829\t5185007.147 0.000" - args: +init=IGNF:RGF93G +to +init=IGNF:GEOPORTALFXX -f %.3f in: -1d15'48.9240599" 44d05'59.8251878" 0.0 - out: "-96825.467\t4909184.138 0.000" + out: "-96996.446\t4884928.296 0.000" - args: +init=IGNF:RGF93G +to +init=IGNF:GEOPORTALFXX -f %.3f in: 6d50'12.2276489" 44d06'00.0517019" 0.0 - out: "523880.022\t4909191.143 0.000" + out: "524805.116\t4884935.287 0.000" - args: +init=IGNF:RGF93G +to +init=IGNF:GEOPORTALFXX -f %.3f in: 2d20'11.7754730" 42d18'00.0824436" 0.0 - out: "179047.634\t4708817.010 0.000" + out: "179363.805\t4684962.283 0.000" - args: +init=IGNF:RGF93G +to +init=IGNF:GEOPORTALFXX -f %.3f in: 9d32'12.6680218" 41d24'00.3542556" 0.0 - out: "730783.054\t4608637.873 0.000" + out: "732073.508\t4585007.331 0.000" - args: +init=IGNF:RGF93G +to +init=IGNF:MILLER -f %.3f in: 2d20'11.4239243" 50d23'59.7718445" 0.0 out: "260098.730\t6140682.441 0.000" diff --git a/test/gie/builtins.gie b/test/gie/builtins.gie index a301cce7af..f85bf4e06a 100644 --- a/test/gie/builtins.gie +++ b/test/gie/builtins.gie @@ -1629,11 +1629,12 @@ expect -0.002029979 -0.000789630 =============================================================================== # Equidistant Cylindrical (Plate Carree) -# Cyl, Sph +# Cyl, Sph&Ell # lat_ts=[, lat_0=0] =============================================================================== ------------------------------------------------------------------------------- +# Spherical case operation +proj=eqc +a=6400000 ------------------------------------------------------------------------------- tolerance 0.1 mm @@ -1656,6 +1657,143 @@ expect -0.001790493 0.000895247 accept -200 -100 expect -0.001790493 -0.000895247 +------------------------------------------------------------------------------- +# Ellipsoidal case (EPSG:1028) +# Test values from IOGP Guidance Note 7-2, Section 3.2.5 +# WGS84 ellipsoid: a=6378137, 1/f=298.257223563 +# Standard parallel: 0° (lat_ts=0) +# Input: lat=55°, lon=10° +# Expected: E=1113194.91, N=6097230.31 +operation +proj=eqc +ellps=WGS84 +lat_ts=0 +------------------------------------------------------------------------------- +tolerance 0.01 m +accept 10 55 +expect 1113194.91 6097230.31 + +direction inverse +accept 1113194.91 6097230.31 +expect 10 55 + +------------------------------------------------------------------------------- +# Ellipsoidal case - high latitude test +# WGS84 ellipsoid +# Input: lat=71.9993230521°, lon=-143.9999611505° +# Expected: +# E=-16030002.350, N=7992053.817 +operation +proj=eqc +ellps=WGS84 +lat_ts=0 +------------------------------------------------------------------------------- +tolerance 0.1 m +accept -143.9999611505 71.9993230521 +expect -16030002.350 7992053.817 + +direction inverse +accept -16030002.350 7992053.817 +expect -143.9999611505 71.9993230521 + +------------------------------------------------------------------------------- +# Ellipsoidal case - Edge cases +# WGS84 ellipsoid, lat_ts=0 +# Tests for origin, hemispheres, date line, near-pole +operation +proj=eqc +ellps=WGS84 +lat_ts=0 +------------------------------------------------------------------------------- +tolerance 0.001 m + +# Origin (0, 0) +accept 0 0 +expect 0.0 0.0 + +# Southern hemisphere (10, -45) +accept 10 -45 +expect 1113194.90793 -4984944.37798 + +# Date line (180, 30) +accept 180 30 +expect 20037508.34279 3320113.39794 + +# Near north pole (0, 89) +accept 0 89 +expect 0.0 9890271.86440 + +# San Francisco (-122.4194, 37.7749) +accept -122.4194 37.7749 +expect -13627665.27122 4182513.19136 + +direction inverse +# Inverse tests for edge cases +accept 0.0 0.0 +expect 0 0 + +accept 1113194.90793 -4984944.37798 +expect 10 -45 + +accept 20037508.34279 3320113.39794 +expect 180 30 + +accept 0.0 9890271.86440 +expect 0 89 + +accept -13627665.27122 4182513.19136 +expect -122.4194 37.7749 + +------------------------------------------------------------------------------- +# Ellipsoidal case with non-zero standard parallel (lat_ts=45) +# WGS84 ellipsoid +# Tests the ν₁ cos(φ₁) scaling factor for easting +operation +proj=eqc +ellps=WGS84 +lat_ts=45 +------------------------------------------------------------------------------- +tolerance 0.01 m + +# Paris region (2, 49) +accept 2 49 +expect 157693.670 5429627.632 + +# Origin (0, 0) - northing should still be 0 +accept 0 0 +expect 0.0 0.0 + +# High latitude (10, 70) +accept 10 70 +expect 788468.351 7768980.728 + +direction inverse +accept 157693.670 5429627.632 +expect 2 49 + +accept 0.0 0.0 +expect 0 0 + +accept 788468.351 7768980.728 +expect 10 70 + +------------------------------------------------------------------------------- +# Ellipsoidal case with lat_ts=30 and lat_0=45 (non-zero origin) +# Tests meridional arc offset M₀ +operation +proj=eqc +ellps=WGS84 +lat_ts=30 +lat_0=45 +------------------------------------------------------------------------------- +tolerance 0.01 m + +# At the origin latitude, northing should be 0 +accept 0 45 +expect 0.0 0.0 + +# Above origin +accept 0 60 +expect 0.0 1669128.442 + +# Below origin +accept 0 30 +expect 0.0 -1664830.980 + +direction inverse +accept 0.0 0.0 +expect 0 45 + +accept 0.0 1669128.442 +expect 0 60 + +accept 0.0 -1664830.980 +expect 0 30 + =============================================================================== # Equidistant Conic