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          10 : GDALClipCommon::LoadGeometry()
      35             : {
      36          10 :     auto poDS = m_likeDataset.GetDatasetRef();
      37          10 :     OGRLayer *poLyr = nullptr;
      38          10 :     if (!m_likeSQL.empty())
      39           1 :         poLyr = poDS->ExecuteSQL(m_likeSQL.c_str(), nullptr, nullptr);
      40           9 :     else if (!m_likeLayer.empty())
      41           2 :         poLyr = poDS->GetLayerByName(m_likeLayer.c_str());
      42             :     else
      43           7 :         poLyr = poDS->GetLayer(0);
      44             : 
      45          10 :     if (poLyr == nullptr)
      46             :     {
      47             :         return {nullptr,
      48           1 :                 "Failed to identify source layer from clipping dataset."};
      49             :     }
      50             : 
      51           9 :     if (!m_likeWhere.empty())
      52           2 :         poLyr->SetAttributeFilter(m_likeWhere.c_str());
      53             : 
      54          18 :     OGRGeometryCollection oGC;
      55           9 :     oGC.assignSpatialReference(poLyr->GetSpatialRef());
      56             : 
      57          34 :     for (auto &poFeat : poLyr)
      58             :     {
      59          26 :         auto poSrcGeom = std::unique_ptr<OGRGeometry>(poFeat->StealGeometry());
      60          26 :         if (poSrcGeom)
      61             :         {
      62             :             // Only take into account areal geometries.
      63          26 :             if (poSrcGeom->getDimension() == 2)
      64             :             {
      65          26 :                 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          25 :                     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           8 :     if (!m_likeSQL.empty())
      89           1 :         poDS->ReleaseResultSet(poLyr);
      90             : 
      91           8 :     if (oGC.IsEmpty())
      92             :     {
      93           1 :         return {nullptr, "No clipping geometry found"};
      94             :     }
      95             : 
      96           7 :     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          54 : GDALClipCommon::GetClipGeometry()
     105             : {
     106             : 
     107          54 :     std::unique_ptr<OGRGeometry> poClipGeom;
     108             : 
     109          54 :     if (!m_bbox.empty())
     110             :     {
     111          36 :         poClipGeom = std::make_unique<OGRPolygon>(m_bbox[0], m_bbox[1],
     112          36 :                                                   m_bbox[2], m_bbox[3]);
     113             : 
     114          18 :         if (!m_bboxCrs.empty())
     115             :         {
     116           2 :             auto poSRS = new OGRSpatialReference();
     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);
     120           2 :             poSRS->Release();
     121             :         }
     122             :     }
     123          36 :     else if (!m_geometry.empty())
     124             :     {
     125             :         {
     126          30 :             CPLErrorStateBackuper oErrorStateBackuper(CPLQuietErrorHandler);
     127          15 :             auto [poGeom, eErr] =
     128          30 :                 OGRGeometryFactory::createFromWkt(m_geometry.c_str());
     129          15 :             if (eErr == OGRERR_NONE)
     130             :             {
     131          11 :                 poClipGeom = std::move(poGeom);
     132             :             }
     133             :             else
     134             :             {
     135           4 :                 poClipGeom.reset(
     136             :                     OGRGeometryFactory::createFromGeoJson(m_geometry.c_str()));
     137           4 :                 if (poClipGeom && poClipGeom->getSpatialReference() == nullptr)
     138             :                 {
     139           0 :                     auto poSRS = new OGRSpatialReference();
     140           0 :                     poSRS->SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
     141           0 :                     CPL_IGNORE_RET_VAL(poSRS->SetFromUserInput("WGS84"));
     142           0 :                     poClipGeom->assignSpatialReference(poSRS);
     143           0 :                     poSRS->Release();
     144             :                 }
     145             :             }
     146             :         }
     147          15 :         if (!poClipGeom)
     148             :         {
     149             :             return {
     150             :                 nullptr,
     151           2 :                 "Clipping geometry is neither a valid WKT or GeoJSON geometry"};
     152             :         }
     153             : 
     154          13 :         if (!m_geometryCrs.empty())
     155             :         {
     156           3 :             auto poSRS = new OGRSpatialReference();
     157           3 :             poSRS->SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
     158             :             // Validity of CRS already checked by GDALAlgorithm
     159           3 :             CPL_IGNORE_RET_VAL(poSRS->SetFromUserInput(m_geometryCrs.c_str()));
     160           3 :             poClipGeom->assignSpatialReference(poSRS);
     161           3 :             poSRS->Release();
     162             :         }
     163             :     }
     164          21 :     else if (auto poLikeDS = m_likeDataset.GetDatasetRef())
     165             :     {
     166          21 :         if (poLikeDS->GetLayerCount() > 1 && m_likeLayer.empty() &&
     167           2 :             m_likeSQL.empty())
     168             :         {
     169             :             return {
     170             :                 nullptr,
     171             :                 "Only single layer dataset can be specified with --like when "
     172           1 :                 "neither --like-layer or --like-sql have been specified"};
     173             :         }
     174          18 :         else if (poLikeDS->GetLayerCount() > 0)
     175             :         {
     176          10 :             std::string errMsg;
     177          10 :             std::tie(poClipGeom, errMsg) = LoadGeometry();
     178          10 :             if (!poClipGeom)
     179           3 :                 return {nullptr, errMsg};
     180             :         }
     181           8 :         else if (poLikeDS->GetRasterCount() > 0)
     182             :         {
     183             :             double adfGT[6];
     184           7 :             if (poLikeDS->GetGeoTransform(adfGT) != CE_None)
     185             :             {
     186             :                 return {
     187             :                     nullptr,
     188           2 :                     CPLSPrintf(
     189             :                         "Dataset '%s' has no geotransform matrix. Its bounds "
     190             :                         "cannot be established.",
     191           6 :                         poLikeDS->GetDescription())};
     192             :             }
     193           5 :             auto poLikeSRS = poLikeDS->GetSpatialRef();
     194           5 :             const double dfTLX = adfGT[0];
     195           5 :             const double dfTLY = adfGT[3];
     196             : 
     197           5 :             double dfTRX = 0;
     198           5 :             double dfTRY = 0;
     199           5 :             GDALApplyGeoTransform(adfGT, poLikeDS->GetRasterXSize(), 0, &dfTRX,
     200             :                                   &dfTRY);
     201             : 
     202           5 :             double dfBLX = 0;
     203           5 :             double dfBLY = 0;
     204           5 :             GDALApplyGeoTransform(adfGT, 0, poLikeDS->GetRasterYSize(), &dfBLX,
     205             :                                   &dfBLY);
     206             : 
     207           5 :             double dfBRX = 0;
     208           5 :             double dfBRY = 0;
     209           5 :             GDALApplyGeoTransform(adfGT, poLikeDS->GetRasterXSize(),
     210           5 :                                   poLikeDS->GetRasterYSize(), &dfBRX, &dfBRY);
     211             : 
     212          10 :             auto poPoly = std::make_unique<OGRPolygon>();
     213          10 :             auto poLR = std::make_unique<OGRLinearRing>();
     214           5 :             poLR->addPoint(dfTLX, dfTLY);
     215           5 :             poLR->addPoint(dfTRX, dfTRY);
     216           5 :             poLR->addPoint(dfBRX, dfBRY);
     217           5 :             poLR->addPoint(dfBLX, dfBLY);
     218           5 :             poLR->addPoint(dfTLX, dfTLY);
     219           5 :             poPoly->addRingDirectly(poLR.release());
     220           5 :             poPoly->assignSpatialReference(poLikeSRS);
     221           5 :             poClipGeom = std::move(poPoly);
     222             :         }
     223             :         else
     224             :         {
     225           1 :             return {nullptr, "Cannot get extent from clip dataset"};
     226             :         }
     227             :     }
     228             :     else
     229             :     {
     230           2 :         return {nullptr, "--bbox, --geometry or --like must be specified"};
     231             :     }
     232             : 
     233          43 :     return {std::move(poClipGeom), std::string()};
     234             : }
     235             : 
     236             : //! @endcond