Line data    Source code

       1             : /******************************************************************************
       2             : *
       3             :  * Project:  GDAL
       4             :  * Purpose:  "gdal vector combine" subcommand
       5             :  * Author:   Daniel Baston
       6             :  *
       7             :  ******************************************************************************
       8             :  * Copyright (c) 2025-2026, ISciences LLC
       9             :  *
      10             :  * SPDX-License-Identifier: MIT
      11             :  ****************************************************************************/
      12             : 
      13             : #include "gdalalg_vector_combine.h"
      14             : 
      15             : #include "cpl_error.h"
      16             : #include "gdal_priv.h"
      17             : #include "gdalalg_vector_geom.h"
      18             : #include "ogr_geometry.h"
      19             : 
      20             : #include <algorithm>
      21             : #include <cinttypes>
      22             : #include <optional>
      23             : 
      24             : #ifndef _
      25             : #define _(x) (x)
      26             : #endif
      27             : 
      28             : //! @cond Doxygen_Suppress
      29             : 
      30          62 : GDALVectorCombineAlgorithm::GDALVectorCombineAlgorithm(bool standaloneStep)
      31             :     : GDALVectorPipelineStepAlgorithm(NAME, DESCRIPTION, HELP_URL,
      32          62 :                                       standaloneStep)
      33             : {
      34             :     AddArg("group-by", 0,
      35         124 :            _("Names of field(s) by which inputs should be grouped"), &m_groupBy)
      36             :         .AddValidationAction(
      37          15 :             [this]()
      38             :             {
      39          30 :                 auto fields = m_groupBy;
      40             : 
      41          15 :                 std::sort(fields.begin(), fields.end());
      42          15 :                 if (std::adjacent_find(fields.begin(), fields.end()) !=
      43          30 :                     fields.end())
      44             :                 {
      45           1 :                     CPLError(
      46             :                         CE_Failure, CPLE_AppDefined,
      47             :                         "--group-by must be a list of unique field names.");
      48           1 :                     return false;
      49             :                 }
      50          14 :                 return true;
      51          62 :             });
      52             : 
      53             :     AddArg("keep-nested", 0,
      54             :            _("Avoid combining the components of multipart geometries"),
      55          62 :            &m_keepNested);
      56          62 : }
      57             : 
      58             : namespace
      59             : {
      60             : class GDALVectorCombineOutputLayer final
      61             :     : public GDALVectorNonStreamingAlgorithmLayer
      62             : {
      63             :   public:
      64          19 :     explicit GDALVectorCombineOutputLayer(
      65             :         OGRLayer &srcLayer, int geomFieldIndex,
      66             :         const std::vector<std::string> &groupBy, bool keepNested)
      67          19 :         : GDALVectorNonStreamingAlgorithmLayer(srcLayer, geomFieldIndex),
      68          19 :           m_groupBy(groupBy), m_defn(OGRFeatureDefn::CreateFeatureDefn(
      69          19 :                                   srcLayer.GetLayerDefn()->GetName())),
      70          19 :           m_keepNested(keepNested)
      71             :     {
      72          19 :         m_defn->Reference();
      73             : 
      74          19 :         const OGRFeatureDefn *srcDefn = m_srcLayer.GetLayerDefn();
      75             : 
      76             :         // Copy field definitions for attribute fields used in
      77             :         // --group-by. All other attributes are discarded.
      78          25 :         for (const auto &fieldName : m_groupBy)
      79             :         {
      80             :             // RunStep already checked that the field exists
      81           6 :             const auto iField = srcDefn->GetFieldIndex(fieldName.c_str());
      82           6 :             CPLAssert(iField >= 0);
      83             : 
      84           6 :             m_srcFieldIndices.push_back(iField);
      85           6 :             m_defn->AddFieldDefn(srcDefn->GetFieldDefn(iField));
      86             :         }
      87             : 
      88             :         // Create a new geometry field corresponding to each input geometry
      89             :         // field. An appropriate type is worked out below.
      90          19 :         m_defn->SetGeomType(wkbNone);  // Remove default geometry field
      91          39 :         for (const OGRGeomFieldDefn *srcGeomDefn : srcDefn->GetGeomFields())
      92             :         {
      93          20 :             const auto eSrcGeomType = srcGeomDefn->GetType();
      94          20 :             const bool bHasZ = OGR_GT_HasZ(eSrcGeomType);
      95          20 :             const bool bHasM = OGR_GT_HasM(eSrcGeomType);
      96             : 
      97             :             OGRwkbGeometryType eDstGeomType =
      98          20 :                 OGR_GT_SetModifier(wkbGeometryCollection, bHasZ, bHasM);
      99             : 
     100             :             // If the layer claims to have single-part geometries, choose a more
     101             :             // specific output type like "MultiPoint" rather than "GeometryCollection"
     102          36 :             if (wkbFlatten(eSrcGeomType) != wkbUnknown &&
     103          16 :                 !OGR_GT_IsSubClassOf(wkbFlatten(eSrcGeomType),
     104             :                                      wkbGeometryCollection))
     105             :             {
     106          16 :                 eDstGeomType = OGR_GT_GetCollection(eSrcGeomType);
     107             :             }
     108             : 
     109             :             auto dstGeomDefn = std::make_unique<OGRGeomFieldDefn>(
     110          40 :                 srcGeomDefn->GetNameRef(), eDstGeomType);
     111          20 :             dstGeomDefn->SetSpatialRef(srcGeomDefn->GetSpatialRef());
     112          20 :             m_defn->AddGeomFieldDefn(std::move(dstGeomDefn));
     113             :         }
     114          19 :     }
     115             : 
     116          38 :     ~GDALVectorCombineOutputLayer() override
     117          19 :     {
     118          19 :         m_defn->Release();
     119          38 :     }
     120             : 
     121          19 :     GIntBig GetFeatureCount(int bForce) override
     122             :     {
     123          19 :         if (m_poAttrQuery == nullptr && m_poFilterGeom == nullptr)
     124             :         {
     125          13 :             return static_cast<GIntBig>(m_features.size());
     126             :         }
     127             : 
     128           6 :         return OGRLayer::GetFeatureCount(bForce);
     129             :     }
     130             : 
     131         203 :     const OGRFeatureDefn *GetLayerDefn() const override
     132             :     {
     133         203 :         return m_defn;
     134             :     }
     135             : 
     136           4 :     OGRErr IGetExtent(int iGeomField, OGREnvelope *psExtent,
     137             :                       bool bForce) override
     138             :     {
     139           4 :         return m_srcLayer.GetExtent(iGeomField, psExtent, bForce);
     140             :     }
     141             : 
     142           0 :     OGRErr IGetExtent3D(int iGeomField, OGREnvelope3D *psExtent,
     143             :                         bool bForce) override
     144             :     {
     145           0 :         return m_srcLayer.GetExtent3D(iGeomField, psExtent, bForce);
     146             :     }
     147             : 
     148         160 :     std::unique_ptr<OGRFeature> GetNextProcessedFeature() override
     149             :     {
     150         160 :         if (!m_itFeature)
     151             :         {
     152          59 :             m_itFeature = m_features.begin();
     153             :         }
     154             : 
     155         160 :         if (m_itFeature.value() == m_features.end())
     156             :         {
     157          35 :             return nullptr;
     158             :         }
     159             : 
     160             :         std::unique_ptr<OGRFeature> feature(
     161         250 :             m_itFeature.value()->second->Clone());
     162         125 :         feature->SetFID(m_nProcessedFeaturesRead++);
     163         125 :         ++m_itFeature.value();
     164         125 :         return feature;
     165             :     }
     166             : 
     167          19 :     bool Process(GDALProgressFunc pfnProgress, void *pProgressData) override
     168             :     {
     169          19 :         const int nGeomFields = m_srcLayer.GetLayerDefn()->GetGeomFieldCount();
     170             : 
     171             :         const GIntBig nLayerFeatures =
     172          19 :             m_srcLayer.TestCapability(OLCFastFeatureCount)
     173          19 :                 ? m_srcLayer.GetFeatureCount(false)
     174          19 :                 : -1;
     175             :         const double dfInvLayerFeatures =
     176          19 :             1.0 / std::max(1.0, static_cast<double>(nLayerFeatures));
     177             : 
     178          19 :         GIntBig nFeaturesRead = 0;
     179             : 
     180          38 :         std::vector<std::string> fieldValues(m_srcFieldIndices.size());
     181             :         const auto srcDstFieldMap =
     182          38 :             m_defn->ComputeMapForSetFrom(m_srcLayer.GetLayerDefn(), true);
     183          98 :         for (const auto &srcFeature : m_srcLayer)
     184             :         {
     185         127 :             for (size_t iDstField = 0; iDstField < m_srcFieldIndices.size();
     186             :                  iDstField++)
     187             :             {
     188          48 :                 const int iSrcField = m_srcFieldIndices[iDstField];
     189          48 :                 fieldValues[iDstField] =
     190          48 :                     srcFeature->GetFieldAsString(iSrcField);
     191             :             }
     192             : 
     193             :             OGRFeature *dstFeature;
     194             : 
     195          79 :             if (auto it = m_features.find(fieldValues); it == m_features.end())
     196             :             {
     197          31 :                 it = m_features
     198          62 :                          .insert(std::pair(
     199          93 :                              fieldValues, std::make_unique<OGRFeature>(m_defn)))
     200             :                          .first;
     201          31 :                 dstFeature = it->second.get();
     202             : 
     203          31 :                 dstFeature->SetFrom(srcFeature.get(), srcDstFieldMap.data(),
     204             :                                     false);
     205             : 
     206          63 :                 for (int iGeomField = 0; iGeomField < nGeomFields; iGeomField++)
     207             :                 {
     208             :                     OGRGeomFieldDefn *poGeomDefn =
     209          32 :                         m_defn->GetGeomFieldDefn(iGeomField);
     210          32 :                     const auto eGeomType = poGeomDefn->GetType();
     211             : 
     212             :                     std::unique_ptr<OGRGeometry> poGeom(
     213          32 :                         OGRGeometryFactory::createGeometry(eGeomType));
     214          32 :                     poGeom->assignSpatialReference(poGeomDefn->GetSpatialRef());
     215             : 
     216          32 :                     dstFeature->SetGeomField(iGeomField, std::move(poGeom));
     217             :                 }
     218             :             }
     219             :             else
     220             :             {
     221          48 :                 dstFeature = it->second.get();
     222             :             }
     223             : 
     224         160 :             for (int iGeomField = 0; iGeomField < nGeomFields; iGeomField++)
     225             :             {
     226             :                 OGRGeomFieldDefn *poGeomFieldDefn =
     227          81 :                     m_defn->GetGeomFieldDefn(iGeomField);
     228             : 
     229             :                 std::unique_ptr<OGRGeometry> poSrcGeom(
     230          81 :                     srcFeature->StealGeometry(iGeomField));
     231          81 :                 if (poSrcGeom != nullptr && !poSrcGeom->IsEmpty())
     232             :                 {
     233          79 :                     const auto eSrcType = poSrcGeom->getGeometryType();
     234          79 :                     const auto bSrcIsCollection = OGR_GT_IsSubClassOf(
     235             :                         wkbFlatten(eSrcType), wkbGeometryCollection);
     236             :                     const auto bDstIsUntypedCollection =
     237          79 :                         wkbFlatten(poGeomFieldDefn->GetType()) ==
     238          79 :                         wkbGeometryCollection;
     239             : 
     240             :                     // Did this geometry unexpectedly have Z?
     241          79 :                     if (OGR_GT_HasZ(eSrcType) !=
     242          79 :                         OGR_GT_HasZ(poGeomFieldDefn->GetType()))
     243             :                     {
     244           4 :                         AddZ(iGeomField);
     245             :                     }
     246             : 
     247             :                     // Did this geometry unexpectedly have M?
     248          79 :                     if (OGR_GT_HasM(eSrcType) !=
     249          79 :                         OGR_GT_HasM(poGeomFieldDefn->GetType()))
     250             :                     {
     251          10 :                         AddM(iGeomField);
     252             :                     }
     253             : 
     254             :                     // Do we need to change the output from a typed collection
     255             :                     // like MultiPolygon to a generic GeometryCollection?
     256          79 :                     if (m_keepNested && bSrcIsCollection &&
     257           3 :                         !bDstIsUntypedCollection)
     258             :                     {
     259           2 :                         SetTypeGeometryCollection(iGeomField);
     260             :                     }
     261             : 
     262             :                     OGRGeometryCollection *poDstGeom =
     263          79 :                         cpl::down_cast<OGRGeometryCollection *>(
     264             :                             dstFeature->GetGeomFieldRef(iGeomField));
     265             : 
     266          79 :                     if (m_keepNested || !bSrcIsCollection)
     267             :                     {
     268          76 :                         if (poDstGeom->addGeometry(std::move(poSrcGeom)) !=
     269             :                             OGRERR_NONE)
     270             :                         {
     271           0 :                             CPLError(
     272             :                                 CE_Failure, CPLE_AppDefined,
     273             :                                 "Failed to add geometry of type %s to output "
     274             :                                 "feature of type %s",
     275             :                                 OGRGeometryTypeToName(eSrcType),
     276             :                                 OGRGeometryTypeToName(
     277           0 :                                     poDstGeom->getGeometryType()));
     278           0 :                             return false;
     279             :                         }
     280             :                     }
     281             :                     else
     282             :                     {
     283             :                         std::unique_ptr<OGRGeometryCollection>
     284             :                             poSrcGeomCollection(
     285           3 :                                 poSrcGeom.release()->toGeometryCollection());
     286           6 :                         if (poDstGeom->addGeometryComponents(
     287           6 :                                 std::move(poSrcGeomCollection)) != OGRERR_NONE)
     288             :                         {
     289           0 :                             CPLError(CE_Failure, CPLE_AppDefined,
     290             :                                      "Failed to add components from geometry "
     291             :                                      "of type %s to output "
     292             :                                      "feature of type %s",
     293             :                                      OGRGeometryTypeToName(eSrcType),
     294             :                                      OGRGeometryTypeToName(
     295           0 :                                          poDstGeom->getGeometryType()));
     296           0 :                             return false;
     297             :                         }
     298             :                     }
     299             :                 }
     300             :             }
     301             : 
     302          79 :             if (pfnProgress && nLayerFeatures > 0 &&
     303           0 :                 !pfnProgress(static_cast<double>(++nFeaturesRead) *
     304             :                                  dfInvLayerFeatures,
     305             :                              "", pProgressData))
     306             :             {
     307           0 :                 CPLError(CE_Failure, CPLE_UserInterrupt, "Interrupted by user");
     308           0 :                 return false;
     309             :             }
     310             :         }
     311             : 
     312          19 :         if (pfnProgress)
     313             :         {
     314           0 :             pfnProgress(1.0, "", pProgressData);
     315             :         }
     316             : 
     317          19 :         return true;
     318             :     }
     319             : 
     320          37 :     int TestCapability(const char *pszCap) const override
     321             :     {
     322          37 :         if (EQUAL(pszCap, OLCFastFeatureCount))
     323             :         {
     324           0 :             return true;
     325             :         }
     326             : 
     327          37 :         if (EQUAL(pszCap, OLCStringsAsUTF8) ||
     328          24 :             EQUAL(pszCap, OLCFastGetExtent) ||
     329          22 :             EQUAL(pszCap, OLCFastGetExtent3D) ||
     330          22 :             EQUAL(pszCap, OLCCurveGeometries) ||
     331          18 :             EQUAL(pszCap, OLCMeasuredGeometries) ||
     332          14 :             EQUAL(pszCap, OLCZGeometries))
     333             :         {
     334          26 :             return m_srcLayer.TestCapability(pszCap);
     335             :         }
     336             : 
     337          11 :         return false;
     338             :     }
     339             : 
     340          97 :     void ResetReading() override
     341             :     {
     342          97 :         m_itFeature.reset();
     343          97 :         m_nProcessedFeaturesRead = 0;
     344          97 :     }
     345             : 
     346             :     CPL_DISALLOW_COPY_ASSIGN(GDALVectorCombineOutputLayer)
     347             : 
     348             :   private:
     349          10 :     void AddM(int iGeomField)
     350             :     {
     351             :         OGRGeomFieldDefn *poGeomFieldDefn =
     352          10 :             m_defn->GetGeomFieldDefn(iGeomField);
     353          20 :         whileUnsealing(poGeomFieldDefn)
     354          10 :             ->SetType(OGR_GT_SetM(poGeomFieldDefn->GetType()));
     355             : 
     356          20 :         for (auto &[_, poFeature] : m_features)
     357             :         {
     358          10 :             poFeature->GetGeomFieldRef(iGeomField)->setMeasured(true);
     359             :         }
     360          10 :     }
     361             : 
     362           4 :     void AddZ(int iGeomField)
     363             :     {
     364             :         OGRGeomFieldDefn *poGeomFieldDefn =
     365           4 :             m_defn->GetGeomFieldDefn(iGeomField);
     366           8 :         whileUnsealing(poGeomFieldDefn)
     367           4 :             ->SetType(OGR_GT_SetZ(poGeomFieldDefn->GetType()));
     368             : 
     369           8 :         for (auto &[_, poFeature] : m_features)
     370             :         {
     371           4 :             poFeature->GetGeomFieldRef(iGeomField)->set3D(true);
     372             :         }
     373           4 :     }
     374             : 
     375           2 :     void SetTypeGeometryCollection(int iGeomField)
     376             :     {
     377             :         OGRGeomFieldDefn *poGeomFieldDefn =
     378           2 :             m_defn->GetGeomFieldDefn(iGeomField);
     379           2 :         const bool hasZ = OGR_GT_HasZ(poGeomFieldDefn->GetType());
     380           2 :         const bool hasM = OGR_GT_HasM(poGeomFieldDefn->GetType());
     381             : 
     382           4 :         whileUnsealing(poGeomFieldDefn)
     383           2 :             ->SetType(OGR_GT_SetModifier(wkbGeometryCollection, hasZ, hasM));
     384             : 
     385           4 :         for (auto &[_, poFeature] : m_features)
     386             :         {
     387             :             std::unique_ptr<OGRGeometry> poTmpGeom(
     388           2 :                 poFeature->StealGeometry(iGeomField));
     389           4 :             poTmpGeom = OGRGeometryFactory::forceTo(std::move(poTmpGeom),
     390           2 :                                                     poGeomFieldDefn->GetType());
     391           2 :             CPLAssert(poTmpGeom);
     392           2 :             poFeature->SetGeomField(iGeomField, std::move(poTmpGeom));
     393             :         }
     394           2 :     }
     395             : 
     396             :     const std::vector<std::string> m_groupBy{};
     397             :     std::vector<int> m_srcFieldIndices{};
     398             :     std::map<std::vector<std::string>, std::unique_ptr<OGRFeature>>
     399             :         m_features{};
     400             :     std::optional<decltype(m_features)::const_iterator> m_itFeature{};
     401             :     OGRFeatureDefn *const m_defn;
     402             :     GIntBig m_nProcessedFeaturesRead = 0;
     403             :     const bool m_keepNested;
     404             : };
     405             : }  // namespace
     406             : 
     407          20 : bool GDALVectorCombineAlgorithm::RunStep(GDALPipelineStepRunContext &ctxt)
     408             : {
     409          20 :     auto poSrcDS = m_inputDataset[0].GetDatasetRef();
     410          40 :     auto poDstDS = std::make_unique<GDALVectorNonStreamingAlgorithmDataset>();
     411             : 
     412          40 :     GDALVectorAlgorithmLayerProgressHelper progressHelper(ctxt);
     413             : 
     414          39 :     for (auto &&poSrcLayer : poSrcDS->GetLayers())
     415             :     {
     416          20 :         if (m_inputLayerNames.empty() ||
     417           0 :             std::find(m_inputLayerNames.begin(), m_inputLayerNames.end(),
     418          20 :                       poSrcLayer->GetDescription()) != m_inputLayerNames.end())
     419             :         {
     420          20 :             const auto poSrcLayerDefn = poSrcLayer->GetLayerDefn();
     421          20 :             if (poSrcLayerDefn->GetGeomFieldCount() == 0)
     422             :             {
     423           0 :                 if (m_inputLayerNames.empty())
     424           0 :                     continue;
     425           0 :                 ReportError(CE_Failure, CPLE_AppDefined,
     426             :                             "Specified layer '%s' has no geometry field",
     427           0 :                             poSrcLayer->GetDescription());
     428           0 :                 return false;
     429             :             }
     430             : 
     431             :             // Check that all attributes exist
     432          26 :             for (const auto &fieldName : m_groupBy)
     433             :             {
     434             :                 const int iSrcFieldIndex =
     435           7 :                     poSrcLayerDefn->GetFieldIndex(fieldName.c_str());
     436           7 :                 if (iSrcFieldIndex == -1)
     437             :                 {
     438           1 :                     ReportError(CE_Failure, CPLE_AppDefined,
     439             :                                 "Specified attribute field '%s' does not exist "
     440             :                                 "in layer '%s'",
     441             :                                 fieldName.c_str(),
     442           1 :                                 poSrcLayer->GetDescription());
     443           1 :                     return false;
     444             :                 }
     445             :             }
     446             : 
     447          19 :             progressHelper.AddProcessedLayer(*poSrcLayer);
     448             :         }
     449             :     }
     450             : 
     451          38 :     for ([[maybe_unused]] auto [poSrcLayer, bProcessed, layerProgressFunc,
     452          76 :                                 layerProgressData] : progressHelper)
     453             :     {
     454             :         auto poLayer = std::make_unique<GDALVectorCombineOutputLayer>(
     455          19 :             *poSrcLayer, -1, m_groupBy, m_keepNested);
     456             : 
     457          19 :         if (!poDstDS->AddProcessedLayer(std::move(poLayer), layerProgressFunc,
     458             :                                         layerProgressData.get()))
     459             :         {
     460           0 :             return false;
     461             :         }
     462             :     }
     463             : 
     464          19 :     m_outputDataset.Set(std::move(poDstDS));
     465             : 
     466          19 :     return true;
     467             : }
     468             : 
     469             : GDALVectorCombineAlgorithmStandalone::~GDALVectorCombineAlgorithmStandalone() =
     470             :     default;
     471             : 
     472             : //! @endcond