Line data    Source code

       1             : /******************************************************************************
       2             :  *
       3             :  * Project:  GDAL
       4             :  * Purpose:  "clip" step of "raster pipeline", or "gdal raster 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_raster_clip.h"
      14             : 
      15             : #include "gdal_priv.h"
      16             : #include "gdal_utils.h"
      17             : 
      18             : #include <algorithm>
      19             : #include <cmath>
      20             : 
      21             : //! @cond Doxygen_Suppress
      22             : 
      23             : #ifndef _
      24             : #define _(x) (x)
      25             : #endif
      26             : 
      27             : /************************************************************************/
      28             : /*          GDALRasterClipAlgorithm::GDALRasterClipAlgorithm()          */
      29             : /************************************************************************/
      30             : 
      31          49 : GDALRasterClipAlgorithm::GDALRasterClipAlgorithm(bool standaloneStep)
      32             :     : GDALRasterPipelineStepAlgorithm(NAME, DESCRIPTION, HELP_URL,
      33          49 :                                       standaloneStep)
      34             : {
      35          49 :     constexpr const char *EXCLUSION_GROUP = "bbox-window-geometry-like";
      36          49 :     AddBBOXArg(&m_bbox, _("Clipping bounding box as xmin,ymin,xmax,ymax"))
      37          49 :         .SetMutualExclusionGroup(EXCLUSION_GROUP);
      38          98 :     AddArg("bbox-crs", 0, _("CRS of clipping bounding box"), &m_bboxCrs)
      39          98 :         .SetIsCRSArg()
      40          49 :         .AddHiddenAlias("bbox_srs");
      41             : 
      42             :     AddArg("window", 0, _("Raster window as col,line,width,height in pixels"),
      43          98 :            &m_window)
      44          49 :         .SetRepeatedArgAllowed(false)
      45          49 :         .SetMinCount(4)
      46          49 :         .SetMaxCount(4)
      47          49 :         .SetDisplayHintAboutRepetition(false)
      48          98 :         .SetMutualExclusionGroup(EXCLUSION_GROUP)
      49             :         .AddValidationAction(
      50          20 :             [this]()
      51             :             {
      52           7 :                 CPLAssert(m_window.size() == 4);
      53           7 :                 if (m_window[2] <= 0 || m_window[3] <= 0)
      54             :                 {
      55           2 :                     CPLError(CE_Failure, CPLE_AppDefined,
      56             :                              "Value of 'window' should be "
      57             :                              "col,line,width,height with "
      58             :                              "width > 0 and height > 0");
      59           2 :                     return false;
      60             :                 }
      61           5 :                 return true;
      62          49 :             });
      63             : 
      64          98 :     AddArg("geometry", 0, _("Clipping geometry (WKT or GeoJSON)"), &m_geometry)
      65          49 :         .SetMutualExclusionGroup(EXCLUSION_GROUP);
      66          98 :     AddArg("geometry-crs", 0, _("CRS of clipping geometry"), &m_geometryCrs)
      67          98 :         .SetIsCRSArg()
      68          49 :         .AddHiddenAlias("geometry_srs");
      69             :     AddArg("like", 0, _("Dataset to use as a template for bounds"),
      70          98 :            &m_likeDataset, GDAL_OF_RASTER | GDAL_OF_VECTOR)
      71          98 :         .SetMetaVar("DATASET")
      72          49 :         .SetMutualExclusionGroup(EXCLUSION_GROUP);
      73             :     AddArg("like-sql", 0, ("SELECT statement to run on the 'like' dataset"),
      74          98 :            &m_likeSQL)
      75          98 :         .SetMetaVar("SELECT-STATEMENT")
      76          49 :         .SetMutualExclusionGroup("sql-where");
      77             :     AddArg("like-layer", 0, ("Name of the layer of the 'like' dataset"),
      78          98 :            &m_likeLayer)
      79          49 :         .SetMetaVar("LAYER-NAME");
      80             :     AddArg("like-where", 0, ("WHERE SQL clause to run on the 'like' dataset"),
      81          98 :            &m_likeWhere)
      82          98 :         .SetMetaVar("WHERE-EXPRESSION")
      83          49 :         .SetMutualExclusionGroup("sql-where");
      84             :     AddArg("only-bbox", 0,
      85             :            _("For 'geometry' and 'like', only consider their bounding box"),
      86          49 :            &m_onlyBBOX);
      87             :     AddArg("allow-bbox-outside-source", 0,
      88             :            _("Allow clipping box to include pixels outside input dataset"),
      89          49 :            &m_allowExtentOutsideSource);
      90             :     AddArg("add-alpha", 0,
      91             :            _("Adds an alpha mask band to the destination when the source "
      92             :              "raster have none."),
      93          49 :            &m_addAlpha);
      94          49 : }
      95             : 
      96             : /************************************************************************/
      97             : /*                 GDALRasterClipAlgorithm::RunStep()                   */
      98             : /************************************************************************/
      99             : 
     100          24 : bool GDALRasterClipAlgorithm::RunStep(GDALPipelineStepRunContext &)
     101             : {
     102          24 :     auto poSrcDS = m_inputDataset[0].GetDatasetRef();
     103          24 :     CPLAssert(poSrcDS);
     104          24 :     CPLAssert(m_outputDataset.GetName().empty());
     105          24 :     CPLAssert(!m_outputDataset.GetDatasetRef());
     106             : 
     107          24 :     if (!m_window.empty())
     108             :     {
     109           2 :         if (m_addAlpha)
     110             :         {
     111           1 :             ReportError(CE_Failure, CPLE_NotSupported,
     112             :                         "'alpha' argument is not supported with 'window'");
     113           1 :             return false;
     114             :         }
     115             : 
     116           1 :         CPLStringList aosOptions;
     117           1 :         aosOptions.AddString("-of");
     118           1 :         aosOptions.AddString("VRT");
     119             : 
     120           1 :         aosOptions.AddString("-srcwin");
     121           1 :         aosOptions.AddString(CPLSPrintf("%d", m_window[0]));
     122           1 :         aosOptions.AddString(CPLSPrintf("%d", m_window[1]));
     123           1 :         aosOptions.AddString(CPLSPrintf("%d", m_window[2]));
     124           1 :         aosOptions.AddString(CPLSPrintf("%d", m_window[3]));
     125             : 
     126           1 :         if (!m_allowExtentOutsideSource)
     127             :         {
     128             :             // Unless we've specifically allowed the bounding box to extend beyond
     129             :             // the source raster, raise an error.
     130           1 :             aosOptions.AddString("-epo");
     131             :         }
     132             : 
     133             :         GDALTranslateOptions *psOptions =
     134           1 :             GDALTranslateOptionsNew(aosOptions.List(), nullptr);
     135             : 
     136           1 :         GDALDatasetH hSrcDS = GDALDataset::ToHandle(poSrcDS);
     137           1 :         auto poRetDS = GDALDataset::FromHandle(
     138             :             GDALTranslate("", hSrcDS, psOptions, nullptr));
     139           1 :         GDALTranslateOptionsFree(psOptions);
     140             : 
     141           1 :         const bool bOK = poRetDS != nullptr;
     142           1 :         if (bOK)
     143             :         {
     144           1 :             m_outputDataset.Set(std::unique_ptr<GDALDataset>(poRetDS));
     145             :         }
     146             : 
     147           1 :         return bOK;
     148             :     }
     149             : 
     150          22 :     GDALGeoTransform gt;
     151          22 :     if (poSrcDS->GetGeoTransform(gt) != CE_None)
     152             :     {
     153           1 :         ReportError(
     154             :             CE_Failure, CPLE_NotSupported,
     155             :             "Clipping is not supported on a raster without a geotransform");
     156           1 :         return false;
     157             :     }
     158          21 :     if (gt[2] != 0 && gt[4] != 0)
     159             :     {
     160           1 :         ReportError(CE_Failure, CPLE_NotSupported,
     161             :                     "Clipping is not supported on a raster whose geotransform "
     162             :                     "has rotation terms");
     163           1 :         return false;
     164             :     }
     165             : 
     166          40 :     auto [poClipGeom, errMsg] = GetClipGeometry();
     167          20 :     if (!poClipGeom)
     168             :     {
     169           3 :         ReportError(CE_Failure, CPLE_AppDefined, "%s", errMsg.c_str());
     170           3 :         return false;
     171             :     }
     172             : 
     173          17 :     auto poLikeDS = m_likeDataset.GetDatasetRef();
     174          18 :     if (!poClipGeom->getSpatialReference() && poLikeDS &&
     175           1 :         poLikeDS->GetLayerCount() == 0)
     176             :     {
     177           1 :         ReportError(CE_Failure, CPLE_AppDefined,
     178             :                     "Dataset '%s' has no CRS. Its bounds cannot be used.",
     179           1 :                     poLikeDS->GetDescription());
     180           1 :         return false;
     181             :     }
     182             : 
     183          32 :     CPLStringList aosOptions;
     184          16 :     aosOptions.AddString("-of");
     185          16 :     aosOptions.AddString("VRT");
     186             : 
     187          16 :     OGREnvelope env;
     188          16 :     poClipGeom->getEnvelope(&env);
     189             : 
     190          16 :     if (m_onlyBBOX)
     191             :     {
     192           2 :         auto poPoly = std::make_unique<OGRPolygon>(env);
     193           1 :         poPoly->assignSpatialReference(poClipGeom->getSpatialReference());
     194           1 :         poClipGeom = std::move(poPoly);
     195             :     }
     196             : 
     197          16 :     const bool bBottomUpRaster = gt[5] > 0;
     198             : 
     199          16 :     if (poClipGeom->IsRectangle() && !m_addAlpha && !bBottomUpRaster)
     200             :     {
     201           9 :         aosOptions.AddString("-projwin");
     202           9 :         aosOptions.AddString(CPLSPrintf("%.17g", env.MinX));
     203           9 :         aosOptions.AddString(CPLSPrintf("%.17g", env.MaxY));
     204           9 :         aosOptions.AddString(CPLSPrintf("%.17g", env.MaxX));
     205           9 :         aosOptions.AddString(CPLSPrintf("%.17g", env.MinY));
     206             : 
     207           9 :         auto poClipGeomSRS = poClipGeom->getSpatialReference();
     208           9 :         if (poClipGeomSRS)
     209             :         {
     210           2 :             const char *const apszOptions[] = {"FORMAT=WKT2", nullptr};
     211           4 :             const std::string osWKT = poClipGeomSRS->exportToWkt(apszOptions);
     212           2 :             aosOptions.AddString("-projwin_srs");
     213           2 :             aosOptions.AddString(osWKT.c_str());
     214             :         }
     215             : 
     216           9 :         if (!m_allowExtentOutsideSource)
     217             :         {
     218             :             // Unless we've specifically allowed the bounding box to extend beyond
     219             :             // the source raster, raise an error.
     220           6 :             aosOptions.AddString("-epo");
     221             :         }
     222             : 
     223             :         GDALTranslateOptions *psOptions =
     224           9 :             GDALTranslateOptionsNew(aosOptions.List(), nullptr);
     225             : 
     226           9 :         GDALDatasetH hSrcDS = GDALDataset::ToHandle(poSrcDS);
     227           9 :         auto poRetDS = GDALDataset::FromHandle(
     228             :             GDALTranslate("", hSrcDS, psOptions, nullptr));
     229           9 :         GDALTranslateOptionsFree(psOptions);
     230             : 
     231           9 :         const bool bOK = poRetDS != nullptr;
     232           9 :         if (bOK)
     233             :         {
     234           7 :             m_outputDataset.Set(std::unique_ptr<GDALDataset>(poRetDS));
     235             :         }
     236             : 
     237           9 :         return bOK;
     238             :     }
     239             :     else
     240             :     {
     241           7 :         if (bBottomUpRaster)
     242             :         {
     243           1 :             gt[3] += gt[5] * poSrcDS->GetRasterYSize();
     244           1 :             gt[5] = -gt[5];
     245             :         }
     246             : 
     247             :         {
     248             :             auto poClipGeomInSrcSRS =
     249          14 :                 std::unique_ptr<OGRGeometry>(poClipGeom->clone());
     250           7 :             if (poClipGeom->getSpatialReference() && poSrcDS->GetSpatialRef())
     251           1 :                 poClipGeomInSrcSRS->transformTo(poSrcDS->GetSpatialRef());
     252           7 :             poClipGeomInSrcSRS->getEnvelope(&env);
     253             :         }
     254             : 
     255           7 :         OGREnvelope rasterEnv;
     256           7 :         poSrcDS->GetExtent(&rasterEnv, nullptr);
     257           7 :         if (!m_allowExtentOutsideSource && !rasterEnv.Contains(env))
     258             :         {
     259           1 :             ReportError(CE_Failure, CPLE_AppDefined,
     260             :                         "Clipping geometry is partially or totally outside the "
     261             :                         "extent of the raster. You can set the "
     262             :                         "'allow-bbox-outside-source' argument to proceed.");
     263           1 :             return false;
     264             :         }
     265             : 
     266           6 :         if (m_addAlpha)
     267             :         {
     268           2 :             aosOptions.AddString("-dstalpha");
     269             :         }
     270             : 
     271           6 :         aosOptions.AddString("-cutline");
     272           6 :         aosOptions.AddString(poClipGeom->exportToWkt());
     273             : 
     274           6 :         aosOptions.AddString("-wo");
     275           6 :         aosOptions.AddString("CUTLINE_ALL_TOUCHED=YES");
     276             : 
     277           6 :         auto poClipGeomSRS = poClipGeom->getSpatialReference();
     278           6 :         if (poClipGeomSRS)
     279             :         {
     280           1 :             const char *const apszOptions[] = {"FORMAT=WKT2", nullptr};
     281           2 :             const std::string osWKT = poClipGeomSRS->exportToWkt(apszOptions);
     282           1 :             aosOptions.AddString("-cutline_srs");
     283           1 :             aosOptions.AddString(osWKT.c_str());
     284             :         }
     285             : 
     286           6 :         constexpr double REL_EPS_PIXEL = 1e-3;
     287             :         const double dfMinX =
     288           6 :             gt[0] + floor((env.MinX - gt[0]) / gt[1] + REL_EPS_PIXEL) * gt[1];
     289             :         const double dfMinY =
     290           6 :             gt[3] + ceil((env.MinY - gt[3]) / gt[5] - REL_EPS_PIXEL) * gt[5];
     291             :         const double dfMaxX =
     292           6 :             gt[0] + ceil((env.MaxX - gt[0]) / gt[1] - REL_EPS_PIXEL) * gt[1];
     293             :         const double dfMaxY =
     294           6 :             gt[3] + floor((env.MaxY - gt[3]) / gt[5] + REL_EPS_PIXEL) * gt[5];
     295             : 
     296           6 :         aosOptions.AddString("-te");
     297           6 :         aosOptions.AddString(CPLSPrintf("%.17g", dfMinX));
     298             :         aosOptions.AddString(
     299           6 :             CPLSPrintf("%.17g", bBottomUpRaster ? dfMaxY : dfMinY));
     300           6 :         aosOptions.AddString(CPLSPrintf("%.17g", dfMaxX));
     301             :         aosOptions.AddString(
     302           6 :             CPLSPrintf("%.17g", bBottomUpRaster ? dfMinY : dfMaxY));
     303             : 
     304           6 :         aosOptions.AddString("-tr");
     305           6 :         aosOptions.AddString(CPLSPrintf("%.17g", gt[1]));
     306           6 :         aosOptions.AddString(CPLSPrintf("%.17g", std::fabs(gt[5])));
     307             : 
     308             :         GDALWarpAppOptions *psOptions =
     309           6 :             GDALWarpAppOptionsNew(aosOptions.List(), nullptr);
     310             : 
     311           6 :         GDALDatasetH hSrcDS = GDALDataset::ToHandle(poSrcDS);
     312           6 :         auto poRetDS = GDALDataset::FromHandle(
     313             :             GDALWarp("", nullptr, 1, &hSrcDS, psOptions, nullptr));
     314           6 :         GDALWarpAppOptionsFree(psOptions);
     315             : 
     316           6 :         const bool bOK = poRetDS != nullptr;
     317           6 :         if (bOK)
     318             :         {
     319           6 :             m_outputDataset.Set(std::unique_ptr<GDALDataset>(poRetDS));
     320             :         }
     321             : 
     322           6 :         return bOK;
     323             :     }
     324             : }
     325             : 
     326             : GDALRasterClipAlgorithmStandalone::~GDALRasterClipAlgorithmStandalone() =
     327             :     default;
     328             : 
     329             : //! @endcond