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