Merge pull request #3535 from rouault/fix_3533

proj_create_crs_to_crs_from_pj(): add a ONLY_BEST=YES option (fixes #3533)

data/proj.ini

@@ -15,6 +15,16 @@ cache_size_MB = 300
 
 cache_ttl_sec = 86400
 
+; Can be set to on so that by default the lack of a known resource files needed
+; for the best transformation PROJ would normally use causes an error, or off
+; to accept missing resource files without errors or warnings.
+; This default value itself is overriden by the PROJ_ONLY_BEST_DEFAULT environment
+; variable if set, and then by the ONLY_BEST setting that can be
+; passed to the proj_create_crs_to_crs() method, or with the --only-best
+; option of the cs2cs program.
+; (added in PROJ 9.2)
+; only_best_default = on
+
 ; Filename of the Certificate Authority (CA) bundle.
 ; Can be overriden with the PROJ_CURL_CA_BUNDLE / CURL_CA_BUNDLE environment variable.
 ; (added in PROJ 9.0)

docs/source/apps/cs2cs.rst

@@ -13,7 +13,8 @@ Synopsis
 
     | **cs2cs** [**-eEfIlrstvwW** [args]]
     |           [[--area <name_or_code>] | [--bbox <west_long,south_lat,east_long,north_lat>]]
-    |           [--authority <name>] [--no-ballpark] [--accuracy <accuracy>] [--3d]
+    |           [--authority <name>] [--3d]
+    |           [--accuracy <accuracy>] [--only-best[=yes|=no]] [--no-ballpark]
     |           ([*+opt[=arg]* ...] [+to *+opt[=arg]* ...] | {source_crs} {target_crs})
     |           file ...
 
@@ -166,6 +167,23 @@ The following control parameters can appear in any order:
     `south_lat` and `north_lat` in the [-90,90]. `west_long` is generally lower than
     `east_long`, except in the case where the area of interest crosses the antimeridian.
 
+.. option:: --only-best[=yes|=no]
+
+    .. versionadded:: 9.2.0
+
+    Force `cs2cs` to only use the best transformation known by PROJ.
+    `cs2cs` will return an error if a grid needed for the best transformation is missing.
+
+    Best transformation should be understood as the most accurate transformation
+    available among all relevant for the point to transform, and if all known
+    grids required to perform such transformation were accessible (either locally
+    or through network).
+
+    Note that the default value for this option can be also set with the
+    :envvar:`PROJ_ONLY_BEST_DEFAULT` environment variable, or with the
+    ``only_best_default`` setting of :ref:`proj-ini` (:option:`--only-best`
+    when specified overrides such default value).
+
 .. option:: --no-ballpark
 
     .. versionadded:: 8.0.0
@@ -284,12 +302,42 @@ The x-y output data will appear as three lines of:
 
 ::
 
-    1402285.98  5076292.42 0.00
+    1402285.93  5076292.58 0.00
 
-.. note::
 
-    To get those exact values, you have need to have all current grids installed
-    locally or use networking capabilities mentioned above.
+To get those exact values, you have need to have all current grids installed
+(in that instance the NADCON5 :file:`us_noaa_nadcon5_nad27_nad83_1986_conus.tif` grid)
+locally or use networking capabilities mentioned above.
+
+To make sure you will get the optimal result, you may add :option:`--only-best`.
+Assuming the above mentionned grid is *not* available,
+
+::
+
+    echo -111.5 45.25919444444 | cs2cs --only-best +proj=latlong +datum=NAD83 +to +proj=utm +zone=10 +datum=NAD27
+
+would return:
+
+::
+
+    Attempt to use coordinate operation axis order change (2D) + Inverse of NAD27 to NAD83 (7) + axis order change (2D) + UTM zone 10N failed. Grid us_noaa_nadcon5_nad27_nad83_1986_conus.tif is not available. Consult https://proj.org/resource_files.html for guidance.
+    *   * inf
+
+Otherwise, if you don't have the grid available and you don't specify :option:`--only-best`:
+
+::
+
+    echo -111.5 45.25919444444 | cs2cs --only-best +proj=latlong +datum=NAD83 +to +proj=utm +zone=10 +datum=NAD27
+
+would return:
+
+::
+
+    1402224.57  5076275.42 0.00
+
+which is the result when the NAD27 and NAD83 datums are dealt as identical,
+which is an approximation at a level of several tens of metres.
+
 
 Using EPSG CRS codes
 --------------------

docs/source/development/reference/functions.rst

@@ -201,6 +201,18 @@ paragraph for more details.
     - ALLOW_BALLPARK=YES/NO: can be set to NO to disallow the use of
       :term:`Ballpark transformation` in the candidate coordinate operations.
 
+    - ONLY_BEST=YES/NO: (PROJ >= 9.2)
+      Can be set to YES to cause PROJ to error out if the best
+      transformation, known of PROJ, and usable by PROJ if all grids known and
+      usable by PROJ were accessible, cannot be used. Best transformation should
+      be understood as the transformation returned by
+      :cpp:func:`proj_get_suggested_operation` if all known grids were
+      accessible (either locally or through network).
+      Note that the default value for this option can be also set with the
+      :envvar:`PROJ_ONLY_BEST_DEFAULT` environment variable, or with the
+      ``only_best_default`` setting of :ref:`proj-ini` (the ONLY_BEST option
+      when specified overrides such default value).
+
     - FORCE_OVER=YES/NO: can be set to YES to force the ``+over`` flag on the transformation
       returned by this function. See :ref:`longitude_wrapping`
 

docs/source/resource_files.rst

@@ -117,37 +117,7 @@ network related parameters.
 
 Its default content is:
 
-::
-
-    [general]
-    ; Lines starting by ; are commented lines.
-    ;
-
-    ; Network capabilities disabled by default.
-    ; Can be overridden with the PROJ_NETWORK=ON environment variable.
-    ; network = on
-
-    ; Can be overridden with the PROJ_NETWORK_ENDPOINT environment variable.
-    cdn_endpoint = https://cdn.proj.org
-
-    cache_enabled = on
-
-    cache_size_MB = 300
-
-    cache_ttl_sec = 86400
-
-    ; Filename of the Certificate Authority (CA) bundle.
-    ; Can be overriden with the PROJ_CURL_CA_BUNDLE / CURL_CA_BUNDLE environment variable.
-    ; (added in PROJ 9.0)
-    ; ca_bundle_path = /path/to/cabundle.pem
-
-    ; Transverse Mercator (and UTM) default algorithm: auto, evenden_snyder or poder_engsager
-    ; * evenden_snyder is the fastest, but less accurate far from central meridian
-    ; * poder_engsager is slower, but more accurate far from central meridian
-    ; * default will auto-select between the two above depending on the coordinate
-    ;   to transform and will use evenden_snyder if the error in doing so is below
-    ;   0.1 mm (for an ellipsoid of the size of Earth)
-    tmerc_default_algo = poder_engsager
+.. literalinclude:: ../../data/proj.ini
 
 
 Transformation grids

src/4D_api.cpp

@@ -273,6 +273,42 @@ int pj_get_suggested_operation(PJ_CONTEXT*,
     return iBest;
 }
 
+/**************************************************************************************/
+static void warnAboutMissingGrid(PJ* P)
+/**************************************************************************************/
+{
+    std::string msg("Attempt to use coordinate operation ");
+    msg += proj_get_name(P);
+    msg += " failed.";
+    int gridUsed = proj_coordoperation_get_grid_used_count(P->ctx, P);
+    for( int i = 0; i < gridUsed; ++i )
+    {
+        const char* gridName = "";
+        int available = FALSE;
+        if( proj_coordoperation_get_grid_used(
+                P->ctx, P, i, &gridName, nullptr, nullptr,
+                nullptr, nullptr, nullptr, &available) &&
+            !available )
+        {
+            msg += " Grid ";
+            msg += gridName;
+            msg += " is not available. "
+                   "Consult https://proj.org/resource_files.html for guidance.";
+        }
+    }
+    if( !P->errorIfBestTransformationNotAvailable &&
+        P->warnIfBestTransformationNotAvailable )
+    {
+        msg += " This might become an error in a future PROJ major release. "
+               "Set the ONLY_BEST option to YES or NO. "
+               "This warning will no longer be emitted (for the current transformation instance).";
+        P->warnIfBestTransformationNotAvailable = false;
+    }
+    pj_log(P->ctx,
+           P->errorIfBestTransformationNotAvailable ? PJ_LOG_ERROR : PJ_LOG_DEBUG,
+           msg.c_str());
+}
+
 /**************************************************************************************/
 PJ_COORD proj_trans (PJ *P, PJ_DIRECTION direction, PJ_COORD coord) {
 /***************************************************************************************
@@ -348,6 +384,12 @@ similarly, but prefers the 2D resp. 3D interfaces if available.
             if( res.xyzt.x != HUGE_VAL ) {
                 return res;
             }
+            else if( P->errorIfBestTransformationNotAvailable ||
+                     P->warnIfBestTransformationNotAvailable ) {
+                warnAboutMissingGrid(alt.pj);
+                if( P->errorIfBestTransformationNotAvailable )
+                    return res;
+            }
             if( iRetry == N_MAX_RETRY ) {
                 break;
             }
@@ -1894,11 +1936,14 @@ PJ  *proj_create_crs_to_crs_from_pj (PJ_CONTEXT *ctx, const PJ *source_crs, cons
     if( !ctx ) {
         ctx = pj_get_default_ctx();
     }
+    pj_load_ini(ctx); // to set ctx->errorIfBestTransformationNotAvailableDefault
 
     const char* authority = nullptr;
     double accuracy = -1;
     bool allowBallparkTransformations = true;
     bool forceOver = false;
+    bool warnIfBestTransformationNotAvailable = ctx->warnIfBestTransformationNotAvailableDefault;
+    bool errorIfBestTransformationNotAvailable = ctx->errorIfBestTransformationNotAvailableDefault;
     for (auto iter = options; iter && iter[0]; ++iter) {
         const char *value;
         if ((value = getOptionValue(*iter, "AUTHORITY="))) {
@@ -1915,6 +1960,17 @@ PJ  *proj_create_crs_to_crs_from_pj (PJ_CONTEXT *ctx, const PJ *source_crs, cons
                             "Invalid value for ALLOW_BALLPARK option.");
                 return nullptr;
             }
+        } else if ((value = getOptionValue(*iter, "ONLY_BEST="))) {
+            warnIfBestTransformationNotAvailable = false;
+            if( ci_equal(value, "yes") )
+                errorIfBestTransformationNotAvailable = true;
+            else if( ci_equal(value, "no") )
+                errorIfBestTransformationNotAvailable = false;
+            else {
+                ctx->logger(ctx->logger_app_data, PJ_LOG_ERROR,
+                            "Invalid value for ONLY_BEST option.");
+                return nullptr;
+            }
         }
         else if ((value = getOptionValue(*iter, "FORCE_OVER="))) {
             if (ci_equal(value, "yes")) {
@@ -1963,7 +2019,9 @@ PJ  *proj_create_crs_to_crs_from_pj (PJ_CONTEXT *ctx, const PJ *source_crs, cons
         ctx, operation_ctx, PROJ_SPATIAL_CRITERION_PARTIAL_INTERSECTION);
     proj_operation_factory_context_set_grid_availability_use(
         ctx, operation_ctx,
-        proj_context_is_network_enabled(ctx) ?
+        (errorIfBestTransformationNotAvailable ||
+         warnIfBestTransformationNotAvailable ||
+         proj_context_is_network_enabled(ctx)) ?
             PROJ_GRID_AVAILABILITY_KNOWN_AVAILABLE:
             PROJ_GRID_AVAILABILITY_DISCARD_OPERATION_IF_MISSING_GRID);
 
@@ -1984,19 +2042,47 @@ PJ  *proj_create_crs_to_crs_from_pj (PJ_CONTEXT *ctx, const PJ *source_crs, cons
 
     ctx->forceOver = forceOver;
 
+    const int old_debug_level = ctx->debug_level;
+    if( errorIfBestTransformationNotAvailable || warnIfBestTransformationNotAvailable )
+        ctx->debug_level = PJ_LOG_NONE;
     PJ* P = proj_list_get(ctx, op_list, 0);
+    ctx->debug_level = old_debug_level;
     assert(P);
 
+    if( P != nullptr ) {
+        P->errorIfBestTransformationNotAvailable = errorIfBestTransformationNotAvailable;
+        P->warnIfBestTransformationNotAvailable = warnIfBestTransformationNotAvailable;
+    }
+
     if( P == nullptr || op_count == 1 ||
         proj_get_type(source_crs) == PJ_TYPE_GEOCENTRIC_CRS ||
         proj_get_type(target_crs) == PJ_TYPE_GEOCENTRIC_CRS ) {
         proj_list_destroy(op_list);
         ctx->forceOver = false;
+
+        if( P != nullptr &&
+            (errorIfBestTransformationNotAvailable ||
+             warnIfBestTransformationNotAvailable) &&
+            !proj_coordoperation_is_instantiable(ctx, P) )
+        {
+            warnAboutMissingGrid(P);
+            if( errorIfBestTransformationNotAvailable ) {
+                proj_destroy(P);
+                return nullptr;
+            }
+        }
+
+        if( P != nullptr ) {
+            P->over = forceOver;
+        }
         return P;
     }
 
+    if( errorIfBestTransformationNotAvailable || warnIfBestTransformationNotAvailable )
+        ctx->debug_level = PJ_LOG_NONE;
     auto preparedOpList = pj_create_prepared_operations(ctx, source_crs, target_crs,
                                                    op_list);
+    ctx->debug_level = old_debug_level;
 
     ctx->forceOver = false;
     proj_list_destroy(op_list);
@@ -2007,6 +2093,12 @@ PJ  *proj_create_crs_to_crs_from_pj (PJ_CONTEXT *ctx, const PJ *source_crs, cons
         return nullptr;
     }
 
+    for( auto& op: preparedOpList ) {
+        op.pj->over = forceOver;
+        op.pj->errorIfBestTransformationNotAvailable = errorIfBestTransformationNotAvailable;
+        op.pj->warnIfBestTransformationNotAvailable = warnIfBestTransformationNotAvailable;
+    }
+
     // If there's finally juste a single result, return it directly
     if( preparedOpList.size() == 1 )
     {
@@ -2019,6 +2111,7 @@ PJ  *proj_create_crs_to_crs_from_pj (PJ_CONTEXT *ctx, const PJ *source_crs, cons
     P->alternativeCoordinateOperations = std::move(preparedOpList);
     // The returned P is rather dummy
     P->descr = "Set of coordinate operations";
+    P->over = forceOver;
     P->iso_obj = nullptr;
     P->fwd = nullptr;
     P->inv = nullptr;

src/apps/cs2cs.cpp

@@ -78,7 +78,8 @@ static const char *oterr = "*\t*"; /* output line for unprojectable input */
 static const char *usage =
     "%s\nusage: %s [-dDeEfIlrstvwW [args]]\n"
     "              [[--area name_or_code] | [--bbox west_long,south_lat,east_long,north_lat]]\n"
-    "              [--authority {name}] [--accuracy {accuracy}] [--no-ballpark] [--3d]\n"
+    "              [--authority {name}] [--3d]\n"
+    "              [--accuracy {accuracy}] [--only-best[=yes|=no]] [--no-ballpark]\n"
     "              [+opt[=arg] ...] [+to +opt[=arg] ...] [file ...]\n";
 
 static double (*informat)(const char *,
@@ -418,6 +419,8 @@ int main(int argc, char **argv) {
     const char* authority = nullptr;
     double accuracy = -1;
     bool allowBallpark = true;
+    bool onlyBestSet = false;
+    bool errorIfBestTransformationNotAvailable = false;
     bool promoteTo3D = false;
 
     /* process run line arguments */
@@ -484,6 +487,15 @@ int main(int argc, char **argv) {
         else if (strcmp(*argv, "--no-ballpark") == 0 ) {
             allowBallpark = false;
         }
+        else if (strcmp(*argv, "--only-best") == 0 ||
+                 strcmp(*argv, "--only-best=yes") == 0 ) {
+            onlyBestSet = true;
+            errorIfBestTransformationNotAvailable = true;
+        }
+        else if (strcmp(*argv, "--only-best=no") == 0 ) {
+            onlyBestSet = true;
+            errorIfBestTransformationNotAvailable = false;
+        }
         else if (strcmp(*argv, "--3d") == 0 ) {
             promoteTo3D = true;
         }
@@ -893,6 +905,14 @@ int main(int argc, char **argv) {
     if( !allowBallpark ) {
         options.push_back("ALLOW_BALLPARK=NO");
     }
+    if( onlyBestSet ) {
+        if( errorIfBestTransformationNotAvailable ) {
+            options.push_back("ONLY_BEST=YES");
+        }
+        else {
+            options.push_back("ONLY_BEST=NO");
+        }
+    }
     options.push_back(nullptr);
     transformation = proj_create_crs_to_crs_from_pj(nullptr, src, dst,
                                                     pj_area, options.data());

src/ctx.cpp

@@ -150,6 +150,8 @@ pj_ctx::pj_ctx(const pj_ctx& other) :
     lastFullErrorMessage(std::string()),
     last_errno(0),
     debug_level(other.debug_level),
+    errorIfBestTransformationNotAvailableDefault(other.errorIfBestTransformationNotAvailableDefault),
+    warnIfBestTransformationNotAvailableDefault(other.warnIfBestTransformationNotAvailableDefault),
     logger(other.logger),
     logger_app_data(other.logger_app_data),
     cpp_context(other.cpp_context ? other.cpp_context->clone(this) : nullptr),