Fix build with -DPROJ_INTERNAL_CPP_NAMESPACE

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 }

docs/source/apps/geod.rst

@@ -4,6 +4,12 @@
 geod
 ================================================================================
 
+.. _invgeod:
+
+================================================================================
+invgeod
+================================================================================
+
 Synopsis
 ********
 

docs/source/apps/proj.rst

@@ -4,6 +4,12 @@
 proj
 ================================================================================
 
+.. _invproj:
+
+================================================================================
+invproj
+================================================================================
+
 .. only:: html
 
     Cartographic projection filter.

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 <https://gdal.org/programs/gdalsrsinfo.html>`__ utility that offers
 partly similar services.
 
 

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 <https://clang-analyzer.llvm.org/scan-build.html>`__ 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
 

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.

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 <https://gdal.org/>`_ as a frontend to PROJ and transform your data with the
-    :program:`ogr2ogr` application.
+    `ogr2ogr <https://gdal.org/programs/ogr2ogr.html>`__ 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 <projinfo>` 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?
 --------------------------------------------------------------------------------

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
 

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                                                       |

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

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:
 

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:
 
 ::
 

src/apps/cct.cpp

@@ -78,6 +78,7 @@ Thomas Knudsen, [email protected], 2016-05-25/2017-10-26
 #include <string.h>
 #include <stdarg.h>
 
+#include <algorithm>
 #include <fstream> // std::ifstream
 #include <iostream>
 
@@ -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<PJ_LOG_LEVEL>(verbose));
     }

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 );

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;

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<std::pair<std::string, std::string>> 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);

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<const crs::CompoundCRS *>(horiz);

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, */