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