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