aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--docs/source/operations/projections/ortho.rst10
-rw-r--r--include/proj/internal/coordinateoperation_constants.hpp3
-rw-r--r--include/proj/internal/esri_projection_mappings.hpp19
-rw-r--r--scripts/build_esri_projection_mapping.py12
-rw-r--r--src/iso19111/coordinateoperation.cpp2
-rw-r--r--src/iso19111/io.cpp37
-rw-r--r--src/proj_constants.h3
-rw-r--r--src/projections/ortho.cpp87
-rw-r--r--test/gie/builtins.gie68
-rw-r--r--test/unit/test_io.cpp68
10 files changed, 295 insertions, 14 deletions
diff --git a/docs/source/operations/projections/ortho.rst b/docs/source/operations/projections/ortho.rst
index 2e8de9ca..2d0e08a8 100644
--- a/docs/source/operations/projections/ortho.rst
+++ b/docs/source/operations/projections/ortho.rst
@@ -10,7 +10,7 @@ around a given latitude and longitude.
+---------------------+--------------------------------------------------------+
| **Classification** | Azimuthal |
+---------------------+--------------------------------------------------------+
-| **Available forms** | Forward and inverse, spherical projection |
+| **Available forms** | Forward and inverse, spherical and ellipsoidal |
+---------------------+--------------------------------------------------------+
| **Defined area** | Global, although only one hemisphere can be seen at a |
| | time |
@@ -31,6 +31,12 @@ around a given latitude and longitude.
proj-string: ``+proj=ortho``
+
+.. note:: Before PROJ 7.2, only the spherical formulation was implemented. If
+ wanting to replicate PROJ < 7.2 results with newer versions, the
+ ellipsoid must be forced to a sphere, for example by adding a ``+f=0``
+ parameter.
+
Parameters
################################################################################
@@ -40,6 +46,8 @@ Parameters
.. include:: ../options/lat_0.rst
+.. include:: ../options/ellps.rst
+
.. include:: ../options/R.rst
.. include:: ../options/x_0.rst
diff --git a/include/proj/internal/coordinateoperation_constants.hpp b/include/proj/internal/coordinateoperation_constants.hpp
index 79b47a6a..7038ba9f 100644
--- a/include/proj/internal/coordinateoperation_constants.hpp
+++ b/include/proj/internal/coordinateoperation_constants.hpp
@@ -729,6 +729,9 @@ static const MethodMapping projectionMethodMappings[] = {
{EPSG_NAME_METHOD_ORTHOGRAPHIC, EPSG_CODE_METHOD_ORTHOGRAPHIC,
"Orthographic", "ortho", nullptr, paramsNatOrigin},
+ {PROJ_WKT2_NAME_ORTHOGRAPHIC_SPHERICAL, 0, "Orthographic", "ortho", "f=0",
+ paramsNatOrigin},
+
{PROJ_WKT2_NAME_METHOD_PATTERSON, 0, "Patterson", "patterson", nullptr,
paramsLonNatOrigin},
diff --git a/include/proj/internal/esri_projection_mappings.hpp b/include/proj/internal/esri_projection_mappings.hpp
index 752c03ad..c230f26b 100644
--- a/include/proj/internal/esri_projection_mappings.hpp
+++ b/include/proj/internal/esri_projection_mappings.hpp
@@ -600,6 +600,19 @@ static const ESRIParamMapping paramsESRI_Orthographic[] = {
EPSG_CODE_PARAMETER_LATITUDE_OF_NATURAL_ORIGIN, "0.0", false},
{nullptr, nullptr, 0, "0.0", false}};
+static const ESRIParamMapping paramsESRI_Local[] = {
+ {"False_Easting", EPSG_NAME_PARAMETER_FALSE_EASTING,
+ EPSG_CODE_PARAMETER_FALSE_EASTING, "0.0", false},
+ {"False_Northing", EPSG_NAME_PARAMETER_FALSE_NORTHING,
+ EPSG_CODE_PARAMETER_FALSE_NORTHING, "0.0", false},
+ {"Scale_Factor", nullptr, 0, "1.0", false},
+ {"Azimuth", nullptr, 0, "0.0", false},
+ {"Longitude_Of_Center", EPSG_NAME_PARAMETER_LONGITUDE_OF_NATURAL_ORIGIN,
+ EPSG_CODE_PARAMETER_LONGITUDE_OF_NATURAL_ORIGIN, "0.0", false},
+ {"Latitude_Of_Center", EPSG_NAME_PARAMETER_LATITUDE_OF_NATURAL_ORIGIN,
+ EPSG_CODE_PARAMETER_LATITUDE_OF_NATURAL_ORIGIN, "0.0", false},
+ {nullptr, nullptr, 0, "0.0", false}};
+
static const ESRIParamMapping paramsESRI_Winkel_Tripel[] = {
{"False_Easting", EPSG_NAME_PARAMETER_FALSE_EASTING,
EPSG_CODE_PARAMETER_FALSE_EASTING, "0.0", false},
@@ -1004,8 +1017,10 @@ static const ESRIMethodMapping esriMappings[] = {
EPSG_CODE_METHOD_KROVAK_NORTH_ORIENTED, paramsESRI_Krovak_alt2},
{"New_Zealand_Map_Grid", EPSG_NAME_METHOD_NZMG, EPSG_CODE_METHOD_NZMG,
paramsESRI_New_Zealand_Map_Grid},
- {"Orthographic", EPSG_NAME_METHOD_ORTHOGRAPHIC,
- EPSG_CODE_METHOD_ORTHOGRAPHIC, paramsESRI_Orthographic},
+ {"Orthographic", PROJ_WKT2_NAME_ORTHOGRAPHIC_SPHERICAL, 0,
+ paramsESRI_Orthographic},
+ {"Local", EPSG_NAME_METHOD_ORTHOGRAPHIC, EPSG_CODE_METHOD_ORTHOGRAPHIC,
+ paramsESRI_Local},
{"Winkel_Tripel", "Winkel Tripel", 0, paramsESRI_Winkel_Tripel},
{"Aitoff", "Aitoff", 0, paramsESRI_Aitoff},
{"Flat_Polar_Quartic", PROJ_WKT2_NAME_METHOD_FLAT_POLAR_QUARTIC, 0,
diff --git a/scripts/build_esri_projection_mapping.py b/scripts/build_esri_projection_mapping.py
index 5206e191..12eeb852 100644
--- a/scripts/build_esri_projection_mapping.py
+++ b/scripts/build_esri_projection_mapping.py
@@ -436,11 +436,23 @@ config_str = """
- Longitude_Of_Origin: EPSG_NAME_PARAMETER_LONGITUDE_OF_NATURAL_ORIGIN
- Latitude_Of_Origin: EPSG_NAME_PARAMETER_LATITUDE_OF_NATURAL_ORIGIN
+# ESRI's Orthographic is a spherical-only formulation. The ellipsoidal capable
+# name is Local. See below
- Orthographic:
+ WKT2_name: PROJ_WKT2_NAME_ORTHOGRAPHIC_SPHERICAL
+ Params:
+ - False_Easting: EPSG_NAME_PARAMETER_FALSE_EASTING
+ - False_Northing: EPSG_NAME_PARAMETER_FALSE_NORTHING
+ - Longitude_Of_Center: EPSG_NAME_PARAMETER_LONGITUDE_OF_NATURAL_ORIGIN
+ - Latitude_Of_Center: EPSG_NAME_PARAMETER_LATITUDE_OF_NATURAL_ORIGIN
+
+- Local:
WKT2_name: EPSG_NAME_METHOD_ORTHOGRAPHIC
Params:
- False_Easting: EPSG_NAME_PARAMETER_FALSE_EASTING
- False_Northing: EPSG_NAME_PARAMETER_FALSE_NORTHING
+ - Scale_Factor: 1.0
+ - Azimuth: 0.0
- Longitude_Of_Center: EPSG_NAME_PARAMETER_LONGITUDE_OF_NATURAL_ORIGIN
- Latitude_Of_Center: EPSG_NAME_PARAMETER_LATITUDE_OF_NATURAL_ORIGIN
diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp
index 65ef2b77..ef262cc8 100644
--- a/src/iso19111/coordinateoperation.cpp
+++ b/src/iso19111/coordinateoperation.cpp
@@ -4251,7 +4251,7 @@ ConversionNNPtr Conversion::createObliqueStereographic(
* This method is defined as [EPSG:9840]
* (https://www.epsg-registry.org/export.htm?gml=urn:ogc:def:method:EPSG::9840)
*
- * \note At the time of writing, PROJ only implements the spherical formulation
+ * \note Before PROJ 7.2, only the spherical formulation was implemented.
*
* @param properties See \ref general_properties of the conversion. If the name
* is not provided, it is automatically set.
diff --git a/src/iso19111/io.cpp b/src/iso19111/io.cpp
index 6680ab00..16ab22f7 100644
--- a/src/iso19111/io.cpp
+++ b/src/iso19111/io.cpp
@@ -3290,7 +3290,7 @@ ConversionNNPtr WKTParser::Private::buildProjectionFromESRI(
}
// Compare parameters present with the ones expected in the mapping
- const ESRIMethodMapping *esriMapping = esriMappings[0];
+ const ESRIMethodMapping *esriMapping = nullptr;
int bestMatchCount = -1;
for (const auto &mapping : esriMappings) {
int matchCount = 0;
@@ -3298,12 +3298,20 @@ ConversionNNPtr WKTParser::Private::buildProjectionFromESRI(
auto iter = mapParamNameToValue.find(param->esri_name);
if (iter != mapParamNameToValue.end()) {
if (param->wkt2_name == nullptr) {
+ bool ok = true;
try {
if (io::asDouble(param->fixed_value) ==
io::asDouble(iter->second)) {
matchCount++;
+ } else {
+ ok = false;
}
} catch (const std::exception &) {
+ ok = false;
+ }
+ if (!ok) {
+ matchCount = -1;
+ break;
}
} else {
matchCount++;
@@ -3317,6 +3325,10 @@ ConversionNNPtr WKTParser::Private::buildProjectionFromESRI(
bestMatchCount = matchCount;
}
}
+ if (esriMapping == nullptr) {
+ return buildProjectionStandard(baseGeodCRS, projCRSNode, projectionNode,
+ defaultLinearUnit, defaultAngularUnit);
+ }
std::map<std::string, const char *> mapWKT2NameToESRIName;
for (const auto *param = esriMapping->params; param->esri_name; ++param) {
@@ -8898,8 +8910,29 @@ static bool is_in_stringlist(const std::string &str, const char *stringlist) {
CRSNNPtr PROJStringParser::Private::buildProjectedCRS(
int iStep, GeographicCRSNNPtr geogCRS, int iUnitConvert, int iAxisSwap) {
auto &step = steps_[iStep];
- auto mappings = getMappingsFromPROJName(step.name);
+ const auto mappings = getMappingsFromPROJName(step.name);
const MethodMapping *mapping = mappings.empty() ? nullptr : mappings[0];
+
+ if (mappings.size() >= 2) {
+ // To distinguish for example +ortho from +ortho +f=0
+ for (const auto *mappingIter : mappings) {
+ if (mappingIter->proj_name_aux != nullptr &&
+ strchr(mappingIter->proj_name_aux, '=') == nullptr &&
+ hasParamValue(step, mappingIter->proj_name_aux)) {
+ mapping = mappingIter;
+ break;
+ } else if (mappingIter->proj_name_aux != nullptr &&
+ strchr(mappingIter->proj_name_aux, '=') != nullptr) {
+ const auto tokens = split(mappingIter->proj_name_aux, '=');
+ if (tokens.size() == 2 &&
+ getParamValue(step, tokens[0]) == tokens[1]) {
+ mapping = mappingIter;
+ break;
+ }
+ }
+ }
+ }
+
if (mapping) {
mapping = selectSphericalOrEllipsoidal(mapping, geogCRS);
}
diff --git a/src/proj_constants.h b/src/proj_constants.h
index 5b626223..f64bf496 100644
--- a/src/proj_constants.h
+++ b/src/proj_constants.h
@@ -183,6 +183,9 @@
#define EPSG_NAME_METHOD_ORTHOGRAPHIC "Orthographic"
#define EPSG_CODE_METHOD_ORTHOGRAPHIC 9840
+#define PROJ_WKT2_NAME_ORTHOGRAPHIC_SPHERICAL \
+ "Orthographic (Spherical)"
+
#define PROJ_WKT2_NAME_METHOD_PATTERSON "Patterson"
#define EPSG_NAME_METHOD_AMERICAN_POLYCONIC "American Polyconic"
diff --git a/src/projections/ortho.cpp b/src/projections/ortho.cpp
index 5f9366c3..ff9d9341 100644
--- a/src/projections/ortho.cpp
+++ b/src/projections/ortho.cpp
@@ -4,7 +4,7 @@
#include "proj_internal.h"
#include <math.h>
-PROJ_HEAD(ortho, "Orthographic") "\n\tAzi, Sph";
+PROJ_HEAD(ortho, "Orthographic") "\n\tAzi, Sph&Ell";
namespace { // anonymous namespace
enum Mode {
@@ -19,6 +19,7 @@ namespace { // anonymous namespace
struct pj_opaque {
double sinph0;
double cosph0;
+ double nu0;
enum Mode mode;
};
} // anonymous namespace
@@ -121,6 +122,72 @@ static PJ_LP ortho_s_inverse (PJ_XY xy, PJ *P) { /* Spheroidal, invers
}
+static PJ_XY ortho_e_forward (PJ_LP lp, PJ *P) { /* Ellipsoidal, forward */
+ PJ_XY xy;
+ struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque);
+
+ // From EPSG guidance note 7.2
+ const double cosphi = cos(lp.phi);
+ const double sinphi = sin(lp.phi);
+ const double coslam = cos(lp.lam);
+ const double sinlam = sin(lp.lam);
+ const double nu = 1.0 / sqrt(1.0 - P->es * sinphi * sinphi);
+ xy.x = nu * cosphi * sinlam;
+ xy.y = nu * (sinphi * Q->cosph0 - cosphi * Q->sinph0 * coslam) +
+ P->es * (Q->nu0 * Q->sinph0 - nu * sinphi) * Q->cosph0;
+
+ return xy;
+}
+
+
+
+static PJ_LP ortho_e_inverse (PJ_XY xy, PJ *P) { /* Ellipsoidal, inverse */
+ PJ_LP lp;
+ struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque);
+
+ // From EPSG guidance note 7.2
+
+ // It suggests as initial guess:
+ // lp.lam = 0;
+ // lp.phi = P->phi0;
+ // But for poles, this will not converge well. Better use:
+ lp = ortho_s_inverse(xy, P);
+
+ for( int i = 0; i < 20; i++ )
+ {
+ const double cosphi = cos(lp.phi);
+ const double sinphi = sin(lp.phi);
+ const double coslam = cos(lp.lam);
+ const double sinlam = sin(lp.lam);
+ const double one_minus_es_sinphi2 = 1.0 - P->es * sinphi * sinphi;
+ const double nu = 1.0 / sqrt(one_minus_es_sinphi2);
+ PJ_XY xy_new;
+ xy_new.x = nu * cosphi * sinlam;
+ xy_new.y = nu * (sinphi * Q->cosph0 - cosphi * Q->sinph0 * coslam) +
+ P->es * (Q->nu0 * Q->sinph0 - nu * sinphi) * Q->cosph0;
+ const double rho = (1.0 - P->es) * nu / one_minus_es_sinphi2;
+ const double J11 = -rho * sinphi * sinlam;
+ const double J12 = nu * cosphi * coslam;
+ const double J21 = rho * (cosphi * Q->cosph0 + sinphi * Q->sinph0 * coslam);
+ const double J22 = nu * Q->sinph0 * Q->cosph0 * sinlam;
+ const double D = J11 * J22 - J12 * J21;
+ const double dx = xy.x - xy_new.x;
+ const double dy = xy.y - xy_new.y;
+ const double dphi = (J22 * dx - J12 * dy) / D;
+ const double dlam = (-J21 * dx + J11 * dy) / D;
+ lp.phi += dphi;
+ if( lp.phi > M_PI_2) lp.phi = M_PI_2;
+ else if( lp.phi < -M_PI_2) lp.phi = -M_PI_2;
+ lp.lam += dlam;
+ if( fabs(dphi) < 1e-12 && fabs(dlam) < 1e-12 )
+ {
+ return lp;
+ }
+ }
+ pj_ctx_set_errno(P->ctx, PJD_ERR_NON_CONVERGENT);
+ return lp;
+}
+
PJ *PROJECTION(ortho) {
struct pj_opaque *Q = static_cast<struct pj_opaque*>(pj_calloc (1, sizeof (struct pj_opaque)));
@@ -128,17 +195,25 @@ PJ *PROJECTION(ortho) {
return pj_default_destructor(P, ENOMEM);
P->opaque = Q;
+ Q->sinph0 = sin(P->phi0);
+ Q->cosph0 = cos(P->phi0);
if (fabs(fabs(P->phi0) - M_HALFPI) <= EPS10)
Q->mode = P->phi0 < 0. ? S_POLE : N_POLE;
else if (fabs(P->phi0) > EPS10) {
Q->mode = OBLIQ;
- Q->sinph0 = sin(P->phi0);
- Q->cosph0 = cos(P->phi0);
} else
Q->mode = EQUIT;
- P->inv = ortho_s_inverse;
- P->fwd = ortho_s_forward;
- P->es = 0.;
+ if( P->es == 0 )
+ {
+ P->inv = ortho_s_inverse;
+ P->fwd = ortho_s_forward;
+ }
+ else
+ {
+ Q->nu0 = 1.0 / sqrt(1.0 - P->es * Q->sinph0 * Q->sinph0);
+ P->inv = ortho_e_inverse;
+ P->fwd = ortho_e_forward;
+ }
return P;
}
diff --git a/test/gie/builtins.gie b/test/gie/builtins.gie
index f755526e..48afc00a 100644
--- a/test/gie/builtins.gie
+++ b/test/gie/builtins.gie
@@ -4266,7 +4266,7 @@ expect -223374.577355253 -111701.072127637
===============================================================================
# Orthographic
-# Azi, Sph.
+# Azi, Sph&Ell.
===============================================================================
-------------------------------------------------------------------------------
@@ -4421,6 +4421,72 @@ expect failure errno tolerance_condition
accept 0.70710678118 0.7071067812
expect 45 0
+
+-------------------------------------------------------------------------------
+# Test the ellipsoidal formulation
+
+# Test data from Guidance Note 7 part 2, March 2020, p. 90
+-------------------------------------------------------------------------------
+operation +proj=ortho +ellps=WGS84 +lat_0=55 +lon_0=5
+-------------------------------------------------------------------------------
+tolerance 1 mm
+accept 2.12955 53.80939444444444
+expect -189011.711 -128640.567
+roundtrip 1
+
+# Equatorial
+operation +proj=ortho +ellps=WGS84
+tolerance 0.1 mm
+accept 0 0
+expect 0 0
+roundtrip 1
+
+accept 0 89.99
+expect 0 6356752.2167
+roundtrip 1
+
+accept 0 -89.99
+expect 0 -6356752.2167
+roundtrip 1
+
+# Consistant with WGS84 semi-minor axis
+# The inverse transformation doesn't converge due to properties of the projection
+accept 0 90
+expect 0 6356752.3142
+
+# North pole tests
+operation +proj=ortho +ellps=WGS84 +lat_0=90
+tolerance 0.1 mm
+accept 0 90
+expect 0 0
+roundtrip 1
+
+accept 1 89.999
+expect 1.9493 -111.6770
+roundtrip 1
+
+# Consistant with WGS84 semi-major axis
+# The inverse transformation doesn't converge due to properties of the projection
+accept 0 0
+expect 0 -6378137
+
+# South pole tests
+operation +proj=ortho +ellps=WGS84 +lat_0=-90
+tolerance 0.1 mm
+accept 0 -90
+expect 0 0
+roundtrip 1
+
+accept 1 -89.999
+expect 1.9493 111.6770
+roundtrip 1
+
+# Consistant with WGS84 semi-major axis
+# The inverse transformation doesn't converge due to properties of the projection
+accept 0 0
+expect 0 6378137
+
+
===============================================================================
# Perspective Conic
# Conic, Sph
diff --git a/test/unit/test_io.cpp b/test/unit/test_io.cpp
index f2464e5c..34396fcf 100644
--- a/test/unit/test_io.cpp
+++ b/test/unit/test_io.cpp
@@ -5225,6 +5225,21 @@ static const struct {
{"False_Northing", 2},
{"Longitude_Of_Center", 3},
{"Latitude_Of_Center", 4}},
+ "Orthographic (Spherical)",
+ {
+ {"Latitude of natural origin", 4},
+ {"Longitude of natural origin", 3},
+ {"False easting", 1},
+ {"False northing", 2},
+ }},
+
+ {"Local",
+ {{"False_Easting", 1},
+ {"False_Northing", 2},
+ {"Scale_Factor", 1},
+ {"Azimuth", 0},
+ {"Longitude_Of_Center", 3},
+ {"Latitude_Of_Center", 4}},
"Orthographic",
{
{"Latitude of natural origin", 4},
@@ -5233,6 +5248,24 @@ static const struct {
{"False northing", 2},
}},
+ // Local with unsupported value for Azimuth
+ {"Local",
+ {{"False_Easting", 1},
+ {"False_Northing", 2},
+ {"Scale_Factor", 1},
+ {"Azimuth", 123},
+ {"Longitude_Of_Center", 3},
+ {"Latitude_Of_Center", 4}},
+ "Local",
+ {
+ {"False_Easting", 1},
+ {"False_Northing", 2},
+ {"Scale_Factor", 1},
+ {"Azimuth", 123},
+ {"Longitude_Of_Center", 3},
+ {"Latitude_Of_Center", 4},
+ }},
+
{"Winkel_Tripel",
{{"False_Easting", 1},
{"False_Northing", 2},
@@ -5543,7 +5576,6 @@ static const struct {
{"Longitude_Of_Origin", 3},
{"Latitude_Of_Origin", 4}},
},
-
};
TEST(wkt_parse, esri_projcs) {
@@ -9201,6 +9233,40 @@ TEST(io, projparse_non_earth_ellipsoid) {
// ---------------------------------------------------------------------------
+TEST(io, projparse_ortho_ellipsoidal) {
+ std::string input("+proj=ortho +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 "
+ "+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_ORTHOGRAPHIC);
+ EXPECT_EQ(
+ crs->exportToPROJString(
+ PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4)
+ .get()),
+ input);
+}
+
+// ---------------------------------------------------------------------------
+
+TEST(io, projparse_ortho_spherical) {
+ std::string input("+proj=ortho +f=0 +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 "
+ "+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()->nameStr(),
+ PROJ_WKT2_NAME_ORTHOGRAPHIC_SPHERICAL);
+ EXPECT_EQ(
+ crs->exportToPROJString(
+ PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4)
+ .get()),
+ input);
+}
+
+// ---------------------------------------------------------------------------
+
TEST(io, projparse_axisswap_unitconvert_longlat_proj) {
std::string input =
"+type=crs +proj=pipeline +step +proj=axisswap +order=2,1 +step "
generated by cgit v1.2.3 (git 2.25.1) at 2025-12-15 18:17:38 +0000