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 139 : GDALRasterReprojectAlgorithm::GDALRasterReprojectAlgorithm(bool standaloneStep)
34 : : GDALRasterPipelineStepAlgorithm(NAME, DESCRIPTION, HELP_URL,
35 139 : standaloneStep)
36 : {
37 :
38 278 : AddArg(GDAL_ARG_NAME_INPUT_CRS, 's', _("Input CRS"), &m_srcCrs)
39 278 : .SetIsCRSArg()
40 278 : .AddHiddenAlias("s_srs")
41 139 : .AddHiddenAlias("src-crs");
42 :
43 : AddArg("like", 0,
44 : _("Dataset to use as a template for target bounds, CRS, size and "
45 : "nodata"),
46 278 : &m_likeDataset, GDAL_OF_RASTER)
47 139 : .SetMetaVar("DATASET");
48 :
49 278 : AddArg(GDAL_ARG_NAME_OUTPUT_CRS, 'd', _("Output CRS"), &m_dstCrs)
50 278 : .SetIsCRSArg()
51 278 : .AddHiddenAlias("t_srs")
52 139 : .AddHiddenAlias("dst-crs");
53 :
54 139 : GDALRasterReprojectUtils::AddResamplingArg(this, m_resampling);
55 :
56 : AddArg("resolution", 0, _("Target resolution (in destination CRS units)"),
57 278 : &m_resolution)
58 139 : .SetMinCount(2)
59 139 : .SetMaxCount(2)
60 139 : .SetMinValueExcluded(0)
61 139 : .SetRepeatedArgAllowed(false)
62 139 : .SetDisplayHintAboutRepetition(false)
63 278 : .SetMetaVar("<xres>,<yres>")
64 139 : .SetMutualExclusionGroup("resolution-size");
65 :
66 278 : AddArg("size", 0, _("Target size in pixels"), &m_size)
67 139 : .SetMinCount(2)
68 139 : .SetMaxCount(2)
69 139 : .SetMinValueIncluded(0)
70 139 : .SetRepeatedArgAllowed(false)
71 139 : .SetDisplayHintAboutRepetition(false)
72 278 : .SetMetaVar("<width>,<height>")
73 139 : .SetMutualExclusionGroup("resolution-size");
74 :
75 : auto &arg = AddBBOXArg(&m_bbox,
76 139 : _("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 139 : });
94 :
95 278 : AddArg("bbox-crs", 0, _("CRS of target bounding box"), &m_bboxCrs)
96 278 : .SetIsCRSArg()
97 139 : .AddHiddenAlias("bbox_srs");
98 :
99 : AddArg("target-aligned-pixels", 0,
100 : _("Round target extent to target resolution"),
101 278 : &m_targetAlignedPixels)
102 278 : .AddHiddenAlias("tap")
103 139 : .SetCategory(GAAC_ADVANCED);
104 : AddArg("input-nodata", 0,
105 : _("Set nodata values for input bands ('None' to unset)."),
106 278 : &m_srcNoData)
107 139 : .SetMinCount(1)
108 278 : .AddHiddenAlias("src-nodata")
109 139 : .SetRepeatedArgAllowed(false)
110 139 : .SetCategory(GAAC_ADVANCED);
111 : AddArg("output-nodata", 0,
112 : _("Set nodata values for output bands ('None' to unset)."),
113 278 : &m_dstNoData)
114 139 : .SetMinCount(1)
115 278 : .AddHiddenAlias("dst-nodata")
116 139 : .SetRepeatedArgAllowed(false)
117 139 : .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 278 : &m_addAlpha)
122 139 : .SetCategory(GAAC_ADVANCED);
123 :
124 139 : GDALRasterReprojectUtils::AddWarpOptTransformOptErrorThresholdArg(
125 139 : this, m_warpOptions, m_transformOptions, m_errorThreshold);
126 :
127 139 : AddNumThreadsArg(&m_numThreads, &m_numThreadsStr);
128 139 : }
129 :
130 : /************************************************************************/
131 : /* GDALRasterReprojectUtils::AddResamplingArg() */
132 : /************************************************************************/
133 :
134 : /*static */ void
135 230 : GDALRasterReprojectUtils::AddResamplingArg(GDALAlgorithm *alg,
136 : std::string &resampling)
137 : {
138 460 : alg->AddArg("resampling", 'r', _("Resampling method"), &resampling)
139 : .SetChoices("nearest", "bilinear", "cubic", "cubicspline", "lanczos",
140 : "average", "rms", "mode", "min", "max", "med", "q1", "q3",
141 230 : "sum")
142 230 : .SetDefault("nearest")
143 230 : .SetHiddenChoices("near");
144 230 : }
145 :
146 : /************************************************************************/
147 : /* AddWarpOptTransformOptErrorThresholdArg() */
148 : /************************************************************************/
149 :
150 : /* static */
151 230 : 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 460 : alg->AddArg("warp-option", 0, _("Warping option(s)"), &warpOptions)
158 460 : .AddAlias("wo")
159 460 : .SetMetaVar("<NAME>=<VALUE>")
160 460 : .SetCategory(GAAC_ADVANCED)
161 230 : .SetPackedValuesAllowed(false);
162 10 : arg.AddValidationAction([alg, &arg]()
163 240 : { 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 230 : });
172 : }
173 : {
174 : auto &arg = alg->AddArg("transform-option", 0, _("Transform option(s)"),
175 460 : &transformOptions)
176 460 : .AddAlias("to")
177 460 : .SetMetaVar("<NAME>=<VALUE>")
178 460 : .SetCategory(GAAC_ADVANCED)
179 230 : .SetPackedValuesAllowed(false);
180 4 : arg.AddValidationAction([alg, &arg]()
181 234 : { 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 230 : });
191 : }
192 460 : alg->AddArg("error-threshold", 0, _("Error threshold"), &errorThreshold)
193 460 : .AddAlias("et")
194 230 : .SetMinValueIncluded(0)
195 230 : .SetCategory(GAAC_ADVANCED);
196 230 : }
197 :
198 : /************************************************************************/
199 : /* GDALRasterReprojectAlgorithm::CanHandleNextStep() */
200 : /************************************************************************/
201 :
202 40 : bool GDALRasterReprojectAlgorithm::CanHandleNextStep(
203 : GDALPipelineStepAlgorithm *poNextStep) const
204 : {
205 78 : return poNextStep->GetName() == GDALRasterWriteAlgorithm::NAME &&
206 78 : poNextStep->GetOutputFormat() != "stream";
207 : }
208 :
209 : /************************************************************************/
210 : /* GDALRasterReprojectAlgorithm::RunStep() */
211 : /************************************************************************/
212 :
213 23 : bool GDALRasterReprojectAlgorithm::RunStep(GDALPipelineStepRunContext &ctxt)
214 : {
215 23 : auto poSrcDS = m_inputDataset[0].GetDatasetRef();
216 23 : CPLAssert(poSrcDS);
217 23 : CPLAssert(m_outputDataset.GetName().empty());
218 23 : CPLAssert(!m_outputDataset.GetDatasetRef());
219 :
220 46 : CPLStringList aosOptions;
221 46 : std::string outputFilename;
222 :
223 : // --like provide defaults: override if not explicitly set
224 23 : 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 23 : if (ctxt.m_poNextUsableStep)
274 : {
275 19 : CPLAssert(CanHandleNextStep(ctxt.m_poNextUsableStep));
276 19 : outputFilename = ctxt.m_poNextUsableStep->GetOutputDataset().GetName();
277 19 : const auto &format = ctxt.m_poNextUsableStep->GetOutputFormat();
278 19 : if (!format.empty())
279 : {
280 7 : aosOptions.AddString("-of");
281 7 : aosOptions.AddString(format.c_str());
282 : }
283 :
284 19 : bool bFoundNumThreads = false;
285 2 : for (const std::string &co :
286 23 : 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 38 : if (!bFoundNumThreads && m_numThreads > 1 &&
297 37 : (EQUAL(format.c_str(), "GTIFF") || EQUAL(format.c_str(), "COG") ||
298 18 : (format.empty() &&
299 31 : 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 23 : if (!m_srcCrs.empty())
312 : {
313 9 : aosOptions.AddString("-s_srs");
314 9 : aosOptions.AddString(m_srcCrs.c_str());
315 : }
316 23 : if (!m_dstCrs.empty())
317 : {
318 16 : aosOptions.AddString("-t_srs");
319 16 : aosOptions.AddString(m_dstCrs.c_str());
320 : }
321 23 : if (!m_resampling.empty())
322 : {
323 23 : aosOptions.AddString("-r");
324 23 : aosOptions.AddString(m_resampling.c_str());
325 : }
326 23 : 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 23 : 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 23 : 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 23 : if (!m_bboxCrs.empty())
347 : {
348 2 : aosOptions.AddString("-te_srs");
349 2 : aosOptions.AddString(m_bboxCrs.c_str());
350 : }
351 23 : if (m_targetAlignedPixels)
352 : {
353 1 : aosOptions.AddString("-tap");
354 : }
355 23 : 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 23 : 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 23 : if (m_addAlpha)
380 : {
381 1 : aosOptions.AddString("-dstalpha");
382 : }
383 :
384 23 : bool bFoundNumThreads = false;
385 26 : 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 23 : 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 21 : aosOptions.AddString("-wo");
405 21 : aosOptions.AddString(CPLSPrintf("NUM_THREADS=%d", m_numThreads));
406 : }
407 :
408 23 : for (const std::string &opt : m_transformOptions)
409 : {
410 1 : aosOptions.AddString("-to");
411 1 : aosOptions.AddString(opt.c_str());
412 : }
413 22 : if (std::isfinite(m_errorThreshold))
414 : {
415 0 : aosOptions.AddString("-et");
416 0 : aosOptions.AddString(CPLSPrintf("%.17g", m_errorThreshold));
417 : }
418 :
419 22 : bool bOK = false;
420 : GDALWarpAppOptions *psOptions =
421 22 : GDALWarpAppOptionsNew(aosOptions.List(), nullptr);
422 22 : if (psOptions)
423 : {
424 22 : if (ctxt.m_poNextUsableStep)
425 : {
426 18 : GDALWarpAppOptionsSetProgress(psOptions, ctxt.m_pfnProgress,
427 : ctxt.m_pProgressData);
428 : }
429 22 : GDALDatasetH hSrcDS = GDALDataset::ToHandle(poSrcDS);
430 22 : auto poRetDS = GDALDataset::FromHandle(GDALWarp(
431 : outputFilename.c_str(), nullptr, 1, &hSrcDS, psOptions, nullptr));
432 22 : GDALWarpAppOptionsFree(psOptions);
433 22 : bOK = poRetDS != nullptr;
434 22 : if (bOK)
435 : {
436 19 : m_outputDataset.Set(std::unique_ptr<GDALDataset>(poRetDS));
437 : }
438 : }
439 22 : return bOK;
440 : }
441 :
442 : GDALRasterReprojectAlgorithmStandalone::
443 : ~GDALRasterReprojectAlgorithmStandalone() = default;
444 :
445 : //! @endcond
|
