Line data    Source code

       1             : /******************************************************************************
       2             :  *
       3             :  * Project:  GDAL
       4             :  * Purpose:  Base classes for some geometry-related vector algorithms
       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_geom.h"
      14             : #include "cpl_enumerate.h"
      15             : 
      16             : #include <algorithm>
      17             : #include <cinttypes>
      18             : 
      19             : //! @cond Doxygen_Suppress
      20             : 
      21             : #ifndef _
      22             : #define _(x) (x)
      23             : #endif
      24             : 
      25             : /************************************************************************/
      26             : /*                  GDALVectorGeomAbstractAlgorithm()                   */
      27             : /************************************************************************/
      28             : 
      29         854 : GDALVectorGeomAbstractAlgorithm::GDALVectorGeomAbstractAlgorithm(
      30             :     const std::string &name, const std::string &description,
      31         854 :     const std::string &helpURL, bool standaloneStep, OptionsBase &opts)
      32             :     : GDALVectorPipelineStepAlgorithm(name, description, helpURL,
      33             :                                       standaloneStep),
      34         854 :       m_activeLayer(opts.m_activeLayer)
      35             : {
      36         854 :     AddActiveLayerArg(&opts.m_activeLayer);
      37             :     AddArg("active-geometry", 0,
      38             :            _("Geometry field name to which to restrict the processing (if not "
      39             :              "specified, all)"),
      40         854 :            &opts.m_geomField);
      41         854 : }
      42             : 
      43             : /************************************************************************/
      44             : /*              GDALVectorGeomAbstractAlgorithm::RunStep()              */
      45             : /************************************************************************/
      46             : 
      47          71 : bool GDALVectorGeomAbstractAlgorithm::RunStep(GDALPipelineStepRunContext &)
      48             : {
      49          71 :     auto poSrcDS = m_inputDataset[0].GetDatasetRef();
      50          71 :     CPLAssert(poSrcDS);
      51          71 :     CPLAssert(m_outputDataset.GetName().empty());
      52          71 :     CPLAssert(!m_outputDataset.GetDatasetRef());
      53             : 
      54          71 :     auto outDS = std::make_unique<GDALVectorPipelineOutputDataset>(*poSrcDS);
      55             : 
      56         145 :     for (auto &&poSrcLayer : poSrcDS->GetLayers())
      57             :     {
      58          80 :         if (m_activeLayer.empty() ||
      59           6 :             m_activeLayer == poSrcLayer->GetDescription())
      60             :         {
      61          71 :             outDS->AddLayer(*poSrcLayer, CreateAlgLayer(*poSrcLayer));
      62             :         }
      63             :         else
      64             :         {
      65           6 :             outDS->AddLayer(
      66             :                 *poSrcLayer,
      67           6 :                 std::make_unique<GDALVectorPipelinePassthroughLayer>(
      68             :                     *poSrcLayer));
      69             :         }
      70             :     }
      71             : 
      72          71 :     m_outputDataset.Set(std::move(outDS));
      73             : 
      74         142 :     return true;
      75             : }
      76             : 
      77             : #ifdef HAVE_GEOS
      78             : 
      79             : // GEOSGeom_releaseCollection allows us to take ownership of the contents of
      80             : // a GeometryCollection. We can then incrementally free the geometries as
      81             : // we write them to features. It requires GEOS >= 3.12.
      82             : #if GEOS_VERSION_MAJOR > 3 ||                                                  \
      83             :     (GEOS_VERSION_MAJOR == 3 && GEOS_VERSION_MINOR >= 12)
      84             : #define GDAL_GEOS_NON_STREAMING_ALGORITHM_DATASET_INCREMENTAL
      85             : #endif
      86             : 
      87             : /************************************************************************/
      88             : /*                  GDALGeosNonStreamingAlgorithmLayer                  */
      89             : /************************************************************************/
      90             : 
      91          34 : GDALGeosNonStreamingAlgorithmLayer::GDALGeosNonStreamingAlgorithmLayer(
      92          34 :     OGRLayer &srcLayer, int geomFieldIndex)
      93             :     : GDALVectorNonStreamingAlgorithmLayer(srcLayer, geomFieldIndex),
      94          34 :       m_poGeosContext{OGRGeometry::createGEOSContext()}
      95             : {
      96          34 : }
      97             : 
      98          34 : GDALGeosNonStreamingAlgorithmLayer::~GDALGeosNonStreamingAlgorithmLayer()
      99             : {
     100          34 :     Cleanup();
     101          34 :     if (m_poGeosContext != nullptr)
     102             :     {
     103          34 :         finishGEOS_r(m_poGeosContext);
     104             :     }
     105          34 : }
     106             : 
     107          68 : void GDALGeosNonStreamingAlgorithmLayer::Cleanup()
     108             : {
     109          68 :     m_readPos = 0;
     110          68 :     m_apoFeatures.clear();
     111             : 
     112          68 :     if (m_poGeosContext != nullptr)
     113             :     {
     114          68 :         for (auto &poGeom : m_apoGeosInputs)
     115             :         {
     116           0 :             GEOSGeom_destroy_r(m_poGeosContext, poGeom);
     117             :         }
     118          68 :         m_apoGeosInputs.clear();
     119             : 
     120          68 :         if (m_poGeosResultAsCollection != nullptr)
     121             :         {
     122           0 :             GEOSGeom_destroy_r(m_poGeosContext, m_poGeosResultAsCollection);
     123           0 :             m_poGeosResultAsCollection = nullptr;
     124             :         }
     125             : 
     126          68 :         if (m_papoGeosResults != nullptr)
     127             :         {
     128         196 :             for (size_t i = 0; i < m_nGeosResultSize; i++)
     129             :             {
     130         168 :                 if (m_papoGeosResults[i] != nullptr)
     131             :                 {
     132          37 :                     GEOSGeom_destroy_r(m_poGeosContext, m_papoGeosResults[i]);
     133             :                 }
     134             :             }
     135             : 
     136          28 :             GEOSFree_r(m_poGeosContext, m_papoGeosResults);
     137          28 :             m_nGeosResultSize = 0;
     138          28 :             m_papoGeosResults = nullptr;
     139             :         }
     140             :     }
     141          68 : }
     142             : 
     143          34 : bool GDALGeosNonStreamingAlgorithmLayer::ConvertInputsToGeos(
     144             :     OGRLayer &srcLayer, int geomFieldIndex, GDALProgressFunc pfnProgress,
     145             :     void *pProgressData)
     146             : {
     147          34 :     const GIntBig nLayerFeatures = srcLayer.TestCapability(OLCFastFeatureCount)
     148          34 :                                        ? srcLayer.GetFeatureCount(false)
     149          34 :                                        : -1;
     150             :     const double dfInvLayerFeatures =
     151          34 :         1.0 / std::max(1.0, static_cast<double>(nLayerFeatures));
     152          34 :     const double dfProgressRatio = dfInvLayerFeatures * 0.5;
     153             : 
     154             :     const bool sameDefn =
     155          34 :         CPL_TO_BOOL(GetLayerDefn()->IsSame(srcLayer.GetLayerDefn()));
     156             : 
     157         202 :     for (auto &feature : srcLayer)
     158             :     {
     159         174 :         const OGRGeometry *poSrcGeom = feature->GetGeomFieldRef(geomFieldIndex);
     160             : 
     161         174 :         if (PolygonsOnly())
     162             :         {
     163             :             const auto eFGType = poSrcGeom
     164         174 :                                      ? wkbFlatten(poSrcGeom->getGeometryType())
     165         174 :                                      : wkbUnknown;
     166         174 :             if (eFGType != wkbPolygon && eFGType != wkbMultiPolygon &&
     167           6 :                 eFGType != wkbCurvePolygon && eFGType != wkbMultiSurface)
     168             :             {
     169           6 :                 CPLError(CE_Failure, CPLE_AppDefined,
     170             :                          "Coverage checking can only be performed on "
     171             :                          "polygonal geometries. Feature %" PRId64
     172             :                          " does not have one",
     173           6 :                          static_cast<int64_t>(feature->GetFID()));
     174           6 :                 return false;
     175             :             }
     176             :         }
     177             : 
     178         168 :         if (poSrcGeom)
     179             :         {
     180             :             GEOSGeometry *geosGeom =
     181         168 :                 poSrcGeom->exportToGEOS(m_poGeosContext, false);
     182         168 :             if (!geosGeom)
     183             :             {
     184             :                 // should not happen normally
     185           0 :                 CPLError(CE_Failure, CPLE_AppDefined,
     186             :                          "Geometry of feature %" PRId64
     187             :                          " failed to convert to GEOS",
     188           0 :                          static_cast<int64_t>(feature->GetFID()));
     189           0 :                 return false;
     190             :             }
     191             : 
     192         168 :             m_apoGeosInputs.push_back(geosGeom);
     193             :         }
     194             :         else
     195             :         {
     196           0 :             m_apoGeosInputs.push_back(GEOSGeom_createEmptyCollection_r(
     197             :                 m_poGeosContext, GEOS_GEOMETRYCOLLECTION));
     198             :         }
     199             : 
     200         168 :         feature->SetGeometry(nullptr);  // free some memory
     201             : 
     202         168 :         if (sameDefn)
     203             :         {
     204          84 :             feature->SetFDefnUnsafe(GetLayerDefn());
     205          84 :             m_apoFeatures.push_back(
     206         168 :                 std::unique_ptr<OGRFeature>(feature.release()));
     207             :         }
     208             :         else
     209             :         {
     210         168 :             auto newFeature = std::make_unique<OGRFeature>(GetLayerDefn());
     211          84 :             newFeature->SetFrom(feature.get(), true);
     212          84 :             newFeature->SetFID(feature->GetFID());
     213          84 :             m_apoFeatures.push_back(std::move(newFeature));
     214             :         }
     215             : 
     216         168 :         if (pfnProgress && nLayerFeatures > 0 &&
     217           0 :             !pfnProgress(static_cast<double>(m_apoFeatures.size()) *
     218             :                              dfProgressRatio,
     219             :                          "", pProgressData))
     220             :         {
     221           0 :             CPLError(CE_Failure, CPLE_UserInterrupt, "Interrupted by user");
     222           0 :             return false;
     223             :         }
     224             :     }
     225             : 
     226          28 :     return true;
     227             : }
     228             : 
     229          34 : bool GDALGeosNonStreamingAlgorithmLayer::Process(GDALProgressFunc pfnProgress,
     230             :                                                  void *pProgressData)
     231             : {
     232          34 :     Cleanup();
     233             : 
     234          34 :     if (!ConvertInputsToGeos(m_srcLayer, m_geomFieldIndex, pfnProgress,
     235             :                              pProgressData))
     236             :     {
     237           6 :         return false;
     238             :     }
     239             : 
     240          28 :     if (!ProcessGeos() || m_poGeosResultAsCollection == nullptr)
     241             :     {
     242           0 :         return false;
     243             :     }
     244             : 
     245             : #ifdef GDAL_GEOS_NON_STREAMING_ALGORITHM_DATASET_INCREMENTAL
     246          28 :     m_nGeosResultSize =
     247          28 :         GEOSGetNumGeometries_r(m_poGeosContext, m_poGeosResultAsCollection);
     248          28 :     m_papoGeosResults = GEOSGeom_releaseCollection_r(
     249             :         m_poGeosContext, m_poGeosResultAsCollection, &m_nGeosResultSize);
     250          28 :     GEOSGeom_destroy_r(m_poGeosContext, m_poGeosResultAsCollection);
     251          28 :     m_poGeosResultAsCollection = nullptr;
     252          28 :     CPLAssert(m_apoFeatures.size() == m_nGeosResultSize);
     253             : #endif
     254             : 
     255          28 :     return true;
     256             : }
     257             : 
     258             : std::unique_ptr<OGRFeature>
     259        1149 : GDALGeosNonStreamingAlgorithmLayer::GetNextProcessedFeature()
     260             : {
     261        1149 :     GEOSGeometry *poGeosResult = nullptr;
     262             : 
     263        1324 :     while (poGeosResult == nullptr && m_readPos < m_apoFeatures.size())
     264             :     {
     265             :         // Have we already constructed a result OGRGeometry when previously
     266             :         // accessing this feature?
     267        1079 :         if (m_apoFeatures[m_readPos]->GetGeometryRef() != nullptr)
     268             :         {
     269             :             return std::unique_ptr<OGRFeature>(
     270         904 :                 m_apoFeatures[m_readPos++]->Clone());
     271             :         }
     272             : 
     273             : #ifdef GDAL_GEOS_NON_STREAMING_ALGORITHM_DATASET_INCREMENTAL
     274         175 :         poGeosResult = m_papoGeosResults[m_readPos];
     275             : #else
     276             :         poGeosResult = const_cast<GEOSGeometry *>(GEOSGetGeometryN_r(
     277             :             m_poGeosContext, m_poGeosResultAsCollection, m_readPos));
     278             : #endif
     279             : 
     280         175 :         if (poGeosResult != nullptr)
     281             :         {
     282             :             const bool skipFeature =
     283         131 :                 SkipEmpty() && GEOSisEmpty_r(m_poGeosContext, poGeosResult);
     284             : 
     285         131 :             if (skipFeature)
     286             :             {
     287             : #ifdef GDAL_GEOS_NON_STREAMING_ALGORITHM_DATASET_INCREMENTAL
     288          44 :                 GEOSGeom_destroy_r(m_poGeosContext, poGeosResult);
     289          44 :                 m_papoGeosResults[m_readPos] = nullptr;
     290             : #endif
     291          44 :                 poGeosResult = nullptr;
     292             :             }
     293             :         }
     294             : 
     295         175 :         m_readPos++;
     296             :     }
     297             : 
     298         245 :     if (poGeosResult == nullptr)
     299             :     {
     300             : #ifndef GDAL_GEOS_NON_STREAMING_ALGORITHM_DATASET_INCREMENTAL
     301             :         GEOSGeom_destroy_r(m_poGeosContext, m_poGeosResultAsCollection);
     302             :         m_poGeosResultAsCollection = nullptr;
     303             : #endif
     304         158 :         return nullptr;
     305             :     }
     306             : 
     307          87 :     const auto eLayerGeomType = GetLayerDefn()->GetGeomType();
     308             : 
     309             :     std::unique_ptr<OGRGeometry> poResultGeom(
     310         174 :         OGRGeometryFactory::createFromGEOS(m_poGeosContext, poGeosResult));
     311             : 
     312             : #ifdef GDAL_GEOS_NON_STREAMING_ALGORITHM_DATASET_INCREMENTAL
     313          87 :     GEOSGeom_destroy_r(m_poGeosContext, poGeosResult);
     314          87 :     m_papoGeosResults[m_readPos - 1] = nullptr;
     315             : #endif
     316             : 
     317         174 :     if (poResultGeom && eLayerGeomType != wkbUnknown &&
     318          87 :         wkbFlatten(poResultGeom->getGeometryType()) !=
     319          87 :             wkbFlatten(eLayerGeomType))
     320             :     {
     321          46 :         poResultGeom = OGRGeometryFactory::forceTo(std::move(poResultGeom),
     322          23 :                                                    eLayerGeomType);
     323             :     }
     324             : 
     325          87 :     if (poResultGeom == nullptr)
     326             :     {
     327           0 :         CPLError(CE_Failure, CPLE_AppDefined,
     328             :                  "Failed to convert result from GEOS");
     329           0 :         return nullptr;
     330             :     }
     331             : 
     332         174 :     poResultGeom->assignSpatialReference(
     333          87 :         GetLayerDefn()->GetGeomFieldDefn(0)->GetSpatialRef());
     334             : 
     335          87 :     auto poFeature = m_apoFeatures[m_readPos - 1].get();
     336          87 :     poFeature->SetGeometry(std::move(poResultGeom));
     337             : 
     338          87 :     return std::unique_ptr<OGRFeature>(poFeature->Clone());
     339             : }
     340             : 
     341             : #undef GDAL_GEOS_NON_STREAMING_ALGORITHM_DATASET_INCREMENTAL
     342             : 
     343         349 : void GDALGeosNonStreamingAlgorithmLayer::ResetReading()
     344             : {
     345         349 :     m_readPos = 0;
     346         349 : }
     347             : 
     348             : #endif
     349             : 
     350             : //! @endcond