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