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