Line data    Source code

       1             : /******************************************************************************
       2             :  *
       3             :  * Project:  GDAL
       4             :  * Purpose:  Common code for gdalalg_raster_clip and gdalalg_vector_clip
       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_clip_common.h"
      14             : 
      15             : #include "ogrsf_frmts.h"
      16             : 
      17             : //! @cond Doxygen_Suppress
      18             : 
      19             : #ifndef _
      20             : #define _(x) (x)
      21             : #endif
      22             : 
      23             : /************************************************************************/
      24             : /*                          ~GDALClipCommon()                           */
      25             : /************************************************************************/
      26             : 
      27             : GDALClipCommon::~GDALClipCommon() = default;
      28             : 
      29             : /************************************************************************/
      30             : /*                            LoadGeometry()                            */
      31             : /************************************************************************/
      32             : 
      33             : std::pair<std::unique_ptr<OGRGeometry>, std::string>
      34          17 : GDALClipCommon::LoadGeometry()
      35             : {
      36          17 :     auto poDS = m_likeDataset.GetDatasetRef();
      37          17 :     OGRLayer *poLyr = nullptr;
      38          17 :     if (!m_likeSQL.empty())
      39           1 :         poLyr = poDS->ExecuteSQL(m_likeSQL.c_str(), nullptr, nullptr);
      40          16 :     else if (!m_likeLayer.empty())
      41           2 :         poLyr = poDS->GetLayerByName(m_likeLayer.c_str());
      42             :     else
      43          14 :         poLyr = poDS->GetLayer(0);
      44             : 
      45          17 :     if (poLyr == nullptr)
      46             :     {
      47             :         return {nullptr,
      48           1 :                 "Failed to identify source layer from clipping dataset."};
      49             :     }
      50             : 
      51          16 :     if (!m_likeWhere.empty())
      52           2 :         poLyr->SetAttributeFilter(m_likeWhere.c_str());
      53             : 
      54          32 :     OGRGeometryCollection oGC;
      55          16 :     oGC.assignSpatialReference(poLyr->GetSpatialRef());
      56             : 
      57         382 :     for (auto &poFeat : poLyr)
      58             :     {
      59         367 :         auto poSrcGeom = std::unique_ptr<OGRGeometry>(poFeat->StealGeometry());
      60         367 :         if (poSrcGeom)
      61             :         {
      62             :             // Only take into account areal geometries.
      63         367 :             if (poSrcGeom->getDimension() == 2)
      64             :             {
      65         367 :                 if (!poSrcGeom->IsValid())
      66             :                 {
      67             :                     return {
      68             :                         nullptr,
      69           1 :                         CPLSPrintf("Geometry of feature " CPL_FRMT_GIB " of %s "
      70             :                                    "is invalid. You may be able to correct it "
      71             :                                    "with 'gdal vector geom make-valid'.",
      72           3 :                                    poFeat->GetFID(), poDS->GetDescription())};
      73             :                 }
      74             :                 else
      75             :                 {
      76         366 :                     oGC.addGeometry(std::move(poSrcGeom));
      77             :                 }
      78             :             }
      79             :             else
      80             :             {
      81           0 :                 CPLErrorOnce(CE_Warning, CPLE_AppDefined,
      82             :                              "Non-polygonal geometry encountered in clipping "
      83             :                              "dataset will be ignored.");
      84             :             }
      85             :         }
      86             :     }
      87             : 
      88          15 :     if (!m_likeSQL.empty())
      89           1 :         poDS->ReleaseResultSet(poLyr);
      90             : 
      91          15 :     if (oGC.IsEmpty())
      92             :     {
      93           1 :         return {nullptr, "No clipping geometry found"};
      94             :     }
      95             : 
      96          14 :     return {std::unique_ptr<OGRGeometry>(oGC.UnaryUnion()), std::string()};
      97             : }
      98             : 
      99             : /************************************************************************/
     100             : /*                          GetClipGeometry()                           */
     101             : /************************************************************************/
     102             : 
     103             : std::pair<std::unique_ptr<OGRGeometry>, std::string>
     104          78 : GDALClipCommon::GetClipGeometry()
     105             : {
     106             : 
     107          78 :     std::unique_ptr<OGRGeometry> poClipGeom;
     108             : 
     109          78 :     if (!m_bbox.empty())
     110             :     {
     111          44 :         poClipGeom = std::make_unique<OGRPolygon>(m_bbox[0], m_bbox[1],
     112          44 :                                                   m_bbox[2], m_bbox[3]);
     113             : 
     114          22 :         if (!m_bboxCrs.empty())
     115             :         {
     116           4 :             auto poSRS = OGRSpatialReferenceRefCountedPtr::makeInstance();
     117           2 :             poSRS->SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
     118           2 :             CPL_IGNORE_RET_VAL(poSRS->SetFromUserInput(m_bboxCrs.c_str()));
     119           2 :             poClipGeom->assignSpatialReference(poSRS.get());
     120             :         }
     121             :     }
     122          56 :     else if (!m_geometry.empty())
     123             :     {
     124             :         {
     125          42 :             CPLErrorStateBackuper oErrorStateBackuper(CPLQuietErrorHandler);
     126          21 :             auto [poGeom, eErr] =
     127          42 :                 OGRGeometryFactory::createFromWkt(m_geometry.c_str());
     128          21 :             if (eErr == OGRERR_NONE)
     129             :             {
     130          16 :                 poClipGeom = std::move(poGeom);
     131             :             }
     132             :             else
     133             :             {
     134           5 :                 poClipGeom.reset(
     135             :                     OGRGeometryFactory::createFromGeoJson(m_geometry.c_str()));
     136           5 :                 if (poClipGeom && poClipGeom->getSpatialReference() == nullptr)
     137             :                 {
     138             :                     auto poSRS =
     139           0 :                         OGRSpatialReferenceRefCountedPtr::makeInstance();
     140           0 :                     poSRS->SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
     141           0 :                     CPL_IGNORE_RET_VAL(poSRS->SetFromUserInput("WGS84"));
     142           0 :                     poClipGeom->assignSpatialReference(poSRS.get());
     143             :                 }
     144             :             }
     145             :         }
     146          21 :         if (!poClipGeom)
     147             :         {
     148             :             return {
     149             :                 nullptr,
     150           3 :                 "Clipping geometry is neither a valid WKT or GeoJSON geometry"};
     151             :         }
     152             : 
     153          18 :         if (!m_geometryCrs.empty())
     154             :         {
     155          10 :             auto poSRS = OGRSpatialReferenceRefCountedPtr::makeInstance();
     156           5 :             poSRS->SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
     157             :             // Validity of CRS already checked by GDALAlgorithm
     158           5 :             CPL_IGNORE_RET_VAL(poSRS->SetFromUserInput(m_geometryCrs.c_str()));
     159           5 :             poClipGeom->assignSpatialReference(poSRS.get());
     160             :         }
     161             :     }
     162          35 :     else if (auto poLikeDS = m_likeDataset.GetDatasetRef())
     163             :     {
     164          35 :         if (poLikeDS->GetLayerCount() > 1 && m_likeLayer.empty() &&
     165           2 :             m_likeSQL.empty())
     166             :         {
     167             :             return {
     168             :                 nullptr,
     169             :                 "Only single layer dataset can be specified with --like when "
     170           1 :                 "neither --like-layer or --like-sql have been specified"};
     171             :         }
     172          32 :         else if (poLikeDS->GetLayerCount() > 0)
     173             :         {
     174          17 :             std::string errMsg;
     175          17 :             std::tie(poClipGeom, errMsg) = LoadGeometry();
     176          17 :             if (!poClipGeom)
     177           3 :                 return {nullptr, errMsg};
     178             :         }
     179          15 :         else if (poLikeDS->GetRasterCount() > 0)
     180             :         {
     181          14 :             GDALGeoTransform gt;
     182          14 :             if (poLikeDS->GetGeoTransform(gt) != CE_None)
     183             :             {
     184             :                 return {
     185             :                     nullptr,
     186           2 :                     CPLSPrintf(
     187             :                         "Dataset '%s' has no geotransform matrix. Its bounds "
     188             :                         "cannot be established.",
     189           6 :                         poLikeDS->GetDescription())};
     190             :             }
     191          12 :             auto poLikeSRS = poLikeDS->GetSpatialRef();
     192          12 :             const double dfTLX = gt.xorig;
     193          12 :             const double dfTLY = gt.yorig;
     194             : 
     195          12 :             double dfTRX = 0;
     196          12 :             double dfTRY = 0;
     197          12 :             gt.Apply(poLikeDS->GetRasterXSize(), 0, &dfTRX, &dfTRY);
     198             : 
     199          12 :             double dfBLX = 0;
     200          12 :             double dfBLY = 0;
     201          12 :             gt.Apply(0, poLikeDS->GetRasterYSize(), &dfBLX, &dfBLY);
     202             : 
     203          12 :             double dfBRX = 0;
     204          12 :             double dfBRY = 0;
     205          12 :             gt.Apply(poLikeDS->GetRasterXSize(), poLikeDS->GetRasterYSize(),
     206             :                      &dfBRX, &dfBRY);
     207             : 
     208          24 :             auto poPoly = std::make_unique<OGRPolygon>();
     209          24 :             auto poLR = std::make_unique<OGRLinearRing>();
     210          12 :             poLR->addPoint(dfTLX, dfTLY);
     211          12 :             poLR->addPoint(dfTRX, dfTRY);
     212          12 :             poLR->addPoint(dfBRX, dfBRY);
     213          12 :             poLR->addPoint(dfBLX, dfBLY);
     214          12 :             poLR->addPoint(dfTLX, dfTLY);
     215          12 :             poPoly->addRingDirectly(poLR.release());
     216          12 :             poPoly->assignSpatialReference(poLikeSRS);
     217          12 :             poClipGeom = std::move(poPoly);
     218             :         }
     219             :         else
     220             :         {
     221           1 :             return {nullptr, "Cannot get extent from clip dataset"};
     222             :         }
     223             :     }
     224             :     else
     225             :     {
     226           2 :         return {nullptr, "--bbox, --geometry or --like must be specified"};
     227             :     }
     228             : 
     229          66 :     return {std::move(poClipGeom), std::string()};
     230             : }
     231             : 
     232             : //! @endcond