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         129 : GDALRasterClipAlgorithm::GDALRasterClipAlgorithm(bool standaloneStep)
      32             :     : GDALRasterPipelineStepAlgorithm(NAME, DESCRIPTION, HELP_URL,
      33         129 :                                       standaloneStep)
      34             : {
      35         129 :     constexpr const char *EXCLUSION_GROUP = "bbox-window-geometry-like";
      36         129 :     AddBBOXArg(&m_bbox, _("Clipping bounding box as xmin,ymin,xmax,ymax"))
      37         129 :         .SetMutualExclusionGroup(EXCLUSION_GROUP);
      38         258 :     AddArg("bbox-crs", 0, _("CRS of clipping bounding box"), &m_bboxCrs)
      39         258 :         .SetIsCRSArg()
      40         129 :         .AddHiddenAlias("bbox_srs");
      41             : 
      42             :     AddArg("window", 0, _("Raster window as col,line,width,height in pixels"),
      43         258 :            &m_window)
      44         129 :         .SetRepeatedArgAllowed(false)
      45         129 :         .SetMinCount(4)
      46         129 :         .SetMaxCount(4)
      47         129 :         .SetDisplayHintAboutRepetition(false)
      48         258 :         .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         129 :             });
      63             : 
      64         258 :     AddArg("geometry", 0, _("Clipping geometry (WKT or GeoJSON)"), &m_geometry)
      65         129 :         .SetMutualExclusionGroup(EXCLUSION_GROUP);
      66         258 :     AddArg("geometry-crs", 0, _("CRS of clipping geometry"), &m_geometryCrs)
      67         258 :         .SetIsCRSArg()
      68         129 :         .AddHiddenAlias("geometry_srs");
      69             :     AddArg("like", 0, _("Dataset to use as a template for bounds"),
      70         258 :            &m_likeDataset, GDAL_OF_RASTER | GDAL_OF_VECTOR)
      71         258 :         .SetMetaVar("DATASET")
      72         129 :         .SetMutualExclusionGroup(EXCLUSION_GROUP);
      73             :     AddArg("like-sql", 0, ("SELECT statement to run on the 'like' dataset"),
      74         258 :            &m_likeSQL)
      75         258 :         .SetMetaVar("SELECT-STATEMENT")
      76         129 :         .SetMutualExclusionGroup("sql-where");
      77             :     AddArg("like-layer", 0, ("Name of the layer of the 'like' dataset"),
      78         258 :            &m_likeLayer)
      79         129 :         .SetMetaVar("LAYER-NAME");
      80             :     AddArg("like-where", 0, ("WHERE SQL clause to run on the 'like' dataset"),
      81         258 :            &m_likeWhere)
      82         258 :         .SetMetaVar("WHERE-EXPRESSION")
      83         129 :         .SetMutualExclusionGroup("sql-where");
      84             :     AddArg("only-bbox", 0,
      85             :            _("For 'geometry' and 'like', only consider their bounding box"),
      86         129 :            &m_onlyBBOX);
      87             :     AddArg("allow-bbox-outside-source", 0,
      88             :            _("Allow clipping box to include pixels outside input dataset"),
      89         129 :            &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         129 :            &m_addAlpha);
      96         129 : }
      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             : 
     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          32 :     GDALGeoTransform gt;
     153          32 :     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          31 :     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          60 :     auto [poClipGeom, errMsg] = GetClipGeometry();
     169          30 :     if (!poClipGeom)
     170             :     {
     171           3 :         ReportError(CE_Failure, CPLE_AppDefined, "%s", errMsg.c_str());
     172           3 :         return false;
     173             :     }
     174             : 
     175          27 :     auto poLikeDS = m_likeDataset.GetDatasetRef();
     176          28 :     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          52 :     CPLStringList aosOptions;
     186          26 :     aosOptions.AddString("-of");
     187          26 :     aosOptions.AddString("VRT");
     188             : 
     189          26 :     OGREnvelope env;
     190          26 :     poClipGeom->getEnvelope(&env);
     191             : 
     192          26 :     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          26 :     const bool bBottomUpRaster = gt.yscale > 0;
     200             : 
     201          26 :     if (poClipGeom->IsRectangle() && !m_addAlpha && !bBottomUpRaster)
     202             :     {
     203          13 :         aosOptions.AddString("-projwin");
     204          13 :         aosOptions.AddString(CPLSPrintf("%.17g", env.MinX));
     205          13 :         aosOptions.AddString(CPLSPrintf("%.17g", env.MaxY));
     206          13 :         aosOptions.AddString(CPLSPrintf("%.17g", env.MaxX));
     207          13 :         aosOptions.AddString(CPLSPrintf("%.17g", env.MinY));
     208             : 
     209          13 :         auto poClipGeomSRS = poClipGeom->getSpatialReference();
     210          13 :         if (poClipGeomSRS)
     211             :         {
     212           4 :             const char *const apszOptions[] = {"FORMAT=WKT2", nullptr};
     213           8 :             const std::string osWKT = poClipGeomSRS->exportToWkt(apszOptions);
     214           4 :             aosOptions.AddString("-projwin_srs");
     215           4 :             aosOptions.AddString(osWKT.c_str());
     216             :         }
     217             : 
     218          13 :         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          10 :             aosOptions.AddString("-epo");
     227             :         }
     228             : 
     229             :         GDALTranslateOptions *psOptions =
     230          13 :             GDALTranslateOptionsNew(aosOptions.List(), nullptr);
     231             : 
     232          13 :         GDALDatasetH hSrcDS = GDALDataset::ToHandle(poSrcDS);
     233          13 :         auto poRetDS = GDALDataset::FromHandle(
     234             :             GDALTranslate("", hSrcDS, psOptions, nullptr));
     235          13 :         GDALTranslateOptionsFree(psOptions);
     236             : 
     237          13 :         const bool bOK = poRetDS != nullptr;
     238          13 :         if (bOK)
     239             :         {
     240          11 :             m_outputDataset.Set(std::unique_ptr<GDALDataset>(poRetDS));
     241             :         }
     242             : 
     243          13 :         return bOK;
     244             :     }
     245             :     else
     246             :     {
     247          13 :         if (bBottomUpRaster)
     248             :         {
     249           1 :             gt.yorig += gt.yscale * poSrcDS->GetRasterYSize();
     250           1 :             gt.yscale = -gt.yscale;
     251             :         }
     252             : 
     253             :         {
     254             :             auto poClipGeomInSrcSRS =
     255          26 :                 std::unique_ptr<OGRGeometry>(poClipGeom->clone());
     256          13 :             if (poClipGeom->getSpatialReference() && poSrcDS->GetSpatialRef())
     257           7 :                 poClipGeomInSrcSRS->transformTo(poSrcDS->GetSpatialRef());
     258          13 :             poClipGeomInSrcSRS->getEnvelope(&env);
     259             :         }
     260             : 
     261          13 :         OGREnvelope rasterEnv;
     262          13 :         poSrcDS->GetExtent(&rasterEnv, nullptr);
     263          13 :         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          12 :         if (m_addAlpha)
     273             :         {
     274           2 :             aosOptions.AddString("-dstalpha");
     275             :         }
     276             : 
     277          12 :         aosOptions.AddString("-cutline");
     278          12 :         aosOptions.AddString(poClipGeom->exportToWkt());
     279             : 
     280          12 :         aosOptions.AddString("-wo");
     281          12 :         aosOptions.AddString("CUTLINE_ALL_TOUCHED=YES");
     282             : 
     283          12 :         auto poClipGeomSRS = poClipGeom->getSpatialReference();
     284          12 :         if (poClipGeomSRS)
     285             :         {
     286           7 :             const char *const apszOptions[] = {"FORMAT=WKT2", nullptr};
     287          14 :             const std::string osWKT = poClipGeomSRS->exportToWkt(apszOptions);
     288           7 :             aosOptions.AddString("-cutline_srs");
     289           7 :             aosOptions.AddString(osWKT.c_str());
     290             :         }
     291             : 
     292          12 :         constexpr double REL_EPS_PIXEL = 1e-3;
     293          12 :         const double dfMinX =
     294          12 :             gt.xorig +
     295          12 :             floor((env.MinX - gt.xorig) / gt.xscale + REL_EPS_PIXEL) *
     296          12 :                 gt.xscale;
     297          12 :         const double dfMinY =
     298          12 :             gt.yorig +
     299          12 :             ceil((env.MinY - gt.yorig) / gt.yscale - REL_EPS_PIXEL) * gt.yscale;
     300          12 :         const double dfMaxX =
     301          12 :             gt.xorig +
     302          12 :             ceil((env.MaxX - gt.xorig) / gt.xscale - REL_EPS_PIXEL) * gt.xscale;
     303          12 :         const double dfMaxY =
     304          12 :             gt.yorig +
     305          12 :             floor((env.MaxY - gt.yorig) / gt.yscale + REL_EPS_PIXEL) *
     306          12 :                 gt.yscale;
     307             : 
     308          12 :         aosOptions.AddString("-te");
     309          12 :         aosOptions.AddString(CPLSPrintf("%.17g", dfMinX));
     310             :         aosOptions.AddString(
     311          12 :             CPLSPrintf("%.17g", bBottomUpRaster ? dfMaxY : dfMinY));
     312          12 :         aosOptions.AddString(CPLSPrintf("%.17g", dfMaxX));
     313             :         aosOptions.AddString(
     314          12 :             CPLSPrintf("%.17g", bBottomUpRaster ? dfMinY : dfMaxY));
     315             : 
     316          12 :         aosOptions.AddString("-tr");
     317          12 :         aosOptions.AddString(CPLSPrintf("%.17g", gt.xscale));
     318          12 :         aosOptions.AddString(CPLSPrintf("%.17g", 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