Merge pull request #1231 from rouault/fix_1229

proj_create_crs_to_crs(): defer selection of actual coordinate operation until proj_trans() is called (fixes #1229)

10 files changed, 514 insertions, 74 deletions

diff --git a/docs/source/apps/cs2cs.rst b/docs/source/apps/cs2cs.rst
index 8df67d72..579dc65e 100644
--- a/docs/source/apps/cs2cs.rst
+++ b/docs/source/apps/cs2cs.rst
@@ -154,6 +154,10 @@ When using a WKT definition or a AUTHORITY:CODE, the axis order of the CRS will
 be enforced. So for example if using EPSG:4326, the first value expected (or
 returned) will be a latitude.
 
+Internally, :program:`cs2cs` uses the :c:func:`proj_create_crs_to_crs` function
+to compute the appropriate coordinate operation, so implementation details of
+this function directly impact the results returned by the program.
+
 The environment parameter :envvar:`PROJ_LIB` establishes the
 directory for resource files (database, datum shift grids, etc.)
 
@@ -195,7 +199,7 @@ The x-y output data will appear as three lines of:
 
 ::
 
-    1402285.98  5076292.42 -0.00
+    1402293.44  5076292.68 0.00
 
 
 Using EPSG codes
diff --git a/docs/source/development/reference/functions.rst b/docs/source/development/reference/functions.rst
index 87271117..4052ff82 100644
--- a/docs/source/development/reference/functions.rst
+++ b/docs/source/development/reference/functions.rst
@@ -108,11 +108,20 @@ paragraph for more details.
           When using that syntax, the axis order and unit for geographic CRS will
           be longitude, latitude, and the unit degrees.
 
+        - the name of a CRS as found in the PROJ database, e.g "WGS84", "NAD27", etc.
+
         - more generally any string accepted by :c:func:`proj_create`
 
     An "area of use" can be specified in area. When it is supplied, the more
     accurate transformation between two given systems can be chosen.
 
+    When no area of use is specific and several coordinate operations are possible
+    depending on the area of use, this function will internally store those
+    candidate coordinate operations in the return PJ object. Each subsequent
+    coordinate transformation done with :c:func:`proj_trans` will then select
+    the appropriate coordinate operation by comparing the input coordinates with
+    the area of use of the candidate coordinate operations.
+
     Example call:
 
     .. code-block:: C
diff --git a/include/proj/crs.hpp b/include/proj/crs.hpp
index 1c012666..6a09ea78 100644
--- a/include/proj/crs.hpp
+++ b/include/proj/crs.hpp
@@ -158,8 +158,7 @@ class PROJ_GCC_DLL SingleCRS : public CRS {
 
     PROJ_DLL const datum::DatumPtr &datum() PROJ_PURE_DECL;
     PROJ_DLL const datum::DatumEnsemblePtr &datumEnsemble() PROJ_PURE_DECL;
-    PROJ_DLL const cs::CoordinateSystemNNPtr &
-    coordinateSystem() PROJ_PURE_DECL;
+    PROJ_DLL const cs::CoordinateSystemNNPtr &coordinateSystem() PROJ_PURE_DECL;
 
     PROJ_PRIVATE :
         //! @cond Doxygen_Suppress
diff --git a/src/4D_api.cpp b/src/4D_api.cpp
index 642fbb1f..67460f3a 100644
--- a/src/4D_api.cpp
+++ b/src/4D_api.cpp
@@ -29,6 +29,7 @@
 
 #define FROM_PROJ_CPP
 
+#include <assert.h>
 #include <errno.h>
 #include <stddef.h>
 #include <stdio.h>
@@ -38,10 +39,13 @@
 #include <strings.h>
 #endif
 
+#include <algorithm>
+#include <limits>
+
 #include "proj.h"
+#include "proj_experimental.h"
 #include "proj_internal.h"
 #include "proj_math.h"
-#include "proj_internal.h"
 #include "geodesic.h"
 
 #include "proj/common.hpp"
@@ -182,18 +186,52 @@ available.
 See also pj_approx_2D_trans and pj_approx_3D_trans in pj_internal.c, which work
 similarly, but prefers the 2D resp. 3D interfaces if available.
 ***************************************************************************************/
-    if (nullptr==P)
+    if (nullptr==P || direction == PJ_IDENT)
         return coord;
     if (P->inverted)
         direction = opposite_direction(direction);
 
+    if( !P->alternativeCoordinateOperations.empty() ) {
+        int i = 0;
+        for( const auto &alt: P->alternativeCoordinateOperations ) {
+            if( direction == PJ_FWD ) {
+                if( coord.xyzt.x >= alt.minxSrc &&
+                    coord.xyzt.y >= alt.minySrc &&
+                    coord.xyzt.x <= alt.maxxSrc &&
+                    coord.xyzt.y <= alt.maxySrc ) {
+                    if( P->iCurCoordOp != i ) {
+                        std::string msg("Using coordinate operation ");
+                        msg += alt.name;
+                        pj_log(P->ctx, PJ_LOG_TRACE, msg.c_str());
+                        P->iCurCoordOp = i;
+                    }
+                    return pj_fwd4d( coord, alt.pj );
+                }
+            } else {
+                if( coord.xyzt.x >= alt.minxDst &&
+                    coord.xyzt.y >= alt.minyDst &&
+                    coord.xyzt.x <= alt.maxxDst &&
+                    coord.xyzt.y <= alt.maxyDst ) {
+                    if( P->iCurCoordOp != i ) {
+                        std::string msg("Using coordinate operation ");
+                        msg += alt.name;
+                        pj_log(P->ctx, PJ_LOG_TRACE, msg.c_str());
+                        P->iCurCoordOp = i;
+                    }
+                    return pj_inv4d( coord, alt.pj );
+                }
+            }
+            i ++;
+        }
+        proj_errno_set (P, EINVAL);
+        return proj_coord_error ();
+    }
+
     switch (direction) {
         case PJ_FWD:
             return pj_fwd4d (coord, P);
         case PJ_INV:
             return  pj_inv4d (coord, P);
-        case PJ_IDENT:
-            return coord;
         default:
             break;
     }
@@ -780,118 +818,389 @@ std::string pj_add_type_crs_if_needed(const std::string& str)
 }
 
 /*****************************************************************************/
+static PJ* op_to_pj(PJ_CONTEXT* ctx, PJ* op)
+/*****************************************************************************/
+{
+    auto proj_string = proj_as_proj_string(ctx, op, PJ_PROJ_5, nullptr);
+    if( !proj_string) {
+        return nullptr;
+    }
+
+    if( proj_string[0] == '\0' ) {
+        /* Null transform ? */
+        return proj_create(ctx, "proj=affine");
+    } else {
+        return proj_create(ctx, proj_string);
+    }
+}
+
+/*****************************************************************************/
+static void reproject_bbox(PJ* pjGeogToCrs,
+                          double west_lon, double south_lat,
+                          double east_lon, double north_lat,
+                          double& minx,
+                          double& miny,
+                          double& maxx,
+                          double& maxy) {
+/*****************************************************************************/
+
+    minx = -std::numeric_limits<double>::max();
+    miny = -std::numeric_limits<double>::max();
+    maxx = std::numeric_limits<double>::max();
+    maxy = std::numeric_limits<double>::max();
+
+    if( !(west_lon == -180.0 && east_lon == 180.0 &&
+          south_lat == -90.0 && north_lat == 90.0) )
+    {
+        minx = -minx;
+        miny = -miny;
+        maxx = -maxx;
+        maxy = -maxy;
+
+        double x[21 * 4], y[21 * 4];
+        for( int j = 0; j <= 20; j++ )
+        {
+            x[j] = west_lon + j * (east_lon - west_lon) / 20;
+            y[j] = south_lat;
+            x[21+j] = west_lon + j * (east_lon - west_lon) / 20;
+            y[21+j] = north_lat;
+            x[21*2+j] = west_lon;
+            y[21*2+j] = south_lat + j * (north_lat - south_lat) / 20;
+            x[21*3+j] = east_lon;
+            y[21*3+j] = south_lat + j * (north_lat - south_lat) / 20;
+        }
+        proj_trans_generic (
+            pjGeogToCrs, PJ_FWD,
+                x, sizeof(double), 21 * 4,
+                y, sizeof(double), 21 * 4,
+                nullptr, 0, 0,
+                nullptr, 0, 0);
+        for( int j = 0; j < 21 * 4; j++ )
+        {
+            if( x[j] != HUGE_VAL && y[j] != HUGE_VAL )
+            {
+                minx = std::min(minx, x[j]);
+                miny = std::min(miny, y[j]);
+                maxx = std::max(maxx, x[j]);
+                maxy = std::max(maxy, y[j]);
+            }
+        }
+    }
+}
+
+
+/*****************************************************************************/
+static PJ* add_coord_op_to_list(PJ_CONTEXT* ctx, PJ* op,
+                            double west_lon, double south_lat,
+                            double east_lon, double north_lat,
+                            PJ* pjGeogToSrc,
+                            PJ* pjGeogToDst,
+                            std::vector<PJconsts::CoordOperation>& altCoordOps) {
+/*****************************************************************************/
+
+    double minxSrc;
+    double minySrc;
+    double maxxSrc;
+    double maxySrc;
+    double minxDst;
+    double minyDst;
+    double maxxDst;
+    double maxyDst;
+
+    reproject_bbox(pjGeogToSrc, west_lon, south_lat, east_lon, north_lat,
+                  minxSrc, minySrc, maxxSrc, maxySrc);
+    reproject_bbox(pjGeogToDst, west_lon, south_lat, east_lon, north_lat,
+                  minxDst, minyDst, maxxDst, maxyDst);
+
+    if( minxSrc <= maxxSrc && minxDst <= maxxDst )
+    {
+        const char* c_name = proj_get_name(op);
+        std::string name(c_name ? c_name : "");
+        auto pj = op_to_pj(ctx, op);
+        proj_destroy(op);
+        op = nullptr;
+        if( pj )
+        {
+            altCoordOps.emplace_back(minxSrc, minySrc, maxxSrc, maxySrc,
+                                     minxDst, minyDst, maxxDst, maxyDst,
+                                     pj, name);
+        }
+    }
+    return op;
+};
+
+/*****************************************************************************/
+static PJ* create_operation_to_base_geog_crs(PJ_CONTEXT* ctx, PJ* crs) {
+/*****************************************************************************/
+    // Create a geographic 2D long-lat degrees CRS that is related to the
+    // CRS
+    auto geodetic_crs = proj_crs_get_geodetic_crs(ctx, crs);
+    if( !geodetic_crs ) {
+        proj_context_log_debug(ctx, "Cannot find geodetic CRS matching CRS");
+        return nullptr;
+    }
+
+    auto geodetic_crs_type = proj_get_type(geodetic_crs);
+    if( geodetic_crs_type == PJ_TYPE_GEOCENTRIC_CRS ||
+        geodetic_crs_type == PJ_TYPE_GEOGRAPHIC_2D_CRS ||
+        geodetic_crs_type == PJ_TYPE_GEOGRAPHIC_3D_CRS )
+    {
+        auto datum = proj_crs_get_datum(ctx, geodetic_crs);
+        if( datum )
+        {
+            auto cs = proj_create_ellipsoidal_2D_cs(
+                ctx, PJ_ELLPS2D_LONGITUDE_LATITUDE, nullptr, 0);
+            auto temp = proj_create_geographic_crs_from_datum(
+                ctx,"unnamed", datum, cs);
+            proj_destroy(datum);
+            proj_destroy(cs);
+            proj_destroy(geodetic_crs);
+            geodetic_crs = temp;
+            geodetic_crs_type = proj_get_type(geodetic_crs);
+        }
+    }
+    if( geodetic_crs_type != PJ_TYPE_GEOGRAPHIC_2D_CRS )
+    {
+        // Shouldn't happen
+        proj_context_log_debug(ctx, "Cannot find geographic CRS matching CRS");
+        proj_destroy(geodetic_crs);
+        return nullptr;
+    }
+
+    // Create the transformation from this geographic 2D CRS to the source CRS
+    auto operation_ctx = proj_create_operation_factory_context(ctx, nullptr);
+    proj_operation_factory_context_set_spatial_criterion(
+        ctx, operation_ctx, PROJ_SPATIAL_CRITERION_PARTIAL_INTERSECTION);
+    proj_operation_factory_context_set_grid_availability_use(
+        ctx, operation_ctx, PROJ_GRID_AVAILABILITY_DISCARD_OPERATION_IF_MISSING_GRID);
+    auto op_list_to_geodetic = proj_create_operations(
+        ctx, geodetic_crs, crs, operation_ctx);
+    proj_operation_factory_context_destroy(operation_ctx);
+    proj_destroy(geodetic_crs);
+
+    if( op_list_to_geodetic == nullptr ||
+        proj_list_get_count(op_list_to_geodetic) == 0 )
+    {
+        proj_context_log_debug(ctx, "Cannot compute transformation from geographic CRS to CRS");
+        proj_list_destroy(op_list_to_geodetic);
+        return nullptr;
+    }
+    auto opGeogToCrs = proj_list_get(ctx, op_list_to_geodetic, 0);
+    assert(opGeogToCrs);
+    proj_list_destroy(op_list_to_geodetic);
+    auto P = op_to_pj(ctx, opGeogToCrs);
+    proj_destroy(opGeogToCrs);
+    return P;
+}
+
+/*****************************************************************************/
 PJ  *proj_create_crs_to_crs (PJ_CONTEXT *ctx, const char *source_crs, const char *target_crs, PJ_AREA *area) {
 /******************************************************************************
     Create a transformation pipeline between two known coordinate reference
     systems.
 
-    source_crs and target_crs can be :
-    - a "AUTHORITY:CODE", like EPSG:25832. When using that syntax for a source
-      CRS, the created pipeline will expect that the values passed to proj_trans()
-      respect the axis order and axis unit of the official definition (
-      so for example, for EPSG:4326, with latitude first and longitude next,
-      in degrees). Similarly, when using that syntax for a target CRS, output
-      values will be emitted according to the official definition of this CRS.
-    - a PROJ string, like "+proj=longlat +datum=WGS84".
-      When using that syntax, the axis order and unit for geographic CRS will
-      be longitude, latitude, and the unit degrees.
-    - more generally any string accepted by proj_create()
-
-    An "area of use" can be specified in area. When it is supplied, the more
-    accurate transformation between two given systems can be chosen.
-
-    Example call:
-
-        PJ *P = proj_create_crs_to_crs(0, "EPSG:25832", "EPSG:25833", NULL);
+    See docs/source/development/reference/functions.rst
 
 ******************************************************************************/
-    const char* proj_string;
-
     if( !ctx ) {
         ctx = pj_get_default_ctx();
     }
 
+    PJ* src;
+    PJ* dst;
     try
     {
         std::string source_crs_modified(pj_add_type_crs_if_needed(source_crs));
         std::string target_crs_modified(pj_add_type_crs_if_needed(target_crs));
 
-        auto src = proj_create(ctx, source_crs_modified.c_str());
+        src = proj_create(ctx, source_crs_modified.c_str());
         if( !src ) {
             proj_context_log_debug(ctx, "Cannot instantiate source_crs");
             return nullptr;
         }
 
-        auto dst = proj_create(ctx, target_crs_modified.c_str());
+        dst = proj_create(ctx, target_crs_modified.c_str());
         if( !dst ) {
             proj_context_log_debug(ctx, "Cannot instantiate target_crs");
             proj_destroy(src);
             return nullptr;
         }
+    }
+    catch( const std::exception& )
+    {
+        return nullptr;
+    }
 
-        auto operation_ctx = proj_create_operation_factory_context(ctx, nullptr);
-        if( !operation_ctx ) {
-            proj_destroy(src);
-            proj_destroy(dst);
-            return nullptr;
-        }
+    auto operation_ctx = proj_create_operation_factory_context(ctx, nullptr);
+    if( !operation_ctx ) {
+        proj_destroy(src);
+        proj_destroy(dst);
+        return nullptr;
+    }
 
-        if( area && area->bbox_set ) {
-            proj_operation_factory_context_set_area_of_interest(
-                                                ctx,
-                                                operation_ctx,
-                                                area->west_lon_degree,
-                                                area->south_lat_degree,
-                                                area->east_lon_degree,
-                                                area->north_lat_degree);
-        }
+    if( area && area->bbox_set ) {
+        proj_operation_factory_context_set_area_of_interest(
+                                            ctx,
+                                            operation_ctx,
+                                            area->west_lon_degree,
+                                            area->south_lat_degree,
+                                            area->east_lon_degree,
+                                            area->north_lat_degree);
+    }
 
-        proj_operation_factory_context_set_grid_availability_use(
-            ctx, operation_ctx, PROJ_GRID_AVAILABILITY_DISCARD_OPERATION_IF_MISSING_GRID);
+    proj_operation_factory_context_set_spatial_criterion(
+        ctx, operation_ctx, PROJ_SPATIAL_CRITERION_PARTIAL_INTERSECTION);
+    proj_operation_factory_context_set_grid_availability_use(
+        ctx, operation_ctx, PROJ_GRID_AVAILABILITY_DISCARD_OPERATION_IF_MISSING_GRID);
 
-        auto op_list = proj_create_operations(ctx, src, dst, operation_ctx);
+    auto op_list = proj_create_operations(ctx, src, dst, operation_ctx);
 
+    if( !op_list ) {
         proj_operation_factory_context_destroy(operation_ctx);
         proj_destroy(src);
         proj_destroy(dst);
+        return nullptr;
+    }
 
-        if( !op_list ) {
-            return nullptr;
-        }
-
-        if( proj_list_get_count(op_list) == 0 ) {
-            proj_list_destroy(op_list);
-            proj_context_log_debug(ctx, "No operation found matching criteria");
-            return nullptr;
-        }
+    auto op_count = proj_list_get_count(op_list);
+    if( op_count == 0 ) {
+        proj_list_destroy(op_list);
+        proj_operation_factory_context_destroy(operation_ctx);
+        proj_destroy(src);
+        proj_destroy(dst);
+        proj_context_log_debug(ctx, "No operation found matching criteria");
+        return nullptr;
+    }
 
+    PJ* P;
+    {
         auto op = proj_list_get(ctx, op_list, 0);
+        assert(op);
+        P = op_to_pj(ctx, op);
+        proj_destroy(op);
+    }
+
+    if( P == nullptr || op_count == 1 || (area && area->bbox_set) ||
+        proj_get_type(src) == PJ_TYPE_GEOCENTRIC_CRS ||
+        proj_get_type(dst) == PJ_TYPE_GEOCENTRIC_CRS ) {
         proj_list_destroy(op_list);
-        if( !op ) {
-            return nullptr;
-        }
+        proj_operation_factory_context_destroy(operation_ctx);
+        proj_destroy(src);
+        proj_destroy(dst);
+        return P;
+    }
+
+    auto pjGeogToSrc = create_operation_to_base_geog_crs(ctx, src);
+    if( !pjGeogToSrc )
+    {
+        proj_list_destroy(op_list);
+        proj_operation_factory_context_destroy(operation_ctx);
+        proj_destroy(src);
+        proj_destroy(dst);
+        proj_context_log_debug(ctx,
+            "Cannot create transformation from geographic CRS of source CRS to source CRS");
+        proj_destroy(P);
+        return nullptr;
+    }
+
+    auto pjGeogToDst = create_operation_to_base_geog_crs(ctx, dst);
+    if( !pjGeogToDst )
+    {
+        proj_list_destroy(op_list);
+        proj_operation_factory_context_destroy(operation_ctx);
+        proj_destroy(src);
+        proj_destroy(dst);
+        proj_context_log_debug(ctx,
+            "Cannot create transformation from geographic CRS of target CRS to target CRS");
+        proj_destroy(P);
+        proj_destroy(pjGeogToSrc);
+        return nullptr;
+    }
+
+    try
+    {
+        // Iterate over source->target candidate transformations and reproject
+        // their long-lat bounding box into the source CRS.
+        for( int i = 0; i < op_count; i++ )
+        {
+            auto op = proj_list_get(ctx, op_list, i);
+            assert(op);
+            double west_lon = 0.0;
+            double south_lat = 0.0;
+            double east_lon = 0.0;
+            double north_lat = 0.0;
+
+            const char* name = proj_get_name(op);
+            if( name && (strstr(name, "Null geographic offset") ||
+                         strstr(name, "Null geocentric translation")) )
+            {
+                // Skip default transformations
+            }
+            else if( proj_get_area_of_use(ctx, op,
+                        &west_lon, &south_lat, &east_lon, &north_lat, nullptr) )
+            {
+                if( west_lon <= east_lon )
+                {
+                    op = add_coord_op_to_list(ctx, op,
+                        west_lon, south_lat, east_lon, north_lat,
+                        pjGeogToSrc, pjGeogToDst,
+                        P->alternativeCoordinateOperations);
+                }
+                else
+                {
+                    auto op_clone = proj_clone(ctx, op);
+
+                    op = add_coord_op_to_list(ctx, op,
+                        west_lon, south_lat, 180, north_lat,
+                        pjGeogToSrc, pjGeogToDst,
+                        P->alternativeCoordinateOperations);
+                    op_clone = add_coord_op_to_list(ctx, op_clone,
+                        -180, south_lat, east_lon, north_lat,
+                        pjGeogToSrc, pjGeogToDst,
+                        P->alternativeCoordinateOperations);
+                    proj_destroy(op_clone);
+                }
+            }
 
-        proj_string = proj_as_proj_string(ctx, op, PJ_PROJ_5, nullptr);
-        if( !proj_string) {
             proj_destroy(op);
-            proj_context_log_debug(ctx, "Cannot export operation as a PROJ string");
-            return nullptr;
         }
 
-        PJ* P;
-        if( proj_string[0] == '\0' ) {
-            /* Null transform ? */
-            P = proj_create(ctx, "proj=affine");
+        proj_list_destroy(op_list);
+
+        proj_operation_factory_context_destroy(operation_ctx);
+        proj_destroy(src);
+        proj_destroy(dst);
+        proj_destroy(pjGeogToSrc);
+        proj_destroy(pjGeogToDst);
+
+        // If there's finally juste a single result, return it directly
+        if( P->alternativeCoordinateOperations.size() == 1 ) {
+            auto retP = P->alternativeCoordinateOperations[0].pj;
+            P->alternativeCoordinateOperations[0].pj = nullptr;
+            proj_destroy(P);
+            P = retP;
         } else {
-            P = proj_create(ctx, proj_string);
+            // The returned P is rather dummy
+            P->iso_obj = nullptr;
+            P->fwd = nullptr;
+            P->inv = nullptr;
+            P->fwd3d = nullptr;
+            P->inv3d = nullptr;
+            P->fwd4d = nullptr;
+            P->inv4d = nullptr;
         }
 
-        proj_destroy(op);
-
         return P;
     }
     catch( const std::exception& )
     {
+        proj_list_destroy(op_list);
+        proj_operation_factory_context_destroy(operation_ctx);
+        proj_destroy(src);
+        proj_destroy(dst);
+        proj_destroy(pjGeogToSrc);
+        proj_destroy(pjGeogToDst);
+        proj_destroy(P);
         return nullptr;
     }
 }
@@ -1132,11 +1441,20 @@ PJ_PROJ_INFO proj_pj_info(PJ *P) {
     if (nullptr==P)
         return pjinfo;
 
+    /* coordinate operation description */
+    if( P->iCurCoordOp >= 0 ) {
+        P = P->alternativeCoordinateOperations[P->iCurCoordOp].pj;
+    } else if( !P->alternativeCoordinateOperations.empty() ) {
+        pjinfo.id = "unknown";
+        pjinfo.description = "unavailable until proj_trans is called";
+        pjinfo.definition = "unavailable until proj_trans is called";
+        return pjinfo;
+    }
+
     /* projection id */
     if (pj_param(P->ctx, P->params, "tproj").i)
         pjinfo.id = pj_param(P->ctx, P->params, "sproj").s;
 
-    /* coordinate operation description */
     if( P->iso_obj ) {
         pjinfo.description = P->iso_obj->nameStr().c_str();
     } else {
diff --git a/src/iso19111/c_api.cpp b/src/iso19111/c_api.cpp
index 06a3c02e..de11f181 100644
--- a/src/iso19111/c_api.cpp
+++ b/src/iso19111/c_api.cpp
@@ -337,7 +337,7 @@ PJ *proj_create(PJ_CONTEXT *ctx, const char *text) {
     assert(text);
 
     // Only connect to proj.db if needed
-    if( strstr(text, "proj=") == nullptr || strstr(text, "init=") != nullptr ) {
+    if (strstr(text, "proj=") == nullptr || strstr(text, "init=") != nullptr) {
         getDBcontextNoException(ctx, __FUNCTION__);
     }
     try {
diff --git a/src/iso19111/crs.cpp b/src/iso19111/crs.cpp
index 74b94a1f..1f03e262 100644
--- a/src/iso19111/crs.cpp
+++ b/src/iso19111/crs.cpp
@@ -1653,8 +1653,7 @@ CRSNNPtr GeographicCRS::_shallowClone() const {
  *
  * @return a EllipsoidalCS.
  */
-const cs::EllipsoidalCSNNPtr &
-GeographicCRS::coordinateSystem() PROJ_PURE_DEFN {
+const cs::EllipsoidalCSNNPtr &GeographicCRS::coordinateSystem() PROJ_PURE_DEFN {
     return d->coordinateSystem_;
 }
 
diff --git a/src/proj_internal.h b/src/proj_internal.h
index b97afdec..f5196939 100644
--- a/src/proj_internal.h
+++ b/src/proj_internal.h
@@ -60,6 +60,7 @@
 
 #include "proj/common.hpp"
 #include "proj/coordinateoperation.hpp"
+
 #include <string>
 #include <vector>
 
@@ -519,6 +520,50 @@ struct PJconsts {
     mutable std::vector<NS_PROJ::operation::GridDescription> gridsNeeded{};
 
     /*************************************************************************************
+     proj_create_crs_to_crs() alternative coordinate operations
+    **************************************************************************************/
+
+    struct CoordOperation
+    {
+        double minxSrc = 0.0;
+        double minySrc = 0.0;
+        double maxxSrc = 0.0;
+        double maxySrc = 0.0;
+        double minxDst = 0.0;
+        double minyDst = 0.0;
+        double maxxDst = 0.0;
+        double maxyDst = 0.0;
+        PJ* pj = nullptr;
+        std::string name{};
+
+        CoordOperation(double minxSrcIn, double minySrcIn, double maxxSrcIn, double maxySrcIn,
+                       double minxDstIn, double minyDstIn, double maxxDstIn, double maxyDstIn,
+                       PJ* pjIn, const std::string& nameIn):
+            minxSrc(minxSrcIn), minySrc(minySrcIn), maxxSrc(maxxSrcIn), maxySrc(maxySrcIn),
+            minxDst(minxDstIn), minyDst(minyDstIn), maxxDst(maxxDstIn), maxyDst(maxyDstIn),
+            pj(pjIn), name(nameIn) {}
+
+        CoordOperation(const CoordOperation&) = delete;
+
+        CoordOperation(CoordOperation&& other):
+            minxSrc(other.minxSrc), minySrc(other.minySrc), maxxSrc(other.maxxSrc), maxySrc(other.maxySrc),
+            minxDst(other.minxDst), minyDst(other.minyDst), maxxDst(other.maxxDst), maxyDst(other.maxyDst),
+            name(std::move(other.name)) {
+                pj = other.pj;
+                other.pj = nullptr;
+            }
+
+        CoordOperation& operator=(const CoordOperation&) = delete;
+
+        ~CoordOperation()
+        {
+            proj_destroy(pj);
+        }
+    };
+    std::vector<CoordOperation> alternativeCoordinateOperations{};
+    int iCurCoordOp = -1;
+
+    /*************************************************************************************
 
                  E N D   O F    G E N E R A L   P A R A M E T E R   S T R U C T
 
diff --git a/test/cli/td_out.dist b/test/cli/td_out.dist
index f6b2a219..cdbe0fa9 100644
--- a/test/cli/td_out.dist
+++ b/test/cli/td_out.dist
@@ -23,3 +23,7 @@ edge or even a wee bit outside (#141).
 -5.5000000000001 52.0000000000001	-5.498893534472	52.000109529717 0.000000000000
 -5.4999 51.9999	-5.498793541695	52.000009529743 0.000000000000
 -5.5001 52.0	-5.500100000000	52.000000000000 0.000000000000
+##############################################################
+NAD27 -> NAD83: 1st through ntv1, 2nd through conus
+44d00'00.000"N 111d00'00.000"W 0.0	43d59'59.732"N	111d0'3.208"W 0.000
+39d00'00.000"N 111d00'00.000"W 0.0	38d59'59.912"N	111d0'2.604"W 0.000
diff --git a/test/cli/testdatumfile b/test/cli/testdatumfile
index d048d8e6..e8995150 100755
--- a/test/cli/testdatumfile
+++ b/test/cli/testdatumfile
@@ -97,7 +97,16 @@ $EXE +proj=latlong +datum=WGS84 \
 -5.4999 51.9999
 -5.5001 52.0
 EOF
+#
+echo "##############################################################" >> ${OUT}
+echo "NAD27 -> NAD83: 1st through ntv1, 2nd through conus" >> ${OUT}
+#
+$EXE NAD27 NAD83 -E >>${OUT} <<EOF
+44d00'00.000"N 111d00'00.000"W 0.0
+39d00'00.000"N 111d00'00.000"W 0.0
+EOF
 
+# Cleanup
 rm -rf "dir with \" space"
 
 #
diff --git a/test/unit/gie_self_tests.cpp b/test/unit/gie_self_tests.cpp
index 1b5ce83b..2d298ea3 100644
--- a/test/unit/gie_self_tests.cpp
+++ b/test/unit/gie_self_tests.cpp
@@ -659,4 +659,57 @@ TEST(gie, horner_selftest) {
     proj_destroy(P);
 }
 
+// ---------------------------------------------------------------------------
+
+TEST(gie, proj_create_crs_to_crs_PULKOVO42_ETRS89) {
+    auto P = proj_create_crs_to_crs(PJ_DEFAULT_CTX, "EPSG:4179", "EPSG:4258",
+                                    nullptr);
+    ASSERT_TRUE(P != nullptr);
+    PJ_COORD c;
+
+    EXPECT_EQ(std::string(proj_pj_info(P).definition),
+              "unavailable until proj_trans is called");
+    EXPECT_EQ(proj_get_name(P), nullptr);
+    EXPECT_EQ(P->fwd, nullptr);
+    EXPECT_EQ(P->fwd3d, nullptr);
+    EXPECT_EQ(P->fwd4d, nullptr);
+
+    // Romania
+    c.xyz.x = 45; // Lat
+    c.xyz.y = 25; // Long
+    c.xyz.z = 0;
+    c = proj_trans(P, PJ_FWD, c);
+    EXPECT_NEAR(c.xy.x, 44.999701238, 1e-9);
+    EXPECT_NEAR(c.xy.y, 24.998474948, 1e-9);
+    EXPECT_EQ(std::string(proj_pj_info(P).definition),
+              "proj=pipeline step proj=axisswap order=2,1 step "
+              "proj=unitconvert xy_in=deg xy_out=rad step proj=cart "
+              "ellps=krass step proj=helmert x=2.3287 y=-147.0425 z=-92.0802 "
+              "rx=0.3092483 ry=-0.32482185 rz=-0.49729934 s=5.68906266 "
+              "convention=coordinate_frame step inv proj=cart ellps=GRS80 step "
+              "proj=unitconvert xy_in=rad xy_out=deg step proj=axisswap "
+              "order=2,1");
+
+    c = proj_trans(P, PJ_INV, c);
+    EXPECT_NEAR(c.xy.x, 45, 1e-9);
+    EXPECT_NEAR(c.xy.y, 25, 1e-9);
+
+    // Poland
+    c.xyz.x = 52; // Lat
+    c.xyz.y = 20; // Long
+    c.xyz.z = 0;
+    c = proj_trans(P, PJ_FWD, c);
+    EXPECT_NEAR(c.xy.x, 51.999714150, 1e-9);
+    EXPECT_NEAR(c.xy.y, 19.998187811, 1e-9);
+    EXPECT_EQ(std::string(proj_pj_info(P).definition),
+              "proj=pipeline step proj=axisswap order=2,1 step "
+              "proj=unitconvert xy_in=deg xy_out=rad step proj=cart "
+              "ellps=krass step proj=helmert x=33.4 y=-146.6 z=-76.3 rx=-0.359 "
+              "ry=-0.053 rz=0.844 s=-0.84 convention=position_vector step inv "
+              "proj=cart ellps=GRS80 step proj=unitconvert xy_in=rad "
+              "xy_out=deg step proj=axisswap order=2,1");
+
+    proj_destroy(P);
+}
+
 } // namespace