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