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