Line data    Source code

       1             : /******************************************************************************
       2             :  *
       3             :  * Project:  GDAL
       4             :  * Purpose:  gdal "vector layer-algebra" 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_layer_algebra.h"
      14             : #include "gdalalg_vector_write.h"
      15             : 
      16             : #include "cpl_conv.h"
      17             : #include "gdal_priv.h"
      18             : #include "gdal_utils.h"
      19             : #include "ogr_api.h"
      20             : #include "ogrsf_frmts.h"
      21             : 
      22             : #include <algorithm>
      23             : 
      24             : //! @cond Doxygen_Suppress
      25             : 
      26             : #ifndef _
      27             : #define _(x) (x)
      28             : #endif
      29             : 
      30             : /************************************************************************/
      31             : /*                  GDALVectorLayerAlgebraAlgorithm()                   */
      32             : /************************************************************************/
      33             : 
      34         126 : GDALVectorLayerAlgebraAlgorithm::GDALVectorLayerAlgebraAlgorithm(
      35         126 :     bool standaloneStep)
      36             :     : GDALVectorPipelineStepAlgorithm(
      37             :           NAME, DESCRIPTION, HELP_URL,
      38           0 :           ConstructorOptions()
      39         126 :               .SetStandaloneStep(standaloneStep)
      40         126 :               .SetInputDatasetMaxCount(1)
      41         126 :               .SetAddInputLayerNameArgument(false)
      42         126 :               .SetAddDefaultArguments(false)
      43         126 :               .SetAddUpsertArgument(false)
      44         126 :               .SetAddSkipErrorsArgument(false)
      45         126 :               .SetOutputLayerNameAvailableInPipelineStep(true)
      46         252 :               .SetOutputFormatCreateCapability(GDAL_DCAP_CREATE))
      47             : {
      48         126 :     if (standaloneStep)
      49             :     {
      50          83 :         AddProgressArg();
      51             :     }
      52             : 
      53             :     auto &opArg =
      54         252 :         AddArg("operation", 0, _("Operation to perform"), &m_operation)
      55             :             .SetChoices("union", "intersection", "sym-difference", "identity",
      56         126 :                         "update", "clip", "erase")
      57         126 :             .SetRequired();
      58         126 :     if (standaloneStep)
      59          83 :         opArg.SetPositional();
      60             : 
      61         126 :     if (standaloneStep)
      62             :     {
      63          83 :         AddVectorInputArgs(false);
      64             :     }
      65             :     else
      66             :     {
      67          43 :         AddVectorHiddenInputDatasetArg();
      68             :     }
      69             : 
      70             :     {
      71             :         auto &arg = AddArg("method", 0, _("Method vector dataset"),
      72         252 :                            &m_methodDataset, GDAL_OF_VECTOR)
      73         126 :                         .SetRequired();
      74         126 :         if (standaloneStep)
      75          83 :             arg.SetPositional();
      76             : 
      77         126 :         SetAutoCompleteFunctionForFilename(arg, GDAL_OF_VECTOR);
      78             :     }
      79             : 
      80         126 :     if (standaloneStep)
      81             :     {
      82          83 :         AddVectorOutputArgs(false, false);
      83             :     }
      84             :     else
      85             :     {
      86          43 :         AddOutputLayerNameArg(/* hiddenForCLI = */ false,
      87             :                               /* shortNameOutputLayerAllowed = */ false);
      88             :     }
      89             : 
      90             :     AddArg(GDAL_ARG_NAME_INPUT_LAYER, 0, _("Input layer name"),
      91         126 :            &m_inputLayerName);
      92             : 
      93         126 :     AddArg("method-layer", 0, _("Method layer name"), &m_methodLayerName);
      94             : 
      95         126 :     AddGeometryTypeArg(&m_geometryType);
      96             : 
      97             :     AddArg("input-prefix", 0,
      98         252 :            _("Prefix for fields corresponding to input layer"), &m_inputPrefix)
      99         126 :         .SetCategory(GAAC_ADVANCED);
     100             :     AddArg("input-field", 0, _("Input field(s) to add to output layer"),
     101         252 :            &m_inputFields)
     102         252 :         .SetCategory(GAAC_ADVANCED)
     103         126 :         .SetMutualExclusionGroup("input-field");
     104             :     AddArg("no-input-field", 0, _("Do not add any input field to output layer"),
     105         252 :            &m_noInputFields)
     106         252 :         .SetCategory(GAAC_ADVANCED)
     107         126 :         .SetMutualExclusionGroup("input-field");
     108             :     AddArg("all-input-field", 0, _("Add all input fields to output layer"),
     109         252 :            &m_allInputFields)
     110         252 :         .SetCategory(GAAC_ADVANCED)
     111         126 :         .SetMutualExclusionGroup("input-field");
     112             : 
     113             :     AddArg("method-prefix", 0,
     114             :            _("Prefix for fields corresponding to method layer"),
     115         252 :            &m_methodPrefix)
     116         126 :         .SetCategory(GAAC_ADVANCED);
     117             :     AddArg("method-field", 0, _("Method field(s) to add to output layer"),
     118         252 :            &m_methodFields)
     119         252 :         .SetCategory(GAAC_ADVANCED)
     120         126 :         .SetMutualExclusionGroup("method-field");
     121             :     AddArg("no-method-field", 0,
     122         252 :            _("Do not add any method field to output layer"), &m_noMethodFields)
     123         252 :         .SetCategory(GAAC_ADVANCED)
     124         126 :         .SetMutualExclusionGroup("method-field");
     125             :     AddArg("all-method-field", 0, _("Add all method fields to output layer"),
     126         252 :            &m_allMethodFields)
     127         252 :         .SetCategory(GAAC_ADVANCED)
     128         126 :         .SetMutualExclusionGroup("method-field");
     129         126 : }
     130             : 
     131             : /************************************************************************/
     132             : /*         GDALVectorLayerAlgebraAlgorithm::CanHandleNextStep()         */
     133             : /************************************************************************/
     134             : 
     135           2 : bool GDALVectorLayerAlgebraAlgorithm::CanHandleNextStep(
     136             :     GDALPipelineStepAlgorithm *poNextStep) const
     137             : {
     138           3 :     return poNextStep->GetName() == GDALVectorWriteAlgorithm::NAME &&
     139           3 :            poNextStep->GetOutputFormat() != "stream";
     140             : }
     141             : 
     142             : /************************************************************************/
     143             : /*               GDALRasterPolygonizeAlgorithm::RunImpl()               */
     144             : /************************************************************************/
     145             : 
     146          38 : bool GDALVectorLayerAlgebraAlgorithm::RunImpl(GDALProgressFunc pfnProgress,
     147             :                                               void *pProgressData)
     148             : {
     149          38 :     GDALPipelineStepRunContext stepCtxt;
     150          38 :     stepCtxt.m_pfnProgress = pfnProgress;
     151          38 :     stepCtxt.m_pProgressData = pProgressData;
     152          38 :     return RunPreStepPipelineValidations() && RunStep(stepCtxt);
     153             : }
     154             : 
     155             : /************************************************************************/
     156             : /*              GDALVectorLayerAlgebraAlgorithm::RunImpl()              */
     157             : /************************************************************************/
     158             : 
     159          41 : bool GDALVectorLayerAlgebraAlgorithm::RunStep(GDALPipelineStepRunContext &ctxt)
     160             : {
     161             : #ifdef HAVE_GEOS
     162          41 :     auto poSrcDS = m_inputDataset[0].GetDatasetRef();
     163          41 :     CPLAssert(poSrcDS);
     164          41 :     auto poMethodDS = m_methodDataset.GetDatasetRef();
     165          41 :     CPLAssert(poMethodDS);
     166             : 
     167          41 :     if (poSrcDS == poMethodDS)
     168             :     {
     169           1 :         ReportError(CE_Failure, CPLE_NotSupported,
     170             :                     "Input and method datasets must be different");
     171           1 :         return false;
     172             :     }
     173             : 
     174          40 :     auto poWriteStep = ctxt.m_poNextUsableStep ? ctxt.m_poNextUsableStep : this;
     175             : 
     176          40 :     GDALDataset *poDstDS = nullptr;
     177          40 :     bool bTemporaryFile = false;
     178          40 :     std::unique_ptr<GDALDataset> poNewRetDS;
     179          80 :     std::string outputLayerName;
     180          40 :     OGRLayer *poDstLayer = nullptr;
     181          40 :     if (!CreateDatasetSingleOutputLayerIfNeeded(ctxt, "output", poDstDS,
     182             :                                                 bTemporaryFile, poNewRetDS,
     183             :                                                 outputLayerName, poDstLayer))
     184             :     {
     185          11 :         return false;
     186             :     }
     187             : 
     188             :     OGRLayer *poInputLayer;
     189          29 :     if (m_inputLayerName.empty() && poSrcDS->GetLayerCount() == 1)
     190          28 :         poInputLayer = poSrcDS->GetLayer(0);
     191             :     else
     192           1 :         poInputLayer = poSrcDS->GetLayerByName(m_inputLayerName.c_str());
     193          29 :     if (!poInputLayer)
     194             :     {
     195           1 :         ReportError(CE_Failure, CPLE_AppDefined, "Cannot get input layer '%s'",
     196             :                     m_inputLayerName.c_str());
     197           1 :         return false;
     198             :     }
     199             : 
     200             :     OGRLayer *poMethodLayer;
     201          28 :     if (m_methodLayerName.empty() && poMethodDS->GetLayerCount() == 1)
     202          27 :         poMethodLayer = poMethodDS->GetLayer(0);
     203             :     else
     204           1 :         poMethodLayer = poMethodDS->GetLayerByName(m_methodLayerName.c_str());
     205          28 :     if (!poMethodLayer)
     206             :     {
     207           1 :         ReportError(CE_Failure, CPLE_AppDefined, "Cannot get method layer '%s'",
     208             :                     m_methodLayerName.c_str());
     209           1 :         return false;
     210             :     }
     211             : 
     212          27 :     const auto poInputSRS = poInputLayer->GetSpatialRef();
     213          27 :     const auto poMethodSRS = poMethodLayer->GetSpatialRef();
     214          27 :     if (poInputSRS && !poMethodSRS)
     215             :     {
     216           1 :         ReportError(
     217             :             CE_Warning, CPLE_AppDefined,
     218             :             "Input layer has a CRS, but method layer has none. Assuming "
     219             :             "geometries of method layer to be expressed in input layer CRS.");
     220             :     }
     221          26 :     else if (!poInputSRS && poMethodSRS)
     222             :     {
     223           1 :         ReportError(
     224             :             CE_Warning, CPLE_AppDefined,
     225             :             "Method layer has a CRS, but input layer has none. Assuming "
     226             :             "geometries of input layer to be expressed in method layer CRS.");
     227             :     }
     228          25 :     else if (poInputSRS && poMethodSRS && !poInputSRS->IsSame(poMethodSRS))
     229             :     {
     230           1 :         ReportError(
     231             :             CE_Warning, CPLE_AppDefined,
     232             :             "Input and method layer have non-equivalent CRS. No on-the-fly "
     233             :             "reprojection is performed, and thus results may be incorrect.");
     234             :     }
     235             : 
     236          27 :     if (!poDstLayer)
     237             :     {
     238             :         const CPLStringList aosLayerCreationOptions(
     239          24 :             poWriteStep->GetLayerCreationOptions());
     240             : 
     241             :         const OGRwkbGeometryType eType =
     242          24 :             !m_geometryType.empty() ? OGRFromOGCGeomType(m_geometryType.c_str())
     243          22 :                                     : poInputLayer->GetGeomType();
     244          48 :         poDstLayer = poDstDS->CreateLayer(outputLayerName.c_str(),
     245          24 :                                           poInputLayer->GetSpatialRef(), eType,
     246             :                                           aosLayerCreationOptions.List());
     247          24 :         if (!poDstLayer)
     248           1 :             return false;
     249             :     }
     250             : 
     251          52 :     CPLStringList aosOptions;
     252             : 
     253          26 :     if (m_inputFields.empty() && !m_noInputFields)
     254          19 :         m_allInputFields = true;
     255             : 
     256          26 :     if (m_methodFields.empty() && !m_noMethodFields && !m_allMethodFields)
     257             :     {
     258          36 :         if (m_operation == "update" || m_operation == "clip" ||
     259          17 :             m_operation == "erase")
     260           3 :             m_noMethodFields = true;
     261             :         else
     262          16 :             m_allMethodFields = true;
     263             :     }
     264             : 
     265          26 :     if (m_noInputFields && m_noMethodFields)
     266             :     {
     267           6 :         aosOptions.SetNameValue("ADD_INPUT_FIELDS", "NO");
     268           6 :         aosOptions.SetNameValue("ADD_METHOD_FIELDS", "NO");
     269             :     }
     270             :     else
     271             :     {
     272             :         // Copy fields from input or method layer to output layer
     273             :         const auto CopyFields =
     274          36 :             [poDstLayer](OGRLayer *poSrcLayer, const std::string &prefix,
     275          81 :                          const std::vector<std::string> &srcFields)
     276             :         {
     277             :             const auto contains =
     278           4 :                 [](const std::vector<std::string> &v, const std::string &s)
     279           4 :             { return std::find(v.begin(), v.end(), s) != v.end(); };
     280             : 
     281          36 :             const auto poOutFDefn = poDstLayer->GetLayerDefn();
     282          36 :             const auto poFDefn = poSrcLayer->GetLayerDefn();
     283          36 :             const int nCount = poFDefn->GetFieldCount();
     284          87 :             for (int i = 0; i < nCount; ++i)
     285             :             {
     286          53 :                 const auto poSrcFieldDefn = poFDefn->GetFieldDefn(i);
     287          53 :                 const char *pszName = poSrcFieldDefn->GetNameRef();
     288          53 :                 if (srcFields.empty() || contains(srcFields, pszName))
     289             :                 {
     290          51 :                     OGRFieldDefn oField(*poSrcFieldDefn);
     291          51 :                     const std::string outName = prefix + pszName;
     292          51 :                     whileUnsealing(&oField)->SetName(outName.c_str());
     293          96 :                     if (poOutFDefn->GetFieldIndex(outName.c_str()) < 0 &&
     294          45 :                         poDstLayer->CreateField(&oField) != OGRERR_NONE)
     295             :                     {
     296           2 :                         return false;
     297             :                     }
     298             :                 }
     299             :             }
     300          38 :             for (const auto &name : srcFields)
     301             :             {
     302           4 :                 if (poFDefn->GetFieldIndex(name.c_str()) < 0)
     303             :                 {
     304           2 :                     CPLError(CE_Warning, CPLE_AppDefined, "Unknown field '%s'",
     305             :                              name.c_str());
     306             :                 }
     307             :             }
     308          34 :             return true;
     309          20 :         };
     310             : 
     311          20 :         if (!m_noInputFields)
     312             :         {
     313          79 :             if (!GetArg("input-prefix")->IsExplicitlySet() &&
     314          79 :                 m_inputPrefix.empty() && !m_noMethodFields)
     315             :             {
     316          16 :                 m_inputPrefix = "input_";
     317             :             }
     318          20 :             if (!m_inputPrefix.empty())
     319             :             {
     320          16 :                 aosOptions.SetNameValue("INPUT_PREFIX", m_inputPrefix.c_str());
     321             :             }
     322          20 :             if (!CopyFields(poInputLayer, m_inputPrefix, m_inputFields))
     323           2 :                 return false;
     324             :         }
     325             : 
     326          19 :         if (!m_noMethodFields)
     327             :         {
     328          63 :             if (!GetArg("method-prefix")->IsExplicitlySet() &&
     329          63 :                 m_methodPrefix.empty() && !m_noInputFields)
     330             :             {
     331          15 :                 m_methodPrefix = "method_";
     332             :             }
     333          16 :             if (!m_methodPrefix.empty())
     334             :             {
     335             :                 aosOptions.SetNameValue("METHOD_PREFIX",
     336          15 :                                         m_methodPrefix.c_str());
     337             :             }
     338          16 :             if (!CopyFields(poMethodLayer, m_methodPrefix, m_methodFields))
     339           1 :                 return false;
     340             :         }
     341             :     }
     342             : 
     343             :     aosOptions.SetNameValue(
     344             :         "OUTPUT_GEOMETRY_TYPE",
     345          24 :         OGRToOGCGeomType(poDstLayer->GetGeomType(), false, true));
     346             : 
     347             :     const std::map<std::string, decltype(&OGRLayer::Union)>
     348             :         mapOperationToMethod = {
     349             :             {"union", &OGRLayer::Union},
     350             :             {"intersection", &OGRLayer::Intersection},
     351             :             {"sym-difference", &OGRLayer::SymDifference},
     352             :             {"identity", &OGRLayer::Identity},
     353             :             {"update", &OGRLayer::Update},
     354             :             {"clip", &OGRLayer::Clip},
     355             :             {"erase", &OGRLayer::Erase},
     356         216 :         };
     357             : 
     358          24 :     const auto oIter = mapOperationToMethod.find(m_operation);
     359          24 :     CPLAssert(oIter != mapOperationToMethod.end());
     360          24 :     const auto pFunc = oIter->second;
     361          24 :     bool bOK = (poInputLayer->*pFunc)(poMethodLayer, poDstLayer,
     362          24 :                                       aosOptions.List(), ctxt.m_pfnProgress,
     363          24 :                                       ctxt.m_pProgressData) == OGRERR_NONE;
     364          24 :     if (bOK && poNewRetDS)
     365             :     {
     366          18 :         if (bTemporaryFile)
     367             :         {
     368           2 :             bOK = poNewRetDS->FlushCache() == CE_None;
     369             : #if !defined(__APPLE__)
     370             :             // For some unknown reason, unlinking the file on MacOSX
     371             :             // leads to later "disk I/O error". See https://github.com/OSGeo/gdal/issues/13794
     372           2 :             VSIUnlink(poNewRetDS->GetDescription());
     373             : #endif
     374             :         }
     375             : 
     376          18 :         m_outputDataset.Set(std::move(poNewRetDS));
     377             :     }
     378             : 
     379          24 :     return bOK;
     380             : #else
     381             :     (void)ctxt;
     382             :     ReportError(CE_Failure, CPLE_NotSupported,
     383             :                 "This algorithm is only supported for builds against GEOS");
     384             :     return false;
     385             : #endif
     386             : }
     387             : 
     388             : GDALVectorLayerAlgebraAlgorithmStandalone::
     389             :     ~GDALVectorLayerAlgebraAlgorithmStandalone() = default;
     390             : 
     391             : //! @endcond