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          38 : GDALVectorClipAlgorithm::GDALVectorClipAlgorithm(bool standaloneStep)
      32             :     : GDALVectorPipelineStepAlgorithm(NAME, DESCRIPTION, HELP_URL,
      33          38 :                                       standaloneStep)
      34             : {
      35          38 :     AddBBOXArg(&m_bbox, _("Clipping bounding box as xmin,ymin,xmax,ymax"))
      36          38 :         .SetMutualExclusionGroup("bbox-geometry-like");
      37          76 :     AddArg("bbox-crs", 0, _("CRS of clipping bounding box"), &m_bboxCrs)
      38          38 :         .SetIsCRSArg()
      39          38 :         .AddHiddenAlias("bbox_srs");
      40          76 :     AddArg("geometry", 0, _("Clipping geometry (WKT or GeoJSON)"), &m_geometry)
      41          38 :         .SetMutualExclusionGroup("bbox-geometry-like");
      42          76 :     AddArg("geometry-crs", 0, _("CRS of clipping geometry"), &m_geometryCrs)
      43          38 :         .SetIsCRSArg()
      44          38 :         .AddHiddenAlias("geometry_srs");
      45             :     AddArg("like", 0, _("Dataset to use as a template for bounds"),
      46          76 :            &m_likeDataset, GDAL_OF_RASTER | GDAL_OF_VECTOR)
      47          76 :         .SetMetaVar("DATASET")
      48          38 :         .SetMutualExclusionGroup("bbox-geometry-like");
      49             :     AddArg("like-sql", 0, ("SELECT statement to run on the 'like' dataset"),
      50          76 :            &m_likeSQL)
      51          76 :         .SetMetaVar("SELECT-STATEMENT")
      52          38 :         .SetMutualExclusionGroup("sql-where");
      53             :     AddArg("like-layer", 0, ("Name of the layer of the 'like' dataset"),
      54          76 :            &m_likeLayer)
      55          38 :         .SetMetaVar("LAYER-NAME");
      56             :     AddArg("like-where", 0, ("WHERE SQL clause to run on the 'like' dataset"),
      57          76 :            &m_likeWhere)
      58          76 :         .SetMetaVar("WHERE-EXPRESSION")
      59          38 :         .SetMutualExclusionGroup("sql-where");
      60          38 : }
      61             : 
      62             : /************************************************************************/
      63             : /*                   GDALVectorClipAlgorithmDataset                     */
      64             : /************************************************************************/
      65             : 
      66             : namespace
      67             : {
      68             : class GDALVectorClipAlgorithmLayer final : public GDALVectorPipelineOutputLayer
      69             : {
      70             :   public:
      71          28 :     GDALVectorClipAlgorithmLayer(OGRLayer &oSrcLayer,
      72             :                                  std::unique_ptr<OGRGeometry> poClipGeom)
      73          28 :         : GDALVectorPipelineOutputLayer(oSrcLayer),
      74          28 :           m_poClipGeom(std::move(poClipGeom)),
      75          28 :           m_eSrcLayerGeomType(oSrcLayer.GetGeomType()),
      76          56 :           m_eFlattenSrcLayerGeomType(wkbFlatten(m_eSrcLayerGeomType)),
      77          56 :           m_bSrcLayerGeomTypeIsCollection(OGR_GT_IsSubClassOf(
      78          56 :               m_eFlattenSrcLayerGeomType, wkbGeometryCollection))
      79             :     {
      80          28 :         SetDescription(oSrcLayer.GetDescription());
      81          28 :         oSrcLayer.SetSpatialFilter(m_poClipGeom.get());
      82          28 :     }
      83             : 
      84         292 :     OGRFeatureDefn *GetLayerDefn() override
      85             :     {
      86         292 :         return m_srcLayer.GetLayerDefn();
      87             :     }
      88             : 
      89          76 :     void TranslateFeature(
      90             :         std::unique_ptr<OGRFeature> poSrcFeature,
      91             :         std::vector<std::unique_ptr<OGRFeature>> &apoOutFeatures) override
      92             :     {
      93          76 :         auto poGeom = poSrcFeature->GetGeometryRef();
      94          76 :         if (!poGeom)
      95           0 :             return;
      96             : 
      97             :         auto poIntersection = std::unique_ptr<OGRGeometry>(
      98          76 :             poGeom->Intersection(m_poClipGeom.get()));
      99          76 :         if (!poIntersection)
     100           0 :             return;
     101             : 
     102             :         const auto eFeatGeomType =
     103          76 :             wkbFlatten(poIntersection->getGeometryType());
     104          76 :         if (m_eFlattenSrcLayerGeomType != wkbUnknown &&
     105          62 :             m_eFlattenSrcLayerGeomType != eFeatGeomType)
     106             :         {
     107             :             // If the intersection is a collection of geometry and the
     108             :             // layer geometry type is of non-collection type, create
     109             :             // one feature per element of the collection.
     110          10 :             if (!m_bSrcLayerGeomTypeIsCollection &&
     111           2 :                 OGR_GT_IsSubClassOf(eFeatGeomType, wkbGeometryCollection))
     112             :             {
     113             :                 auto poGeomColl = std::unique_ptr<OGRGeometryCollection>(
     114           2 :                     poIntersection.release()->toGeometryCollection());
     115           3 :                 for (const auto *poSubGeom : poGeomColl.get())
     116             :                 {
     117             :                     auto poDstFeature =
     118           4 :                         std::unique_ptr<OGRFeature>(poSrcFeature->Clone());
     119           2 :                     poDstFeature->SetGeometry(poSubGeom);
     120           2 :                     apoOutFeatures.push_back(std::move(poDstFeature));
     121             :                 }
     122             :             }
     123           7 :             else if (OGR_GT_GetCollection(eFeatGeomType) ==
     124           7 :                      m_eFlattenSrcLayerGeomType)
     125             :             {
     126           5 :                 poIntersection.reset(OGRGeometryFactory::forceTo(
     127           5 :                     poIntersection.release(), m_eSrcLayerGeomType));
     128           5 :                 poSrcFeature->SetGeometryDirectly(poIntersection.release());
     129           5 :                 apoOutFeatures.push_back(std::move(poSrcFeature));
     130             :             }
     131           2 :             else if (m_eFlattenSrcLayerGeomType == wkbGeometryCollection)
     132             :             {
     133           2 :                 auto poGeomColl = std::make_unique<OGRGeometryCollection>();
     134           1 :                 poGeomColl->addGeometry(std::move(poIntersection));
     135           1 :                 poSrcFeature->SetGeometryDirectly(poGeomColl.release());
     136           1 :                 apoOutFeatures.push_back(std::move(poSrcFeature));
     137           8 :             }
     138             :             // else discard geometries of incompatible type with the
     139             :             // layer geometry type
     140             :         }
     141             :         else
     142             :         {
     143          68 :             poSrcFeature->SetGeometryDirectly(poIntersection.release());
     144          68 :             apoOutFeatures.push_back(std::move(poSrcFeature));
     145             :         }
     146             :     }
     147             : 
     148          22 :     int TestCapability(const char *pszCap) override
     149             :     {
     150          22 :         if (EQUAL(pszCap, OLCStringsAsUTF8) ||
     151          22 :             EQUAL(pszCap, OLCCurveGeometries) || EQUAL(pszCap, OLCZGeometries))
     152           0 :             return m_srcLayer.TestCapability(pszCap);
     153          22 :         return false;
     154             :     }
     155             : 
     156             :   private:
     157             :     std::unique_ptr<OGRGeometry> m_poClipGeom{};
     158             :     const OGRwkbGeometryType m_eSrcLayerGeomType;
     159             :     const OGRwkbGeometryType m_eFlattenSrcLayerGeomType;
     160             :     const bool m_bSrcLayerGeomTypeIsCollection;
     161             : 
     162             :     CPL_DISALLOW_COPY_ASSIGN(GDALVectorClipAlgorithmLayer)
     163             : };
     164             : 
     165             : }  // namespace
     166             : 
     167             : /************************************************************************/
     168             : /*                           LoadGeometry()                             */
     169             : /************************************************************************/
     170             : 
     171           8 : static std::unique_ptr<OGRGeometry> LoadGeometry(GDALDataset *poDS,
     172             :                                                  const std::string &osSQL,
     173             :                                                  const std::string &osLyr,
     174             :                                                  const std::string &osWhere)
     175             : {
     176           8 :     OGRLayer *poLyr = nullptr;
     177           8 :     if (!osSQL.empty())
     178           1 :         poLyr = poDS->ExecuteSQL(osSQL.c_str(), nullptr, nullptr);
     179           7 :     else if (!osLyr.empty())
     180           2 :         poLyr = poDS->GetLayerByName(osLyr.c_str());
     181             :     else
     182           5 :         poLyr = poDS->GetLayer(0);
     183             : 
     184           8 :     if (poLyr == nullptr)
     185             :     {
     186           1 :         CPLError(CE_Failure, CPLE_AppDefined,
     187             :                  "Failed to identify source layer from clipping dataset.");
     188           1 :         return nullptr;
     189             :     }
     190             : 
     191           7 :     if (!osWhere.empty())
     192           2 :         poLyr->SetAttributeFilter(osWhere.c_str());
     193             : 
     194          14 :     OGRGeometryCollection oGC;
     195             : 
     196           7 :     const auto poSRSSrc = poLyr->GetSpatialRef();
     197           7 :     if (poSRSSrc)
     198             :     {
     199           1 :         auto poSRSClone = poSRSSrc->Clone();
     200           1 :         oGC.assignSpatialReference(poSRSClone);
     201           1 :         poSRSClone->Release();
     202             :     }
     203             : 
     204          12 :     for (auto &poFeat : poLyr)
     205             :     {
     206           6 :         auto poSrcGeom = std::unique_ptr<OGRGeometry>(poFeat->StealGeometry());
     207           6 :         if (poSrcGeom)
     208             :         {
     209             :             // Only take into account areal geometries.
     210           6 :             if (poSrcGeom->getDimension() == 2)
     211             :             {
     212           6 :                 if (!poSrcGeom->IsValid())
     213             :                 {
     214           1 :                     CPLError(CE_Failure, CPLE_AppDefined,
     215             :                              "Geometry of feature " CPL_FRMT_GIB " of %s "
     216             :                              "is invalid.",
     217           1 :                              poFeat->GetFID(), poDS->GetDescription());
     218           1 :                     return nullptr;
     219             :                 }
     220             :                 else
     221             :                 {
     222           5 :                     oGC.addGeometry(std::move(poSrcGeom));
     223             :                 }
     224             :             }
     225             :         }
     226             :     }
     227             : 
     228           6 :     if (!osSQL.empty())
     229           1 :         poDS->ReleaseResultSet(poLyr);
     230             : 
     231           6 :     if (oGC.IsEmpty())
     232             :     {
     233           1 :         CPLError(CE_Failure, CPLE_AppDefined, "No clipping geometry found");
     234           1 :         return nullptr;
     235             :     }
     236             : 
     237           5 :     return std::unique_ptr<OGRGeometry>(oGC.UnaryUnion());
     238             : }
     239             : 
     240             : /************************************************************************/
     241             : /*                 GDALVectorClipAlgorithm::RunStep()                   */
     242             : /************************************************************************/
     243             : 
     244          32 : bool GDALVectorClipAlgorithm::RunStep(GDALProgressFunc, void *)
     245             : {
     246          32 :     CPLAssert(m_inputDataset.GetDatasetRef());
     247          32 :     CPLAssert(m_outputDataset.GetName().empty());
     248          32 :     CPLAssert(!m_outputDataset.GetDatasetRef());
     249             : 
     250          32 :     auto poSrcDS = m_inputDataset.GetDatasetRef();
     251             : 
     252          32 :     std::unique_ptr<OGRGeometry> poClipGeom;
     253             : 
     254          32 :     const int nLayerCount = poSrcDS->GetLayerCount();
     255          32 :     bool bSrcLayerHasSRS = false;
     256          49 :     for (int i = 0; i < nLayerCount; ++i)
     257             :     {
     258          30 :         auto poSrcLayer = poSrcDS->GetLayer(i);
     259          30 :         if (poSrcLayer && poSrcLayer->GetSpatialRef())
     260             :         {
     261          13 :             bSrcLayerHasSRS = true;
     262          13 :             break;
     263             :         }
     264             :     }
     265             : 
     266          32 :     if (!m_bbox.empty())
     267             :     {
     268          20 :         poClipGeom = std::make_unique<OGRPolygon>(m_bbox[0], m_bbox[1],
     269          20 :                                                   m_bbox[2], m_bbox[3]);
     270             : 
     271          10 :         if (!m_bboxCrs.empty())
     272             :         {
     273           1 :             auto poSRS = new OGRSpatialReference();
     274           1 :             poSRS->SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
     275           1 :             CPL_IGNORE_RET_VAL(poSRS->SetFromUserInput(m_bboxCrs.c_str()));
     276           1 :             poClipGeom->assignSpatialReference(poSRS);
     277           1 :             poSRS->Release();
     278             :         }
     279             :     }
     280          22 :     else if (!m_geometry.empty())
     281             :     {
     282             :         {
     283          16 :             CPLErrorStateBackuper oErrorStateBackuper(CPLQuietErrorHandler);
     284           8 :             auto [poGeom, eErr] =
     285          16 :                 OGRGeometryFactory::createFromWkt(m_geometry.c_str());
     286           8 :             if (eErr == OGRERR_NONE)
     287             :             {
     288           5 :                 poClipGeom = std::move(poGeom);
     289             :             }
     290             :             else
     291             :             {
     292           3 :                 poClipGeom.reset(
     293             :                     OGRGeometryFactory::createFromGeoJson(m_geometry.c_str()));
     294           3 :                 if (poClipGeom)
     295             :                 {
     296           2 :                     auto poSRS = new OGRSpatialReference();
     297           2 :                     poSRS->SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
     298           2 :                     CPL_IGNORE_RET_VAL(poSRS->SetFromUserInput("WGS84"));
     299           2 :                     poClipGeom->assignSpatialReference(poSRS);
     300           2 :                     poSRS->Release();
     301             :                 }
     302             :             }
     303             :         }
     304           8 :         if (!poClipGeom)
     305             :         {
     306           1 :             ReportError(
     307             :                 CE_Failure, CPLE_AppDefined,
     308             :                 "Clipping geometry is neither a valid WKT or GeoJSON geometry");
     309           1 :             return false;
     310             :         }
     311             : 
     312           7 :         if (!m_geometryCrs.empty())
     313             :         {
     314           2 :             auto poSRS = new OGRSpatialReference();
     315           2 :             poSRS->SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
     316           2 :             CPL_IGNORE_RET_VAL(poSRS->SetFromUserInput(m_geometryCrs.c_str()));
     317           2 :             poClipGeom->assignSpatialReference(poSRS);
     318           2 :             poSRS->Release();
     319             :         }
     320             :     }
     321          14 :     else if (auto poLikeDS = m_likeDataset.GetDatasetRef())
     322             :     {
     323          15 :         if (poLikeDS->GetLayerCount() > 1 && m_likeLayer.empty() &&
     324           2 :             m_likeSQL.empty())
     325             :         {
     326           1 :             ReportError(
     327             :                 CE_Failure, CPLE_AppDefined,
     328             :                 "Only single layer dataset can be specified with --like when "
     329             :                 "neither --like-layer or --like-sql have been specified");
     330           1 :             return false;
     331             :         }
     332          12 :         else if (poLikeDS->GetLayerCount() > 0)
     333             :         {
     334             :             poClipGeom =
     335           8 :                 LoadGeometry(poLikeDS, m_likeSQL, m_likeLayer, m_likeWhere);
     336           8 :             if (!poClipGeom)
     337           3 :                 return false;
     338             :         }
     339           4 :         else if (poLikeDS->GetRasterCount() > 0)
     340             :         {
     341             :             double adfGT[6];
     342           3 :             if (poLikeDS->GetGeoTransform(adfGT) != CE_None)
     343             :             {
     344           1 :                 ReportError(
     345             :                     CE_Failure, CPLE_AppDefined,
     346             :                     "Dataset '%s' has no geotransform matrix. Its bounds "
     347             :                     "cannot be established.",
     348           1 :                     poLikeDS->GetDescription());
     349           1 :                 return false;
     350             :             }
     351           2 :             auto poLikeSRS = poLikeDS->GetSpatialRef();
     352           2 :             if (bSrcLayerHasSRS && !poLikeSRS)
     353             :             {
     354           0 :                 ReportError(CE_Warning, CPLE_AppDefined,
     355             :                             "Dataset '%s' has no SRS. Assuming its SRS is the "
     356             :                             "same as the input vector.",
     357           0 :                             poLikeDS->GetDescription());
     358             :             }
     359           2 :             const double dfTLX = adfGT[0];
     360           2 :             const double dfTLY = adfGT[3];
     361             :             const double dfTRX =
     362           2 :                 adfGT[0] + poLikeDS->GetRasterXSize() * adfGT[1];
     363             :             const double dfTRY =
     364           2 :                 adfGT[3] + poLikeDS->GetRasterXSize() * adfGT[4];
     365             :             const double dfBLX =
     366           2 :                 adfGT[0] + poLikeDS->GetRasterYSize() * adfGT[2];
     367             :             const double dfBLY =
     368           2 :                 adfGT[3] + poLikeDS->GetRasterYSize() * adfGT[5];
     369           4 :             const double dfBRX = adfGT[0] +
     370           2 :                                  poLikeDS->GetRasterXSize() * adfGT[1] +
     371           2 :                                  poLikeDS->GetRasterYSize() * adfGT[2];
     372           4 :             const double dfBRY = adfGT[3] +
     373           2 :                                  poLikeDS->GetRasterXSize() * adfGT[4] +
     374           2 :                                  poLikeDS->GetRasterYSize() * adfGT[5];
     375             : 
     376           4 :             auto poPoly = std::make_unique<OGRPolygon>();
     377           4 :             auto poLR = std::make_unique<OGRLinearRing>();
     378           2 :             poLR->addPoint(dfTLX, dfTLY);
     379           2 :             poLR->addPoint(dfTRX, dfTRY);
     380           2 :             poLR->addPoint(dfBRX, dfBRY);
     381           2 :             poLR->addPoint(dfBLX, dfBLY);
     382           2 :             poLR->addPoint(dfTLX, dfTLY);
     383           2 :             poPoly->addRingDirectly(poLR.release());
     384           2 :             poPoly->assignSpatialReference(poLikeSRS);
     385           2 :             poClipGeom = std::move(poPoly);
     386             :         }
     387             :         else
     388             :         {
     389           1 :             ReportError(CE_Failure, CPLE_AppDefined,
     390             :                         "Cannot get extent from clip dataset");
     391           1 :             return false;
     392             :         }
     393             :     }
     394             :     else
     395             :     {
     396           1 :         ReportError(CE_Failure, CPLE_AppDefined,
     397             :                     "--bbox, --geometry or --like must be specified");
     398           1 :         return false;
     399             :     }
     400             : 
     401          24 :     auto outDS = std::make_unique<GDALVectorPipelineOutputDataset>(*poSrcDS);
     402             : 
     403          24 :     bool ret = true;
     404          52 :     for (int i = 0; ret && i < nLayerCount; ++i)
     405             :     {
     406          28 :         auto poSrcLayer = poSrcDS->GetLayer(i);
     407          28 :         ret = (poSrcLayer != nullptr);
     408          28 :         if (ret)
     409             :         {
     410             :             auto poClipGeomForLayer =
     411          56 :                 std::unique_ptr<OGRGeometry>(poClipGeom->clone());
     412          34 :             if (poClipGeomForLayer->getSpatialReference() &&
     413           6 :                 poSrcLayer->GetSpatialRef())
     414             :             {
     415          10 :                 ret = poClipGeomForLayer->transformTo(
     416           5 :                           poSrcLayer->GetSpatialRef()) == OGRERR_NONE;
     417             :             }
     418          28 :             if (ret)
     419             :             {
     420          56 :                 outDS->AddLayer(
     421             :                     *poSrcLayer,
     422          56 :                     std::make_unique<GDALVectorClipAlgorithmLayer>(
     423          28 :                         *poSrcLayer, std::move(poClipGeomForLayer)));
     424             :             }
     425             :         }
     426             :     }
     427             : 
     428          24 :     if (ret)
     429          24 :         m_outputDataset.Set(std::move(outDS));
     430             : 
     431          24 :     return ret;
     432             : }
     433             : 
     434             : //! @endcond