Merge branch 'master' into update-grid-alternatives.sql-for-Icealnd

9 files changed, 144 insertions, 27 deletions

diff --git a/src/4D_api.cpp b/src/4D_api.cpp
index 26e9456e..f37594b5 100644
--- a/src/4D_api.cpp
+++ b/src/4D_api.cpp
@@ -1557,7 +1557,7 @@ PJ_GRID_INFO proj_grid_info(const char *gridname) {
     memset(&grinfo, 0, sizeof(PJ_GRID_INFO));
 
     /* in case the grid wasn't found */
-    if (gridinfo->filename == nullptr) {
+    if (gridinfo->filename == nullptr || gridinfo->ct == nullptr) {
         pj_gridinfo_free(ctx, gridinfo);
         strcpy(grinfo.format, "missing");
         return grinfo;
diff --git a/src/Makefile.am b/src/Makefile.am
index 0a209f9b..2e14d1dd 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -11,7 +11,7 @@ AM_CPPFLAGS =	-DPROJ_LIB=\"$(pkgdatadir)\" \
 AM_CXXFLAGS =    @CXX_WFLAGS@ @FLTO_FLAG@
 
 include_HEADERS = proj.h proj_experimental.h proj_constants.h proj_api.h geodesic.h \
-	org_proj4_PJ.h proj_symbol_rename.h proj_json_streaming_writer.hpp
+	org_proj4_PJ.h proj_symbol_rename.h
 
 EXTRA_DIST = bin_cct.cmake bin_gie.cmake bin_cs2cs.cmake \
 	bin_geod.cmake bin_proj.cmake bin_projinfo.cmake \
@@ -210,6 +210,7 @@ libproj_la_SOURCES = \
 	wkt2_parser.h wkt2_parser.cpp \
 	wkt2_generated_parser.h wkt2_generated_parser.c \
 	\
+	proj_json_streaming_writer.hpp \
 	proj_json_streaming_writer.cpp \
 	\
 	tracing.cpp
diff --git a/src/geodesic.h b/src/geodesic.h
index 5d230531..e2265c89 100644
--- a/src/geodesic.h
+++ b/src/geodesic.h
@@ -158,7 +158,7 @@
                       GEODESIC_VERSION_PATCH)
 
 #if !defined(GEOD_DLL)
-#if defined(_MSC_VER)
+#if defined(_MSC_VER) && defined(PROJ_MSVC_DLL_EXPORT)
 #define GEOD_DLL __declspec(dllexport)
 #elif defined(__GNUC__)
 #define GEOD_DLL __attribute__ ((visibility("default")))
diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp
index 27f11100..0469cae1 100644
--- a/src/iso19111/coordinateoperation.cpp
+++ b/src/iso19111/coordinateoperation.cpp
@@ -11301,9 +11301,17 @@ struct FilterResults {
                         setOfSetOfGrids.end()) {
                         continue;
                     }
+
+                    const bool sameNameOrEmptyName =
+                        ((!curExtent && !lastExtent) ||
+                         (curExtent && lastExtent &&
+                          !curExtent->description()->empty() &&
+                          *(curExtent->description()) ==
+                              *(lastExtent->description())));
+
                     // If we have already found a operation without grids for
                     // that extent, no need to add any lower accuracy operation
-                    if (!lastHasGrids) {
+                    if (!lastHasGrids && sameNameOrEmptyName) {
                         continue;
                     }
                     // If we had only operations involving grids, but one
@@ -11580,11 +11588,13 @@ CoordinateOperationFactory::Private::findOpsInRegistryDirect(
 
             const auto authorities(getCandidateAuthorities(
                 authFactory, srcAuthName, targetAuthName));
+            std::vector<CoordinateOperationNNPtr> res;
             for (const auto &authority : authorities) {
+                const auto authName =
+                    authority == "any" ? std::string() : authority;
                 const auto tmpAuthFactory = io::AuthorityFactory::create(
-                    authFactory->databaseContext(),
-                    authority == "any" ? std::string() : authority);
-                auto res =
+                    authFactory->databaseContext(), authName);
+                auto resTmp =
                     tmpAuthFactory->createFromCoordinateReferenceSystemCodes(
                         srcAuthName, srcCode, targetAuthName, targetCode,
                         context.context->getUsePROJAlternativeGridNames(),
@@ -11593,6 +11603,10 @@ CoordinateOperationFactory::Private::findOpsInRegistryDirect(
                                 DISCARD_OPERATION_IF_MISSING_GRID,
                         context.context->getDiscardSuperseded(), true, false,
                         context.extent1, context.extent2);
+                res.insert(res.end(), resTmp.begin(), resTmp.end());
+                if (authName == "PROJ") {
+                    continue;
+                }
                 if (!res.empty()) {
                     resNonEmptyBeforeFiltering = true;
                     auto resFiltered =
@@ -12566,6 +12580,18 @@ void CoordinateOperationFactory::Private::createOperationsWithDatumPivot(
 
     // Start in priority with candidates that have exactly the same name as
     // the sourcCRS and targetCRS. Typically for the case of init=IGNF:XXXX
+
+    // Transformation from IGNF:NTFP to IGNF:RGF93G,
+    // using
+    // NTF geographiques Paris (gr) vers NTF GEOGRAPHIQUES GREENWICH (DMS) +
+    // NOUVELLE TRIANGULATION DE LA FRANCE (NTF) vers RGF93 (ETRS89)
+    // that is using ntf_r93.gsb, is horribly dependent
+    // of IGNF:RGF93G being returned before IGNF:RGF93GEO in candidatesDstGeod.
+    // If RGF93GEO is returned before then we go through WGS84 and use
+    // instead a Helmert transformation.
+    // The below logic is thus quite fragile, and attempts at changing it
+    // result in degraded results for other use cases...
+
     for (const auto &candidateSrcGeod : candidatesSrcGeod) {
         if (candidateSrcGeod->nameStr() == sourceCRS->nameStr()) {
             for (const auto &candidateDstGeod : candidatesDstGeod) {
@@ -12619,7 +12645,7 @@ void CoordinateOperationFactory::Private::createOperationsWithDatumPivot(
 #endif
             createTransformations(candidateSrcGeod, candidateDstGeod,
                                   opsFirst[0], isNullFirst);
-            if (!res.empty() && hasResultSetOnlyResultsWithPROJStep(res)) {
+            if (!res.empty() && !hasResultSetOnlyResultsWithPROJStep(res)) {
                 return;
             }
         }
@@ -12917,9 +12943,40 @@ bool CoordinateOperationFactory::Private::createOperationsFromDatabase(
     bool sameGeodeticDatum = false;
 
     if (vertSrc || vertDst) {
-        createOperationsFromDatabaseWithVertCRS(sourceCRS, targetCRS, context,
-                                                geogSrc, geogDst, vertSrc,
-                                                vertDst, res);
+        if (res.empty()) {
+            if (geogSrc &&
+                geogSrc->coordinateSystem()->axisList().size() == 2 &&
+                vertDst) {
+                auto dbContext =
+                    context.context->getAuthorityFactory()->databaseContext();
+                auto resTmp = findOpsInRegistryDirect(
+                    sourceCRS->promoteTo3D(std::string(), dbContext), targetCRS,
+                    context, resFindDirectNonEmptyBeforeFiltering);
+                for (auto &op : resTmp) {
+                    auto newOp = op->shallowClone();
+                    setCRSs(newOp.get(), sourceCRS, targetCRS);
+                    res.emplace_back(newOp);
+                }
+            } else if (geogDst &&
+                       geogDst->coordinateSystem()->axisList().size() == 2 &&
+                       vertSrc) {
+                auto dbContext =
+                    context.context->getAuthorityFactory()->databaseContext();
+                auto resTmp = findOpsInRegistryDirect(
+                    sourceCRS, targetCRS->promoteTo3D(std::string(), dbContext),
+                    context, resFindDirectNonEmptyBeforeFiltering);
+                for (auto &op : resTmp) {
+                    auto newOp = op->shallowClone();
+                    setCRSs(newOp.get(), sourceCRS, targetCRS);
+                    res.emplace_back(newOp);
+                }
+            }
+        }
+        if (res.empty()) {
+            createOperationsFromDatabaseWithVertCRS(sourceCRS, targetCRS,
+                                                    context, geogSrc, geogDst,
+                                                    vertSrc, vertDst, res);
+        }
     } else if (geodSrc && geodDst) {
 
         const auto &srcDatum = geodSrc->datum();
@@ -13400,7 +13457,7 @@ void CoordinateOperationFactory::Private::createOperationsGeodToGeod(
                 util::IComparable::Criterion::EQUIVALENT)) {
             res.emplace_back(
                 Conversion::createGeographicGeocentric(sourceCRS, targetCRS));
-        } else if (isSrcGeocentric) {
+        } else if (isSrcGeocentric && geogDst) {
             std::string interm_crs_name(geogDst->nameStr());
             interm_crs_name += " (geocentric)";
             auto interm_crs =
@@ -13533,7 +13590,7 @@ void CoordinateOperationFactory::Private::createOperationsBoundToGeog(
             }
         }
         // If the datum are equivalent, this is also fine
-    } else if (geogCRSOfBaseOfBoundSrc && hubSrcGeog->datum() &&
+    } else if (geogCRSOfBaseOfBoundSrc && hubSrcGeog && hubSrcGeog->datum() &&
                geogDst->datum() &&
                hubSrcGeog->datum()->_isEquivalentTo(
                    geogDst->datum().get(),
@@ -13562,7 +13619,7 @@ void CoordinateOperationFactory::Private::createOperationsBoundToGeog(
         // Case of "+proj=latlong +ellps=clrk66
         // +nadgrids=ntv1_can.dat,conus"
         // to "+proj=latlong +datum=NAD83"
-    } else if (geogCRSOfBaseOfBoundSrc && hubSrcGeog->datum() &&
+    } else if (geogCRSOfBaseOfBoundSrc && hubSrcGeog && hubSrcGeog->datum() &&
                geogDst->datum() &&
                geogCRSOfBaseOfBoundSrc->ellipsoid()->_isEquivalentTo(
                    datum::Ellipsoid::CLARKE_1866.get(),
diff --git a/src/iso19111/factory.cpp b/src/iso19111/factory.cpp
index 57850303..f0ad157e 100644
--- a/src/iso19111/factory.cpp
+++ b/src/iso19111/factory.cpp
@@ -1433,9 +1433,17 @@ AuthorityFactory::AuthorityFactory(const DatabaseContextNNPtr &context,
 AuthorityFactoryNNPtr
 AuthorityFactory::create(const DatabaseContextNNPtr &context,
                          const std::string &authorityName) {
-
-    auto factory = AuthorityFactory::nn_make_shared<AuthorityFactory>(
-        context, authorityName);
+    const auto getFactory = [&context, &authorityName]() {
+        for (const auto &knownName : {"EPSG", "ESRI", "PROJ"}) {
+            if (ci_equal(authorityName, knownName)) {
+                return AuthorityFactory::nn_make_shared<AuthorityFactory>(
+                    context, knownName);
+            }
+        }
+        return AuthorityFactory::nn_make_shared<AuthorityFactory>(
+            context, authorityName);
+    };
+    auto factory = getFactory();
     factory->d->setThis(factory);
     return factory;
 }
@@ -5601,6 +5609,7 @@ std::list<crs::GeodeticCRSNNPtr> AuthorityFactory::createGeodeticCRSFromDatum(
         sql += " AND type = ?";
         params.emplace_back(geodetic_crs_type);
     }
+    sql += " ORDER BY auth_name, code";
     auto sqlRes = d->run(sql, params);
     std::list<crs::GeodeticCRSNNPtr> res;
     for (const auto &row : sqlRes) {
diff --git a/src/iso19111/io.cpp b/src/iso19111/io.cpp
index 0d98e2de..e0e6152a 100644
--- a/src/iso19111/io.cpp
+++ b/src/iso19111/io.cpp
@@ -5669,10 +5669,7 @@ static BaseObjectNNPtr createFromUserInput(const std::string &text,
         DatabaseContextNNPtr dbContextNNPtr(NN_NO_CHECK(dbContext));
         const auto &authName = tokens[0];
         const auto &code = tokens[1];
-        static const std::string epsg_lowercase("epsg");
-        auto factory = AuthorityFactory::create(
-            dbContextNNPtr,
-            authName == epsg_lowercase ? Identifier::EPSG : authName);
+        auto factory = AuthorityFactory::create(dbContextNNPtr, authName);
         try {
             return factory->createCoordinateReferenceSystem(code);
         } catch (...) {
@@ -8582,7 +8579,11 @@ CRSNNPtr PROJStringParser::Private::buildProjectedCRS(
                     break;
                 }
             }
-            if (getNumericValue(getParamValue(step, "a")) == 6378137) {
+            if (getNumericValue(getParamValue(step, "a")) == 6378137 &&
+                getAngularValue(getParamValue(step, "lon_0")) == 0.0 &&
+                getAngularValue(getParamValue(step, "lat_0")) == 0.0 &&
+                getAngularValue(getParamValue(step, "x_0")) == 0.0 &&
+                getAngularValue(getParamValue(step, "y_0")) == 0.0) {
                 bWebMercator = true;
             }
         } else if (hasParamValue(step, "lat_ts")) {
diff --git a/src/lib_proj.cmake b/src/lib_proj.cmake
index c6a6e111..bc27bafe 100644
--- a/src/lib_proj.cmake
+++ b/src/lib_proj.cmake
@@ -283,6 +283,7 @@ set(SRC_LIBPROJ_CORE
   wkt_parser.cpp
   wkt_parser.hpp
   zpoly1.cpp
+  proj_json_streaming_writer.hpp
   proj_json_streaming_writer.cpp
   tracing.cpp
   ${CMAKE_CURRENT_BINARY_DIR}/proj_config.h
@@ -293,7 +294,6 @@ set(HEADERS_LIBPROJ
   proj.h
   proj_experimental.h
   proj_constants.h
-  proj_json_streaming_writer.hpp
   geodesic.h
 )
 
@@ -443,7 +443,7 @@ endif()
 include_directories(${SQLITE3_INCLUDE_DIR})
 target_link_libraries(${PROJ_CORE_TARGET} ${SQLITE3_LIBRARY})
 
-if(MSVC)
+if(MSVC AND BUILD_LIBPROJ_SHARED)
   target_compile_definitions(${PROJ_CORE_TARGET}
     PRIVATE PROJ_MSVC_DLL_EXPORT=1)
 endif()
diff --git a/src/open_lib.cpp b/src/open_lib.cpp
index 7fedb28b..24c31033 100644
--- a/src/open_lib.cpp
+++ b/src/open_lib.cpp
@@ -151,6 +151,51 @@ void pj_set_searchpath ( int count, const char **path )
     proj_context_set_search_paths( nullptr, count, const_cast<const char* const*>(path) );
 }
 
+#ifdef _WIN32
+#include <windows.h>
+#include <sys/stat.h>
+static const char *get_path_from_win32_projlib(const char *name, std::string& out) {
+    /* Check if proj.db lieves in a share/proj dir parallel to bin/proj.dll */
+    /* Based in https://stackoverflow.com/questions/9112893/how-to-get-path-to-executable-in-c-running-on-windows */
+
+    DWORD path_size = 1024;
+
+    for (;;) {
+        out.resize(path_size);
+        memset(&out[0], 0, path_size);
+        DWORD result     = GetModuleFileNameA(nullptr, &out[0], path_size - 1);
+        DWORD last_error = GetLastError();
+
+        if (result == 0) {
+            return nullptr;
+        }
+        else if (result == path_size - 1) {
+            if (ERROR_INSUFFICIENT_BUFFER != last_error) {
+                return nullptr;
+            }
+            path_size = path_size * 2;
+        }
+        else {
+            break;
+        }
+    }
+    // Now remove the program's name. It was (example) "C:\programs\gmt6\bin\gdal_translate.exe"
+    size_t k = strlen(out.c_str());
+    while (k > 0 && out[--k] != '\\') {}
+    out.resize(k);
+
+    out += "/../share/proj/";
+    out += name;
+
+    struct stat fileInfo; 
+    if (stat(out.c_str(), &fileInfo) == 0)	// Check if file exists (probably there are simpler ways)
+        return out.c_str();
+    else {
+        return nullptr;
+    }
+}
+#endif
+
 /************************************************************************/
 /*                          pj_open_lib_ex()                            */
 /************************************************************************/
@@ -229,6 +274,10 @@ pj_open_lib_ex(projCtx ctx, const char *name, const char *mode,
                 if( fid )
                     break;
             }
+#ifdef _WIN32
+        /* check if it lives in a ../share/proj dir of the proj dll */
+        } else if ((sysname = get_path_from_win32_projlib(name, fname)) != nullptr) {
+#endif
         /* or hardcoded path */
         } else if ((sysname = proj_lib_name) != nullptr) {
             fname = sysname;
diff --git a/src/tests/test228.cpp b/src/tests/test228.cpp
index b45fb286..8ae17c87 100644
--- a/src/tests/test228.cpp
+++ b/src/tests/test228.cpp
@@ -66,7 +66,7 @@ static void* thread_main(void* unused)
     p_WGS84_proj=pj_init_plus_ctx(p_proj_ctxt,"+proj=longlat "
             "+ellps=WGS84 +datum=WGS84");
     p_OSGB36_proj=pj_init_plus_ctx(p_proj_ctxt,
-            "+proj=longlat +ellps=airy +datum=OSGB36 +nadgrids=OSTN02_NTv2.gsb");
+            "+proj=longlat +ellps=airy +datum=OSGB36 +nadgrids=OSTN15_NTv2_OSGBtoETRS.gsb");
 
     while(run)
     {
@@ -81,8 +81,8 @@ static void* thread_main(void* unused)
         y *= RAD_TO_DEG;
         /*printf("%.18f %.18f\n", x, y); */
         assert(proj_ret == 0);
-        assert(fabs(x - -5.198965360936369962) < 1e-15);
-        assert(fabs(y - 49.999396034285531698) < 1e-15);
+        assert(fabs(x - -5.198965207267856492) < 1e-15);
+        assert(fabs(y - 49.999396074140378232) < 1e-15);
     }
 
     pj_free (p_OSGB36_proj);