Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Implement EPSG:1026 Mercator (Spherical) method #3741

Merged
merged 2 commits into from
May 31, 2023
Merged

Implement EPSG:1026 Mercator (Spherical) method #3741

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
d069819
Add +R_C as an ellipsoid spherification parameter to use the radius o…
rouault May 25, 2023
2ba6637
Map EPSG:1026 Mercator (Spherical) method to +proj=merc +R_C
rouault May 25, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
16 changes: 11 additions & 5 deletions docs/source/usage/ellipsoids.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -65,23 +65,29 @@ the ellipsoid into a sphere with features defined by the ellipsoid.
Ellipsoid spherification parameters
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

.. option:: +R_A=<value>
.. option:: +R_A

A sphere with the same surface area as the ellipsoid.

.. option:: +R_V=<value>
.. option:: +R_V

A sphere with the same volume as the ellipsoid.

.. option:: +R_a=<value>
.. option:: +R_C

.. versionadded:: 9.3.0

A sphere whose radius is the radius of the conformal sphere at :math:`\phi_0`.

.. option:: +R_a

A sphere with :math:`R = (a + b)/2` (arithmetic mean).

.. option:: +R_g=<value>
.. option:: +R_g

A sphere with :math:`R = \sqrt{ab}` (geometric mean).

.. option:: +R_h=<value>
.. option:: +R_h

A sphere with :math:`R = 2ab/(a+b)` (harmonic mean).

Expand Down
5 changes: 5 additions & 0 deletions include/proj/coordinateoperation.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1292,6 +1292,11 @@ class PROJ_GCC_DLL Conversion : public SingleOperation {
const common::Angle &centerLong, const common::Length &falseEasting,
const common::Length &falseNorthing);

PROJ_DLL static ConversionNNPtr createMercatorSpherical(
const util::PropertyMap &properties, const common::Angle &centerLat,
const common::Angle &centerLong, const common::Length &falseEasting,
const common::Length &falseNorthing);

PROJ_DLL static ConversionNNPtr
createMollweide(const util::PropertyMap &properties,
const common::Angle &centerLong,
Expand Down
1 change: 1 addition & 0 deletions scripts/reference_exported_symbols.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -589,6 +589,7 @@ osgeo::proj::operation::Conversion::createLambertConicConformal_2SP_Michigan(osg
osgeo::proj::operation::Conversion::createLambertConicConformal_2SP(osgeo::proj::util::PropertyMap const&, osgeo::proj::common::Angle const&, osgeo::proj::common::Angle const&, osgeo::proj::common::Angle const&, osgeo::proj::common::Angle const&, osgeo::proj::common::Length const&, osgeo::proj::common::Length const&)
osgeo::proj::operation::Conversion::createLambertCylindricalEqualArea(osgeo::proj::util::PropertyMap const&, osgeo::proj::common::Angle const&, osgeo::proj::common::Angle const&, osgeo::proj::common::Length const&, osgeo::proj::common::Length const&)
osgeo::proj::operation::Conversion::createLambertCylindricalEqualAreaSpherical(osgeo::proj::util::PropertyMap const&, osgeo::proj::common::Angle const&, osgeo::proj::common::Angle const&, osgeo::proj::common::Length const&, osgeo::proj::common::Length const&)
osgeo::proj::operation::Conversion::createMercatorSpherical(osgeo::proj::util::PropertyMap const&, osgeo::proj::common::Angle const&, osgeo::proj::common::Angle const&, osgeo::proj::common::Length const&, osgeo::proj::common::Length const&)
osgeo::proj::operation::Conversion::createMercatorVariantA(osgeo::proj::util::PropertyMap const&, osgeo::proj::common::Angle const&, osgeo::proj::common::Angle const&, osgeo::proj::common::Scale const&, osgeo::proj::common::Length const&, osgeo::proj::common::Length const&)
osgeo::proj::operation::Conversion::createMercatorVariantB(osgeo::proj::util::PropertyMap const&, osgeo::proj::common::Angle const&, osgeo::proj::common::Angle const&, osgeo::proj::common::Length const&, osgeo::proj::common::Length const&)
osgeo::proj::operation::Conversion::createMillerCylindrical(osgeo::proj::util::PropertyMap const&, osgeo::proj::common::Angle const&, osgeo::proj::common::Length const&, osgeo::proj::common::Length const&)
Expand Down
33 changes: 27 additions & 6 deletions src/ell_set.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -49,16 +49,21 @@ int pj_ellipsoid(PJ *P) {

Spherification parameters supported are:
R_A, which gives a sphere with the same surface area as the
ellipsoid R_V, which gives a sphere with the same volume as the ellipsoid
ellipsoid R_V, which gives a sphere with the same volume as the
ellipsoid

R_a, which gives a sphere with R = (a + b)/2 (arithmetic mean)
R_g, which gives a sphere with R = sqrt(a*b) (geometric mean)
R_h, which gives a sphere with R = 2*a*b/(a+b) (harmonic mean)

R_lat_a=phi, which gives a sphere with R being the arithmetic mean
of of the corresponding ellipsoid at latitude phi. R_lat_g=phi, which gives
a sphere with R being the geometric mean of of the corresponding ellipsoid
at latitude phi.
of of the corresponding ellipsoid at latitude phi.
R_lat_g=phi, which gives
a sphere with R being the geometric mean of of the corresponding
ellipsoid at latitude phi.

R_C, which gives a sphere with the radius of the conformal sphere
at phi0.

If R is given as size parameter, any shape and spherification parameters
given are ignored.
Expand Down Expand Up @@ -349,8 +354,8 @@ static const double RV6 = 55 / 1296.;
/***************************************************************************************/
static int ellps_spherification(PJ *P) {
/***************************************************************************************/
const char *keys[] = {"R_A", "R_V", "R_a", "R_g",
"R_h", "R_lat_a", "R_lat_g"};
const char *keys[] = {"R_A", "R_V", "R_a", "R_g",
"R_h", "R_lat_a", "R_lat_g", "R_C"};
size_t len, i;
paralist *par = nullptr;

Expand Down Expand Up @@ -428,6 +433,22 @@ static int ellps_spherification(PJ *P) {
else /* geometric */
P->a *= sqrt(1 - P->es) / t;
break;

/* R_C - a sphere with R = radius of conformal sphere, taken at a
* latitude that is phi0 (note: at least for mercator. for other
* projection methods, this could be phi1)
* Formula from IOGP Publication 373-7-2 – Geomatics Guidance Note number 7,
* part 2 1.1 Ellipsoid parameters
*/
case 7:
t = sin(P->phi0);
t = 1 - P->es * t * t;
if (t == 0.) {
proj_log_error(P, _("Invalid eccentricity"));
return proj_errno_set(P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
}
P->a *= sqrt(1 - P->es) / t;
break;
}

if (P->a <= 0.) {
Expand Down
19 changes: 19 additions & 0 deletions src/iso19111/io.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -11180,7 +11180,26 @@ PROJStringParser::Private::buildProjectedCRS(int iStep,
}
}
} else if (hasParamValue(step, "lat_ts")) {
if (hasParamValue(step, "R_C") &&
!geodCRS->ellipsoid()->isSphere() &&
getAngularValue(getParamValue(step, "lat_ts")) != 0) {
throw ParsingException("lat_ts != 0 not supported for "
"spherical Mercator on an ellipsoid");
}
mapping = getMapping(EPSG_CODE_METHOD_MERCATOR_VARIANT_B);
} else if (hasParamValue(step, "R_C")) {
const auto &k = getParamValueK(step);
if (!k.empty() && getNumericValue(k) != 1.0) {
if (geodCRS->ellipsoid()->isSphere()) {
mapping = getMapping(EPSG_CODE_METHOD_MERCATOR_VARIANT_A);
} else {
throw ParsingException(
"k_0 != 1 not supported for spherical Mercator on an "
"ellipsoid");
}
} else {
mapping = getMapping(EPSG_CODE_METHOD_MERCATOR_SPHERICAL);
}
} else {
mapping = getMapping(EPSG_CODE_METHOD_MERCATOR_VARIANT_A);
}
Expand Down
32 changes: 32 additions & 0 deletions src/iso19111/operation/conversion.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1751,6 +1751,38 @@ ConversionNNPtr Conversion::createPopularVisualisationPseudoMercator(

// ---------------------------------------------------------------------------

/** \brief Instantiate a conversion based on the
* <a href="../../../operations/projections/merc.html">
* Mercator</a> projection method, using its spherical formulation
*
* When used with an ellipsoid, the radius used is the radius of the conformal
* sphere at centerLat.
*
* This method is defined as
* <a
* href="https://epsg.org/coord-operation-method_1026/Mercator-Spherical.html">
* EPSG:1026</a>.
*
* @param properties See \ref general_properties of the conversion. If the name
* is not provided, it is automatically set.
* @param centerLat See \ref center_latitude . Usually 0
* @param centerLong See \ref center_longitude . Usually 0
* @param falseEasting See \ref false_easting . Usually 0
* @param falseNorthing See \ref false_northing . Usually 0
* @return a new Conversion.
* @since 9.3
*/
ConversionNNPtr Conversion::createMercatorSpherical(
const util::PropertyMap &properties, const common::Angle &centerLat,
const common::Angle &centerLong, const common::Length &falseEasting,
const common::Length &falseNorthing) {
return create(
properties, EPSG_CODE_METHOD_MERCATOR_SPHERICAL,
createParams(centerLat, centerLong, falseEasting, falseNorthing));
}

// ---------------------------------------------------------------------------

/** \brief Instantiate a conversion based on the
* <a href="../../../operations/projections/moll.html">
* Mollweide</a> projection method.
Expand Down
4 changes: 4 additions & 0 deletions src/iso19111/operation/parammappings.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -770,6 +770,9 @@ static const MethodMapping projectionMethodMappings[] = {
// handled manually
"webmerc", nullptr, paramsNatOrigin},

{EPSG_NAME_METHOD_MERCATOR_SPHERICAL, EPSG_CODE_METHOD_MERCATOR_SPHERICAL,
nullptr, "merc", "R_C", paramsNatOrigin},

{PROJ_WKT2_NAME_METHOD_MOLLWEIDE, 0, "Mollweide", "moll", nullptr,
paramsLongNatOrigin},

Expand Down Expand Up @@ -938,6 +941,7 @@ const struct MethodNameCode methodNameCodes[] = {
METHOD_NAME_CODE(KROVAK),
METHOD_NAME_CODE(LAMBERT_AZIMUTHAL_EQUAL_AREA),
METHOD_NAME_CODE(POPULAR_VISUALISATION_PSEUDO_MERCATOR),
METHOD_NAME_CODE(MERCATOR_SPHERICAL),
METHOD_NAME_CODE(MERCATOR_VARIANT_A),
METHOD_NAME_CODE(MERCATOR_VARIANT_B),
METHOD_NAME_CODE(OBLIQUE_STEREOGRAPHIC),
Expand Down
3 changes: 3 additions & 0 deletions src/proj_constants.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -187,6 +187,9 @@
"Popular Visualisation Pseudo Mercator"
#define EPSG_CODE_METHOD_POPULAR_VISUALISATION_PSEUDO_MERCATOR 1024

#define EPSG_NAME_METHOD_MERCATOR_SPHERICAL "Mercator (Spherical)"
#define EPSG_CODE_METHOD_MERCATOR_SPHERICAL 1026

#define PROJ_WKT2_NAME_METHOD_MOLLWEIDE "Mollweide"

#define PROJ_WKT2_NAME_METHOD_NATURAL_EARTH "Natural Earth"
Expand Down
16 changes: 16 additions & 0 deletions test/gie/builtins.gie
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3749,6 +3749,22 @@ direction inverse
accept 0 1e-15
expect 0 57.295779513e-15

-------------------------------------------------------------------------------
# Mercator Spherical
operation +proj=merc +R_C +ellps=WGS84 +lat_0=45
-------------------------------------------------------------------------------
accept 2 49
expect 221892.515234695253 6253822.971804938279

-------------------------------------------------------------------------------
# Mercator Spherical
operation +proj=merc +R_C +a=6371007 +b=6371007
-------------------------------------------------------------------------------
# Test point from IOGP guidance note 7.2
tolerance 1cm
accept -100.33333333333333 24.381786944444446
expect -11156569.90 2796869.94


===============================================================================
# Miller Oblated Stereographic
Expand Down
17 changes: 17 additions & 0 deletions test/unit/test_io.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10618,6 +10618,23 @@ TEST(io, projparse_cea_ellipsoidal_with_k_0) {

// ---------------------------------------------------------------------------

TEST(io, projparse_merc_spherical_on_ellipsoid) {
std::string input("+proj=merc +R_C +lat_0=1 +lon_0=2 +x_0=3 +y_0=4 "
"+ellps=WGS84 +units=m +no_defs +type=crs");
auto obj = PROJStringParser().createFromPROJString(input);
auto crs = nn_dynamic_pointer_cast<ProjectedCRS>(obj);
ASSERT_TRUE(crs != nullptr);
EXPECT_EQ(crs->derivingConversion()->method()->getEPSGCode(),
EPSG_CODE_METHOD_MERCATOR_SPHERICAL);
EXPECT_EQ(
crs->exportToPROJString(
PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4)
.get()),
input);
}

// ---------------------------------------------------------------------------

TEST(io, projparse_geos_sweep_x) {
auto obj = PROJStringParser().createFromPROJString(
"+proj=geos +sweep=x +h=1 +type=crs");
Expand Down
79 changes: 79 additions & 0 deletions test/unit/test_operation.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3376,6 +3376,85 @@ TEST(operation, webmerc_import_from_broken_esri_WGS_84_Pseudo_Mercator) {

// ---------------------------------------------------------------------------

TEST(operation, mercator_spherical_export) {
auto conv = Conversion::createMercatorSpherical(
PropertyMap(), Angle(1), Angle(2), Length(3), Length(4));
EXPECT_TRUE(conv->validateParameters().empty());

EXPECT_EQ(conv->exportToPROJString(PROJStringFormatter::create().get()),
"+proj=merc +R_C +lat_0=1 +lon_0=2 +x_0=3 +y_0=4");

EXPECT_EQ(conv->exportToWKT(WKTFormatter::create().get()),
"CONVERSION[\"Mercator (Spherical)\",\n"
" METHOD[\"Mercator (Spherical)\",\n"
" ID[\"EPSG\",1026]],\n"
" PARAMETER[\"Latitude of natural origin\",1,\n"
" ANGLEUNIT[\"degree\",0.0174532925199433],\n"
" ID[\"EPSG\",8801]],\n"
" PARAMETER[\"Longitude of natural origin\",2,\n"
" ANGLEUNIT[\"degree\",0.0174532925199433],\n"
" ID[\"EPSG\",8802]],\n"
" PARAMETER[\"False easting\",3,\n"
" LENGTHUNIT[\"metre\",1],\n"
" ID[\"EPSG\",8806]],\n"
" PARAMETER[\"False northing\",4,\n"
" LENGTHUNIT[\"metre\",1],\n"
" ID[\"EPSG\",8807]]]");
}

// ---------------------------------------------------------------------------

TEST(operation, mercator_spherical_import) {
auto wkt2 = "PROJCRS[\"unknown\",\n"
" BASEGEOGCRS[\"unknown\",\n"
" DATUM[\"World Geodetic System 1984\",\n"
" ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n"
" LENGTHUNIT[\"metre\",1]],\n"
" ID[\"EPSG\",6326]],\n"
" PRIMEM[\"Greenwich\",0,\n"
" ANGLEUNIT[\"degree\",0.0174532925199433],\n"
" ID[\"EPSG\",8901]]],\n"
" CONVERSION[\"unknown\",\n"
" METHOD[\"Mercator (Spherical)\",\n"
" ID[\"EPSG\",1026]],\n"
" PARAMETER[\"Latitude of natural origin\",0,\n"
" ANGLEUNIT[\"degree\",0.0174532925199433],\n"
" ID[\"EPSG\",8801]],\n"
" PARAMETER[\"Longitude of natural origin\",0,\n"
" ANGLEUNIT[\"degree\",0.0174532925199433],\n"
" ID[\"EPSG\",8802]],\n"
" PARAMETER[\"False easting\",0,\n"
" LENGTHUNIT[\"metre\",1],\n"
" ID[\"EPSG\",8806]],\n"
" PARAMETER[\"False northing\",0,\n"
" LENGTHUNIT[\"metre\",1],\n"
" ID[\"EPSG\",8807]]],\n"
" CS[Cartesian,2],\n"
" AXIS[\"(E)\",east,\n"
" ORDER[1],\n"
" LENGTHUNIT[\"metre\",1,\n"
" ID[\"EPSG\",9001]]],\n"
" AXIS[\"(N)\",north,\n"
" ORDER[2],\n"
" LENGTHUNIT[\"metre\",1,\n"
" ID[\"EPSG\",9001]]]]";

auto obj = WKTParser().createFromWKT(wkt2);
auto crs = nn_dynamic_pointer_cast<ProjectedCRS>(obj);
ASSERT_TRUE(crs != nullptr);

EXPECT_EQ(crs->exportToPROJString(PROJStringFormatter::create().get()),
"+proj=merc +R_C +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 "
"+units=m +no_defs +type=crs");

EXPECT_EQ(
crs->exportToWKT(
WKTFormatter::create(WKTFormatter::Convention::WKT2_2019).get()),
wkt2);
}

// ---------------------------------------------------------------------------

TEST(operation, mollweide_export) {

auto conv = Conversion::createMollweide(PropertyMap(), Angle(1), Length(2),
Expand Down