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