Line data    Source code

       1             : /******************************************************************************
       2             :  *
       3             :  * Project:  GDAL
       4             :  * Purpose:  gdal "vector grid" subcommand
       5             :  * Author:   Even Rouault <even dot rouault at spatialys.com>
       6             :  *
       7             :  ******************************************************************************
       8             :  * Copyright (c) 2025, Even Rouault <even dot rouault at spatialys.com>
       9             :  *
      10             :  * SPDX-License-Identifier: MIT
      11             :  ****************************************************************************/
      12             : 
      13             : #include "gdalalg_vector_grid.h"
      14             : #include "gdalalg_vector_grid_average.h"
      15             : #include "gdalalg_vector_grid_data_metrics.h"
      16             : #include "gdalalg_vector_grid_invdist.h"
      17             : #include "gdalalg_vector_grid_invdistnn.h"
      18             : #include "gdalalg_vector_grid_linear.h"
      19             : #include "gdalalg_vector_grid_nearest.h"
      20             : 
      21             : #include "cpl_conv.h"
      22             : #include "gdal_priv.h"
      23             : #include "gdal_utils.h"
      24             : #include "ogrsf_frmts.h"
      25             : 
      26             : //! @cond Doxygen_Suppress
      27             : 
      28             : #ifndef _
      29             : #define _(x) (x)
      30             : #endif
      31             : 
      32             : /************************************************************************/
      33             : /*          GDALVectorGridAlgorithm::GDALVectorGridAlgorithm()          */
      34             : /************************************************************************/
      35             : 
      36         112 : GDALVectorGridAlgorithm::GDALVectorGridAlgorithm(bool standaloneStep)
      37             :     : GDALPipelineStepAlgorithm(
      38             :           NAME, DESCRIPTION, HELP_URL,
      39           0 :           ConstructorOptions()
      40         112 :               .SetStandaloneStep(standaloneStep)
      41         224 :               .SetOutputFormatCreateCapability(GDAL_DCAP_CREATE))
      42             : {
      43         112 :     if (standaloneStep)
      44             :     {
      45          98 :         RegisterSubAlgorithm<GDALVectorGridAverageAlgorithmStandalone>();
      46          98 :         RegisterSubAlgorithm<GDALVectorGridInvdistAlgorithmStandalone>();
      47          98 :         RegisterSubAlgorithm<GDALVectorGridInvdistNNAlgorithmStandalone>();
      48          98 :         RegisterSubAlgorithm<GDALVectorGridLinearAlgorithmStandalone>();
      49          98 :         RegisterSubAlgorithm<GDALVectorGridNearestAlgorithmStandalone>();
      50          98 :         RegisterSubAlgorithm<GDALVectorGridMinimumAlgorithmStandalone>();
      51          98 :         RegisterSubAlgorithm<GDALVectorGridMaximumAlgorithmStandalone>();
      52          98 :         RegisterSubAlgorithm<GDALVectorGridRangeAlgorithmStandalone>();
      53          98 :         RegisterSubAlgorithm<GDALVectorGridCountAlgorithmStandalone>();
      54             :         RegisterSubAlgorithm<
      55          98 :             GDALVectorGridAverageDistanceAlgorithmStandalone>();
      56             :         RegisterSubAlgorithm<
      57          98 :             GDALVectorGridAverageDistancePointsAlgorithmStandalone>();
      58             :     }
      59             :     else
      60             :     {
      61          14 :         RegisterSubAlgorithm<GDALVectorGridAverageAlgorithm>();
      62          14 :         RegisterSubAlgorithm<GDALVectorGridInvdistAlgorithm>();
      63          14 :         RegisterSubAlgorithm<GDALVectorGridInvdistNNAlgorithm>();
      64          14 :         RegisterSubAlgorithm<GDALVectorGridLinearAlgorithm>();
      65          14 :         RegisterSubAlgorithm<GDALVectorGridNearestAlgorithm>();
      66          14 :         RegisterSubAlgorithm<GDALVectorGridMinimumAlgorithm>();
      67          14 :         RegisterSubAlgorithm<GDALVectorGridMaximumAlgorithm>();
      68          14 :         RegisterSubAlgorithm<GDALVectorGridRangeAlgorithm>();
      69          14 :         RegisterSubAlgorithm<GDALVectorGridCountAlgorithm>();
      70          14 :         RegisterSubAlgorithm<GDALVectorGridAverageDistanceAlgorithm>();
      71          14 :         RegisterSubAlgorithm<GDALVectorGridAverageDistancePointsAlgorithm>();
      72             :     }
      73         112 : }
      74             : 
      75             : /************************************************************************/
      76             : /*                  GDALVectorGridAlgorithm::RunStep()                  */
      77             : /************************************************************************/
      78             : 
      79           0 : bool GDALVectorGridAlgorithm::RunStep(GDALPipelineStepRunContext &)
      80             : {
      81           0 :     CPLError(CE_Failure, CPLE_AppDefined,
      82             :              "The Run() method should not be called directly on the \"gdal "
      83             :              "vector grid\" program.");
      84           0 :     return false;
      85             : }
      86             : 
      87             : /************************************************************************/
      88             : /*                  GDALVectorGeomAlgorithm::RunImpl()                  */
      89             : /************************************************************************/
      90             : 
      91           1 : bool GDALVectorGridAlgorithm::RunImpl(GDALProgressFunc, void *)
      92             : {
      93           1 :     CPLError(CE_Failure, CPLE_AppDefined,
      94             :              "The Run() method should not be called directly on the \"gdal "
      95             :              "vector grid\" program.");
      96           1 :     return false;
      97             : }
      98             : 
      99             : /************************************************************************/
     100             : /*                  GDALVectorGridAbstractAlgorithm()                   */
     101             : /************************************************************************/
     102             : 
     103         201 : GDALVectorGridAbstractAlgorithm::GDALVectorGridAbstractAlgorithm(
     104             :     const std::string &name, const std::string &description,
     105         201 :     const std::string &helpURL, bool standaloneStep)
     106             :     : GDALPipelineStepAlgorithm(
     107             :           name, description, helpURL,
     108           0 :           ConstructorOptions()
     109         201 :               .SetStandaloneStep(standaloneStep)
     110         402 :               .SetOutputFormatCreateCapability(GDAL_DCAP_CREATE))
     111             : {
     112         201 :     AddProgressArg();
     113         201 :     if (standaloneStep)
     114             :     {
     115         178 :         AddOutputFormatArg(&m_format).AddMetadataItem(
     116         534 :             GAAMDI_REQUIRED_CAPABILITIES, {GDAL_DCAP_RASTER, GDAL_DCAP_CREATE});
     117         178 :         AddOpenOptionsArg(&m_openOptions);
     118         178 :         AddInputFormatsArg(&m_inputFormats)
     119         356 :             .AddMetadataItem(GAAMDI_REQUIRED_CAPABILITIES, {GDAL_DCAP_VECTOR});
     120         178 :         AddInputDatasetArg(&m_inputDataset, GDAL_OF_VECTOR);
     121         178 :         AddOutputDatasetArg(&m_outputDataset, GDAL_OF_RASTER);
     122         178 :         AddCreationOptionsArg(&m_creationOptions);
     123             :     }
     124             :     else
     125             :     {
     126          23 :         AddVectorHiddenInputDatasetArg();
     127             :     }
     128             :     AddArg("extent", 0, _("Set the target georeferenced extent"),
     129         402 :            &m_targetExtent)
     130         201 :         .SetMinCount(4)
     131         201 :         .SetMaxCount(4)
     132         201 :         .SetRepeatedArgAllowed(false)
     133         201 :         .SetMetaVar("<xmin>,<ymin>,<xmax>,<ymax>");
     134         402 :     AddArg("resolution", 0, _("Set the target resolution"), &m_targetResolution)
     135         201 :         .SetMinCount(2)
     136         201 :         .SetMaxCount(2)
     137         201 :         .SetRepeatedArgAllowed(false)
     138         402 :         .SetMetaVar("<xres>,<yres>")
     139         201 :         .SetMutualExclusionGroup("size-or-resolution");
     140             :     AddArg("size", 0, _("Set the target size in pixels and lines"),
     141         402 :            &m_targetSize)
     142         201 :         .SetMinCount(2)
     143         201 :         .SetMaxCount(2)
     144         201 :         .SetRepeatedArgAllowed(false)
     145         402 :         .SetMetaVar("<xsize>,<ysize>")
     146         201 :         .SetMutualExclusionGroup("size-or-resolution");
     147         201 :     AddOutputDataTypeArg(&m_outputType).SetDefault(m_outputType);
     148         402 :     AddArg("crs", 0, _("Override the projection for the output file"), &m_crs)
     149         402 :         .AddHiddenAlias("srs")
     150         201 :         .SetIsCRSArg(/*noneAllowed=*/false);
     151         201 :     AddOverwriteArg(&m_overwrite);
     152         201 :     AddLayerNameArg(&m_layers)
     153         402 :         .SetMutualExclusionGroup("layer-sql")
     154         201 :         .AddAlias("layer");
     155         402 :     AddArg("sql", 0, _("SQL statement"), &m_sql)
     156         201 :         .SetReadFromFileAtSyntaxAllowed()
     157         402 :         .SetMetaVar("<statement>|@<filename>")
     158         201 :         .SetRemoveSQLCommentsEnabled()
     159         201 :         .SetMutualExclusionGroup("layer-sql");
     160             :     AddBBOXArg(
     161             :         &m_bbox,
     162         201 :         _("Select only points contained within the specified bounding box"));
     163         402 :     AddArg("zfield", 0, _("Field name from which to get Z values."), &m_zField)
     164         201 :         .AddHiddenAlias("z-field");
     165             :     AddArg("zoffset", 0,
     166             :            _("Value to add to the Z field value (applied before zmultiply)"),
     167         402 :            &m_zOffset)
     168         201 :         .SetDefault(m_zOffset)
     169         201 :         .AddHiddenAlias("z-offset");
     170             :     AddArg("zmultiply", 0,
     171             :            _("Multiplication factor for the Z field value (applied after "
     172             :              "zoffset)"),
     173         402 :            &m_zMultiply)
     174         201 :         .SetDefault(m_zMultiply)
     175         201 :         .AddHiddenAlias("z-multiply");
     176             : 
     177         201 :     AddValidationAction(
     178          94 :         [this]()
     179             :         {
     180          91 :             bool ret = true;
     181          91 :             if (!m_targetResolution.empty() && m_targetExtent.empty())
     182             :             {
     183           1 :                 ReportError(CE_Failure, CPLE_AppDefined,
     184             :                             "'resolution' should be defined when 'extent' is.");
     185           1 :                 ret = false;
     186             :             }
     187          91 :             return ret;
     188             :         });
     189         201 : }
     190             : 
     191             : /************************************************************************/
     192             : /*              GDALVectorGridAbstractAlgorithm::RunStep()              */
     193             : /************************************************************************/
     194             : 
     195          78 : bool GDALVectorGridAbstractAlgorithm::RunStep(GDALPipelineStepRunContext &ctxt)
     196             : {
     197          78 :     auto poSrcDS = m_inputDataset[0].GetDatasetRef();
     198          78 :     CPLAssert(poSrcDS);
     199          78 :     CPLAssert(!m_outputDataset.GetDatasetRef());
     200             : 
     201         156 :     CPLStringList aosOptions;
     202             : 
     203         156 :     std::string outputFilename;
     204          78 :     if (m_standaloneStep)
     205             :     {
     206          77 :         outputFilename = m_outputDataset.GetName();
     207          77 :         if (!m_format.empty())
     208             :         {
     209          73 :             aosOptions.AddString("-of");
     210          73 :             aosOptions.AddString(m_format.c_str());
     211             :         }
     212             : 
     213          78 :         for (const std::string &co : m_creationOptions)
     214             :         {
     215           1 :             aosOptions.AddString("-co");
     216           1 :             aosOptions.AddString(co.c_str());
     217             :         }
     218             :     }
     219             :     else
     220             :     {
     221           1 :         outputFilename = CPLGenerateTempFilenameSafe("_grid.tif");
     222             : 
     223           1 :         aosOptions.AddString("-of");
     224           1 :         aosOptions.AddString("GTiff");
     225             : 
     226           1 :         aosOptions.AddString("-co");
     227           1 :         aosOptions.AddString("TILED=YES");
     228             :     }
     229             : 
     230          78 :     if (!m_targetExtent.empty())
     231             :     {
     232           6 :         aosOptions.AddString("-txe");
     233           6 :         aosOptions.AddString(CPLSPrintf("%.17g", m_targetExtent[0]));
     234           6 :         aosOptions.AddString(CPLSPrintf("%.17g", m_targetExtent[2]));
     235           6 :         aosOptions.AddString("-tye");
     236           6 :         aosOptions.AddString(CPLSPrintf("%.17g", m_targetExtent[1]));
     237           6 :         aosOptions.AddString(CPLSPrintf("%.17g", m_targetExtent[3]));
     238             :     }
     239             : 
     240          78 :     if (!m_bbox.empty())
     241             :     {
     242           2 :         aosOptions.AddString("-clipsrc");
     243          10 :         for (double v : m_bbox)
     244             :         {
     245           8 :             aosOptions.AddString(CPLSPrintf("%.17g", v));
     246             :         }
     247             :     }
     248             : 
     249          78 :     if (!m_targetResolution.empty())
     250             :     {
     251           1 :         aosOptions.AddString("-tr");
     252           3 :         for (double targetResolution : m_targetResolution)
     253             :         {
     254           2 :             aosOptions.AddString(CPLSPrintf("%.17g", targetResolution));
     255             :         }
     256             :     }
     257             : 
     258          78 :     if (!m_targetSize.empty())
     259             :     {
     260           2 :         aosOptions.AddString("-outsize");
     261           6 :         for (int targetSize : m_targetSize)
     262             :         {
     263           4 :             aosOptions.AddString(CPLSPrintf("%d", targetSize));
     264             :         }
     265             :     }
     266             : 
     267          78 :     if (!m_outputType.empty())
     268             :     {
     269          78 :         aosOptions.AddString("-ot");
     270          78 :         aosOptions.AddString(m_outputType.c_str());
     271             :     }
     272             : 
     273          78 :     if (!m_crs.empty())
     274             :     {
     275           1 :         aosOptions.AddString("-a_srs");
     276           1 :         aosOptions.AddString(m_crs.c_str());
     277             :     }
     278             : 
     279          78 :     if (m_sql.empty())
     280             :     {
     281          79 :         for (const std::string &layer : m_layers)
     282             :         {
     283           3 :             aosOptions.AddString("-l");
     284           3 :             aosOptions.AddString(layer.c_str());
     285             :         }
     286             :     }
     287             :     else
     288             :     {
     289           2 :         aosOptions.AddString("-sql");
     290           2 :         aosOptions.AddString(m_sql.c_str());
     291             :     }
     292             : 
     293          78 :     if (m_zOffset != 0)
     294             :     {
     295           1 :         aosOptions.AddString("-z_increase");
     296           1 :         aosOptions.AddString(CPLSPrintf("%.17g", m_zOffset));
     297             :     }
     298             : 
     299          78 :     if (m_zMultiply != 0)
     300             :     {
     301          78 :         aosOptions.AddString("-z_multiply");
     302          78 :         aosOptions.AddString(CPLSPrintf("%.17g", m_zMultiply));
     303             :     }
     304             : 
     305          78 :     if (!m_zField.empty())
     306             :     {
     307           2 :         aosOptions.AddString("-zfield");
     308           2 :         aosOptions.AddString(m_zField.c_str());
     309             :     }
     310          76 :     else if (m_sql.empty())
     311             :     {
     312          75 :         const auto CheckLayer = [this](OGRLayer *poLayer)
     313             :         {
     314             :             auto poFeat =
     315         146 :                 std::unique_ptr<OGRFeature>(poLayer->GetNextFeature());
     316          73 :             poLayer->ResetReading();
     317          73 :             if (poFeat)
     318             :             {
     319          73 :                 const auto poGeom = poFeat->GetGeometryRef();
     320          73 :                 if (!poGeom || !poGeom->Is3D())
     321             :                 {
     322           2 :                     ReportError(
     323             :                         CE_Warning, CPLE_AppDefined,
     324             :                         "At least one geometry of layer '%s' lacks a Z "
     325             :                         "component. You may need to set the 'zfield' argument",
     326           2 :                         poLayer->GetName());
     327           2 :                     return false;
     328             :                 }
     329             :             }
     330          71 :             return true;
     331          74 :         };
     332             : 
     333          74 :         if (m_layers.empty())
     334             :         {
     335         142 :             for (auto poLayer : poSrcDS->GetLayers())
     336             :             {
     337          71 :                 if (!CheckLayer(poLayer))
     338           1 :                     break;
     339             :             }
     340             :         }
     341             :         else
     342             :         {
     343           5 :             for (const std::string &layerName : m_layers)
     344             :             {
     345           3 :                 auto poLayer = poSrcDS->GetLayerByName(layerName.c_str());
     346           3 :                 if (poLayer && !CheckLayer(poLayer))
     347           1 :                     break;
     348             :             }
     349             :         }
     350             :     }
     351             : 
     352          78 :     aosOptions.AddString("-a");
     353          78 :     aosOptions.AddString(GetGridAlgorithm().c_str());
     354             : 
     355          78 :     bool bOK = false;
     356             :     std::unique_ptr<GDALGridOptions, decltype(&GDALGridOptionsFree)> psOptions{
     357          78 :         GDALGridOptionsNew(aosOptions.List(), nullptr), GDALGridOptionsFree};
     358             : 
     359          78 :     if (psOptions)
     360             :     {
     361          78 :         GDALGridOptionsSetProgress(psOptions.get(), ctxt.m_pfnProgress,
     362             :                                    ctxt.m_pProgressData);
     363             : 
     364             :         auto poRetDS = std::unique_ptr<GDALDataset>(GDALDataset::FromHandle(
     365             :             GDALGrid(outputFilename.c_str(), GDALDataset::ToHandle(poSrcDS),
     366         156 :                      psOptions.get(), nullptr)));
     367             : 
     368          78 :         if (poRetDS)
     369             :         {
     370          76 :             bOK = true;
     371             : 
     372          76 :             if (!m_standaloneStep)
     373             :             {
     374           1 :                 VSIUnlink(outputFilename.c_str());
     375           1 :                 poRetDS->MarkSuppressOnClose();
     376             :             }
     377             : 
     378          76 :             m_outputDataset.Set(std::move(poRetDS));
     379             :         }
     380             :     }
     381             : 
     382         156 :     return bOK;
     383             : }
     384             : 
     385             : /************************************************************************/
     386             : /*              GDALVectorGridAbstractAlgorithm::RunImpl()              */
     387             : /************************************************************************/
     388             : 
     389          77 : bool GDALVectorGridAbstractAlgorithm::RunImpl(GDALProgressFunc pfnProgress,
     390             :                                               void *pProgressData)
     391             : {
     392          77 :     GDALPipelineStepRunContext stepCtxt;
     393          77 :     stepCtxt.m_pfnProgress = pfnProgress;
     394          77 :     stepCtxt.m_pProgressData = pProgressData;
     395          77 :     return RunPreStepPipelineValidations() && RunStep(stepCtxt);
     396             : }
     397             : 
     398             : /************************************************************************/
     399             : /*           GDALVectorGridAbstractAlgorithm::AddRadiusArg()            */
     400             : /************************************************************************/
     401             : 
     402         201 : GDALInConstructionAlgorithmArg &GDALVectorGridAbstractAlgorithm::AddRadiusArg()
     403             : {
     404         402 :     return AddArg("radius", 0, _("Radius of the search circle"), &m_radius)
     405         402 :         .SetMutualExclusionGroup("radius");
     406             : }
     407             : 
     408             : /************************************************************************/
     409             : /*      GDALVectorGridAbstractAlgorithm::AddRadius1AndRadius2Arg()      */
     410             : /************************************************************************/
     411             : 
     412         167 : void GDALVectorGridAbstractAlgorithm::AddRadius1AndRadius2Arg()
     413             : {
     414         334 :     AddArg("radius1", 0, _("First axis of the search ellipse"), &m_radius1)
     415         167 :         .SetMutualExclusionGroup("radius");
     416         167 :     AddArg("radius2", 0, _("Second axis of the search ellipse"), &m_radius2);
     417             : 
     418         167 :     AddValidationAction(
     419         182 :         [this]()
     420             :         {
     421          77 :             bool ret = true;
     422          77 :             if (m_radius1 > 0 && m_radius2 == 0)
     423             :             {
     424           2 :                 ReportError(CE_Failure, CPLE_AppDefined,
     425             :                             "'radius2' should be defined when 'radius1' is.");
     426           2 :                 ret = false;
     427             :             }
     428          75 :             else if (m_radius2 > 0 && m_radius1 == 0)
     429             :             {
     430           2 :                 ReportError(CE_Failure, CPLE_AppDefined,
     431             :                             "'radius1' should be defined when 'radius2' is.");
     432           2 :                 ret = false;
     433             :             }
     434          77 :             return ret;
     435             :         });
     436         167 : }
     437             : 
     438             : /************************************************************************/
     439             : /*            GDALVectorGridAbstractAlgorithm::AddAngleArg()            */
     440             : /************************************************************************/
     441             : 
     442         167 : GDALInConstructionAlgorithmArg &GDALVectorGridAbstractAlgorithm::AddAngleArg()
     443             : {
     444             :     return AddArg("angle", 0,
     445             :                   _("Angle of search ellipse rotation in degrees (counter "
     446             :                     "clockwise)"),
     447         334 :                   &m_angle)
     448         334 :         .SetDefault(m_angle);
     449             : }
     450             : 
     451             : /************************************************************************/
     452             : /*          GDALVectorGridAbstractAlgorithm::AddMinPointsArg()          */
     453             : /************************************************************************/
     454             : 
     455             : GDALInConstructionAlgorithmArg &
     456         171 : GDALVectorGridAbstractAlgorithm::AddMinPointsArg()
     457             : {
     458             :     return AddArg("min-points", 0, _("Minimum number of data points to use"),
     459         342 :                   &m_minPoints)
     460         342 :         .SetDefault(m_minPoints);
     461             : }
     462             : 
     463             : /************************************************************************/
     464             : /*          GDALVectorGridAbstractAlgorithm::AddMaxPointsArg()          */
     465             : /************************************************************************/
     466             : 
     467             : GDALInConstructionAlgorithmArg &
     468          99 : GDALVectorGridAbstractAlgorithm::AddMaxPointsArg()
     469             : {
     470             :     return AddArg("max-points", 0, _("Maximum number of data points to use"),
     471         198 :                   &m_maxPoints)
     472         198 :         .SetDefault(m_maxPoints);
     473             : }
     474             : 
     475             : /************************************************************************/
     476             : /*   GDALVectorGridAbstractAlgorithm::AddMinMaxPointsPerQuadrantArg()   */
     477             : /************************************************************************/
     478             : 
     479         171 : void GDALVectorGridAbstractAlgorithm::AddMinMaxPointsPerQuadrantArg()
     480             : {
     481             :     AddArg("min-points-per-quadrant", 0,
     482             :            _("Minimum number of data points to use per quadrant"),
     483         342 :            &m_minPointsPerQuadrant)
     484         171 :         .SetDefault(m_minPointsPerQuadrant);
     485             :     AddArg("max-points-per-quadrant", 0,
     486             :            _("Maximum number of data points to use per quadrant"),
     487         342 :            &m_maxPointsPerQuadrant)
     488         171 :         .SetDefault(m_maxPointsPerQuadrant);
     489         171 : }
     490             : 
     491             : /************************************************************************/
     492             : /*           GDALVectorGridAbstractAlgorithm::AddNodataArg()            */
     493             : /************************************************************************/
     494             : 
     495         201 : GDALInConstructionAlgorithmArg &GDALVectorGridAbstractAlgorithm::AddNodataArg()
     496             : {
     497         402 :     return AddArg("nodata", 0, _("Target nodata value"), &m_nodata)
     498         402 :         .SetDefault(m_nodata);
     499             : }
     500             : 
     501             : GDALVectorGridAlgorithmStandalone::~GDALVectorGridAlgorithmStandalone() =
     502             :     default;
     503             : 
     504             : //! @endcond