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