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          66 : GDALRasterClipAlgorithm::GDALRasterClipAlgorithm(bool standaloneStep)
      32             :     : GDALRasterPipelineStepAlgorithm(NAME, DESCRIPTION, HELP_URL,
      33          66 :                                       standaloneStep)
      34             : {
      35          66 :     constexpr const char *EXCLUSION_GROUP = "bbox-window-geometry-like";
      36          66 :     AddBBOXArg(&m_bbox, _("Clipping bounding box as xmin,ymin,xmax,ymax"))
      37          66 :         .SetMutualExclusionGroup(EXCLUSION_GROUP);
      38         132 :     AddArg("bbox-crs", 0, _("CRS of clipping bounding box"), &m_bboxCrs)
      39         132 :         .SetIsCRSArg()
      40          66 :         .AddHiddenAlias("bbox_srs");
      41             : 
      42             :     AddArg("window", 0, _("Raster window as col,line,width,height in pixels"),
      43         132 :            &m_window)
      44          66 :         .SetRepeatedArgAllowed(false)
      45          66 :         .SetMinCount(4)
      46          66 :         .SetMaxCount(4)
      47          66 :         .SetDisplayHintAboutRepetition(false)
      48         132 :         .SetMutualExclusionGroup(EXCLUSION_GROUP)
      49             :         .AddValidationAction(
      50          32 :             [this]()
      51             :             {
      52          11 :                 CPLAssert(m_window.size() == 4);
      53          11 :                 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           9 :                 return true;
      62          66 :             });
      63             : 
      64         132 :     AddArg("geometry", 0, _("Clipping geometry (WKT or GeoJSON)"), &m_geometry)
      65          66 :         .SetMutualExclusionGroup(EXCLUSION_GROUP);
      66         132 :     AddArg("geometry-crs", 0, _("CRS of clipping geometry"), &m_geometryCrs)
      67         132 :         .SetIsCRSArg()
      68          66 :         .AddHiddenAlias("geometry_srs");
      69             :     AddArg("like", 0, _("Dataset to use as a template for bounds"),
      70         132 :            &m_likeDataset, GDAL_OF_RASTER | GDAL_OF_VECTOR)
      71         132 :         .SetMetaVar("DATASET")
      72          66 :         .SetMutualExclusionGroup(EXCLUSION_GROUP);
      73             :     AddArg("like-sql", 0, ("SELECT statement to run on the 'like' dataset"),
      74         132 :            &m_likeSQL)
      75         132 :         .SetMetaVar("SELECT-STATEMENT")
      76          66 :         .SetMutualExclusionGroup("sql-where");
      77             :     AddArg("like-layer", 0, ("Name of the layer of the 'like' dataset"),
      78         132 :            &m_likeLayer)
      79          66 :         .SetMetaVar("LAYER-NAME");
      80             :     AddArg("like-where", 0, ("WHERE SQL clause to run on the 'like' dataset"),
      81         132 :            &m_likeWhere)
      82         132 :         .SetMetaVar("WHERE-EXPRESSION")
      83          66 :         .SetMutualExclusionGroup("sql-where");
      84             :     AddArg("only-bbox", 0,
      85             :            _("For 'geometry' and 'like', only consider their bounding box"),
      86          66 :            &m_onlyBBOX);
      87             :     AddArg("allow-bbox-outside-source", 0,
      88             :            _("Allow clipping box to include pixels outside input dataset"),
      89          66 :            &m_allowExtentOutsideSource);
      90             :     AddArg("add-alpha", 0,
      91             :            _("Adds an alpha mask band to the destination when the source "
      92             :              "raster have none."),
      93          66 :            &m_addAlpha);
      94          66 : }
      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           3 :             aosOptions.AddString("--no-warn-about-outside-window");
     219             :         }
     220             :         else
     221             :         {
     222             :             // Unless we've specifically allowed the bounding box to extend beyond
     223             :             // the source raster, raise an error.
     224           6 :             aosOptions.AddString("-epo");
     225             :         }
     226             : 
     227             :         GDALTranslateOptions *psOptions =
     228           9 :             GDALTranslateOptionsNew(aosOptions.List(), nullptr);
     229             : 
     230           9 :         GDALDatasetH hSrcDS = GDALDataset::ToHandle(poSrcDS);
     231           9 :         auto poRetDS = GDALDataset::FromHandle(
     232             :             GDALTranslate("", hSrcDS, psOptions, nullptr));
     233           9 :         GDALTranslateOptionsFree(psOptions);
     234             : 
     235           9 :         const bool bOK = poRetDS != nullptr;
     236           9 :         if (bOK)
     237             :         {
     238           7 :             m_outputDataset.Set(std::unique_ptr<GDALDataset>(poRetDS));
     239             :         }
     240             : 
     241           9 :         return bOK;
     242             :     }
     243             :     else
     244             :     {
     245           7 :         if (bBottomUpRaster)
     246             :         {
     247           1 :             gt[3] += gt[5] * poSrcDS->GetRasterYSize();
     248           1 :             gt[5] = -gt[5];
     249             :         }
     250             : 
     251             :         {
     252             :             auto poClipGeomInSrcSRS =
     253          14 :                 std::unique_ptr<OGRGeometry>(poClipGeom->clone());
     254           7 :             if (poClipGeom->getSpatialReference() && poSrcDS->GetSpatialRef())
     255           1 :                 poClipGeomInSrcSRS->transformTo(poSrcDS->GetSpatialRef());
     256           7 :             poClipGeomInSrcSRS->getEnvelope(&env);
     257             :         }
     258             : 
     259           7 :         OGREnvelope rasterEnv;
     260           7 :         poSrcDS->GetExtent(&rasterEnv, nullptr);
     261           7 :         if (!m_allowExtentOutsideSource && !rasterEnv.Contains(env))
     262             :         {
     263           1 :             ReportError(CE_Failure, CPLE_AppDefined,
     264             :                         "Clipping geometry is partially or totally outside the "
     265             :                         "extent of the raster. You can set the "
     266             :                         "'allow-bbox-outside-source' argument to proceed.");
     267           1 :             return false;
     268             :         }
     269             : 
     270           6 :         if (m_addAlpha)
     271             :         {
     272           2 :             aosOptions.AddString("-dstalpha");
     273             :         }
     274             : 
     275           6 :         aosOptions.AddString("-cutline");
     276           6 :         aosOptions.AddString(poClipGeom->exportToWkt());
     277             : 
     278           6 :         aosOptions.AddString("-wo");
     279           6 :         aosOptions.AddString("CUTLINE_ALL_TOUCHED=YES");
     280             : 
     281           6 :         auto poClipGeomSRS = poClipGeom->getSpatialReference();
     282           6 :         if (poClipGeomSRS)
     283             :         {
     284           1 :             const char *const apszOptions[] = {"FORMAT=WKT2", nullptr};
     285           2 :             const std::string osWKT = poClipGeomSRS->exportToWkt(apszOptions);
     286           1 :             aosOptions.AddString("-cutline_srs");
     287           1 :             aosOptions.AddString(osWKT.c_str());
     288             :         }
     289             : 
     290           6 :         constexpr double REL_EPS_PIXEL = 1e-3;
     291             :         const double dfMinX =
     292           6 :             gt[0] + floor((env.MinX - gt[0]) / gt[1] + REL_EPS_PIXEL) * gt[1];
     293             :         const double dfMinY =
     294           6 :             gt[3] + ceil((env.MinY - gt[3]) / gt[5] - REL_EPS_PIXEL) * gt[5];
     295             :         const double dfMaxX =
     296           6 :             gt[0] + ceil((env.MaxX - gt[0]) / gt[1] - REL_EPS_PIXEL) * gt[1];
     297             :         const double dfMaxY =
     298           6 :             gt[3] + floor((env.MaxY - gt[3]) / gt[5] + REL_EPS_PIXEL) * gt[5];
     299             : 
     300           6 :         aosOptions.AddString("-te");
     301           6 :         aosOptions.AddString(CPLSPrintf("%.17g", dfMinX));
     302             :         aosOptions.AddString(
     303           6 :             CPLSPrintf("%.17g", bBottomUpRaster ? dfMaxY : dfMinY));
     304           6 :         aosOptions.AddString(CPLSPrintf("%.17g", dfMaxX));
     305             :         aosOptions.AddString(
     306           6 :             CPLSPrintf("%.17g", bBottomUpRaster ? dfMinY : dfMaxY));
     307             : 
     308           6 :         aosOptions.AddString("-tr");
     309           6 :         aosOptions.AddString(CPLSPrintf("%.17g", gt[1]));
     310           6 :         aosOptions.AddString(CPLSPrintf("%.17g", std::fabs(gt[5])));
     311             : 
     312             :         GDALWarpAppOptions *psOptions =
     313           6 :             GDALWarpAppOptionsNew(aosOptions.List(), nullptr);
     314             : 
     315           6 :         GDALDatasetH hSrcDS = GDALDataset::ToHandle(poSrcDS);
     316           6 :         auto poRetDS = GDALDataset::FromHandle(
     317             :             GDALWarp("", nullptr, 1, &hSrcDS, psOptions, nullptr));
     318           6 :         GDALWarpAppOptionsFree(psOptions);
     319             : 
     320           6 :         const bool bOK = poRetDS != nullptr;
     321           6 :         if (bOK)
     322             :         {
     323           6 :             m_outputDataset.Set(std::unique_ptr<GDALDataset>(poRetDS));
     324             :         }
     325             : 
     326           6 :         return bOK;
     327             :     }
     328             : }
     329             : 
     330             : GDALRasterClipAlgorithmStandalone::~GDALRasterClipAlgorithmStandalone() =
     331             :     default;
     332             : 
     333             : //! @endcond