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