diff --git a/docs/source/operations/projections/ortho.rst b/docs/source/operations/projections/ortho.rst index 2e8de9caf4..2d0e08a8bf 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 79b47a6a08..7038ba9fe5 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 752c03ad94..c230f26b6c 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 5206e19156..12eeb8525c 100644 --- a/scripts/build_esri_projection_mapping.py +++ b/scripts/build_esri_projection_mapping.py @@ -436,11 +436,23 @@ - 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 65ef2b7776..ef262cc8e7 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 6680ab00d0..16ab22f751 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 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 5b62622337..f64bf4965d 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 5f9366c3fb..ff9d934151 100644 --- a/src/projections/ortho.cpp +++ b/src/projections/ortho.cpp @@ -4,7 +4,7 @@ #include "proj_internal.h" #include -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(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(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(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 f755526eff..48afc00ad7 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 f2464e5c8a..34396fcfaf 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(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(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 "