Line data    Source code

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