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