From 2e104e092578347de11208e9a3a80a3bf711265d Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Fri, 25 Sep 2020 22:53:38 +0200 Subject: Implement ellipsoidal formulation of +proj=ortho (fixes #397) - Map ESRI 'Local' to +proj=ortho when Scale_Factor = 1 and Azimuth = 0 - Map ESRI 'Orthographic' to a PROJ WKT2 'Orthographic (Spherical)' which maps to +proj=ortho +f=0 to froce spherical evaluation --- docs/source/operations/projections/ortho.rst | 10 ++- .../internal/coordinateoperation_constants.hpp | 3 + include/proj/internal/esri_projection_mappings.hpp | 19 ++++- scripts/build_esri_projection_mapping.py | 12 +++ src/iso19111/coordinateoperation.cpp | 2 +- src/iso19111/io.cpp | 37 ++++++++- src/proj_constants.h | 3 + src/projections/ortho.cpp | 87 ++++++++++++++++++++-- test/gie/builtins.gie | 68 ++++++++++++++++- test/unit/test_io.cpp | 68 ++++++++++++++++- 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 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 -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 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(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 " -- cgit v1.2.3 From 3d88b6fc89f95803bfd4d59c47eef1c05bda710c Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Sat, 26 Sep 2020 14:40:29 +0200 Subject: Ortho: add visibility condition for ellipsoidal case. Credits to @cffk --- src/projections/ortho.cpp | 15 +++++++++++++++ test/gie/builtins.gie | 41 +++++++++++++++++++++++++++++++++++++++++ 2 files changed, 56 insertions(+) diff --git a/src/projections/ortho.cpp b/src/projections/ortho.cpp index ff9d9341..0594205b 100644 --- a/src/projections/ortho.cpp +++ b/src/projections/ortho.cpp @@ -50,6 +50,13 @@ static PJ_XY ortho_s_forward (PJ_LP lp, PJ *P) { /* Spheroidal, forwar break; case OBLIQ: sinphi = sin(lp.phi); + + // Is the point visible from the projection plane ? + // From https://lists.osgeo.org/pipermail/proj/2020-September/009831.html + // this is the dot product of the normal of the ellipsoid at the center of + // the projection and at the point considered for projection. + // [cos(phi)*cos(lambda), cos(phi)*sin(lambda), sin(phi)] + // Also from Snyder's Map Projection - A working manual, equation (5-3), page 149 if (Q->sinph0 * sinphi + Q->cosph0 * cosphi * coslam < - EPS10) return forward_error(P, lp, xy); xy.y = Q->cosph0 * sinphi - Q->sinph0 * cosphi * coslam; @@ -131,6 +138,14 @@ static PJ_XY ortho_e_forward (PJ_LP lp, PJ *P) { /* Ellipsoidal, forwa const double sinphi = sin(lp.phi); const double coslam = cos(lp.lam); const double sinlam = sin(lp.lam); + + // Is the point visible from the projection plane ? + // Same condition as in spherical case + if (Q->sinph0 * sinphi + Q->cosph0 * cosphi * coslam < - EPS10) { + xy.x = HUGE_VAL; xy.y = HUGE_VAL; + return forward_error(P, lp, xy); + } + 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) + diff --git a/test/gie/builtins.gie b/test/gie/builtins.gie index 48afc00a..9d5b6644 100644 --- a/test/gie/builtins.gie +++ b/test/gie/builtins.gie @@ -4434,13 +4434,23 @@ 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 89.99 0 +expect 6378136.9029 0 +roundtrip 1 + +accept -89.99 0 +expect -6378136.9029 0 +roundtrip 1 + accept 0 89.99 expect 0 6356752.2167 roundtrip 1 @@ -4449,12 +4459,33 @@ accept 0 -89.99 expect 0 -6356752.2167 roundtrip 1 +# Point not visible from the projection plane +accept 90.00001 0 +expect failure errno tolerance_condition + +# Point not visible from the projection plane +accept -90.00001 0 +expect failure errno tolerance_condition + +# Consistant with WGS84 semi-major axis +# The inverse transformation doesn't converge due to properties of the projection +accept 90 0 +expect 6378137 0 + +accept -90 0 +expect -6378137 0 + # 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 +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 @@ -4465,12 +4496,18 @@ accept 1 89.999 expect 1.9493 -111.6770 roundtrip 1 +# Point not visible from the projection plane +accept 0 -0.0000001 +expect failure errno tolerance_condition + # 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 @@ -4481,6 +4518,10 @@ accept 1 -89.999 expect 1.9493 111.6770 roundtrip 1 +# Point not visible from the projection plane +accept 0 0.0000001 +expect failure errno tolerance_condition + # Consistant with WGS84 semi-major axis # The inverse transformation doesn't converge due to properties of the projection accept 0 0 -- cgit v1.2.3 From b04894819ea4e4d9d93e03015f0c7c9aa84642fe Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Sat, 26 Sep 2020 18:45:25 +0200 Subject: Ortho ellipsoidal inverse: add non iterative implementations for polar and equatorial --- src/projections/ortho.cpp | 60 +++++++++++++++++++++++++++++++++++++++++++++++ test/gie/builtins.gie | 59 +++++++++++++++++++++++++++++++++++++++------- 2 files changed, 111 insertions(+), 8 deletions(-) diff --git a/src/projections/ortho.cpp b/src/projections/ortho.cpp index 0594205b..ac54d88a 100644 --- a/src/projections/ortho.cpp +++ b/src/projections/ortho.cpp @@ -160,6 +160,66 @@ static PJ_LP ortho_e_inverse (PJ_XY xy, PJ *P) { /* Ellipsoidal, inver PJ_LP lp; struct pj_opaque *Q = static_cast(P->opaque); + if( Q->mode == N_POLE || Q->mode == S_POLE ) + { + // Polar case. Forward case equations can be simplified as: + // xy.x = nu * cosphi * sinlam + // xy.y = nu * -cosphi * coslam * sign(phi0) + // ==> lam = atan2(xy.x, -xy.y * sign(phi0)) + // ==> xy.x^2 + xy.y^2 = nu^2 * cosphi^2 + // rh^2 = cosphi^2 / (1 - es * sinphi^2) + // ==> sinphi^2 = (1 - rh^2) / (1 - es * rh^2) + + const double rh2 = xy.x * xy.x + xy.y * xy.y; + if (rh2 >= 1. - 1e-15) { + if ((rh2 - 1.) > EPS10) { + proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); + proj_log_trace(P, "Point (%.3f, %.3f) is outside the projection boundary"); + lp.lam = HUGE_VAL; lp.phi = HUGE_VAL; + return lp; + } + lp.phi = 0; + } + else { + lp.phi = asin(sqrt((1 - rh2) / (1 - P->es * rh2))) * (Q->mode == N_POLE ? 1 : -1); + } + lp.lam = atan2(xy.x, xy.y * (Q->mode == N_POLE ? -1 : 1)); + return lp; + } + + if( Q->mode == EQUIT ) + { + // Equatorial case. Forward case equations can be simplified as: + // xy.x = nu * cosphi * sinlam + // xy.y = nu * sinphi * (1 - P->es) + // x^2 * (1 - es * sinphi^2) = (1 - sinphi^2) * sinlam^2 + // y^2 / ((1 - es)^2 + y^2 * es) = sinphi^2 + + const auto SQ = [](double x) { return x*x; }; + + // Equation of the ellipse + if( SQ(xy.x) + SQ(xy.y * (P->a / P->b)) > 1 + 1e-11 ) { + proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); + proj_log_trace(P, "Point (%.3f, %.3f) is outside the projection boundary"); + lp.lam = HUGE_VAL; lp.phi = HUGE_VAL; + return lp; + } + + const double sinphi2 = xy.y == 0 ? 0 : 1.0 / (SQ((1 - P->es) / xy.y) + P->es); + if( sinphi2 > 1 - 1e-11 ) { + lp.phi = M_PI_2 * (xy.y > 0 ? 1 : -1); + lp.lam = 0; + return lp; + } + lp.phi = asin(sqrt(sinphi2)) * (xy.y > 0 ? 1 : -1); + const double sinlam = xy.x * sqrt((1 - P->es * sinphi2) / (1 - sinphi2)); + if( fabs(sinlam) - 1 > -1e-15 ) + lp.lam = M_PI_2 * (xy.x > 0 ? 1: -1); + else + lp.lam = asin(sinlam); + return lp; + } + // From EPSG guidance note 7.2 // It suggests as initial guess: diff --git a/test/gie/builtins.gie b/test/gie/builtins.gie index 9d5b6644..b63b7902 100644 --- a/test/gie/builtins.gie +++ b/test/gie/builtins.gie @@ -4270,7 +4270,7 @@ expect -223374.577355253 -111701.072127637 =============================================================================== ------------------------------------------------------------------------------- -# Test the equatorial aspect of the Orthopgraphic projection. +# Test the equatorial aspect of the Orthographic projection. # Test data from Snyder (1987), table 22, p. 151. ------------------------------------------------------------------------------- @@ -4322,7 +4322,7 @@ expect failure errno tolerance_condition ------------------------------------------------------------------------------- -# Test the oblique aspect of the Orthopgraphic projection. +# Test the oblique aspect of the Orthographic projection. # Test data from Snyder (1987), table 23, pp. 152-153. ------------------------------------------------------------------------------- @@ -4356,7 +4356,7 @@ expect failure errno tolerance_condition ------------------------------------------------------------------------------- -# Test the north polar aspect of the Orthopgraphic projection. +# Test the north polar aspect of the Orthographic projection. ------------------------------------------------------------------------------- operation +proj=ortho +R=1 +lat_0=90 +lon_0=0 ------------------------------------------------------------------------------- @@ -4386,7 +4386,7 @@ accept 2 2 expect failure errno tolerance_condition ------------------------------------------------------------------------------- -# Test the south polar aspect of the Orthopgraphic projection. +# Test the south polar aspect of the Orthographic projection. ------------------------------------------------------------------------------- operation +proj=ortho +R=1 +lat_0=-90 +lon_0=0 ------------------------------------------------------------------------------- @@ -4443,6 +4443,22 @@ accept 0 0 expect 0 0 roundtrip 1 +accept 1 1 +expect 111296.9991 110568.7748 +roundtrip 1 + +accept 1 -1 +expect 111296.9991 -110568.7748 +roundtrip 1 + +accept -1 1 +expect -111296.9991 110568.7748 +roundtrip 1 + +accept -1 -1 +expect -111296.9991 -110568.7748 +roundtrip 1 + accept 89.99 0 expect 6378136.9029 0 roundtrip 1 @@ -4468,20 +4484,37 @@ accept -90.00001 0 expect failure errno tolerance_condition # Consistant with WGS84 semi-major axis -# The inverse transformation doesn't converge due to properties of the projection accept 90 0 expect 6378137 0 +roundtrip 1 accept -90 0 expect -6378137 0 +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 +roundtrip 1 accept 0 -90 expect 0 -6356752.3142 +roundtrip 1 + +# Point not visible from the projection plane +direction inverse +accept 0 6356752.3143 +expect failure errno tolerance_condition + +# Point not visible from the projection plane +direction inverse +accept 1000 6356752.314 +expect failure errno tolerance_condition + +# Point not visible from the projection plane +direction inverse +accept 6378137.0001 0 +expect failure errno tolerance_condition ------------------------------------------------------------------------------- # North pole tests @@ -4501,9 +4534,14 @@ accept 0 -0.0000001 expect failure errno tolerance_condition # Consistant with WGS84 semi-major axis -# The inverse transformation doesn't converge due to properties of the projection accept 0 0 expect 0 -6378137 +roundtrip 1 + +# Point not visible from the projection plane +direction inverse +accept 0 -6378137.1 +expect failure errno tolerance_condition ------------------------------------------------------------------------------- # South pole tests @@ -4523,9 +4561,14 @@ accept 0 0.0000001 expect failure errno tolerance_condition # Consistant with WGS84 semi-major axis -# The inverse transformation doesn't converge due to properties of the projection accept 0 0 expect 0 6378137 +roundtrip 1 + +# Point not visible from the projection plane +direction inverse +accept 0 6378137.1 +expect failure errno tolerance_condition =============================================================================== -- cgit v1.2.3 From fedeeec68ff6a65126da35ae54ec75a719ff40ce Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Sun, 27 Sep 2020 13:27:01 +0200 Subject: Ortho ellipsoidal inverse: add domain check for oblique case, and slighly improve initial guessing --- src/projections/ortho.cpp | 26 +++++++++++++++++++++----- test/gie/builtins.gie | 40 ++++++++++++++++++++++++++++++++++++++++ 2 files changed, 61 insertions(+), 5 deletions(-) diff --git a/src/projections/ortho.cpp b/src/projections/ortho.cpp index ac54d88a..b4ecff7f 100644 --- a/src/projections/ortho.cpp +++ b/src/projections/ortho.cpp @@ -20,6 +20,8 @@ struct pj_opaque { double sinph0; double cosph0; double nu0; + double y_shift; + double y_scale; enum Mode mode; }; } // anonymous namespace @@ -155,11 +157,12 @@ static PJ_XY ortho_e_forward (PJ_LP lp, PJ *P) { /* Ellipsoidal, forwa } - static PJ_LP ortho_e_inverse (PJ_XY xy, PJ *P) { /* Ellipsoidal, inverse */ PJ_LP lp; struct pj_opaque *Q = static_cast(P->opaque); + const auto SQ = [](double x) { return x*x; }; + if( Q->mode == N_POLE || Q->mode == S_POLE ) { // Polar case. Forward case equations can be simplified as: @@ -170,7 +173,7 @@ static PJ_LP ortho_e_inverse (PJ_XY xy, PJ *P) { /* Ellipsoidal, inver // rh^2 = cosphi^2 / (1 - es * sinphi^2) // ==> sinphi^2 = (1 - rh^2) / (1 - es * rh^2) - const double rh2 = xy.x * xy.x + xy.y * xy.y; + const double rh2 = SQ(xy.x) + SQ(xy.y); if (rh2 >= 1. - 1e-15) { if ((rh2 - 1.) > EPS10) { proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); @@ -195,8 +198,6 @@ static PJ_LP ortho_e_inverse (PJ_XY xy, PJ *P) { /* Ellipsoidal, inver // x^2 * (1 - es * sinphi^2) = (1 - sinphi^2) * sinlam^2 // y^2 / ((1 - es)^2 + y^2 * es) = sinphi^2 - const auto SQ = [](double x) { return x*x; }; - // Equation of the ellipse if( SQ(xy.x) + SQ(xy.y * (P->a / P->b)) > 1 + 1e-11 ) { proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); @@ -220,13 +221,26 @@ static PJ_LP ortho_e_inverse (PJ_XY xy, PJ *P) { /* Ellipsoidal, inver return lp; } + // Using Q->sinph0 * sinphi + Q->cosph0 * cosphi * coslam == 0 (visibity + // condition of the forward case) in the forward equations, and a lot of + // substition games... + PJ_XY xy_recentered; + xy_recentered.x = xy.x; + xy_recentered.y = (xy.y - Q->y_shift) / Q->y_scale; + if( SQ(xy.x) + SQ(xy_recentered.y) > 1 + 1e-11 ) { + proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); + proj_log_trace(P, "Point (%.3f, %.3f) is outside the projection boundary"); + lp.lam = HUGE_VAL; lp.phi = HUGE_VAL; + return lp; + } + // 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); + lp = ortho_s_inverse(xy_recentered, P); for( int i = 0; i < 20; i++ ) { @@ -286,6 +300,8 @@ PJ *PROJECTION(ortho) { else { Q->nu0 = 1.0 / sqrt(1.0 - P->es * Q->sinph0 * Q->sinph0); + Q->y_shift = P->es * Q->nu0 * Q->sinph0 * Q->cosph0; + Q->y_scale = 1.0 / sqrt(1.0 - P->es * Q->cosph0 * Q->cosph0); P->inv = ortho_e_inverse; P->fwd = ortho_e_forward; } diff --git a/test/gie/builtins.gie b/test/gie/builtins.gie index b63b7902..71622c7a 100644 --- a/test/gie/builtins.gie +++ b/test/gie/builtins.gie @@ -4434,6 +4434,46 @@ accept 2.12955 53.80939444444444 expect -189011.711 -128640.567 roundtrip 1 +------------------------------------------------------------------------------- +# Oblique +------------------------------------------------------------------------------- +operation +proj=ortho +ellps=WGS84 +lat_0=30 + +# On boundary of visibility domain. +direction forward +tolerance 0.1 mm +accept -90 0 +expect -6378137 18504.1253 + +# This test is fragile. Note the slighly important tolerance +# direction inverse +# tolerance 100 mm +# accept -6378137 18504.125313223721605027 +# expect -90 0 + +# Slightly outside +direction inverse +accept -6378137.001 18504.1253 +expect failure errno tolerance_condition + +# On boundary of visibility domain +direction forward +tolerance 0.1 mm +accept 0 -60 +expect 0 -6343601.0991 + +# Just on it, but fails to converge. This test might be fragile +direction inverse +accept 0 -6343601.099075031466782093 +expect failure errno non_convergent + +# Slightly inside +direction inverse +tolerance 0.1 mm +accept 0 -6343600 +expect 0 -59.966377950099655436 + + ------------------------------------------------------------------------------- # Equatorial ------------------------------------------------------------------------------- -- cgit v1.2.3 From 5e7d58ec69764f7ffa8d09991aa0b0102d2f02ab Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Sun, 27 Sep 2020 15:58:54 +0200 Subject: Ortho ellipsoidal inverse: improve accuracy in polar case with (x,y) close to (0,0) --- src/projections/ortho.cpp | 4 ++-- test/gie/builtins.gie | 12 ++++++++---- 2 files changed, 10 insertions(+), 6 deletions(-) diff --git a/src/projections/ortho.cpp b/src/projections/ortho.cpp index b4ecff7f..555f3f71 100644 --- a/src/projections/ortho.cpp +++ b/src/projections/ortho.cpp @@ -171,7 +171,7 @@ static PJ_LP ortho_e_inverse (PJ_XY xy, PJ *P) { /* Ellipsoidal, inver // ==> lam = atan2(xy.x, -xy.y * sign(phi0)) // ==> xy.x^2 + xy.y^2 = nu^2 * cosphi^2 // rh^2 = cosphi^2 / (1 - es * sinphi^2) - // ==> sinphi^2 = (1 - rh^2) / (1 - es * rh^2) + // ==> cosphi^2 = rh^2 * (1 - es) / (1 - es * rh^2) const double rh2 = SQ(xy.x) + SQ(xy.y); if (rh2 >= 1. - 1e-15) { @@ -184,7 +184,7 @@ static PJ_LP ortho_e_inverse (PJ_XY xy, PJ *P) { /* Ellipsoidal, inver lp.phi = 0; } else { - lp.phi = asin(sqrt((1 - rh2) / (1 - P->es * rh2))) * (Q->mode == N_POLE ? 1 : -1); + lp.phi = acos(sqrt(rh2 * P->one_es / (1 - P->es * rh2))) * (Q->mode == N_POLE ? 1 : -1); } lp.lam = atan2(xy.x, xy.y * (Q->mode == N_POLE ? -1 : 1)); return lp; diff --git a/test/gie/builtins.gie b/test/gie/builtins.gie index 71622c7a..02e6766e 100644 --- a/test/gie/builtins.gie +++ b/test/gie/builtins.gie @@ -4565,8 +4565,12 @@ accept 0 90 expect 0 0 roundtrip 1 -accept 1 89.999 -expect 1.9493 -111.6770 +accept 30 45 +expect 2258795.4394 -3912348.4650 +roundtrip 1 + +accept 135 89.999999873385 +expect 0.01 0.01 roundtrip 1 # Point not visible from the projection plane @@ -4592,8 +4596,8 @@ accept 0 -90 expect 0 0 roundtrip 1 -accept 1 -89.999 -expect 1.9493 111.6770 +accept 135 -89.999999873385 +expect 0.01 -0.01 roundtrip 1 # Point not visible from the projection plane -- cgit v1.2.3 From 016e8e60747a2def2dfd7ffc7f6ad2e6aa8ba009 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Wed, 30 Sep 2020 11:17:13 +0200 Subject: ortho.cpp: more precise reference to guidance note --- src/projections/ortho.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/projections/ortho.cpp b/src/projections/ortho.cpp index 555f3f71..75b8199d 100644 --- a/src/projections/ortho.cpp +++ b/src/projections/ortho.cpp @@ -135,7 +135,7 @@ static PJ_XY ortho_e_forward (PJ_LP lp, PJ *P) { /* Ellipsoidal, forwa PJ_XY xy; struct pj_opaque *Q = static_cast(P->opaque); - // From EPSG guidance note 7.2 + // From EPSG guidance note 7.2, March 2020, §3.3.5 Orthographic const double cosphi = cos(lp.phi); const double sinphi = sin(lp.phi); const double coslam = cos(lp.lam); @@ -234,7 +234,7 @@ static PJ_LP ortho_e_inverse (PJ_XY xy, PJ *P) { /* Ellipsoidal, inver return lp; } - // From EPSG guidance note 7.2 + // From EPSG guidance note 7.2, March 2020, §3.3.5 Orthographic // It suggests as initial guess: // lp.lam = 0; -- cgit v1.2.3