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             : 
      20             : //! @cond Doxygen_Suppress
      21             : 
      22             : #ifndef _
      23             : #define _(x) (x)
      24             : #endif
      25             : 
      26             : /************************************************************************/
      27             : /*          GDALRasterClipAlgorithm::GDALRasterClipAlgorithm()          */
      28             : /************************************************************************/
      29             : 
      30          14 : GDALRasterClipAlgorithm::GDALRasterClipAlgorithm(bool standaloneStep)
      31             :     : GDALRasterPipelineStepAlgorithm(NAME, DESCRIPTION, HELP_URL,
      32          14 :                                       standaloneStep)
      33             : {
      34          14 :     AddBBOXArg(&m_bbox, _("Clipping bounding box as xmin,ymin,xmax,ymax"))
      35          14 :         .SetMutualExclusionGroup("exclusion-group");
      36          28 :     AddArg("bbox-crs", 0, _("CRS of clipping bounding box"), &m_bboxCrs)
      37          14 :         .SetIsCRSArg()
      38          14 :         .AddHiddenAlias("bbox_srs");
      39             :     AddArg("like", 0, _("Raster dataset to use as a template for bounds"),
      40          28 :            &m_likeDataset, GDAL_OF_RASTER)
      41          28 :         .SetMetaVar("DATASET")
      42          14 :         .SetMutualExclusionGroup("exclusion-group");
      43          14 : }
      44             : 
      45             : /************************************************************************/
      46             : /*                 GDALRasterClipAlgorithm::RunStep()                   */
      47             : /************************************************************************/
      48             : 
      49           6 : bool GDALRasterClipAlgorithm::RunStep(GDALProgressFunc, void *)
      50             : {
      51           6 :     CPLAssert(m_inputDataset.GetDatasetRef());
      52           6 :     CPLAssert(m_outputDataset.GetName().empty());
      53           6 :     CPLAssert(!m_outputDataset.GetDatasetRef());
      54             : 
      55          12 :     CPLStringList aosOptions;
      56           6 :     aosOptions.AddString("-of");
      57           6 :     aosOptions.AddString("VRT");
      58           6 :     if (!m_bbox.empty())
      59             :     {
      60           2 :         aosOptions.AddString("-projwin");
      61           2 :         aosOptions.AddString(CPLSPrintf("%.17g", m_bbox[0]));  // minx
      62           2 :         aosOptions.AddString(CPLSPrintf("%.17g", m_bbox[3]));  // maxy
      63           2 :         aosOptions.AddString(CPLSPrintf("%.17g", m_bbox[2]));  // maxx
      64           2 :         aosOptions.AddString(CPLSPrintf("%.17g", m_bbox[1]));  // miny
      65             : 
      66           2 :         if (!m_bboxCrs.empty())
      67             :         {
      68           1 :             aosOptions.AddString("-projwin_srs");
      69           1 :             aosOptions.AddString(m_bboxCrs.c_str());
      70             :         }
      71             :     }
      72           4 :     else if (auto poLikeDS = m_likeDataset.GetDatasetRef())
      73             :     {
      74             :         double adfGT[6];
      75           3 :         if (poLikeDS->GetGeoTransform(adfGT) != CE_None)
      76             :         {
      77           1 :             ReportError(CE_Failure, CPLE_AppDefined,
      78             :                         "Dataset '%s' has no geotransform matrix. Its bounds "
      79             :                         "cannot be established.",
      80           1 :                         poLikeDS->GetDescription());
      81           2 :             return false;
      82             :         }
      83           2 :         auto poLikeSRS = poLikeDS->GetSpatialRef();
      84           2 :         if (!poLikeSRS)
      85             :         {
      86           1 :             ReportError(CE_Failure, CPLE_AppDefined,
      87             :                         "Dataset '%s' has no SRS. Its bounds cannot be used.",
      88           1 :                         poLikeDS->GetDescription());
      89           1 :             return false;
      90             :         }
      91           1 :         const double dfTLX = adfGT[0];
      92           1 :         const double dfTLY = adfGT[3];
      93           1 :         const double dfTRX = adfGT[0] + poLikeDS->GetRasterXSize() * adfGT[1];
      94           1 :         const double dfTRY = adfGT[3] + poLikeDS->GetRasterXSize() * adfGT[4];
      95           1 :         const double dfBLX = adfGT[0] + poLikeDS->GetRasterYSize() * adfGT[2];
      96           1 :         const double dfBLY = adfGT[3] + poLikeDS->GetRasterYSize() * adfGT[5];
      97           1 :         const double dfBRX = adfGT[0] + poLikeDS->GetRasterXSize() * adfGT[1] +
      98           1 :                              poLikeDS->GetRasterYSize() * adfGT[2];
      99           1 :         const double dfBRY = adfGT[3] + poLikeDS->GetRasterXSize() * adfGT[4] +
     100           1 :                              poLikeDS->GetRasterYSize() * adfGT[5];
     101             :         const double dfMinX =
     102           1 :             std::min(std::min(dfTLX, dfTRX), std::min(dfBLX, dfBRX));
     103             :         const double dfMinY =
     104           1 :             std::min(std::min(dfTLY, dfTRY), std::min(dfBLY, dfBRY));
     105             :         const double dfMaxX =
     106           1 :             std::max(std::max(dfTLX, dfTRX), std::max(dfBLX, dfBRX));
     107             :         const double dfMaxY =
     108           1 :             std::max(std::max(dfTLY, dfTRY), std::max(dfBLY, dfBRY));
     109             : 
     110           1 :         aosOptions.AddString("-projwin");
     111           1 :         aosOptions.AddString(CPLSPrintf("%.17g", dfMinX));
     112           1 :         aosOptions.AddString(CPLSPrintf("%.17g", dfMaxY));
     113           1 :         aosOptions.AddString(CPLSPrintf("%.17g", dfMaxX));
     114           1 :         aosOptions.AddString(CPLSPrintf("%.17g", dfMinY));
     115             : 
     116           1 :         const char *const apszOptions[] = {"FORMAT=WKT2", nullptr};
     117           2 :         const std::string osWKT = poLikeSRS->exportToWkt(apszOptions);
     118           1 :         aosOptions.AddString("-projwin_srs");
     119           1 :         aosOptions.AddString(osWKT.c_str());
     120             :     }
     121             :     else
     122             :     {
     123           1 :         ReportError(CE_Failure, CPLE_AppDefined,
     124             :                     "Either --bbox or --like must be specified");
     125           1 :         return false;
     126             :     }
     127             : 
     128             :     GDALTranslateOptions *psOptions =
     129           3 :         GDALTranslateOptionsNew(aosOptions.List(), nullptr);
     130             : 
     131           3 :     GDALDatasetH hSrcDS = GDALDataset::ToHandle(m_inputDataset.GetDatasetRef());
     132             :     auto poRetDS =
     133           3 :         GDALDataset::FromHandle(GDALTranslate("", hSrcDS, psOptions, nullptr));
     134           3 :     GDALTranslateOptionsFree(psOptions);
     135             : 
     136           3 :     const bool bOK = poRetDS != nullptr;
     137           3 :     if (bOK)
     138             :     {
     139           3 :         m_outputDataset.Set(std::unique_ptr<GDALDataset>(poRetDS));
     140             :     }
     141             : 
     142           3 :     return bOK;
     143             : }
     144             : 
     145             : //! @endcond