Line data    Source code

       1             : /******************************************************************************
       2             :  *
       3             :  * Project:  GDAL
       4             :  * Purpose:  "clip" step of "vector pipeline", or "gdal vector clip" 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_clip.h"
      14             : 
      15             : #include "gdal_priv.h"
      16             : #include "gdal_utils.h"
      17             : #include "ogrsf_frmts.h"
      18             : 
      19             : #include <algorithm>
      20             : 
      21             : //! @cond Doxygen_Suppress
      22             : 
      23             : #ifndef _
      24             : #define _(x) (x)
      25             : #endif
      26             : 
      27             : /************************************************************************/
      28             : /*          GDALVectorClipAlgorithm::GDALVectorClipAlgorithm()          */
      29             : /************************************************************************/
      30             : 
      31          74 : GDALVectorClipAlgorithm::GDALVectorClipAlgorithm(bool standaloneStep)
      32             :     : GDALVectorPipelineStepAlgorithm(NAME, DESCRIPTION, HELP_URL,
      33          74 :                                       standaloneStep)
      34             : {
      35          74 :     AddActiveLayerArg(&m_activeLayer);
      36          74 :     AddBBOXArg(&m_bbox, _("Clipping bounding box as xmin,ymin,xmax,ymax"))
      37          74 :         .SetMutualExclusionGroup("bbox-geometry-like");
      38         148 :     AddArg("bbox-crs", 0, _("CRS of clipping bounding box"), &m_bboxCrs)
      39         148 :         .SetIsCRSArg()
      40          74 :         .AddHiddenAlias("bbox_srs");
      41         148 :     AddArg("geometry", 0, _("Clipping geometry (WKT or GeoJSON)"), &m_geometry)
      42          74 :         .SetMutualExclusionGroup("bbox-geometry-like");
      43         148 :     AddArg("geometry-crs", 0, _("CRS of clipping geometry"), &m_geometryCrs)
      44         148 :         .SetIsCRSArg()
      45          74 :         .AddHiddenAlias("geometry_srs");
      46             :     AddArg("like", 0, _("Dataset to use as a template for bounds"),
      47         148 :            &m_likeDataset, GDAL_OF_RASTER | GDAL_OF_VECTOR)
      48         148 :         .SetMetaVar("DATASET")
      49          74 :         .SetMutualExclusionGroup("bbox-geometry-like");
      50             :     AddArg("like-sql", 0, ("SELECT statement to run on the 'like' dataset"),
      51         148 :            &m_likeSQL)
      52         148 :         .SetMetaVar("SELECT-STATEMENT")
      53          74 :         .SetMutualExclusionGroup("sql-where");
      54             :     AddArg("like-layer", 0, ("Name of the layer of the 'like' dataset"),
      55         148 :            &m_likeLayer)
      56          74 :         .SetMetaVar("LAYER-NAME");
      57             :     AddArg("like-where", 0, ("WHERE SQL clause to run on the 'like' dataset"),
      58         148 :            &m_likeWhere)
      59         148 :         .SetMetaVar("WHERE-EXPRESSION")
      60          74 :         .SetMutualExclusionGroup("sql-where");
      61          74 : }
      62             : 
      63             : /************************************************************************/
      64             : /*                   GDALVectorClipAlgorithmLayer                       */
      65             : /************************************************************************/
      66             : 
      67             : namespace
      68             : {
      69             : class GDALVectorClipAlgorithmLayer final : public GDALVectorPipelineOutputLayer
      70             : {
      71             :   public:
      72          32 :     GDALVectorClipAlgorithmLayer(OGRLayer &oSrcLayer,
      73             :                                  std::unique_ptr<OGRGeometry> poClipGeom)
      74          32 :         : GDALVectorPipelineOutputLayer(oSrcLayer),
      75          32 :           m_poClipGeom(std::move(poClipGeom)),
      76          32 :           m_eSrcLayerGeomType(oSrcLayer.GetGeomType()),
      77          64 :           m_eFlattenSrcLayerGeomType(wkbFlatten(m_eSrcLayerGeomType)),
      78          64 :           m_bSrcLayerGeomTypeIsCollection(OGR_GT_IsSubClassOf(
      79          32 :               m_eFlattenSrcLayerGeomType, wkbGeometryCollection)),
      80          64 :           m_poFeatureDefn(oSrcLayer.GetLayerDefn()->Clone())
      81             :     {
      82          32 :         SetDescription(oSrcLayer.GetDescription());
      83          32 :         SetMetadata(oSrcLayer.GetMetadata());
      84          32 :         oSrcLayer.SetSpatialFilter(m_poClipGeom.get());
      85          32 :         m_poFeatureDefn->Reference();
      86          32 :     }
      87             : 
      88          64 :     ~GDALVectorClipAlgorithmLayer() override
      89          32 :     {
      90          32 :         m_poFeatureDefn->Release();
      91          64 :     }
      92             : 
      93         426 :     const OGRFeatureDefn *GetLayerDefn() const override
      94             :     {
      95         426 :         return m_poFeatureDefn;
      96             :     }
      97             : 
      98         500 :     void TranslateFeature(
      99             :         std::unique_ptr<OGRFeature> poSrcFeature,
     100             :         std::vector<std::unique_ptr<OGRFeature>> &apoOutFeatures) override
     101             :     {
     102           0 :         std::unique_ptr<OGRGeometry> poIntersection;
     103         500 :         auto poGeom = poSrcFeature->GetGeometryRef();
     104             : 
     105         500 :         if (poGeom == nullptr)
     106           0 :             return;
     107             : 
     108         500 :         poIntersection.reset(poGeom->Intersection(m_poClipGeom.get()));
     109         500 :         if (!poIntersection)
     110           0 :             return;
     111        1000 :         poIntersection->assignSpatialReference(
     112         500 :             m_poFeatureDefn->GetGeomFieldDefn(0)->GetSpatialRef());
     113             : 
     114         500 :         poSrcFeature->SetFDefnUnsafe(m_poFeatureDefn);
     115             : 
     116         500 :         const auto eSrcGeomType = wkbFlatten(poGeom->getGeometryType());
     117             :         const auto eFeatGeomType =
     118         500 :             wkbFlatten(poIntersection->getGeometryType());
     119         500 :         if (eFeatGeomType != eSrcGeomType &&
     120           8 :             m_eFlattenSrcLayerGeomType != wkbUnknown &&
     121           8 :             m_eFlattenSrcLayerGeomType != eFeatGeomType)
     122             :         {
     123             :             // If the intersection is a collection of geometry and the
     124             :             // layer geometry type is of non-collection type, create
     125             :             // one feature per element of the collection.
     126          10 :             if (!m_bSrcLayerGeomTypeIsCollection &&
     127           2 :                 OGR_GT_IsSubClassOf(eFeatGeomType, wkbGeometryCollection))
     128             :             {
     129             :                 auto poGeomColl = std::unique_ptr<OGRGeometryCollection>(
     130           2 :                     poIntersection.release()->toGeometryCollection());
     131           3 :                 for (const auto *poSubGeom : poGeomColl.get())
     132             :                 {
     133             :                     auto poDstFeature =
     134           4 :                         std::unique_ptr<OGRFeature>(poSrcFeature->Clone());
     135           2 :                     poDstFeature->SetGeometry(poSubGeom);
     136           2 :                     apoOutFeatures.push_back(std::move(poDstFeature));
     137             :                 }
     138             :             }
     139           7 :             else if (OGR_GT_GetCollection(eFeatGeomType) ==
     140           7 :                      m_eFlattenSrcLayerGeomType)
     141             :             {
     142           5 :                 poIntersection.reset(OGRGeometryFactory::forceTo(
     143           5 :                     poIntersection.release(), m_eSrcLayerGeomType));
     144           5 :                 poSrcFeature->SetGeometryDirectly(poIntersection.release());
     145           5 :                 apoOutFeatures.push_back(std::move(poSrcFeature));
     146             :             }
     147           2 :             else if (m_eFlattenSrcLayerGeomType == wkbGeometryCollection)
     148             :             {
     149           2 :                 auto poGeomColl = std::make_unique<OGRGeometryCollection>();
     150           1 :                 poGeomColl->addGeometry(std::move(poIntersection));
     151           1 :                 poSrcFeature->SetGeometryDirectly(poGeomColl.release());
     152           1 :                 apoOutFeatures.push_back(std::move(poSrcFeature));
     153           8 :             }
     154             :             // else discard geometries of incompatible type with the
     155             :             // layer geometry type
     156             :         }
     157             :         else
     158             :         {
     159         492 :             poSrcFeature->SetGeometryDirectly(poIntersection.release());
     160         492 :             apoOutFeatures.push_back(std::move(poSrcFeature));
     161             :         }
     162             :     }
     163             : 
     164          62 :     int TestCapability(const char *pszCap) const override
     165             :     {
     166          62 :         if (EQUAL(pszCap, OLCStringsAsUTF8) ||
     167          49 :             EQUAL(pszCap, OLCCurveGeometries) || EQUAL(pszCap, OLCZGeometries))
     168          20 :             return m_srcLayer.TestCapability(pszCap);
     169          42 :         return false;
     170             :     }
     171             : 
     172             :   private:
     173             :     std::unique_ptr<OGRGeometry> const m_poClipGeom{};
     174             :     const OGRwkbGeometryType m_eSrcLayerGeomType;
     175             :     const OGRwkbGeometryType m_eFlattenSrcLayerGeomType;
     176             :     const bool m_bSrcLayerGeomTypeIsCollection;
     177             :     OGRFeatureDefn *const m_poFeatureDefn;
     178             : 
     179             :     CPL_DISALLOW_COPY_ASSIGN(GDALVectorClipAlgorithmLayer)
     180             : };
     181             : 
     182             : }  // namespace
     183             : 
     184             : /************************************************************************/
     185             : /*                 GDALVectorClipAlgorithm::RunStep()                   */
     186             : /************************************************************************/
     187             : 
     188          36 : bool GDALVectorClipAlgorithm::RunStep(GDALPipelineStepRunContext &)
     189             : {
     190          36 :     auto poSrcDS = m_inputDataset[0].GetDatasetRef();
     191          36 :     CPLAssert(poSrcDS);
     192             : 
     193          36 :     CPLAssert(m_outputDataset.GetName().empty());
     194          36 :     CPLAssert(!m_outputDataset.GetDatasetRef());
     195             : 
     196          36 :     const int nLayerCount = poSrcDS->GetLayerCount();
     197          36 :     bool bSrcLayerHasSRS = false;
     198          56 :     for (int i = 0; i < nLayerCount; ++i)
     199             :     {
     200          35 :         auto poSrcLayer = poSrcDS->GetLayer(i);
     201          35 :         if (poSrcLayer &&
     202          37 :             (m_activeLayer.empty() ||
     203          72 :              m_activeLayer == poSrcLayer->GetDescription()) &&
     204          34 :             poSrcLayer->GetSpatialRef())
     205             :         {
     206          15 :             bSrcLayerHasSRS = true;
     207          15 :             break;
     208             :         }
     209             :     }
     210             : 
     211          72 :     auto [poClipGeom, errMsg] = GetClipGeometry();
     212          36 :     if (!poClipGeom)
     213             :     {
     214           8 :         ReportError(CE_Failure, CPLE_AppDefined, "%s", errMsg.c_str());
     215           8 :         return false;
     216             :     }
     217             : 
     218          28 :     auto poLikeDS = m_likeDataset.GetDatasetRef();
     219          29 :     if (bSrcLayerHasSRS && !poClipGeom->getSpatialReference() && poLikeDS &&
     220           1 :         poLikeDS->GetLayerCount() == 0)
     221             :     {
     222           1 :         ReportError(CE_Warning, CPLE_AppDefined,
     223             :                     "Dataset '%s' has no CRS. Assuming its CRS is the "
     224             :                     "same as the input vector.",
     225           1 :                     poLikeDS->GetDescription());
     226             :     }
     227             : 
     228          28 :     auto outDS = std::make_unique<GDALVectorPipelineOutputDataset>(*poSrcDS);
     229             : 
     230          28 :     bool ret = true;
     231          61 :     for (int i = 0; ret && i < nLayerCount; ++i)
     232             :     {
     233          33 :         auto poSrcLayer = poSrcDS->GetLayer(i);
     234          33 :         ret = (poSrcLayer != nullptr);
     235          33 :         if (ret)
     236             :         {
     237          35 :             if (m_activeLayer.empty() ||
     238           2 :                 m_activeLayer == poSrcLayer->GetDescription())
     239             :             {
     240             :                 auto poClipGeomForLayer =
     241          64 :                     std::unique_ptr<OGRGeometry>(poClipGeom->clone());
     242          40 :                 if (poClipGeomForLayer->getSpatialReference() &&
     243           8 :                     poSrcLayer->GetSpatialRef())
     244             :                 {
     245          14 :                     ret = poClipGeomForLayer->transformTo(
     246           7 :                               poSrcLayer->GetSpatialRef()) == OGRERR_NONE;
     247             :                 }
     248          32 :                 if (ret)
     249             :                 {
     250          64 :                     outDS->AddLayer(
     251             :                         *poSrcLayer,
     252          64 :                         std::make_unique<GDALVectorClipAlgorithmLayer>(
     253          32 :                             *poSrcLayer, std::move(poClipGeomForLayer)));
     254             :                 }
     255             :             }
     256             :             else
     257             :             {
     258           2 :                 outDS->AddLayer(
     259             :                     *poSrcLayer,
     260           2 :                     std::make_unique<GDALVectorPipelinePassthroughLayer>(
     261             :                         *poSrcLayer));
     262             :             }
     263             :         }
     264             :     }
     265             : 
     266          28 :     if (ret)
     267          28 :         m_outputDataset.Set(std::move(outDS));
     268             : 
     269          28 :     return ret;
     270             : }
     271             : 
     272             : GDALVectorClipAlgorithmStandalone::~GDALVectorClipAlgorithmStandalone() =
     273             :     default;
     274             : 
     275             : //! @endcond