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          29 : GDALRasterClipAlgorithm::GDALRasterClipAlgorithm(bool standaloneStep)
      32             :     : GDALRasterPipelineStepAlgorithm(NAME, DESCRIPTION, HELP_URL,
      33          29 :                                       standaloneStep)
      34             : {
      35          29 :     AddBBOXArg(&m_bbox, _("Clipping bounding box as xmin,ymin,xmax,ymax"))
      36          29 :         .SetMutualExclusionGroup("bbox-geometry-like");
      37          58 :     AddArg("bbox-crs", 0, _("CRS of clipping bounding box"), &m_bboxCrs)
      38          58 :         .SetIsCRSArg()
      39          29 :         .AddHiddenAlias("bbox_srs");
      40          58 :     AddArg("geometry", 0, _("Clipping geometry (WKT or GeoJSON)"), &m_geometry)
      41          29 :         .SetMutualExclusionGroup("bbox-geometry-like");
      42          58 :     AddArg("geometry-crs", 0, _("CRS of clipping geometry"), &m_geometryCrs)
      43          58 :         .SetIsCRSArg()
      44          29 :         .AddHiddenAlias("geometry_srs");
      45             :     AddArg("like", 0, _("Dataset to use as a template for bounds"),
      46          58 :            &m_likeDataset, GDAL_OF_RASTER | GDAL_OF_VECTOR)
      47          58 :         .SetMetaVar("DATASET")
      48          29 :         .SetMutualExclusionGroup("bbox-geometry-like");
      49             :     AddArg("like-sql", 0, ("SELECT statement to run on the 'like' dataset"),
      50          58 :            &m_likeSQL)
      51          58 :         .SetMetaVar("SELECT-STATEMENT")
      52          29 :         .SetMutualExclusionGroup("sql-where");
      53             :     AddArg("like-layer", 0, ("Name of the layer of the 'like' dataset"),
      54          58 :            &m_likeLayer)
      55          29 :         .SetMetaVar("LAYER-NAME");
      56             :     AddArg("like-where", 0, ("WHERE SQL clause to run on the 'like' dataset"),
      57          58 :            &m_likeWhere)
      58          58 :         .SetMetaVar("WHERE-EXPRESSION")
      59          29 :         .SetMutualExclusionGroup("sql-where");
      60             :     AddArg("only-bbox", 0,
      61             :            _("For 'geometry' and 'like', only consider their bounding box"),
      62          29 :            &m_onlyBBOX);
      63             :     AddArg("allow-bbox-outside-source", 0,
      64             :            _("Allow clipping box to include pixels outside input dataset"),
      65          29 :            &m_allowExtentOutsideSource);
      66             :     AddArg("add-alpha", 0,
      67             :            _("Adds an alpha mask band to the destination when the source "
      68             :              "raster have none."),
      69          29 :            &m_addAlpha);
      70          29 : }
      71             : 
      72             : /************************************************************************/
      73             : /*                 GDALRasterClipAlgorithm::RunStep()                   */
      74             : /************************************************************************/
      75             : 
      76          21 : bool GDALRasterClipAlgorithm::RunStep(GDALProgressFunc, void *)
      77             : {
      78          21 :     auto poSrcDS = m_inputDataset.GetDatasetRef();
      79          21 :     CPLAssert(poSrcDS);
      80          21 :     CPLAssert(m_outputDataset.GetName().empty());
      81          21 :     CPLAssert(!m_outputDataset.GetDatasetRef());
      82             : 
      83             :     double adfGT[6];
      84          21 :     if (poSrcDS->GetGeoTransform(adfGT) != CE_None)
      85             :     {
      86           1 :         ReportError(
      87             :             CE_Failure, CPLE_NotSupported,
      88             :             "Clipping is not supported on a raster without a geotransform");
      89           1 :         return false;
      90             :     }
      91          20 :     if (adfGT[2] != 0 && adfGT[4] != 0)
      92             :     {
      93           1 :         ReportError(CE_Failure, CPLE_NotSupported,
      94             :                     "Clipping is not supported on a raster whose geotransform "
      95             :                     "has rotation terms");
      96           1 :         return false;
      97             :     }
      98             : 
      99          38 :     auto [poClipGeom, errMsg] = GetClipGeometry();
     100          19 :     if (!poClipGeom)
     101             :     {
     102           3 :         ReportError(CE_Failure, CPLE_AppDefined, "%s", errMsg.c_str());
     103           3 :         return false;
     104             :     }
     105             : 
     106          16 :     auto poLikeDS = m_likeDataset.GetDatasetRef();
     107          17 :     if (!poClipGeom->getSpatialReference() && poLikeDS &&
     108           1 :         poLikeDS->GetLayerCount() == 0)
     109             :     {
     110           1 :         ReportError(CE_Failure, CPLE_AppDefined,
     111             :                     "Dataset '%s' has no CRS. Its bounds cannot be used.",
     112           1 :                     poLikeDS->GetDescription());
     113           1 :         return false;
     114             :     }
     115             : 
     116          30 :     CPLStringList aosOptions;
     117          15 :     aosOptions.AddString("-of");
     118          15 :     aosOptions.AddString("VRT");
     119             : 
     120          15 :     OGREnvelope env;
     121          15 :     poClipGeom->getEnvelope(&env);
     122             : 
     123          15 :     if (m_onlyBBOX)
     124             :     {
     125           2 :         auto poPoly = std::make_unique<OGRPolygon>(env);
     126           1 :         poPoly->assignSpatialReference(poClipGeom->getSpatialReference());
     127           1 :         poClipGeom = std::move(poPoly);
     128             :     }
     129             : 
     130          15 :     const bool bBottomUpRaster = adfGT[5] > 0;
     131             : 
     132          15 :     if (poClipGeom->IsRectangle() && !m_addAlpha && !bBottomUpRaster)
     133             :     {
     134           9 :         aosOptions.AddString("-projwin");
     135           9 :         aosOptions.AddString(CPLSPrintf("%.17g", env.MinX));
     136           9 :         aosOptions.AddString(CPLSPrintf("%.17g", env.MaxY));
     137           9 :         aosOptions.AddString(CPLSPrintf("%.17g", env.MaxX));
     138           9 :         aosOptions.AddString(CPLSPrintf("%.17g", env.MinY));
     139             : 
     140           9 :         auto poClipGeomSRS = poClipGeom->getSpatialReference();
     141           9 :         if (poClipGeomSRS)
     142             :         {
     143           2 :             const char *const apszOptions[] = {"FORMAT=WKT2", nullptr};
     144           4 :             const std::string osWKT = poClipGeomSRS->exportToWkt(apszOptions);
     145           2 :             aosOptions.AddString("-projwin_srs");
     146           2 :             aosOptions.AddString(osWKT.c_str());
     147             :         }
     148             : 
     149           9 :         if (!m_allowExtentOutsideSource)
     150             :         {
     151             :             // Unless we've specifically allowed the bounding box to extend beyond
     152             :             // the source raster, raise an error.
     153           6 :             aosOptions.AddString("-epo");
     154             :         }
     155             : 
     156             :         GDALTranslateOptions *psOptions =
     157           9 :             GDALTranslateOptionsNew(aosOptions.List(), nullptr);
     158             : 
     159           9 :         GDALDatasetH hSrcDS = GDALDataset::ToHandle(poSrcDS);
     160           9 :         auto poRetDS = GDALDataset::FromHandle(
     161             :             GDALTranslate("", hSrcDS, psOptions, nullptr));
     162           9 :         GDALTranslateOptionsFree(psOptions);
     163             : 
     164           9 :         const bool bOK = poRetDS != nullptr;
     165           9 :         if (bOK)
     166             :         {
     167           7 :             m_outputDataset.Set(std::unique_ptr<GDALDataset>(poRetDS));
     168             :         }
     169             : 
     170           9 :         return bOK;
     171             :     }
     172             :     else
     173             :     {
     174           6 :         if (bBottomUpRaster)
     175             :         {
     176           1 :             adfGT[3] += adfGT[5] * poSrcDS->GetRasterYSize();
     177           1 :             adfGT[5] = -adfGT[5];
     178             :         }
     179             : 
     180             :         {
     181             :             auto poClipGeomInSrcSRS =
     182          12 :                 std::unique_ptr<OGRGeometry>(poClipGeom->clone());
     183           6 :             if (poClipGeom->getSpatialReference() && poSrcDS->GetSpatialRef())
     184           1 :                 poClipGeomInSrcSRS->transformTo(poSrcDS->GetSpatialRef());
     185           6 :             poClipGeomInSrcSRS->getEnvelope(&env);
     186             :         }
     187             : 
     188          10 :         if (!m_allowExtentOutsideSource &&
     189           4 :             !(env.MinX >= adfGT[0] &&
     190           3 :               env.MaxX <= adfGT[0] + adfGT[1] * poSrcDS->GetRasterXSize() &&
     191           3 :               env.MaxY >= adfGT[3] &&
     192           3 :               env.MinY <= adfGT[3] + adfGT[5] * poSrcDS->GetRasterYSize()))
     193             :         {
     194           1 :             ReportError(CE_Failure, CPLE_AppDefined,
     195             :                         "Clipping geometry is partially or totally outside the "
     196             :                         "extent of the raster. You can set the "
     197             :                         "'allow-bbox-outside-source' argument to proceed.");
     198           1 :             return false;
     199             :         }
     200             : 
     201           5 :         if (m_addAlpha)
     202             :         {
     203           2 :             aosOptions.AddString("-dstalpha");
     204             :         }
     205             : 
     206           5 :         aosOptions.AddString("-cutline");
     207           5 :         aosOptions.AddString(poClipGeom->exportToWkt());
     208             : 
     209           5 :         aosOptions.AddString("-wo");
     210           5 :         aosOptions.AddString("CUTLINE_ALL_TOUCHED=YES");
     211             : 
     212           5 :         auto poClipGeomSRS = poClipGeom->getSpatialReference();
     213           5 :         if (poClipGeomSRS)
     214             :         {
     215           1 :             const char *const apszOptions[] = {"FORMAT=WKT2", nullptr};
     216           2 :             const std::string osWKT = poClipGeomSRS->exportToWkt(apszOptions);
     217           1 :             aosOptions.AddString("-cutline_srs");
     218           1 :             aosOptions.AddString(osWKT.c_str());
     219             :         }
     220             : 
     221           5 :         constexpr double REL_EPS_PIXEL = 1e-3;
     222           5 :         const double dfMinX =
     223           5 :             adfGT[0] +
     224           5 :             floor((env.MinX - adfGT[0]) / adfGT[1] + REL_EPS_PIXEL) * adfGT[1];
     225           5 :         const double dfMinY =
     226           5 :             adfGT[3] +
     227           5 :             ceil((env.MinY - adfGT[3]) / adfGT[5] - REL_EPS_PIXEL) * adfGT[5];
     228           5 :         const double dfMaxX =
     229           5 :             adfGT[0] +
     230           5 :             ceil((env.MaxX - adfGT[0]) / adfGT[1] - REL_EPS_PIXEL) * adfGT[1];
     231           5 :         const double dfMaxY =
     232           5 :             adfGT[3] +
     233           5 :             floor((env.MaxY - adfGT[3]) / adfGT[5] + REL_EPS_PIXEL) * adfGT[5];
     234             : 
     235           5 :         aosOptions.AddString("-te");
     236           5 :         aosOptions.AddString(CPLSPrintf("%.17g", dfMinX));
     237             :         aosOptions.AddString(
     238           5 :             CPLSPrintf("%.17g", bBottomUpRaster ? dfMaxY : dfMinY));
     239           5 :         aosOptions.AddString(CPLSPrintf("%.17g", dfMaxX));
     240             :         aosOptions.AddString(
     241           5 :             CPLSPrintf("%.17g", bBottomUpRaster ? dfMinY : dfMaxY));
     242             : 
     243           5 :         aosOptions.AddString("-tr");
     244           5 :         aosOptions.AddString(CPLSPrintf("%.17g", adfGT[1]));
     245           5 :         aosOptions.AddString(CPLSPrintf("%.17g", std::fabs(adfGT[5])));
     246             : 
     247             :         GDALWarpAppOptions *psOptions =
     248           5 :             GDALWarpAppOptionsNew(aosOptions.List(), nullptr);
     249             : 
     250           5 :         GDALDatasetH hSrcDS = GDALDataset::ToHandle(poSrcDS);
     251           5 :         auto poRetDS = GDALDataset::FromHandle(
     252             :             GDALWarp("", nullptr, 1, &hSrcDS, psOptions, nullptr));
     253           5 :         GDALWarpAppOptionsFree(psOptions);
     254             : 
     255           5 :         const bool bOK = poRetDS != nullptr;
     256           5 :         if (bOK)
     257             :         {
     258           5 :             m_outputDataset.Set(std::unique_ptr<GDALDataset>(poRetDS));
     259             :         }
     260             : 
     261           5 :         return bOK;
     262             :     }
     263             : }
     264             : 
     265             : //! @endcond