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