Line data    Source code

       1             : /******************************************************************************
       2             :  *
       3             :  * Project:  GDAL
       4             :  * Purpose:  "reproject" step of "raster pipeline"
       5             :  * Author:   Even Rouault <even dot rouault at spatialys.com>
       6             :  *
       7             :  ******************************************************************************
       8             :  * Copyright (c) 2024, Even Rouault <even dot rouault at spatialys.com>
       9             :  *
      10             :  * SPDX-License-Identifier: MIT
      11             :  ****************************************************************************/
      12             : 
      13             : #include "gdalalg_raster_reproject.h"
      14             : #include "gdalalg_raster_write.h"
      15             : 
      16             : #include "gdal_alg.h"
      17             : #include "gdal_priv.h"
      18             : #include "gdal_utils.h"
      19             : #include "gdalwarper.h"
      20             : 
      21             : #include <cmath>
      22             : 
      23             : //! @cond Doxygen_Suppress
      24             : 
      25             : #ifndef _
      26             : #define _(x) (x)
      27             : #endif
      28             : 
      29             : /************************************************************************/
      30             : /*     GDALRasterReprojectAlgorithm::GDALRasterReprojectAlgorithm()     */
      31             : /************************************************************************/
      32             : 
      33         150 : GDALRasterReprojectAlgorithm::GDALRasterReprojectAlgorithm(bool standaloneStep)
      34             :     : GDALRasterPipelineStepAlgorithm(NAME, DESCRIPTION, HELP_URL,
      35         150 :                                       standaloneStep)
      36             : {
      37             : 
      38         300 :     AddArg(GDAL_ARG_NAME_INPUT_CRS, 's', _("Input CRS"), &m_srcCrs)
      39         300 :         .SetIsCRSArg()
      40         300 :         .AddHiddenAlias("s_srs")
      41         150 :         .AddHiddenAlias("src-crs");
      42             : 
      43             :     AddArg("like", 0,
      44             :            _("Dataset to use as a template for target bounds, CRS, size and "
      45             :              "nodata"),
      46         300 :            &m_likeDataset, GDAL_OF_RASTER)
      47         150 :         .SetMetaVar("DATASET");
      48             : 
      49         300 :     AddArg(GDAL_ARG_NAME_OUTPUT_CRS, 'd', _("Output CRS"), &m_dstCrs)
      50         300 :         .SetIsCRSArg()
      51         300 :         .AddHiddenAlias("t_srs")
      52         150 :         .AddHiddenAlias("dst-crs");
      53             : 
      54         150 :     GDALRasterReprojectUtils::AddResamplingArg(this, m_resampling);
      55             : 
      56             :     AddArg("resolution", 0, _("Target resolution (in destination CRS units)"),
      57         300 :            &m_resolution)
      58         150 :         .SetMinCount(2)
      59         150 :         .SetMaxCount(2)
      60         150 :         .SetMinValueExcluded(0)
      61         150 :         .SetRepeatedArgAllowed(false)
      62         150 :         .SetDisplayHintAboutRepetition(false)
      63         300 :         .SetMetaVar("<xres>,<yres>")
      64         150 :         .SetMutualExclusionGroup("resolution-size");
      65             : 
      66         300 :     AddArg("size", 0, _("Target size in pixels"), &m_size)
      67         150 :         .SetMinCount(2)
      68         150 :         .SetMaxCount(2)
      69         150 :         .SetMinValueIncluded(0)
      70         150 :         .SetRepeatedArgAllowed(false)
      71         150 :         .SetDisplayHintAboutRepetition(false)
      72         300 :         .SetMetaVar("<width>,<height>")
      73         150 :         .SetMutualExclusionGroup("resolution-size");
      74             : 
      75             :     auto &arg = AddBBOXArg(&m_bbox,
      76         150 :                            _("Target bounding box (in destination CRS units)"));
      77             : 
      78             :     arg.AddValidationAction(
      79          10 :         [this, &arg]()
      80             :         {
      81             :             // Validate it's not empty
      82           8 :             const std::vector<double> &bbox = arg.Get<std::vector<double>>();
      83           8 :             if ((bbox[0] >= bbox[2]) || (bbox[1] >= bbox[3]))
      84             :             {
      85           2 :                 ReportError(CE_Failure, CPLE_AppDefined,
      86             :                             "Invalid bounding box specified");
      87           2 :                 return false;
      88             :             }
      89             :             else
      90             :             {
      91           6 :                 return true;
      92             :             }
      93         150 :         });
      94             : 
      95         300 :     AddArg("bbox-crs", 0, _("CRS of target bounding box"), &m_bboxCrs)
      96         300 :         .SetIsCRSArg()
      97         150 :         .AddHiddenAlias("bbox_srs");
      98             : 
      99             :     AddArg("target-aligned-pixels", 0,
     100             :            _("Round target extent to target resolution"),
     101         300 :            &m_targetAlignedPixels)
     102         300 :         .AddHiddenAlias("tap")
     103         150 :         .SetCategory(GAAC_ADVANCED);
     104             :     AddArg("input-nodata", 0,
     105             :            _("Set nodata values for input bands ('None' to unset)."),
     106         300 :            &m_srcNoData)
     107         150 :         .SetMinCount(1)
     108         300 :         .AddHiddenAlias("src-nodata")
     109         150 :         .SetRepeatedArgAllowed(false)
     110         150 :         .SetCategory(GAAC_ADVANCED);
     111             :     AddArg("output-nodata", 0,
     112             :            _("Set nodata values for output bands ('None' to unset)."),
     113         300 :            &m_dstNoData)
     114         150 :         .SetMinCount(1)
     115         300 :         .AddHiddenAlias("dst-nodata")
     116         150 :         .SetRepeatedArgAllowed(false)
     117         150 :         .SetCategory(GAAC_ADVANCED);
     118             :     AddArg("add-alpha", 0,
     119             :            _("Adds an alpha mask band to the destination when the source "
     120             :              "raster have none."),
     121         300 :            &m_addAlpha)
     122         150 :         .SetCategory(GAAC_ADVANCED);
     123             : 
     124         150 :     GDALRasterReprojectUtils::AddWarpOptTransformOptErrorThresholdArg(
     125         150 :         this, m_warpOptions, m_transformOptions, m_errorThreshold);
     126             : 
     127         150 :     AddNumThreadsArg(&m_numThreads, &m_numThreadsStr);
     128         150 : }
     129             : 
     130             : /************************************************************************/
     131             : /*             GDALRasterReprojectUtils::AddResamplingArg()             */
     132             : /************************************************************************/
     133             : 
     134             : /*static */ void
     135         290 : GDALRasterReprojectUtils::AddResamplingArg(GDALAlgorithm *alg,
     136             :                                            std::string &resampling)
     137             : {
     138         580 :     alg->AddArg("resampling", 'r', _("Resampling method"), &resampling)
     139             :         .SetChoices("nearest", "bilinear", "cubic", "cubicspline", "lanczos",
     140             :                     "average", "rms", "mode", "min", "max", "med", "q1", "q3",
     141         290 :                     "sum")
     142         290 :         .SetDefault("nearest")
     143         290 :         .SetHiddenChoices("near");
     144         290 : }
     145             : 
     146             : /************************************************************************/
     147             : /*              AddWarpOptTransformOptErrorThresholdArg()               */
     148             : /************************************************************************/
     149             : 
     150             : /* static */
     151         243 : void GDALRasterReprojectUtils::AddWarpOptTransformOptErrorThresholdArg(
     152             :     GDALAlgorithm *alg, std::vector<std::string> &warpOptions,
     153             :     std::vector<std::string> &transformOptions, double &errorThreshold)
     154             : {
     155             :     {
     156             :         auto &arg =
     157         486 :             alg->AddArg("warp-option", 0, _("Warping option(s)"), &warpOptions)
     158         486 :                 .AddAlias("wo")
     159         486 :                 .SetMetaVar("<NAME>=<VALUE>")
     160         486 :                 .SetCategory(GAAC_ADVANCED)
     161         243 :                 .SetPackedValuesAllowed(false);
     162          10 :         arg.AddValidationAction([alg, &arg]()
     163         253 :                                 { return alg->ParseAndValidateKeyValue(arg); });
     164             :         arg.SetAutoCompleteFunction(
     165           0 :             [](const std::string &currentValue)
     166             :             {
     167           0 :                 std::vector<std::string> ret;
     168           0 :                 GDALAlgorithm::AddOptionsSuggestions(GDALWarpGetOptionList(), 0,
     169             :                                                      currentValue, ret);
     170           0 :                 return ret;
     171         243 :             });
     172             :     }
     173             :     {
     174             :         auto &arg = alg->AddArg("transform-option", 0, _("Transform option(s)"),
     175         486 :                                 &transformOptions)
     176         486 :                         .AddAlias("to")
     177         486 :                         .SetMetaVar("<NAME>=<VALUE>")
     178         486 :                         .SetCategory(GAAC_ADVANCED)
     179         243 :                         .SetPackedValuesAllowed(false);
     180           4 :         arg.AddValidationAction([alg, &arg]()
     181         247 :                                 { return alg->ParseAndValidateKeyValue(arg); });
     182             :         arg.SetAutoCompleteFunction(
     183           0 :             [](const std::string &currentValue)
     184             :             {
     185           0 :                 std::vector<std::string> ret;
     186           0 :                 GDALAlgorithm::AddOptionsSuggestions(
     187             :                     GDALGetGenImgProjTranformerOptionList(), 0, currentValue,
     188             :                     ret);
     189           0 :                 return ret;
     190         243 :             });
     191             :     }
     192         486 :     alg->AddArg("error-threshold", 0, _("Error threshold"), &errorThreshold)
     193         486 :         .AddAlias("et")
     194         243 :         .SetMinValueIncluded(0)
     195         243 :         .SetCategory(GAAC_ADVANCED);
     196         243 : }
     197             : 
     198             : /************************************************************************/
     199             : /*          GDALRasterReprojectAlgorithm::CanHandleNextStep()           */
     200             : /************************************************************************/
     201             : 
     202          52 : bool GDALRasterReprojectAlgorithm::CanHandleNextStep(
     203             :     GDALPipelineStepAlgorithm *poNextStep) const
     204             : {
     205         102 :     return poNextStep->GetName() == GDALRasterWriteAlgorithm::NAME &&
     206         102 :            poNextStep->GetOutputFormat() != "stream";
     207             : }
     208             : 
     209             : /************************************************************************/
     210             : /*               GDALRasterReprojectAlgorithm::RunStep()                */
     211             : /************************************************************************/
     212             : 
     213          29 : bool GDALRasterReprojectAlgorithm::RunStep(GDALPipelineStepRunContext &ctxt)
     214             : {
     215          29 :     auto poSrcDS = m_inputDataset[0].GetDatasetRef();
     216          29 :     CPLAssert(poSrcDS);
     217          29 :     CPLAssert(m_outputDataset.GetName().empty());
     218          29 :     CPLAssert(!m_outputDataset.GetDatasetRef());
     219             : 
     220          58 :     CPLStringList aosOptions;
     221          58 :     std::string outputFilename;
     222             : 
     223          29 :     aosOptions.AddString("--invoked-from-gdal-algorithm");
     224             : 
     225             :     // --like provide defaults: override if not explicitly set
     226          29 :     if (auto poLikeDS = m_likeDataset.GetDatasetRef())
     227             :     {
     228           2 :         const auto poSpatialRef = poLikeDS->GetSpatialRef();
     229           2 :         if (poSpatialRef)
     230             :         {
     231           2 :             char *pszWKT = nullptr;
     232           2 :             poSpatialRef->exportToWkt(&pszWKT);
     233           2 :             m_dstCrs = pszWKT;
     234           2 :             CPLFree(pszWKT);
     235           2 :             GDALGeoTransform gt;
     236           2 :             if (poLikeDS->GetGeoTransform(gt) == CE_None)
     237             :             {
     238           2 :                 if (gt.IsAxisAligned())
     239             :                 {
     240           1 :                     if (m_resolution.empty())
     241             :                     {
     242           1 :                         m_resolution = {std::abs(gt[1]), std::abs(gt[5])};
     243             :                     }
     244           1 :                     const int nXSize = poLikeDS->GetRasterXSize();
     245           1 :                     const int nYSize = poLikeDS->GetRasterYSize();
     246           1 :                     if (m_size.empty())
     247             :                     {
     248           1 :                         m_size = {nXSize, nYSize};
     249             :                     }
     250           1 :                     if (m_bbox.empty())
     251             :                     {
     252           1 :                         double minX = gt.xorig;
     253           1 :                         double maxY = gt.yorig;
     254           1 :                         double maxX =
     255           1 :                             gt.xorig + nXSize * gt.xscale + nYSize * gt.xrot;
     256           1 :                         double minY =
     257           1 :                             gt.yorig + nXSize * gt.yrot + nYSize * gt.yscale;
     258           1 :                         if (minY > maxY)
     259           0 :                             std::swap(minY, maxY);
     260           1 :                         m_bbox = {minX, minY, maxX, maxY};
     261           1 :                         m_bboxCrs = m_dstCrs;
     262             :                     }
     263             :                 }
     264             :                 else
     265             :                 {
     266           1 :                     ReportError(
     267             :                         CE_Warning, CPLE_AppDefined,
     268             :                         "Dataset provided with --like has a geotransform "
     269             :                         "with rotation. Ignoring it");
     270             :                 }
     271             :             }
     272             :         }
     273             :     }
     274             : 
     275          29 :     if (ctxt.m_poNextUsableStep)
     276             :     {
     277          25 :         CPLAssert(CanHandleNextStep(ctxt.m_poNextUsableStep));
     278          25 :         outputFilename = ctxt.m_poNextUsableStep->GetOutputDataset().GetName();
     279          25 :         const auto &format = ctxt.m_poNextUsableStep->GetOutputFormat();
     280          25 :         if (!format.empty())
     281             :         {
     282           8 :             aosOptions.AddString("-of");
     283           8 :             aosOptions.AddString(format.c_str());
     284             :         }
     285             : 
     286          25 :         bool bFoundNumThreads = false;
     287           2 :         for (const std::string &co :
     288          29 :              ctxt.m_poNextUsableStep->GetCreationOptions())
     289             :         {
     290           2 :             aosOptions.AddString("-co");
     291           2 :             if (STARTS_WITH_CI(co.c_str(), "NUM_THREADS="))
     292           0 :                 bFoundNumThreads = true;
     293           2 :             aosOptions.AddString(co.c_str());
     294             :         }
     295             : 
     296             :         // Forward m_numThreads to GeoTIFF driver if --co NUM_THREADS not
     297             :         // specified
     298          50 :         if (!bFoundNumThreads && m_numThreads > 1 &&
     299          49 :             (EQUAL(format.c_str(), "GTIFF") || EQUAL(format.c_str(), "COG") ||
     300          24 :              (format.empty() &&
     301          42 :               EQUAL(CPLGetExtensionSafe(outputFilename.c_str()).c_str(),
     302             :                     "tif"))))
     303             :         {
     304          17 :             aosOptions.AddString("-co");
     305          17 :             aosOptions.AddString(CPLSPrintf("NUM_THREADS=%d", m_numThreads));
     306             :         }
     307             :     }
     308             :     else
     309             :     {
     310           4 :         aosOptions.AddString("-of");
     311           4 :         aosOptions.AddString("VRT");
     312             :     }
     313          29 :     if (!m_srcCrs.empty())
     314             :     {
     315          10 :         aosOptions.AddString("-s_srs");
     316          10 :         aosOptions.AddString(m_srcCrs.c_str());
     317             :     }
     318          29 :     if (!m_dstCrs.empty())
     319             :     {
     320          22 :         aosOptions.AddString("-t_srs");
     321          22 :         aosOptions.AddString(m_dstCrs.c_str());
     322             :     }
     323          29 :     if (!m_resampling.empty())
     324             :     {
     325          29 :         aosOptions.AddString("-r");
     326          29 :         aosOptions.AddString(m_resampling.c_str());
     327             :     }
     328          29 :     if (!m_resolution.empty())
     329             :     {
     330           2 :         aosOptions.AddString("-tr");
     331           2 :         aosOptions.AddString(CPLSPrintf("%.17g", m_resolution[0]));
     332           2 :         aosOptions.AddString(CPLSPrintf("%.17g", m_resolution[1]));
     333             :     }
     334          29 :     if (!m_size.empty())
     335             :     {
     336           2 :         aosOptions.AddString("-ts");
     337           2 :         aosOptions.AddString(CPLSPrintf("%d", m_size[0]));
     338           2 :         aosOptions.AddString(CPLSPrintf("%d", m_size[1]));
     339             :     }
     340          29 :     if (!m_bbox.empty())
     341             :     {
     342           3 :         aosOptions.AddString("-te");
     343           3 :         aosOptions.AddString(CPLSPrintf("%.17g", m_bbox[0]));
     344           3 :         aosOptions.AddString(CPLSPrintf("%.17g", m_bbox[1]));
     345           3 :         aosOptions.AddString(CPLSPrintf("%.17g", m_bbox[2]));
     346           3 :         aosOptions.AddString(CPLSPrintf("%.17g", m_bbox[3]));
     347             :     }
     348          29 :     if (!m_bboxCrs.empty())
     349             :     {
     350           2 :         aosOptions.AddString("-te_srs");
     351           2 :         aosOptions.AddString(m_bboxCrs.c_str());
     352             :     }
     353          29 :     if (m_targetAlignedPixels)
     354             :     {
     355           1 :         aosOptions.AddString("-tap");
     356             :     }
     357          29 :     if (!m_srcNoData.empty())
     358             :     {
     359           1 :         aosOptions.push_back("-srcnodata");
     360           2 :         std::string s;
     361           2 :         for (const std::string &v : m_srcNoData)
     362             :         {
     363           1 :             if (!s.empty())
     364           0 :                 s += " ";
     365           1 :             s += v;
     366             :         }
     367           1 :         aosOptions.push_back(s);
     368             :     }
     369          29 :     if (!m_dstNoData.empty())
     370             :     {
     371           2 :         aosOptions.push_back("-dstnodata");
     372           4 :         std::string s;
     373           5 :         for (const std::string &v : m_dstNoData)
     374             :         {
     375           3 :             if (!s.empty())
     376           1 :                 s += " ";
     377           3 :             s += v;
     378             :         }
     379           2 :         aosOptions.push_back(s);
     380             :     }
     381          29 :     if (m_addAlpha)
     382             :     {
     383           1 :         aosOptions.AddString("-dstalpha");
     384             :     }
     385             : 
     386          29 :     bool bFoundNumThreads = false;
     387          32 :     for (const std::string &opt : m_warpOptions)
     388             :     {
     389           3 :         aosOptions.AddString("-wo");
     390           3 :         if (STARTS_WITH_CI(opt.c_str(), "NUM_THREADS="))
     391           2 :             bFoundNumThreads = true;
     392           3 :         aosOptions.AddString(opt.c_str());
     393             :     }
     394          29 :     if (bFoundNumThreads)
     395             :     {
     396           2 :         if (GetArg(GDAL_ARG_NAME_NUM_THREADS)->IsExplicitlySet())
     397             :         {
     398           1 :             ReportError(CE_Failure, CPLE_AppDefined,
     399             :                         "--num-threads argument and NUM_THREADS warp options "
     400             :                         "are mutually exclusive.");
     401           1 :             return false;
     402             :         }
     403             :     }
     404             :     else
     405             :     {
     406          27 :         aosOptions.AddString("-wo");
     407          27 :         aosOptions.AddString(CPLSPrintf("NUM_THREADS=%d", m_numThreads));
     408             :     }
     409             : 
     410          29 :     for (const std::string &opt : m_transformOptions)
     411             :     {
     412           1 :         aosOptions.AddString("-to");
     413           1 :         aosOptions.AddString(opt.c_str());
     414             :     }
     415          28 :     if (std::isfinite(m_errorThreshold))
     416             :     {
     417           0 :         aosOptions.AddString("-et");
     418           0 :         aosOptions.AddString(CPLSPrintf("%.17g", m_errorThreshold));
     419             :     }
     420             : 
     421          28 :     bool bOK = false;
     422             :     GDALWarpAppOptions *psOptions =
     423          28 :         GDALWarpAppOptionsNew(aosOptions.List(), nullptr);
     424          28 :     if (psOptions)
     425             :     {
     426          28 :         if (ctxt.m_poNextUsableStep)
     427             :         {
     428          24 :             GDALWarpAppOptionsSetProgress(psOptions, ctxt.m_pfnProgress,
     429             :                                           ctxt.m_pProgressData);
     430             :         }
     431          28 :         GDALDatasetH hSrcDS = GDALDataset::ToHandle(poSrcDS);
     432          28 :         auto poRetDS = GDALDataset::FromHandle(GDALWarp(
     433             :             outputFilename.c_str(), nullptr, 1, &hSrcDS, psOptions, nullptr));
     434          28 :         GDALWarpAppOptionsFree(psOptions);
     435          28 :         bOK = poRetDS != nullptr;
     436          28 :         if (bOK)
     437             :         {
     438          25 :             m_outputDataset.Set(std::unique_ptr<GDALDataset>(poRetDS));
     439             :         }
     440             :     }
     441          28 :     return bOK;
     442             : }
     443             : 
     444             : GDALRasterReprojectAlgorithmStandalone::
     445             :     ~GDALRasterReprojectAlgorithmStandalone() = default;
     446             : 
     447             : //! @endcond