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         895 : GDALVectorGeomAbstractAlgorithm::GDALVectorGeomAbstractAlgorithm(
      30             :     const std::string &name, const std::string &description,
      31         895 :     const std::string &helpURL, bool standaloneStep, OptionsBase &opts)
      32             :     : GDALVectorPipelineStepAlgorithm(name, description, helpURL,
      33             :                                       standaloneStep),
      34         895 :       m_activeLayer(opts.m_activeLayer)
      35             : {
      36         895 :     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         895 :            &opts.m_geomField);
      41         895 : }
      42             : 
      43             : /************************************************************************/
      44             : /*              GDALVectorGeomAbstractAlgorithm::RunStep()              */
      45             : /************************************************************************/
      46             : 
      47          81 : bool GDALVectorGeomAbstractAlgorithm::RunStep(GDALPipelineStepRunContext &)
      48             : {
      49          81 :     auto poSrcDS = m_inputDataset[0].GetDatasetRef();
      50          81 :     CPLAssert(poSrcDS);
      51          81 :     CPLAssert(m_outputDataset.GetName().empty());
      52          81 :     CPLAssert(!m_outputDataset.GetDatasetRef());
      53             : 
      54          81 :     auto outDS = std::make_unique<GDALVectorPipelineOutputDataset>(*poSrcDS);
      55             : 
      56         165 :     for (auto &&poSrcLayer : poSrcDS->GetLayers())
      57             :     {
      58          90 :         if (m_activeLayer.empty() ||
      59           6 :             m_activeLayer == poSrcLayer->GetDescription())
      60             :         {
      61          81 :             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          81 :     m_outputDataset.Set(std::move(outDS));
      73             : 
      74         162 :     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             :     using ScaledProgressData =
     235             :         std::unique_ptr<void, decltype(&GDALDestroyScaledProgress)>;
     236             : 
     237          34 :     constexpr double CONVERSION_WORK_FRACTION = 0.1;
     238             : 
     239             :     ScaledProgressData convertProgress(
     240             :         GDALCreateScaledProgress(0, CONVERSION_WORK_FRACTION, pfnProgress,
     241             :                                  pProgressData),
     242          68 :         GDALDestroyScaledProgress);
     243          68 :     if (!ConvertInputsToGeos(m_srcLayer, m_geomFieldIndex,
     244          34 :                              convertProgress ? GDALScaledProgress : nullptr,
     245             :                              convertProgress.get()))
     246             :     {
     247           6 :         return false;
     248             :     }
     249             : 
     250             :     ScaledProgressData processProgress(
     251             :         GDALCreateScaledProgress(CONVERSION_WORK_FRACTION, 1.0, pfnProgress,
     252             :                                  pProgressData),
     253          56 :         GDALDestroyScaledProgress);
     254          28 :     if (!ProcessGeos(processProgress ? GDALScaledProgress : nullptr,
     255          84 :                      processProgress.get()) ||
     256          28 :         m_poGeosResultAsCollection == nullptr)
     257             :     {
     258           0 :         return false;
     259             :     }
     260             : 
     261             : #ifdef GDAL_GEOS_NON_STREAMING_ALGORITHM_DATASET_INCREMENTAL
     262          28 :     m_nGeosResultSize =
     263          28 :         GEOSGetNumGeometries_r(m_poGeosContext, m_poGeosResultAsCollection);
     264          28 :     m_papoGeosResults = GEOSGeom_releaseCollection_r(
     265             :         m_poGeosContext, m_poGeosResultAsCollection, &m_nGeosResultSize);
     266          28 :     GEOSGeom_destroy_r(m_poGeosContext, m_poGeosResultAsCollection);
     267          28 :     m_poGeosResultAsCollection = nullptr;
     268          28 :     CPLAssert(m_apoFeatures.size() == m_nGeosResultSize);
     269             : #endif
     270             : 
     271          28 :     return true;
     272             : }
     273             : 
     274             : std::unique_ptr<OGRFeature>
     275        1149 : GDALGeosNonStreamingAlgorithmLayer::GetNextProcessedFeature()
     276             : {
     277        1149 :     GEOSGeometry *poGeosResult = nullptr;
     278             : 
     279        1324 :     while (poGeosResult == nullptr && m_readPos < m_apoFeatures.size())
     280             :     {
     281             :         // Have we already constructed a result OGRGeometry when previously
     282             :         // accessing this feature?
     283        1079 :         if (m_apoFeatures[m_readPos]->GetGeometryRef() != nullptr)
     284             :         {
     285             :             return std::unique_ptr<OGRFeature>(
     286         904 :                 m_apoFeatures[m_readPos++]->Clone());
     287             :         }
     288             : 
     289             : #ifdef GDAL_GEOS_NON_STREAMING_ALGORITHM_DATASET_INCREMENTAL
     290         175 :         poGeosResult = m_papoGeosResults[m_readPos];
     291             : #else
     292             :         poGeosResult = const_cast<GEOSGeometry *>(GEOSGetGeometryN_r(
     293             :             m_poGeosContext, m_poGeosResultAsCollection, m_readPos));
     294             : #endif
     295             : 
     296         175 :         if (poGeosResult != nullptr)
     297             :         {
     298             :             const bool skipFeature =
     299         131 :                 SkipEmpty() && GEOSisEmpty_r(m_poGeosContext, poGeosResult);
     300             : 
     301         131 :             if (skipFeature)
     302             :             {
     303             : #ifdef GDAL_GEOS_NON_STREAMING_ALGORITHM_DATASET_INCREMENTAL
     304          44 :                 GEOSGeom_destroy_r(m_poGeosContext, poGeosResult);
     305          44 :                 m_papoGeosResults[m_readPos] = nullptr;
     306             : #endif
     307          44 :                 poGeosResult = nullptr;
     308             :             }
     309             :         }
     310             : 
     311         175 :         m_readPos++;
     312             :     }
     313             : 
     314         245 :     if (poGeosResult == nullptr)
     315             :     {
     316             : #ifndef GDAL_GEOS_NON_STREAMING_ALGORITHM_DATASET_INCREMENTAL
     317             :         GEOSGeom_destroy_r(m_poGeosContext, m_poGeosResultAsCollection);
     318             :         m_poGeosResultAsCollection = nullptr;
     319             : #endif
     320         158 :         return nullptr;
     321             :     }
     322             : 
     323          87 :     const auto eLayerGeomType = GetLayerDefn()->GetGeomType();
     324             : 
     325             :     std::unique_ptr<OGRGeometry> poResultGeom(
     326         174 :         OGRGeometryFactory::createFromGEOS(m_poGeosContext, poGeosResult));
     327             : 
     328             : #ifdef GDAL_GEOS_NON_STREAMING_ALGORITHM_DATASET_INCREMENTAL
     329          87 :     GEOSGeom_destroy_r(m_poGeosContext, poGeosResult);
     330          87 :     m_papoGeosResults[m_readPos - 1] = nullptr;
     331             : #endif
     332             : 
     333         174 :     if (poResultGeom && eLayerGeomType != wkbUnknown &&
     334          87 :         wkbFlatten(poResultGeom->getGeometryType()) !=
     335          87 :             wkbFlatten(eLayerGeomType))
     336             :     {
     337          66 :         poResultGeom = OGRGeometryFactory::forceTo(std::move(poResultGeom),
     338          33 :                                                    eLayerGeomType);
     339             :     }
     340             : 
     341          87 :     if (poResultGeom == nullptr)
     342             :     {
     343           0 :         CPLError(CE_Failure, CPLE_AppDefined,
     344             :                  "Failed to convert result from GEOS");
     345           0 :         return nullptr;
     346             :     }
     347             : 
     348         174 :     poResultGeom->assignSpatialReference(
     349          87 :         GetLayerDefn()->GetGeomFieldDefn(0)->GetSpatialRef());
     350             : 
     351          87 :     auto poFeature = m_apoFeatures[m_readPos - 1].get();
     352          87 :     poFeature->SetGeometry(std::move(poResultGeom));
     353             : 
     354          87 :     return std::unique_ptr<OGRFeature>(poFeature->Clone());
     355             : }
     356             : 
     357             : #undef GDAL_GEOS_NON_STREAMING_ALGORITHM_DATASET_INCREMENTAL
     358             : 
     359         349 : void GDALGeosNonStreamingAlgorithmLayer::ResetReading()
     360             : {
     361         349 :     m_readPos = 0;
     362         349 : }
     363             : 
     364             : #endif
     365             : 
     366             : //! @endcond