From 9df6fd0323634ff8a32eb334628ba429b3d4f0c8 Mon Sep 17 00:00:00 2001 From: Alexandre Barreto Date: Tue, 3 Jan 2023 11:19:22 -0500 Subject: [PATCH] Fix build with -DPROJ_INTERNAL_CPP_NAMESPACE --- docs/source/_extensions/program_with_link.py | 14 +++ docs/source/apps/geod.rst | 6 + docs/source/apps/proj.rst | 6 + docs/source/apps/projinfo.rst | 2 +- docs/source/community/code_contributions.rst | 6 +- docs/source/conf.py | 1 + docs/source/faq.rst | 5 +- docs/source/install.rst | 4 +- docs/source/operations/projections/aeqd.rst | 2 + docs/source/operations/projections/gnom.rst | 23 +++- docs/source/resource_files.rst | 8 +- docs/source/usage/quickstart.rst | 2 +- src/apps/cct.cpp | 3 +- src/conversions/unitconvert.cpp | 22 +++- src/ell_set.cpp | 2 +- src/iso19111/crs.cpp | 18 +-- src/iso19111/operation/conversion.cpp | 6 +- src/pipeline.cpp | 5 +- src/projections/aeqd.cpp | 39 +++---- src/projections/gnom.cpp | 105 ++++++++++++++++-- src/transformations/gridshift.cpp | 4 +- src/transformations/helmert.cpp | 12 +- src/transformations/hgridshift.cpp | 30 ++++- src/transformations/molodensky.cpp | 18 ++- src/transformations/vgridshift.cpp | 30 ++++- test/gie/builtins.gie | 109 ++++++++++++++++++- test/unit/test_c_api.cpp | 18 ++- test/unit/test_crs.cpp | 34 +++++- test/unit/test_io.cpp | 13 ++- 29 files changed, 450 insertions(+), 97 deletions(-) create mode 100644 docs/source/_extensions/program_with_link.py diff --git a/docs/source/_extensions/program_with_link.py b/docs/source/_extensions/program_with_link.py new file mode 100644 index 000000000000..62b31ef06faa --- /dev/null +++ b/docs/source/_extensions/program_with_link.py @@ -0,0 +1,14 @@ +from sphinx.domains.std import StandardDomain + +class MyStandardDomain(StandardDomain): + + def resolve_xref(self, env, fromdocname, builder, typ, target, node, contnode): + if typ == 'program': + typ = 'ref' + return StandardDomain.resolve_xref(self, env, fromdocname, builder, typ, target, node, contnode) + +def setup(app): + # Override the "program" role to be an alias of "ref" + MyStandardDomain.roles["program"] = StandardDomain.roles["ref"] + app.add_domain(MyStandardDomain, override=True) + return { 'parallel_read_safe': True, 'parallel_write_safe': True } diff --git a/docs/source/apps/geod.rst b/docs/source/apps/geod.rst index 4c0a8c2f7191..efda807247b7 100644 --- a/docs/source/apps/geod.rst +++ b/docs/source/apps/geod.rst @@ -4,6 +4,12 @@ geod ================================================================================ +.. _invgeod: + +================================================================================ +invgeod +================================================================================ + Synopsis ******** diff --git a/docs/source/apps/proj.rst b/docs/source/apps/proj.rst index 3229fb1ab652..10361736013c 100644 --- a/docs/source/apps/proj.rst +++ b/docs/source/apps/proj.rst @@ -4,6 +4,12 @@ proj ================================================================================ +.. _invproj: + +================================================================================ +invproj +================================================================================ + .. only:: html Cartographic projection filter. diff --git a/docs/source/apps/projinfo.rst b/docs/source/apps/projinfo.rst index 9805d0e6b8bc..e68cf26e7d9a 100644 --- a/docs/source/apps/projinfo.rst +++ b/docs/source/apps/projinfo.rst @@ -72,7 +72,7 @@ or PROJJSON string). It can also be used to query coordinate operations available between two CRS. -The program is named with some reference to the GDAL :program:`gdalsrsinfo` that offers +The program is named with some reference to the GDAL `gdalsrsinfo `__ utility that offers partly similar services. diff --git a/docs/source/community/code_contributions.rst b/docs/source/community/code_contributions.rst index a21b590b8308..5318578576dc 100644 --- a/docs/source/community/code_contributions.rst +++ b/docs/source/community/code_contributions.rst @@ -125,21 +125,21 @@ Preliminary step: install clang. For example: mv clang+llvm-9.0.0-x86_64-linux-gnu-ubuntu-18.04 clang+llvm-9 export PATH=$PWD/clang+llvm-9/bin:$PATH -Configure PROJ with the :program:`scan-build` utility of clang: +Configure PROJ with the `scan-build `__ utility of clang: :: mkdir csa_build cd csa_build scan-build cmake .. -Build using :program:`scan-build`: +Build using ``scan-build``: :: scan-build make [-j8] If CSA finds errors, they will be emitted during the build. And in which case, -at the end of the build process, :program:`scan-build` will emit a warning message +at the end of the build process, ``scan-build`` will emit a warning message indicating errors have been found and how to display the error report. This is with something like diff --git a/docs/source/conf.py b/docs/source/conf.py index 8980039cdb9b..9ff199fe8c0d 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -38,6 +38,7 @@ 'sphinxcontrib.spelling', 'breathe', 'redirects', + 'program_with_link', ] # Add any paths that contain templates here, relative to this directory. diff --git a/docs/source/faq.rst b/docs/source/faq.rst index 09920c40d50b..fa758b7cf298 100644 --- a/docs/source/faq.rst +++ b/docs/source/faq.rst @@ -23,7 +23,7 @@ separate column. If your data is stored in a common geodata file format chances are that you can use `GDAL `_ as a frontend to PROJ and transform your data with the - :program:`ogr2ogr` application. + `ogr2ogr `__ application. Can I transform from *abc* to *xyz*? -------------------------------------------------------------------------------- @@ -49,8 +49,7 @@ ETRS89/UTM32N (EPSG:25832) and ETRS89/DKTM1 (EPSG:4093): +step +proj=tmerc +lat_0=0 +lon_0=9 +k=0.99998 +x_0=200000 +y_0=-5000000 +ellps=GRS80 -See the :program:`projinfo` :ref:`documentation ` for more info on -how to use it. +See the :program:`projinfo` documentation for more info on how to use it. Coordinate reference system *xyz* is not in the EPSG registry, what do I do? -------------------------------------------------------------------------------- diff --git a/docs/source/install.rst b/docs/source/install.rst index b9cbc3a5f470..edd1d32dfaa2 100644 --- a/docs/source/install.rst +++ b/docs/source/install.rst @@ -102,14 +102,14 @@ On Debian and similar systems (e.g. Ubuntu) the APT package manager is used:: Fedora ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ -On Fedora the :program:`dnf` package manager is used:: +On Fedora the ``dnf`` package manager is used:: sudo dnf install proj Red Hat ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ -On Red Hat based system packages are installed with :program:`yum`:: +On Red Hat based system packages are installed with ``yum``:: sudo yum install proj diff --git a/docs/source/operations/projections/aeqd.rst b/docs/source/operations/projections/aeqd.rst index 3c5566ec96a4..2afadd12fdb3 100644 --- a/docs/source/operations/projections/aeqd.rst +++ b/docs/source/operations/projections/aeqd.rst @@ -8,6 +8,8 @@ Azimuthal Equidistant +---------------------+----------------------------------------------------------+ | **Available forms** | Forward and inverse, spherical and ellipsoidal | +---------------------+----------------------------------------------------------+ +| **Defined area** | Global | ++---------------------+----------------------------------------------------------+ | **Alias** | aeqd | +---------------------+----------------------------------------------------------+ | **Domain** | 2D | diff --git a/docs/source/operations/projections/gnom.rst b/docs/source/operations/projections/gnom.rst index 1bbb5034c0e3..06a9edbf99b7 100644 --- a/docs/source/operations/projections/gnom.rst +++ b/docs/source/operations/projections/gnom.rst @@ -4,12 +4,21 @@ Gnomonic ******************************************************************************** +For a sphere, the gnomonic projection is a projection from the center of +the sphere onto a plane tangent to the center point of the projection. +This projects great circles to straight lines. For an ellipsoid, it is +the limit of a doubly azimuthal projection, a projection where the +azimuths from 2 points are preserved, as the two points merge into the +center point. In this case, geodesics project to approximately straight +lines (these are exactly straight if the geodesic includes the center +point). For details, see Section 8 of :cite:`Karney2013`. + +---------------------+----------------------------------------------------------+ -| **Classification** | Pseudocylindrical | +| **Classification** | Azimuthal | +---------------------+----------------------------------------------------------+ -| **Available forms** | Forward and inverse, spherical projection | +| **Available forms** | Forward and inverse, spherical and ellipsoidal | +---------------------+----------------------------------------------------------+ -| **Defined area** | Global | +| **Defined area** | Within a quarter circumference of the center point | +---------------------+----------------------------------------------------------+ | **Alias** | gnom | +---------------------+----------------------------------------------------------+ @@ -26,7 +35,7 @@ Gnomonic :align: center :alt: Gnomonic - proj-string: ``+proj=gnom +lat_0=90 +lon_0=-50`` + proj-string: ``+proj=gnom +lat_0=90 +lon_0=-50 +R=6.4e6`` Parameters ################################################################################ @@ -37,8 +46,10 @@ Parameters .. include:: ../options/lat_0.rst -.. include:: ../options/R.rst - .. include:: ../options/x_0.rst .. include:: ../options/y_0.rst + +.. include:: ../options/ellps.rst + +.. include:: ../options/R.rst diff --git a/docs/source/resource_files.rst b/docs/source/resource_files.rst index 28797609de09..750c0e5a76d8 100644 --- a/docs/source/resource_files.rst +++ b/docs/source/resource_files.rst @@ -323,7 +323,7 @@ compiler. For Ubuntu something like the following should work. apt-get install gfortran To compile the program do something like the following to produce the binary -:program:`htdp` from the source code. +``htdp`` from the source code. :: @@ -369,13 +369,13 @@ Usage The goal of :file:`crs2crs2grid.py` is to produce a grid shift file for a designated region. The region is defined using the ``-griddef`` switch. When missing a continental US region is used. The script creates a set of sample points for -the grid definition, runs :program:`htdp` against it and then parses the +the grid definition, runs ``htdp`` against it and then parses the resulting points and computes a point by point shift to encode into the final -grid shift file. By default it is assumed that :program:`htdp` is in the +grid shift file. By default it is assumed that ``htdp`` is in the executable path. If not, please provide the path to the executable using the ``-htdp`` switch. -The :program:`htdp` program supports transformations between many CRSes and for each (or +The ``htdp`` program supports transformations between many CRSes and for each (or most?) of them you need to provide a date at which the CRS is fixed. The full set of CRS Ids available in the HTDP program are: diff --git a/docs/source/usage/quickstart.rst b/docs/source/usage/quickstart.rst index 85e7aa95dc8e..8d4e46f79c1f 100644 --- a/docs/source/usage/quickstart.rst +++ b/docs/source/usage/quickstart.rst @@ -35,7 +35,7 @@ utility :program:`proj` we can convert the geodetic coordinates to projected spa If called as above :program:`proj` will be in interactive mode, letting you type the input data manually and getting a response presented on screen. :program:`proj` works as any UNIX filter though, which means that you can also -pipe data to the utility, for instance by using the :program:`echo` command: +pipe data to the utility, for instance by using the ``echo`` command: :: diff --git a/src/apps/cct.cpp b/src/apps/cct.cpp index 62a3714d3f42..f5e385921fa1 100644 --- a/src/apps/cct.cpp +++ b/src/apps/cct.cpp @@ -78,6 +78,7 @@ Thomas Knudsen, thokn@sdfe.dk, 2016-05-25/2017-10-26 #include #include +#include #include // std::ifstream #include @@ -236,7 +237,7 @@ int main(int argc, char **argv) { PJ_DIRECTION direction = opt_given (o, "I")? PJ_INV: PJ_FWD; - verbose = MIN(opt_given (o, "v"), 3); /* log level can't be larger than 3 */ + verbose = std::min(opt_given (o, "v"), 3); /* log level can't be larger than 3 */ if( verbose > 0 ) { proj_log_level (PJ_DEFAULT_CTX, static_cast(verbose)); } diff --git a/src/conversions/unitconvert.cpp b/src/conversions/unitconvert.cpp index 417d1ae4f0a2..eb88b996796e 100644 --- a/src/conversions/unitconvert.cpp +++ b/src/conversions/unitconvert.cpp @@ -323,6 +323,10 @@ static PJ_XYZ forward_3d(PJ_LPZ lpz, PJ *P) { point.lpz = lpz; /* take care of the horizontal components in the 2D function */ + + // Assigning in 2 steps avoids cppcheck warning + // "Overlapping read/write of union is undefined behavior" + // Cf https://github.com/OSGeo/PROJ/pull/3527#pullrequestreview-1233332710 const auto xy = forward_2d(point.lp, P); point.xy = xy; @@ -341,6 +345,10 @@ static PJ_LPZ reverse_3d(PJ_XYZ xyz, PJ *P) { point.xyz = xyz; /* take care of the horizontal components in the 2D function */ + + // Assigning in 2 steps avoids cppcheck warning + // "Overlapping read/write of union is undefined behavior" + // Cf https://github.com/OSGeo/PROJ/pull/3527#pullrequestreview-1233332710 const auto lp = reverse_2d(point.xy, P); point.lp = lp; @@ -358,7 +366,12 @@ static void forward_4d(PJ_COORD& coo, PJ *P) { struct pj_opaque_unitconvert *Q = (struct pj_opaque_unitconvert *) P->opaque; /* delegate unit conversion of physical dimensions to the 3D function */ - coo.xyz = forward_3d(coo.lpz, P); + + // Assigning in 2 steps avoids cppcheck warning + // "Overlapping read/write of union is undefined behavior" + // Cf https://github.com/OSGeo/PROJ/pull/3527#pullrequestreview-1233332710 + const auto xyz = forward_3d(coo.lpz, P); + coo.xyz = xyz; if (Q->t_in_id >= 0) coo.xyzt.t = time_units[Q->t_in_id].t_in( coo.xyzt.t ); @@ -375,7 +388,12 @@ static void reverse_4d(PJ_COORD& coo, PJ *P) { struct pj_opaque_unitconvert *Q = (struct pj_opaque_unitconvert *) P->opaque; /* delegate unit conversion of physical dimensions to the 3D function */ - coo.lpz = reverse_3d(coo.xyz, P); + + // Assigning in 2 steps avoids cppcheck warning + // "Overlapping read/write of union is undefined behavior" + // Cf https://github.com/OSGeo/PROJ/pull/3527#pullrequestreview-1233332710 + const auto lpz = reverse_3d(coo.xyz, P); + coo.lpz = lpz; if (Q->t_out_id >= 0) coo.xyzt.t = time_units[Q->t_out_id].t_in( coo.xyzt.t ); diff --git a/src/ell_set.cpp b/src/ell_set.cpp index f86bf65a28d7..a4c4d7ac1d99 100644 --- a/src/ell_set.cpp +++ b/src/ell_set.cpp @@ -584,7 +584,7 @@ int pj_calc_ellipsoid_params (PJ *P, double a, double es) { /* flattening */ if (0==P->f) P->f = 1 - cos (P->alpha); /* = 1 - sqrt (1 - PIN->es); */ - if (P->f == 1.0) { + if (!(P->f >= 0.0 && P->f < 1.0)) { proj_log_error(P, _("Invalid eccentricity")); proj_errno_set (P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE); return PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE; diff --git a/src/iso19111/crs.cpp b/src/iso19111/crs.cpp index 6dbb8011799a..a63d76ead5c7 100644 --- a/src/iso19111/crs.cpp +++ b/src/iso19111/crs.cpp @@ -4829,7 +4829,8 @@ ProjectedCRS::identify(const io::AuthorityFactoryPtr &authorityFactory) const { } const bool l_implicitCS = hasImplicitCS(); - const auto addCRS = [&](const ProjectedCRSNNPtr &crs, const bool eqName) { + const auto addCRS = [&](const ProjectedCRSNNPtr &crs, const bool eqName, + bool hasNonMatchingId) { const auto &l_unit = cs->axisList()[0]->unit(); if (_isEquivalentTo(crs.get(), util::IComparable::Criterion:: @@ -4849,7 +4850,7 @@ ProjectedCRS::identify(const io::AuthorityFactoryPtr &authorityFactory) const { util::IComparable::Criterion::EQUIVALENT, dbContext))) { if (crs->nameStr() == thisName) { res.clear(); - res.emplace_back(crs, 100); + res.emplace_back(crs, hasNonMatchingId ? 70 : 100); } else { res.emplace_back(crs, eqName ? 90 : 70); } @@ -4890,6 +4891,7 @@ ProjectedCRS::identify(const io::AuthorityFactoryPtr &authorityFactory) const { ci_equal(thisName, "unnamed"); bool foundEquivalentName = false; + bool hasNonMatchingId = false; if (hasCodeCompatibleOfAuthorityFactory(this, authorityFactory)) { // If the CRS has already an id, check in the database for the // official object, and verify that they are equivalent. @@ -4906,11 +4908,14 @@ ProjectedCRS::identify(const io::AuthorityFactoryPtr &authorityFactory) const { EQUIVALENT_EXCEPT_AXIS_ORDER_GEOGCRS, dbContext); res.emplace_back(crs, match ? 100 : 25); - return res; + if (match) { + return res; + } } catch (const std::exception &) { } } } + hasNonMatchingId = true; } else if (!insignificantName) { for (int ipass = 0; ipass < 2; ipass++) { const bool approximateMatch = ipass == 1; @@ -4926,7 +4931,7 @@ ProjectedCRS::identify(const io::AuthorityFactoryPtr &authorityFactory) const { thisName.c_str(), pairObjName.second.c_str()); foundEquivalentName |= eqName; - if (addCRS(crsNN, eqName).second == 100) { + if (addCRS(crsNN, eqName, false).second == 100) { return res; } } @@ -4962,8 +4967,7 @@ ProjectedCRS::identify(const io::AuthorityFactoryPtr &authorityFactory) const { // Sort results res.sort(lambdaSort); - if (!hasCodeCompatibleOfAuthorityFactory(this, authorityFactory) && - !foundEquivalentName && (res.empty() || res.front().second < 50)) { + if (!foundEquivalentName && (res.empty() || res.front().second < 50)) { std::set> alreadyKnown; for (const auto &pair : res) { const auto &ids = pair.first->identifiers(); @@ -4985,7 +4989,7 @@ ProjectedCRS::identify(const io::AuthorityFactoryPtr &authorityFactory) const { continue; } - addCRS(crs, insignificantName); + addCRS(crs, insignificantName, hasNonMatchingId); } res.sort(lambdaSort); diff --git a/src/iso19111/operation/conversion.cpp b/src/iso19111/operation/conversion.cpp index 6a1e19c6378a..0151f5fc9507 100644 --- a/src/iso19111/operation/conversion.cpp +++ b/src/iso19111/operation/conversion.cpp @@ -3985,8 +3985,10 @@ void Conversion::_exportToPROJString( bAxisSpecFound = true; } - // No need to add explicit f=0 if the ellipsoid is a sphere - if (strcmp(mapping->proj_name_aux, "f=0") == 0) { + // No need to add explicit f=0 or R_A if the ellipsoid is a + // sphere + if (strcmp(mapping->proj_name_aux, "f=0") == 0 || + strcmp(mapping->proj_name_aux, "R_A") == 0) { crs::CRS *horiz = l_sourceCRS.get(); const auto compound = dynamic_cast(horiz); diff --git a/src/pipeline.cpp b/src/pipeline.cpp index 8351e7c87128..9f8ac844b493 100644 --- a/src/pipeline.cpp +++ b/src/pipeline.cpp @@ -376,9 +376,8 @@ static void set_ellipsoid(PJ *P) { P->a_orig = P->a; P->es_orig = P->es; - pj_calc_ellipsoid_params (P, P->a, P->es); - - geod_init(P->geod, P->a, P->es / (1 + sqrt(P->one_es))); + if( pj_calc_ellipsoid_params (P, P->a, P->es) == 0 ) + geod_init(P->geod, P->a, P->f); /* Re-attach the dangling list */ /* Note: cur will always be non 0 given argv_sentinel presence, */ diff --git a/src/projections/aeqd.cpp b/src/projections/aeqd.cpp index 92ee1b8676b6..e619840940ac 100644 --- a/src/projections/aeqd.cpp +++ b/src/projections/aeqd.cpp @@ -96,16 +96,16 @@ static PJ_XY aeqd_e_forward (PJ_LP lp, PJ *P) { /* Ellipsoidal, forward struct pj_opaque *Q = static_cast(P->opaque); double coslam, cosphi, sinphi, rho; double azi1, azi2, s12; - double lam1, phi1, lam2, phi2; + double lat1, lon1, lat2, lon2; coslam = cos(lp.lam); - cosphi = cos(lp.phi); - sinphi = sin(lp.phi); switch (Q->mode) { case N_POLE: coslam = - coslam; PROJ_FALLTHROUGH; case S_POLE: + cosphi = cos(lp.phi); + sinphi = sin(lp.phi); rho = fabs(Q->Mp - pj_mlfn(lp.phi, sinphi, cosphi, Q->en)); xy.x = rho * sin(lp.lam); xy.y = rho * coslam; @@ -117,15 +117,15 @@ static PJ_XY aeqd_e_forward (PJ_LP lp, PJ *P) { /* Ellipsoidal, forward break; } - phi1 = P->phi0 / DEG_TO_RAD; - lam1 = P->lam0 / DEG_TO_RAD; - phi2 = lp.phi / DEG_TO_RAD; - lam2 = (lp.lam+P->lam0) / DEG_TO_RAD; + lat1 = P->phi0 / DEG_TO_RAD; + lon1 = 0; + lat2 = lp.phi / DEG_TO_RAD; + lon2 = lp.lam / DEG_TO_RAD; - geod_inverse(&Q->g, phi1, lam1, phi2, lam2, &s12, &azi1, &azi2); + geod_inverse(&Q->g, lat1, lon1, lat2, lon2, &s12, &azi1, &azi2); azi1 *= DEG_TO_RAD; - xy.x = s12 * sin(azi1) / P->a; - xy.y = s12 * cos(azi1) / P->a; + xy.x = s12 * sin(azi1); + xy.y = s12 * cos(azi1); break; } return xy; @@ -231,28 +231,23 @@ static PJ_LP e_guam_inv(PJ_XY xy, PJ *P) { /* Guam elliptical */ static PJ_LP aeqd_e_inverse (PJ_XY xy, PJ *P) { /* Ellipsoidal, inverse */ PJ_LP lp = {0.0,0.0}; struct pj_opaque *Q = static_cast(P->opaque); - double c; - double azi1, azi2, s12, x2, y2, lat1, lon1, lat2, lon2; + double azi1, azi2, s12, lat1, lon1, lat2, lon2; - if ((c = hypot(xy.x, xy.y)) < EPS10) { + if ((s12 = hypot(xy.x, xy.y)) < EPS10) { lp.phi = P->phi0; lp.lam = 0.; return (lp); } if (Q->mode == OBLIQ || Q->mode == EQUIT) { - - x2 = xy.x * P->a; - y2 = xy.y * P->a; lat1 = P->phi0 / DEG_TO_RAD; - lon1 = P->lam0 / DEG_TO_RAD; - azi1 = atan2(x2, y2) / DEG_TO_RAD; - s12 = sqrt(x2 * x2 + y2 * y2); + lon1 = 0; + azi1 = atan2(xy.x, xy.y) / DEG_TO_RAD; // Clockwise from north geod_direct(&Q->g, lat1, lon1, azi1, s12, &lat2, &lon2, &azi2); lp.phi = lat2 * DEG_TO_RAD; lp.lam = lon2 * DEG_TO_RAD; - lp.lam -= P->lam0; } else { /* Polar */ - lp.phi = pj_inv_mlfn(Q->mode == N_POLE ? Q->Mp - c : Q->Mp + c, Q->en); + lp.phi = pj_inv_mlfn(Q->mode == N_POLE ? Q->Mp - s12 : Q->Mp + s12, + Q->en); lp.lam = atan2(xy.x, Q->mode == N_POLE ? -xy.y : xy.y); } return lp; @@ -308,7 +303,7 @@ PJ *PROJECTION(aeqd) { P->opaque = Q; P->destructor = destructor; - geod_init(&Q->g, P->a, P->es / (1 + sqrt(P->one_es))); + geod_init(&Q->g, 1, P->f); if (fabs(fabs(P->phi0) - M_HALFPI) < EPS10) { Q->mode = P->phi0 < 0. ? S_POLE : N_POLE; diff --git a/src/projections/gnom.cpp b/src/projections/gnom.cpp index f56ca2bdde34..d7ff556bfd16 100644 --- a/src/projections/gnom.cpp +++ b/src/projections/gnom.cpp @@ -2,10 +2,11 @@ #include #include +#include +#include "geodesic.h" #include "proj.h" #include "proj_internal.h" -#include PROJ_HEAD(gnom, "Gnomonic") "\n\tAzi, Sph"; @@ -25,6 +26,7 @@ struct pj_opaque { double sinph0; double cosph0; enum Mode mode; + struct geod_geodesic g; }; } // anonymous namespace @@ -122,6 +124,80 @@ static PJ_LP gnom_s_inverse (PJ_XY xy, PJ *P) { /* Spheroidal, inverse return lp; } +static PJ_XY gnom_e_forward (PJ_LP lp, PJ *P) { /* Ellipsoidal, forward */ + PJ_XY xy = {0.0,0.0}; + struct pj_opaque *Q = static_cast(P->opaque); + + double + lat0 = P->phi0 / DEG_TO_RAD, + lon0 = 0, + lat1 = lp.phi / DEG_TO_RAD, + lon1 = lp.lam / DEG_TO_RAD, + azi0, m, M; + + geod_geninverse(&Q->g, lat0, lon0, lat1, lon1, + nullptr, &azi0, nullptr, &m, &M, nullptr, nullptr); + if (M <= 0) { + proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN); + xy.x = xy.y = HUGE_VAL; + } + else { + double rho = m/M; + azi0 *= DEG_TO_RAD; + xy.x = rho * sin(azi0); + xy.y = rho * cos(azi0); + } + return xy; +} + + +static PJ_LP gnom_e_inverse (PJ_XY xy, PJ *P) { /* Ellipsoidal, inverse */ + constexpr int numit_ = 10; + static const double eps_ = 0.01 * sqrt(DBL_EPSILON); + + PJ_LP lp = {0.0,0.0}; + struct pj_opaque *Q = static_cast(P->opaque); + + double + lat0 = P->phi0 / DEG_TO_RAD, + lon0 = 0, + azi0 = atan2(xy.x, xy.y) / DEG_TO_RAD, // Clockwise from north + rho = hypot(xy.x, xy.y), + s = atan(rho); + bool little = rho <= 1; + if (!little) + rho = 1/rho; + struct geod_geodesicline l; + geod_lineinit(&l, &Q->g, lat0, lon0, azi0, + GEOD_LATITUDE | GEOD_LONGITUDE | GEOD_DISTANCE_IN | + GEOD_REDUCEDLENGTH | GEOD_GEODESICSCALE); + + int count = numit_, trip = 0; + double lat1 = 0, lon1 = 0; + while (count--) { + double m, M; + geod_genposition(&l, 0, s, + &lat1, &lon1, nullptr, &s, &m, &M, nullptr, nullptr); + if (trip) + break; + // If little, solve rho(s) = rho with drho(s)/ds = 1/M^2 + // else solve 1/rho(s) = 1/rho with d(1/rho(s))/ds = -1/m^2 + double ds = little ? (m - rho * M) * M : (rho * m - M) * m; + s -= ds; + // Reversed test to allow escape with NaNs + if (!(fabs(ds) >= eps_)) + ++trip; + } + if (trip) { + lp.phi = lat1 * DEG_TO_RAD; + lp.lam = lon1 * DEG_TO_RAD; + } else { + proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN); + lp.phi = lp.lam = HUGE_VAL; + } + return lp; +} + PJ *PROJECTION(gnom) { struct pj_opaque *Q = static_cast(calloc (1, sizeof (struct pj_opaque))); @@ -129,18 +205,25 @@ PJ *PROJECTION(gnom) { return pj_default_destructor (P, PROJ_ERR_OTHER /*ENOMEM*/); P->opaque = Q; - 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 = EQUIT; + if (P->es == 0) { + 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 = EQUIT; + } else { + Q->mode = OBLIQ; + Q->sinph0 = sin(P->phi0); + Q->cosph0 = cos(P->phi0); + } + + P->inv = gnom_s_inverse; + P->fwd = gnom_s_forward; } else { - Q->mode = OBLIQ; - Q->sinph0 = sin(P->phi0); - Q->cosph0 = cos(P->phi0); - } + geod_init(&Q->g, 1, P->f); - P->inv = gnom_s_inverse; - P->fwd = gnom_s_forward; + P->inv = gnom_e_inverse; + P->fwd = gnom_e_forward; + } P->es = 0.; return P; diff --git a/src/transformations/gridshift.cpp b/src/transformations/gridshift.cpp index bce61d23800b..209c380ee41b 100644 --- a/src/transformations/gridshift.cpp +++ b/src/transformations/gridshift.cpp @@ -491,7 +491,7 @@ PJ_LPZ gridshiftData::grid_interpolate(PJ_CONTEXT *ctx, const std::string &type, // --------------------------------------------------------------------------- static PJ_LP normalizeLongitude(const GenericShiftGrid *grid, const PJ_LPZ in, - const osgeo::proj::ExtentAndRes *&extentOut) { + const NS_PROJ::ExtentAndRes *&extentOut) { PJ_LP normalized; normalized.lam = in.lam; normalized.phi = in.phi; @@ -520,7 +520,7 @@ PJ_LPZ gridshiftData::grid_apply_internal( return in; /* normalized longitude of input */ - const osgeo::proj::ExtentAndRes *extent; + const NS_PROJ::ExtentAndRes *extent; PJ_LP normalized_in = normalizeLongitude(grid, in, extent); PJ_LPZ shift = grid_interpolate(ctx, type, normalized_in, grid); diff --git a/src/transformations/helmert.cpp b/src/transformations/helmert.cpp index 027cd4b09756..0f5e8f7de52c 100644 --- a/src/transformations/helmert.cpp +++ b/src/transformations/helmert.cpp @@ -450,7 +450,11 @@ static void helmert_forward_4d (PJ_COORD &point, PJ *P) { build_rot_matrix(P); } - point.xyz = helmert_forward_3d (point.lpz, P); + // Assigning in 2 steps avoids cppcheck warning + // "Overlapping read/write of union is undefined behavior" + // Cf https://github.com/OSGeo/PROJ/pull/3527#pullrequestreview-1233332710 + const auto xyz = helmert_forward_3d (point.lpz, P); + point.xyz = xyz; } @@ -466,7 +470,11 @@ static void helmert_reverse_4d (PJ_COORD& point, PJ *P) { build_rot_matrix(P); } - point.lpz = helmert_reverse_3d (point.xyz, P); + // Assigning in 2 steps avoids cppcheck warning + // "Overlapping read/write of union is undefined behavior" + // Cf https://github.com/OSGeo/PROJ/pull/3527#pullrequestreview-1233332710 + const auto lpz = helmert_reverse_3d (point.xyz, P); + point.lpz = lpz; } /* Arcsecond to radians */ diff --git a/src/transformations/hgridshift.cpp b/src/transformations/hgridshift.cpp index 9ef8f0e1ea06..884cb569822e 100644 --- a/src/transformations/hgridshift.cpp +++ b/src/transformations/hgridshift.cpp @@ -75,13 +75,22 @@ static void forward_4d(PJ_COORD& coo, PJ *P) { /* If transformation is not time restricted, we always call it */ if (Q->t_final==0 || Q->t_epoch==0) { - coo.xyz = forward_3d (coo.lpz, P); + // Assigning in 2 steps avoids cppcheck warning + // "Overlapping read/write of union is undefined behavior" + // Cf https://github.com/OSGeo/PROJ/pull/3527#pullrequestreview-1233332710 + const auto xyz = forward_3d (coo.lpz, P); + coo.xyz = xyz; return; } /* Time restricted - only apply transform if within time bracket */ - if (coo.lpzt.t < Q->t_epoch && Q->t_final > Q->t_epoch) - coo.xyz = forward_3d (coo.lpz, P); + if (coo.lpzt.t < Q->t_epoch && Q->t_final > Q->t_epoch) { + // Assigning in 2 steps avoids cppcheck warning + // "Overlapping read/write of union is undefined behavior" + // Cf https://github.com/OSGeo/PROJ/pull/3527#pullrequestreview-1233332710 + const auto xyz = forward_3d (coo.lpz, P); + coo.xyz = xyz; + } } static void reverse_4d(PJ_COORD& coo, PJ *P) { @@ -89,13 +98,22 @@ static void reverse_4d(PJ_COORD& coo, PJ *P) { /* If transformation is not time restricted, we always call it */ if (Q->t_final==0 || Q->t_epoch==0) { - coo.lpz = reverse_3d (coo.xyz, P); + // Assigning in 2 steps avoids cppcheck warning + // "Overlapping read/write of union is undefined behavior" + // Cf https://github.com/OSGeo/PROJ/pull/3527#pullrequestreview-1233332710 + const auto lpz = reverse_3d (coo.xyz, P); + coo.lpz = lpz; return; } /* Time restricted - only apply transform if within time bracket */ - if (coo.lpzt.t < Q->t_epoch && Q->t_final > Q->t_epoch) - coo.lpz = reverse_3d (coo.xyz, P); + if (coo.lpzt.t < Q->t_epoch && Q->t_final > Q->t_epoch) { + // Assigning in 2 steps avoids cppcheck warning + // "Overlapping read/write of union is undefined behavior" + // Cf https://github.com/OSGeo/PROJ/pull/3527#pullrequestreview-1233332710 + const auto lpz = reverse_3d (coo.xyz, P); + coo.lpz = lpz; + } } static PJ *destructor (PJ *P, int errlev) { diff --git a/src/transformations/molodensky.cpp b/src/transformations/molodensky.cpp index d78a2db20c7c..429ca2aac846 100644 --- a/src/transformations/molodensky.cpp +++ b/src/transformations/molodensky.cpp @@ -215,6 +215,9 @@ static PJ_XY forward_2d(PJ_LP lp, PJ *P) { PJ_COORD point = {{0,0,0,0}}; point.lp = lp; + // Assigning in 2 steps avoids cppcheck warning + // "Overlapping read/write of union is undefined behavior" + // Cf https://github.com/OSGeo/PROJ/pull/3527#pullrequestreview-1233332710 const auto xyz = forward_3d(point.lpz, P); point.xyz = xyz; @@ -227,6 +230,9 @@ static PJ_LP reverse_2d(PJ_XY xy, PJ *P) { point.xy = xy; point.xyz.z = 0; + // Assigning in 2 steps avoids cppcheck warning + // "Overlapping read/write of union is undefined behavior" + // Cf https://github.com/OSGeo/PROJ/pull/3527#pullrequestreview-1233332710 const auto lpz = reverse_3d(point.xyz, P); point.lpz = lpz; @@ -261,7 +267,11 @@ static PJ_XYZ forward_3d(PJ_LPZ lpz, PJ *P) { static void forward_4d(PJ_COORD& obs, PJ *P) { - obs.xyz = forward_3d(obs.lpz, P); + // Assigning in 2 steps avoids cppcheck warning + // "Overlapping read/write of union is undefined behavior" + // Cf https://github.com/OSGeo/PROJ/pull/3527#pullrequestreview-1233332710 + const auto xyz = forward_3d(obs.lpz, P); + obs.xyz = xyz; } @@ -292,7 +302,11 @@ static PJ_LPZ reverse_3d(PJ_XYZ xyz, PJ *P) { static void reverse_4d(PJ_COORD& obs, PJ *P) { - obs.lpz = reverse_3d(obs.xyz, P); + // Assigning in 2 steps avoids cppcheck warning + // "Overlapping read/write of union is undefined behavior" + // Cf https://github.com/OSGeo/PROJ/pull/3527#pullrequestreview-1233332710 + const auto lpz = reverse_3d(obs.xyz, P); + obs.lpz = lpz; } diff --git a/src/transformations/vgridshift.cpp b/src/transformations/vgridshift.cpp index 6218553d93cd..0abac0e67914 100644 --- a/src/transformations/vgridshift.cpp +++ b/src/transformations/vgridshift.cpp @@ -106,13 +106,22 @@ static void forward_4d(PJ_COORD& coo, PJ *P) { /* If transformation is not time restricted, we always call it */ if (Q->t_final==0 || Q->t_epoch==0) { - coo.xyz = forward_3d (coo.lpz, P); + // Assigning in 2 steps avoids cppcheck warning + // "Overlapping read/write of union is undefined behavior" + // Cf https://github.com/OSGeo/PROJ/pull/3527#pullrequestreview-1233332710 + const auto xyz = forward_3d (coo.lpz, P); + coo.xyz = xyz; return; } /* Time restricted - only apply transform if within time bracket */ - if (coo.lpzt.t < Q->t_epoch && Q->t_final > Q->t_epoch) - coo.xyz = forward_3d (coo.lpz, P); + if (coo.lpzt.t < Q->t_epoch && Q->t_final > Q->t_epoch) { + // Assigning in 2 steps avoids cppcheck warning + // "Overlapping read/write of union is undefined behavior" + // Cf https://github.com/OSGeo/PROJ/pull/3527#pullrequestreview-1233332710 + const auto xyz = forward_3d (coo.lpz, P); + coo.xyz = xyz; + } } static void reverse_4d(PJ_COORD& coo, PJ *P) { @@ -120,13 +129,22 @@ static void reverse_4d(PJ_COORD& coo, PJ *P) { /* If transformation is not time restricted, we always call it */ if (Q->t_final==0 || Q->t_epoch==0) { - coo.lpz = reverse_3d (coo.xyz, P); + // Assigning in 2 steps avoids cppcheck warning + // "Overlapping read/write of union is undefined behavior" + // Cf https://github.com/OSGeo/PROJ/pull/3527#pullrequestreview-1233332710 + const auto lpz = reverse_3d (coo.xyz, P); + coo.lpz = lpz; return; } /* Time restricted - only apply transform if within time bracket */ - if (coo.lpzt.t < Q->t_epoch && Q->t_final > Q->t_epoch) - coo.lpz = reverse_3d (coo.xyz, P); + if (coo.lpzt.t < Q->t_epoch && Q->t_final > Q->t_epoch) { + // Assigning in 2 steps avoids cppcheck warning + // "Overlapping read/write of union is undefined behavior" + // Cf https://github.com/OSGeo/PROJ/pull/3527#pullrequestreview-1233332710 + const auto lpz = reverse_3d (coo.xyz, P); + coo.lpz = lpz; + } } static PJ *destructor (PJ *P, int errlev) { diff --git a/test/gie/builtins.gie b/test/gie/builtins.gie index 61f236f42a85..3fe88e7d13d1 100644 --- a/test/gie/builtins.gie +++ b/test/gie/builtins.gie @@ -1767,11 +1767,12 @@ expect -0.001790493 -0.000895247 =============================================================================== # Gnomonic -# Azi, Sph. +# Azi, Sph*Ell =============================================================================== ------------------------------------------------------------------------------- -# Test material from Snyder p. 168, table 26. +# Test material from Snyder p. 168, table 26. Repeat tests with ellispoid of +# flattening 1/200. # Tests the equatorial aspect of the projection. ------------------------------------------------------------------------------- operation +proj=gnom +R=1 @@ -1815,6 +1816,50 @@ direction inverse accept 0 1e8 expect 0 90 +------------------------------------------------------------------------------- +operation +proj=gnom +a=1 +rf=200 +------------------------------------------------------------------------------- +tolerance 0.1 mm +accept 0 0 +expect 0 0 +roundtrip 100 +accept 10 80 +expect 0.1763 5.7232 +roundtrip 100 +accept 20 70 +expect 0.3641 2.9037 +roundtrip 100 +accept 30 60 +expect 0.5778 1.9861 +roundtrip 100 +accept 40 50 +expect 0.8405 1.5459 +roundtrip 100 +accept 50 40 +expect 1.1958 1.2994 +roundtrip 100 +accept 60 30 +expect 1.7435 1.1534 +roundtrip 100 +accept 70 20 +expect 2.7852 1.0711 +roundtrip 100 +accept 80 10 +expect 5.8813 1.0465 +roundtrip 100 +accept 80 80 +expect 5.7134 32.7298 +roundtrip 100 +accept 0 89.99 +expect 0 5700.9222 +roundtrip 100 +accept 180 89.99 +expect failure errno coord_transfm_outside_projection_domain + +# test that extreme northings are mapped to the sphere +direction inverse +accept 0 1e8 +expect 0 90 ------------------------------------------------------------------------------- # Test the northern polar aspect of the gnonomic projection @@ -1833,6 +1878,27 @@ expect failure errno coord_transfm_outside_projection_domain accept 90 0 expect failure errno coord_transfm_outside_projection_domain +------------------------------------------------------------------------------- +operation +proj=gnom +a=1 +rf=200 +lat_0=90 +------------------------------------------------------------------------------- +tolerance 0.1 mm +accept 0 90 +expect 0 0 +roundtrip 100 +accept 45 45 +expect 0.7079 -0.7079 +roundtrip 100 +accept 0 0 +expect 0 -127.4835 +roundtrip 100 +accept 90 0 +expect 127.4835 0 +roundtrip 100 +accept 0 -0.5 +expect failure errno coord_transfm_outside_projection_domain +accept 90 -0.5 +expect failure errno coord_transfm_outside_projection_domain + ------------------------------------------------------------------------------- # Test the southern polar aspect of the gnonomic projection ------------------------------------------------------------------------------- @@ -1850,6 +1916,27 @@ expect failure errno coord_transfm_outside_projection_domain accept 90 0 expect failure errno coord_transfm_outside_projection_domain +------------------------------------------------------------------------------- +operation +proj=gnom +a=1 +rf=200 +lat_0=-90 +------------------------------------------------------------------------------- +tolerance 0.1 mm +accept 0 -90 +expect 0 0 +roundtrip 100 +accept 45 -45 +expect 0.7079 0.7079 +roundtrip 100 +accept 0 0 +expect 0 127.4835 +roundtrip 100 +accept 90 0 +expect 127.4835 0 +roundtrip 100 +accept 0 0.5 +expect failure errno coord_transfm_outside_projection_domain +accept 90 0.5 +expect failure errno coord_transfm_outside_projection_domain + ------------------------------------------------------------------------------- # Test the oblique aspect of the gnonomic projection ------------------------------------------------------------------------------- @@ -1868,6 +1955,24 @@ roundtrip 100 accept 0 -45 expect failure errno coord_transfm_outside_projection_domain +------------------------------------------------------------------------------- +operation +proj=gnom +a=1 +rf=200 +lat_0=45 +------------------------------------------------------------------------------- +tolerance 0.1 mm +accept 0 45 +expect 0 0 +roundtrip 100 +accept 0 0 +expect 0 -0.9897 +roundtrip 100 +accept 0 90 +expect 0 1.0025 +roundtrip 100 +accept 0 -45 +expect 0 -154.8623 +roundtrip 100 +accept 0 -45.5 +expect failure errno coord_transfm_outside_projection_domain =============================================================================== # Goode Homolosine # PCyl, Sph. diff --git a/test/unit/test_c_api.cpp b/test/unit/test_c_api.cpp index 0900941f60f3..9b53468fa6d9 100644 --- a/test/unit/test_c_api.cpp +++ b/test/unit/test_c_api.cpp @@ -6107,12 +6107,26 @@ TEST_F(CApi, proj_trans_bounds_ignore_inf) { double out_bottom; double out_right; double out_top; + // Before the ellipsoidal version of the gnomonic projection was + // implemented the WGS84 ellipsoid was treated as a sphere of radius + // 6378137m and the equator was the "horizon" for the projection with the + // south polar aspect. + // + // The boundary with ndiv = 21 then mapped into a line extending to ymin = + // -89178007.2 which was the projection of lat = -90d/(ndiv+1), lon = 180d. + // + // With the implementation of the ellipsoidal gnonomic projection, the + // horizon is now at lat = +0.3035d. + // + // We move the north edge of the box to lat = -90+4.15*(ndiv+1) = +1.3d. + // The northernmost point on the boundary which is within the horizon is + // now lat = -90+4.15*ndiv = -2.85d, lon = 180d for which y = -116576598.5. int success = - proj_trans_bounds(m_ctxt, P, PJ_FWD, -180.0, -90.0, 180.0, 0.0, + proj_trans_bounds(m_ctxt, P, PJ_FWD, -180.0, -90.0, 180.0, 1.3, &out_left, &out_bottom, &out_right, &out_top, 21); EXPECT_TRUE(success == 1); EXPECT_NEAR(out_left, 0, 1); - EXPECT_NEAR(out_bottom, -89178008, 1); + EXPECT_NEAR(out_bottom, -116576598.5, 1); EXPECT_NEAR(out_right, 0, 1); EXPECT_NEAR(out_top, 0, 1); } diff --git a/test/unit/test_crs.cpp b/test/unit/test_crs.cpp index 053700a66edb..9257902273c7 100644 --- a/test/unit/test_crs.cpp +++ b/test/unit/test_crs.cpp @@ -2617,7 +2617,9 @@ TEST(crs, projectedCRS_identify_db) { sourceCRS->baseCRS(), sourceCRS->derivingConversion(), sourceCRS->coordinateSystem()); auto res = crs->identify(factoryEPSG); - EXPECT_EQ(res.size(), 0U); + EXPECT_EQ(res.size(), 1U); + EXPECT_EQ(res.front().first->getEPSGCode(), 2172); + EXPECT_EQ(res.front().second, 70); } { // Existing code, but not matching content @@ -2632,8 +2634,8 @@ TEST(crs, projectedCRS_identify_db) { sourceCRS->coordinateSystem()); auto res = crs->identify(factoryEPSG); ASSERT_EQ(res.size(), 1U); - EXPECT_EQ(res.front().first->getEPSGCode(), 32631); - EXPECT_EQ(res.front().second, 25); + EXPECT_EQ(res.front().first->getEPSGCode(), 2172); + EXPECT_EQ(res.front().second, 70); } { // Identify by exact name @@ -3252,6 +3254,32 @@ TEST(crs, projectedCRS_identify_db) { EXPECT_EQ(res.front().first->getEPSGCode(), 8353); EXPECT_EQ(res.front().second, 100); } + { + // Identify from a pseudo WKT ESRI with has an AUTHORITY node that + // points to another object. + // Cf + // https://lists.osgeo.org/pipermail/qgis-user/2023-January/052299.html + auto obj = WKTParser().attachDatabaseContext(dbContext).createFromWKT( + "PROJCS[\"ETRS_1989_UTM_Zone_32N_6Stellen\"," + "GEOGCS[\"GCS_ETRS_1989\",DATUM[\"D_ETRS_1989\"," + "SPHEROID[\"GRS_1980\",6378137.0,298.257222101]]," + "PRIMEM[\"Greenwich\",0.0],UNIT[\"Degree\",0.0174532925199433]]," + "PROJECTION[\"Transverse_Mercator\"]," + "PARAMETER[\"False_Easting\",500000.0]," + "PARAMETER[\"False_Northing\",0.0]," + "PARAMETER[\"Central_Meridian\",9.0]," + "PARAMETER[\"Scale_Factor\",0.9996]," + "PARAMETER[\"Latitude_Of_Origin\",0.0]," + "UNIT[\"Meter\",1.0]," + "AUTHORITY[\"ESRI\",\"102328\"]]"); + auto crs = nn_dynamic_pointer_cast(obj); + ASSERT_TRUE(crs != nullptr); + auto allFactory = AuthorityFactory::create(dbContext, std::string()); + auto res = crs->identify(allFactory); + ASSERT_GE(res.size(), 1U); + EXPECT_EQ(res.front().first->getEPSGCode(), 25832); + EXPECT_EQ(res.front().second, 70); + } } // --------------------------------------------------------------------------- diff --git a/test/unit/test_io.cpp b/test/unit/test_io.cpp index da1e12e08429..59ad3d93865d 100644 --- a/test/unit/test_io.cpp +++ b/test/unit/test_io.cpp @@ -10017,12 +10017,19 @@ TEST(io, projparse_aeqd_guam) { // --------------------------------------------------------------------------- TEST(io, projparse_cea_spherical) { - auto obj = PROJStringParser().createFromPROJString( - "+proj=cea +R=6371228 +type=crs"); + const std::string input( + "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +R=6371228 +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_LAMBERT_CYLINDRICAL_EQUAL_AREA_SPHERICAL); + EXPECT_EQ( + crs->exportToPROJString( + PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4) + .get()), + input); auto crs2 = ProjectedCRS::create( PropertyMap(), crs->baseCRS(), @@ -10041,7 +10048,7 @@ TEST(io, projparse_cea_spherical) { // --------------------------------------------------------------------------- TEST(io, projparse_cea_spherical_on_ellipsoid) { - std::string input("+proj=cea +R_A +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 " + std::string input("+proj=cea +R_A +lat_ts=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);