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