Line data    Source code

       1             : /******************************************************************************
       2             :  *
       3             :  * Project:  GDAL
       4             :  * Purpose:  gdal "raster proximity" subcommand
       5             :  * Author:   Alessandro Pasotti <elpaso at itopen dot it>
       6             :  *
       7             :  ******************************************************************************
       8             :  * Copyright (c) 2025, Alessandro Pasotti <elpaso at itopen dot it>
       9             :  *
      10             :  * SPDX-License-Identifier: MIT
      11             :  ****************************************************************************/
      12             : 
      13             : #include "gdalalg_raster_proximity.h"
      14             : 
      15             : #include <cmath>
      16             : #include <limits>
      17             : 
      18             : #include "cpl_conv.h"
      19             : 
      20             : #include "gdal_alg.h"
      21             : #include "gdal_priv.h"
      22             : 
      23             : //! @cond Doxygen_Suppress
      24             : 
      25             : #ifndef _
      26             : #define _(x) (x)
      27             : #endif
      28             : 
      29             : /************************************************************************/
      30             : /*     GDALRasterProximityAlgorithm::GDALRasterProximityAlgorithm()     */
      31             : /************************************************************************/
      32             : 
      33         114 : GDALRasterProximityAlgorithm::GDALRasterProximityAlgorithm(bool standaloneStep)
      34             :     : GDALRasterPipelineNonNativelyStreamingAlgorithm(NAME, DESCRIPTION,
      35         114 :                                                       HELP_URL, standaloneStep)
      36             : {
      37         114 :     AddOutputDataTypeArg(&m_outputDataType)
      38             :         .SetChoices("Int8", "UInt16", "Int16", "UInt32", "Int32", "Float32",
      39         114 :                     "Float64")
      40         114 :         .SetHiddenChoices("Byte")
      41         114 :         .SetDefault(m_outputDataType);
      42         114 :     AddBandArg(&m_inputBand);
      43         114 :     AddArg("target-values", 0, _("Target pixel values"), &m_targetPixelValues);
      44         228 :     AddArg("distance-units", 0, _("Distance units"), &m_distanceUnits)
      45         114 :         .SetChoices("pixel", "geo")
      46         114 :         .SetDefault(m_distanceUnits);
      47             :     AddArg("max-distance", 0,
      48             :            _("Maximum distance. The nodata value will be used for pixels "
      49             :              "beyond this distance"),
      50         228 :            &m_maxDistance)
      51         114 :         .SetDefault(m_maxDistance);
      52             :     AddArg("fixed-value", 0,
      53             :            _("Fixed value for the pixels that are within the "
      54             :              "maximum distance (instead of the actual distance)"),
      55         228 :            &m_fixedBufferValue)
      56         114 :         .SetMinValueIncluded(0)
      57         114 :         .SetDefault(m_fixedBufferValue);
      58             :     AddArg("output-nodata", 0,
      59             :            _("Specify a nodata value to use for pixels that are beyond the "
      60             :              "maximum distance"),
      61         228 :            &m_noDataValue)
      62         114 :         .AddHiddenAlias("nodata");
      63         114 : }
      64             : 
      65             : /************************************************************************/
      66             : /*               GDALRasterProximityAlgorithm::RunStep()                */
      67             : /************************************************************************/
      68             : 
      69          25 : bool GDALRasterProximityAlgorithm::RunStep(GDALPipelineStepRunContext &ctxt)
      70             : {
      71          25 :     auto pfnProgress = ctxt.m_pfnProgress;
      72          25 :     auto pProgressData = ctxt.m_pProgressData;
      73             : 
      74          25 :     auto poSrcDS = m_inputDataset[0].GetDatasetRef();
      75          25 :     CPLAssert(poSrcDS);
      76             : 
      77          25 :     GDALDataType outputType = GDT_Float32;
      78          25 :     if (!m_outputDataType.empty())
      79             :     {
      80          25 :         outputType = GDALGetDataTypeByName(m_outputDataType.c_str());
      81             :     }
      82             : 
      83             :     auto poTmpDS = CreateTemporaryDataset(
      84             :         poSrcDS->GetRasterXSize(), poSrcDS->GetRasterYSize(), 1, outputType,
      85          50 :         /* bTiledIfPossible = */ true, poSrcDS, /* bCopyMetadata = */ false);
      86          25 :     if (!poTmpDS)
      87           1 :         return false;
      88             : 
      89          24 :     const auto srcBand = poSrcDS->GetRasterBand(m_inputBand);
      90          24 :     CPLAssert(srcBand);
      91             : 
      92          24 :     const auto dstBand = poTmpDS->GetRasterBand(1);
      93          24 :     CPLAssert(dstBand);
      94             : 
      95             :     // Build options for GDALComputeProximity
      96          48 :     CPLStringList proximityOptions;
      97             : 
      98          24 :     if (GetArg("max-distance")->IsExplicitlySet())
      99             :     {
     100          22 :         proximityOptions.AddString(CPLSPrintf("MAXDIST=%.17g", m_maxDistance));
     101             :     }
     102             : 
     103          24 :     if (GetArg("distance-units")->IsExplicitlySet())
     104             :     {
     105             :         proximityOptions.AddString(
     106          12 :             CPLSPrintf("DISTUNITS=%s", m_distanceUnits.c_str()));
     107             :     }
     108             : 
     109          24 :     if (GetArg("fixed-value")->IsExplicitlySet())
     110             :     {
     111             :         proximityOptions.AddString(
     112           4 :             CPLSPrintf("FIXED_BUF_VAL=%.17g", m_fixedBufferValue));
     113             :     }
     114             : 
     115          24 :     if (GetArg("output-nodata")->IsExplicitlySet())
     116             :     {
     117             :         // Because GDALComputeProximity does all of its work with 32-bit floats,
     118             :         // we need to check that a manually-specified NoData value can be represented
     119             :         // as a 32-bit float. The default NoData value for some integer types
     120             :         // also cannot be represented as 32-bit floats, but they still round-trip OK,
     121             :         // so we only perform the runtime check for manually-specified values.
     122          11 :         if (!std::isnan(m_noDataValue))
     123             :         {
     124          11 :             const float fNoData = static_cast<float>(m_noDataValue);
     125          11 :             if (static_cast<double>(fNoData) != m_noDataValue)
     126             :             {
     127           2 :                 CPLError(CE_Failure, CPLE_AppDefined,
     128             :                          "gdal raster proximity requires the output NoData "
     129             :                          "value to be representable as a 32-bit floating point "
     130             :                          "number, but the specified value of %g cannot.",
     131             :                          m_noDataValue);
     132           2 :                 return false;
     133             :             }
     134             :         }
     135           9 :         dstBand->SetNoDataValue(m_noDataValue);
     136             :     }
     137          13 :     else if (outputType == GDT_Float16 || outputType == GDT_Float32 ||
     138             :              outputType == GDT_Float64)
     139             :     {
     140           5 :         dstBand->SetNoDataValue(std::numeric_limits<double>::quiet_NaN());
     141             :     }
     142             :     else
     143             :     {
     144             :         double dfNoData;
     145           8 :         if (GDALGetDataTypeMinMaxAsDouble(outputType, nullptr, &dfNoData))
     146             :         {
     147           8 :             dstBand->SetNoDataValue(dfNoData);
     148             :         }
     149             :         else
     150             :         {
     151           0 :             CPLError(CE_Failure, CPLE_AppDefined,
     152             :                      "No default NoData value for output data type %s",
     153             :                      GDALGetDataTypeName(outputType));
     154           0 :             return false;
     155             :         }
     156             :     }
     157             : 
     158             :     // Always set this to YES. Note that this was NOT the
     159             :     // default behavior in the python implementation of the utility.
     160          22 :     proximityOptions.AddString("USE_INPUT_NODATA=YES");
     161             : 
     162          22 :     if (GetArg("target-values")->IsExplicitlySet())
     163             :     {
     164          28 :         std::string targetPixelValues;
     165          29 :         for (const auto &value : m_targetPixelValues)
     166             :         {
     167          15 :             if (!targetPixelValues.empty())
     168           1 :                 targetPixelValues += ",";
     169          15 :             targetPixelValues += CPLSPrintf("%.17g", value);
     170             :         }
     171             :         proximityOptions.AddString(
     172          14 :             CPLSPrintf("VALUES=%s", targetPixelValues.c_str()));
     173             :     }
     174             : 
     175          22 :     const auto error = GDALComputeProximity(srcBand, dstBand, proximityOptions,
     176             :                                             pfnProgress, pProgressData);
     177          22 :     if (error == CE_None)
     178             :     {
     179          22 :         if (pfnProgress)
     180           1 :             pfnProgress(1.0, "", pProgressData);
     181          22 :         m_outputDataset.Set(std::move(poTmpDS));
     182             :     }
     183             : 
     184          22 :     return error == CE_None;
     185             : }
     186             : 
     187             : GDALRasterProximityAlgorithmStandalone::
     188             :     ~GDALRasterProximityAlgorithmStandalone() = default;
     189             : 
     190             : //! @endcond