Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: GDAL
4 : * Purpose: "gdal raster calc" subcommand
5 : * Author: Daniel Baston
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2025, ISciences LLC
9 : *
10 : * SPDX-License-Identifier: MIT
11 : ****************************************************************************/
12 :
13 : #include "gdalalg_raster_calc.h"
14 :
15 : #include "../frmts/vrt/gdal_vrt.h"
16 : #include "../frmts/vrt/vrtdataset.h"
17 :
18 : #include "cpl_float.h"
19 : #include "cpl_vsi_virtual.h"
20 : #include "gdal_priv.h"
21 : #include "gdal_utils.h"
22 : #include "vrtdataset.h"
23 :
24 : #include <algorithm>
25 : #include <cmath>
26 : #include <limits>
27 : #include <optional>
28 :
29 : //! @cond Doxygen_Suppress
30 :
31 : #ifndef _
32 : #define _(x) (x)
33 : #endif
34 :
35 : constexpr const char *DEFAULT_SOURCE_NAME = "X";
36 : constexpr const char *PIPELINE_INPUT_DSN = "";
37 :
38 : struct GDALCalcOptions
39 : {
40 : GDALDataType dstType{GDT_Unknown};
41 : bool checkCRS{true};
42 : bool checkExtent{true};
43 : };
44 :
45 243 : static bool MatchIsCompleteVariableNameWithNoIndex(const std::string &str,
46 : size_t from, size_t to)
47 : {
48 243 : if (to < str.size())
49 : {
50 : // If the character after the end of the match is:
51 : // * alphanumeric or _ : we've matched only part of a variable name
52 : // * [ : we've matched a variable that already has an index
53 : // * ( : we've matched a function name
54 310 : if (std::isalnum(str[to]) || str[to] == '_' || str[to] == '[' ||
55 109 : str[to] == '(')
56 : {
57 93 : return false;
58 : }
59 : }
60 150 : if (from > 0)
61 : {
62 : // If the character before the start of the match is alphanumeric or _,
63 : // we've matched only part of a variable name.
64 91 : if (std::isalnum(str[from - 1]) || str[from - 1] == '_')
65 : {
66 3 : return false;
67 : }
68 : }
69 :
70 147 : return true;
71 : }
72 :
73 : /**
74 : * Add a band subscript to all instances of a specified variable that
75 : * do not already have such a subscript. For example, "X" would be
76 : * replaced with "X[3]" but "X[1]" would be left untouched.
77 : */
78 147 : static std::string SetBandIndices(const std::string &origExpression,
79 : const std::string &variable, int band,
80 : bool &expressionChanged)
81 : {
82 147 : std::string expression = origExpression;
83 147 : expressionChanged = false;
84 :
85 147 : std::string::size_type seekPos = 0;
86 147 : auto pos = expression.find(variable, seekPos);
87 354 : while (pos != std::string::npos)
88 : {
89 207 : auto end = pos + variable.size();
90 :
91 207 : if (MatchIsCompleteVariableNameWithNoIndex(expression, pos, end))
92 : {
93 : // No index specified for variable
94 222 : expression = expression.substr(0, pos + variable.size()) + '[' +
95 333 : std::to_string(band) + ']' + expression.substr(end);
96 111 : expressionChanged = true;
97 : }
98 :
99 207 : seekPos = end;
100 207 : pos = expression.find(variable, seekPos);
101 : }
102 :
103 147 : return expression;
104 : }
105 :
106 72 : static bool PosIsAggregateFunctionArgument(const std::string &expression,
107 : size_t pos)
108 : {
109 : // If this position is a function argument, we should be able to
110 : // scan backwards for a ( and find only variable names, literals or commas.
111 72 : while (pos != 0)
112 : {
113 64 : const char c = expression[pos];
114 64 : if (c == '(')
115 : {
116 24 : pos--;
117 24 : break;
118 : }
119 40 : if (!(isspace(c) || isalnum(c) || c == ',' || c == '.' || c == '[' ||
120 : c == ']' || c == '_'))
121 : {
122 4 : return false;
123 : }
124 36 : pos--;
125 : }
126 :
127 : // Now what we've found the (, the preceding characters should be an
128 : // aggregate function name
129 32 : if (pos < 2)
130 : {
131 8 : return false;
132 : }
133 :
134 24 : if (STARTS_WITH_CI(expression.c_str() + (pos - 2), "avg") ||
135 20 : STARTS_WITH_CI(expression.c_str() + (pos - 2), "sum") ||
136 52 : STARTS_WITH_CI(expression.c_str() + (pos - 2), "min") ||
137 8 : STARTS_WITH_CI(expression.c_str() + (pos - 2), "max"))
138 : {
139 20 : return true;
140 : }
141 :
142 4 : return false;
143 : }
144 :
145 : /**
146 : * Replace X by X[1],X[2],...X[n]
147 : */
148 : static std::string
149 32 : SetBandIndicesFlattenedExpression(const std::string &origExpression,
150 : const std::string &variable, int nBands)
151 : {
152 32 : std::string expression = origExpression;
153 :
154 32 : std::string::size_type seekPos = 0;
155 32 : auto pos = expression.find(variable, seekPos);
156 68 : while (pos != std::string::npos)
157 : {
158 36 : auto end = pos + variable.size();
159 :
160 72 : if (MatchIsCompleteVariableNameWithNoIndex(expression, pos, end) &&
161 36 : PosIsAggregateFunctionArgument(expression, pos))
162 : {
163 20 : std::string newExpr = expression.substr(0, pos);
164 68 : for (int i = 1; i <= nBands; ++i)
165 : {
166 48 : if (i > 1)
167 28 : newExpr += ',';
168 48 : newExpr += variable;
169 48 : newExpr += '[';
170 48 : newExpr += std::to_string(i);
171 48 : newExpr += ']';
172 : }
173 20 : const size_t oldExprSize = expression.size();
174 20 : newExpr += expression.substr(end);
175 20 : expression = std::move(newExpr);
176 20 : end += expression.size() - oldExprSize;
177 : }
178 :
179 36 : seekPos = end;
180 36 : pos = expression.find(variable, seekPos);
181 : }
182 :
183 32 : return expression;
184 : }
185 :
186 : struct SourceProperties
187 : {
188 : int nBands{0};
189 : int nX{0};
190 : int nY{0};
191 : bool hasGT{false};
192 : GDALGeoTransform gt{};
193 : OGRSpatialReferenceRefCountedPtr srs{};
194 : std::vector<std::optional<double>> noData{};
195 : GDALDataType eDT{GDT_Unknown};
196 : };
197 :
198 : static std::optional<SourceProperties>
199 158 : UpdateSourceProperties(SourceProperties &out, GDALDataset *ds,
200 : const GDALCalcOptions &options)
201 : {
202 316 : SourceProperties source;
203 158 : bool srsMismatch = false;
204 158 : bool extentMismatch = false;
205 158 : bool dimensionMismatch = false;
206 :
207 : {
208 158 : source.nBands = ds->GetRasterCount();
209 158 : source.nX = ds->GetRasterXSize();
210 158 : source.nY = ds->GetRasterYSize();
211 158 : source.noData.resize(source.nBands);
212 :
213 158 : if (options.checkExtent)
214 : {
215 152 : ds->GetGeoTransform(source.gt);
216 : }
217 :
218 158 : if (options.checkCRS && out.srs)
219 : {
220 59 : const OGRSpatialReference *srs = ds->GetSpatialRef();
221 59 : srsMismatch = srs && !srs->IsSame(out.srs.get());
222 : }
223 :
224 : // Store the source data type if it is the same for all bands in the source
225 158 : bool bandsHaveSameType = true;
226 418 : for (int i = 1; i <= source.nBands; ++i)
227 : {
228 260 : GDALRasterBand *band = ds->GetRasterBand(i);
229 :
230 260 : if (i == 1)
231 : {
232 158 : source.eDT = band->GetRasterDataType();
233 : }
234 204 : else if (bandsHaveSameType &&
235 102 : source.eDT != band->GetRasterDataType())
236 : {
237 0 : source.eDT = GDT_Unknown;
238 0 : bandsHaveSameType = false;
239 : }
240 :
241 : int success;
242 260 : double noData = band->GetNoDataValue(&success);
243 260 : if (success)
244 : {
245 17 : source.noData[i - 1] = noData;
246 : }
247 : }
248 : }
249 :
250 158 : if (source.nX != out.nX || source.nY != out.nY)
251 : {
252 3 : dimensionMismatch = true;
253 : }
254 :
255 158 : if (source.gt.xorig != out.gt.xorig || source.gt.xrot != out.gt.xrot ||
256 158 : source.gt.yorig != out.gt.yorig || source.gt.yrot != out.gt.yrot)
257 : {
258 6 : extentMismatch = true;
259 : }
260 158 : if (source.gt.xscale != out.gt.xscale || source.gt.yscale != out.gt.yscale)
261 : {
262 : // Resolutions are different. Are the extents the same?
263 9 : double xmaxOut =
264 9 : out.gt.xorig + out.nX * out.gt.xscale + out.nY * out.gt.xrot;
265 9 : double yminOut =
266 9 : out.gt.yorig + out.nX * out.gt.yrot + out.nY * out.gt.yscale;
267 :
268 9 : double xmax = source.gt.xorig + source.nX * source.gt.xscale +
269 9 : source.nY * source.gt.xrot;
270 9 : double ymin = source.gt.yorig + source.nX * source.gt.yrot +
271 9 : source.nY * source.gt.yscale;
272 :
273 : // Max allowable extent misalignment, expressed as fraction of a pixel
274 9 : constexpr double EXTENT_RTOL = 1e-3;
275 :
276 9 : if (std::abs(xmax - xmaxOut) >
277 15 : EXTENT_RTOL * std::abs(source.gt.xscale) ||
278 6 : std::abs(ymin - yminOut) > EXTENT_RTOL * std::abs(source.gt.yscale))
279 : {
280 6 : extentMismatch = true;
281 : }
282 : }
283 :
284 158 : if (options.checkExtent && extentMismatch)
285 : {
286 2 : CPLError(CE_Failure, CPLE_AppDefined,
287 : "Input extents are inconsistent.");
288 2 : return std::nullopt;
289 : }
290 :
291 156 : if (!options.checkExtent && dimensionMismatch)
292 : {
293 1 : CPLError(CE_Failure, CPLE_AppDefined,
294 : "Inputs do not have the same dimensions.");
295 1 : return std::nullopt;
296 : }
297 :
298 : // Find a common resolution
299 155 : if (source.nX > out.nX)
300 : {
301 1 : auto dx = CPLGreatestCommonDivisor(out.gt.xscale, source.gt.xscale);
302 1 : if (std::fabs(dx) < std::numeric_limits<double>::min())
303 : {
304 0 : CPLError(CE_Failure, CPLE_AppDefined,
305 : "Failed to find common resolution for inputs.");
306 0 : return std::nullopt;
307 : }
308 1 : out.nX = static_cast<int>(
309 1 : std::round(static_cast<double>(out.nX) * out.gt.xscale / dx));
310 1 : out.gt.xscale = dx;
311 : }
312 155 : if (source.nY > out.nY)
313 : {
314 1 : auto dy = CPLGreatestCommonDivisor(out.gt.yscale, source.gt.yscale);
315 1 : if (std::fabs(dy) < std::numeric_limits<double>::min())
316 : {
317 0 : CPLError(CE_Failure, CPLE_AppDefined,
318 : "Failed to find common resolution for inputs.");
319 0 : return std::nullopt;
320 : }
321 1 : out.nY = static_cast<int>(
322 1 : std::round(static_cast<double>(out.nY) * out.gt.yscale / dy));
323 1 : out.gt.yscale = dy;
324 : }
325 :
326 155 : if (srsMismatch)
327 : {
328 1 : CPLError(CE_Failure, CPLE_AppDefined,
329 : "Input spatial reference systems are inconsistent.");
330 1 : return std::nullopt;
331 : }
332 :
333 154 : return source;
334 : }
335 :
336 : /** Add one or more derived bands to a VRTDataset, representing the evaluation
337 : * of a single expression
338 : *
339 : * @param poDS VRT dataset
340 : * @param bandType the type of the band(s) to create
341 : * @param expression Expression for which band(s) should be added
342 : * @param dialect Expression dialect
343 : * @param flatten Generate a single band output raster per expression, even if
344 : * input datasets are multiband.
345 : * @param noDataText nodata value to use for the created band, or "none", or ""
346 : * @param pixelFunctionArguments Pixel function arguments.
347 : * @param sources Mapping of source names to DSNs
348 : * @param sourceProps Mapping of source names to properties
349 : * @param fakeSourceFilename If not empty, used instead of real input filenames.
350 : * @param pipelineInputSource A pointer to a dataset representing pipeline input.
351 : * @return true if the band(s) were added, false otherwise
352 : */
353 118 : static bool CreateVRTDerivedBand(
354 : VRTDataset *poDS, GDALDataType bandType, const std::string &expression,
355 : const std::string &dialect, bool flatten, const std::string &noDataText,
356 : const std::vector<std::string> &pixelFunctionArguments,
357 : const std::map<std::string, std::string> &sources,
358 : const std::map<std::string, SourceProperties> &sourceProps,
359 : const std::string &fakeSourceFilename, GDALDataset *pipelineInputSource)
360 : {
361 118 : const char *pszVRTFilename = poDS->GetDescription();
362 :
363 118 : const int nPrevBands = poDS->GetRasterCount();
364 118 : const int nXOut = poDS->GetRasterXSize();
365 118 : const int nYOut = poDS->GetRasterYSize();
366 :
367 236 : CPLStringList papszBandArgs;
368 118 : papszBandArgs.SetNameValue("subclass", "VRTDerivedRasterBand");
369 :
370 118 : int nOutBands = 1; // By default, each expression produces a single output
371 : // band. When processing the expression below, we may
372 : // discover that the expression produces multiple bands,
373 : // in which case this will be updated.
374 :
375 259 : for (int nOutBand = 1; nOutBand <= nOutBands; nOutBand++)
376 : {
377 : // Copy the expression for each output band, because we may modify it
378 : // when adding band indices (e.g., X -> X[1]) to the variables in the
379 : // expression.
380 145 : std::string bandExpression = expression;
381 :
382 145 : if (poDS->AddBand(bandType == GDT_Unknown ? GDT_Float64 : bandType,
383 290 : papszBandArgs) != CE_None)
384 : {
385 0 : return false;
386 : }
387 145 : VRTDerivedRasterBand *poBand = cpl::down_cast<VRTDerivedRasterBand *>(
388 : poDS->GetRasterBand(nPrevBands + nOutBand));
389 :
390 145 : std::optional<double> dstNoData;
391 145 : bool autoSelectNoDataValue = false;
392 145 : if (noDataText.empty())
393 : {
394 140 : autoSelectNoDataValue = true;
395 : }
396 5 : else if (noDataText != "none")
397 : {
398 10 : if (auto parsed = cpl::strict_parse<double>(noDataText);
399 5 : parsed.has_value())
400 : {
401 5 : dstNoData = parsed.value();
402 : }
403 : else
404 : {
405 0 : CPLError(CE_Failure, CPLE_AppDefined,
406 : "Invalid NoData value: %s", noDataText.c_str());
407 0 : return false;
408 : }
409 : }
410 :
411 331 : for (const auto &[source_name, dsn] : sources)
412 : {
413 190 : auto it = sourceProps.find(source_name);
414 190 : CPLAssert(it != sourceProps.end());
415 190 : const auto &props = it->second;
416 :
417 190 : bool expressionAppliedPerBand = false;
418 190 : if (dialect == "builtin")
419 : {
420 43 : expressionAppliedPerBand = !flatten;
421 : }
422 : else
423 : {
424 147 : const int nDefaultInBand = std::min(props.nBands, nOutBand);
425 :
426 147 : if (flatten)
427 : {
428 32 : bandExpression = SetBandIndicesFlattenedExpression(
429 32 : bandExpression, source_name, props.nBands);
430 : }
431 :
432 : bandExpression =
433 294 : SetBandIndices(bandExpression, source_name, nDefaultInBand,
434 147 : expressionAppliedPerBand);
435 : }
436 :
437 190 : if (expressionAppliedPerBand)
438 : {
439 135 : if (nOutBands <= 1)
440 : {
441 94 : nOutBands = props.nBands;
442 : }
443 41 : else if (props.nBands != 1 && props.nBands != nOutBands)
444 : {
445 3 : CPLError(CE_Failure, CPLE_AppDefined,
446 : "Expression cannot operate on all bands of "
447 : "rasters with incompatible numbers of bands "
448 : "(source %s has %d bands but expected to have "
449 : "1 or %d bands).",
450 3 : source_name.c_str(), props.nBands, nOutBands);
451 4 : return false;
452 : }
453 : }
454 :
455 : // Create a source for each input band that is used in
456 : // the expression.
457 519 : for (int nInBand = 1; nInBand <= props.nBands; nInBand++)
458 : {
459 332 : CPLString inBandVariable;
460 332 : if (dialect == "builtin")
461 : {
462 75 : if (!flatten && props.nBands >= 2 && nInBand != nOutBand)
463 11 : continue;
464 : }
465 : else
466 : {
467 : inBandVariable.Printf("%s[%d]", source_name.c_str(),
468 257 : nInBand);
469 257 : if (bandExpression.find(inBandVariable) ==
470 : std::string::npos)
471 : {
472 79 : continue;
473 : }
474 : }
475 :
476 : const std::optional<double> &srcNoData =
477 242 : props.noData[nInBand - 1];
478 :
479 0 : std::unique_ptr<VRTSimpleSource> poSource;
480 242 : if (srcNoData.has_value())
481 : {
482 17 : poSource = std::make_unique<VRTComplexSource>();
483 : }
484 : else
485 : {
486 225 : poSource = std::make_unique<VRTSimpleSource>();
487 : }
488 :
489 242 : if (!inBandVariable.empty())
490 : {
491 178 : poSource->SetName(inBandVariable);
492 : }
493 :
494 242 : if (fakeSourceFilename.empty())
495 : {
496 169 : if (dsn == PIPELINE_INPUT_DSN)
497 : {
498 4 : CPLAssertNotNull(pipelineInputSource);
499 4 : pipelineInputSource->Reference();
500 4 : poSource->SetSrcBand(
501 : pipelineInputSource->GetRasterBand(nInBand));
502 : }
503 : else
504 : {
505 330 : std::string osSourceFilename = dsn;
506 165 : bool bRelativeToVRT = false;
507 165 : if (pszVRTFilename[0])
508 : {
509 0 : std::tie(osSourceFilename, bRelativeToVRT) =
510 : VRTSimpleSource::
511 0 : ComputeSourceNameAndRelativeFlag(
512 0 : CPLGetPathSafe(pszVRTFilename).c_str(),
513 0 : dsn);
514 : }
515 165 : poSource->SetSrcBand(osSourceFilename.c_str(), nInBand);
516 : }
517 : }
518 : else
519 : {
520 73 : poSource->SetSrcBand(fakeSourceFilename.c_str(), nInBand);
521 : }
522 :
523 242 : if (srcNoData.has_value())
524 : {
525 : cpl::down_cast<VRTComplexSource *>(poSource.get())
526 17 : ->SetNoDataValue(srcNoData.value());
527 :
528 17 : if (autoSelectNoDataValue && !dstNoData.has_value())
529 : {
530 8 : dstNoData = srcNoData;
531 : }
532 : }
533 :
534 242 : if (fakeSourceFilename.empty())
535 : {
536 169 : poSource->SetSrcWindow(0, 0, props.nX, props.nY);
537 169 : poSource->SetDstWindow(0, 0, nXOut, nYOut);
538 : }
539 :
540 242 : poBand->AddSource(std::move(poSource));
541 : }
542 :
543 187 : if (dstNoData.has_value())
544 : {
545 17 : if (!GDALIsValueExactAs(dstNoData.value(), bandType))
546 : {
547 1 : CPLError(
548 : CE_Failure, CPLE_AppDefined,
549 : "Band output type %s cannot represent NoData value %g",
550 1 : GDALGetDataTypeName(bandType), dstNoData.value());
551 1 : return false;
552 : }
553 :
554 16 : poBand->SetNoDataValue(dstNoData.value());
555 : }
556 : }
557 :
558 141 : if (dialect == "builtin")
559 : {
560 29 : poBand->SetPixelFunctionName(expression.c_str());
561 : }
562 : else
563 : {
564 112 : poBand->SetPixelFunctionName("expression");
565 112 : poBand->AddPixelFunctionArgument("dialect", "muparser");
566 : // Add the expression as a last step, because we may modify the
567 : // expression as we iterate through the bands.
568 112 : poBand->AddPixelFunctionArgument("expression",
569 : bandExpression.c_str());
570 : }
571 :
572 141 : if (!pixelFunctionArguments.empty())
573 : {
574 16 : const CPLStringList args(pixelFunctionArguments);
575 16 : for (const auto &[key, value] : cpl::IterateNameValue(args))
576 : {
577 8 : poBand->AddPixelFunctionArgument(key, value);
578 : }
579 : }
580 : }
581 :
582 114 : return true;
583 : }
584 :
585 122 : static bool ParseSourceDescriptors(const std::vector<std::string> &inputs,
586 : std::map<std::string, std::string> &datasets,
587 : std::string &firstSourceName,
588 : bool requireSourceNames)
589 : {
590 281 : for (size_t iInput = 0; iInput < inputs.size(); iInput++)
591 : {
592 164 : const std::string &input = inputs[iInput];
593 164 : std::string name;
594 :
595 164 : const auto pos = input.find('=');
596 164 : if (pos == std::string::npos)
597 : {
598 61 : if (requireSourceNames && inputs.size() > 1)
599 : {
600 1 : CPLError(CE_Failure, CPLE_AppDefined,
601 : "Inputs must be named when more than one input is "
602 : "provided.");
603 1 : return false;
604 : }
605 60 : name = DEFAULT_SOURCE_NAME;
606 60 : if (iInput > 0)
607 : {
608 2 : name += std::to_string(iInput);
609 : }
610 : }
611 : else
612 : {
613 103 : name = input.substr(0, pos);
614 : }
615 :
616 : // Check input name is legal
617 347 : for (size_t i = 0; i < name.size(); ++i)
618 : {
619 187 : const char c = name[i];
620 187 : if ((c >= 'a' && c <= 'z') || (c >= 'A' && c <= 'Z'))
621 : {
622 : // ok
623 : }
624 20 : else if (c == '_' || (c >= '0' && c <= '9'))
625 : {
626 19 : if (i == 0)
627 : {
628 : // Reserved constants in MuParser start with an underscore
629 2 : CPLError(
630 : CE_Failure, CPLE_AppDefined,
631 : "Name '%s' is illegal because it starts with a '%c'",
632 : name.c_str(), c);
633 2 : return false;
634 : }
635 : }
636 : else
637 : {
638 1 : CPLError(CE_Failure, CPLE_AppDefined,
639 : "Name '%s' is illegal because character '%c' is not "
640 : "allowed",
641 : name.c_str(), c);
642 1 : return false;
643 : }
644 : }
645 :
646 : std::string dsn =
647 160 : (pos == std::string::npos) ? input : input.substr(pos + 1);
648 :
649 160 : if (!dsn.empty() && dsn.front() == '[' && dsn.back() == ']')
650 : {
651 : dsn = "{\"type\":\"gdal_streamed_alg\", \"command_line\":\"gdal "
652 0 : "raster pipeline " +
653 2 : CPLString(dsn.substr(1, dsn.size() - 2))
654 2 : .replaceAll('\\', "\\\\")
655 2 : .replaceAll('"', "\\\"") +
656 1 : "\"}";
657 : }
658 :
659 160 : if (datasets.find(name) != datasets.end())
660 : {
661 1 : CPLError(CE_Failure, CPLE_AppDefined,
662 : "An input with name '%s' has already been provided",
663 : name.c_str());
664 1 : return false;
665 : }
666 159 : datasets[name] = std::move(dsn);
667 :
668 159 : if (iInput == 0)
669 : {
670 118 : firstSourceName = std::move(name);
671 : }
672 : }
673 :
674 117 : return true;
675 : }
676 :
677 90 : static bool ReadFileLists(const std::vector<GDALArgDatasetValue> &inputDS,
678 : std::vector<std::string> &inputFilenames)
679 : {
680 212 : for (const auto &dsVal : inputDS)
681 : {
682 122 : const auto &input = dsVal.GetName();
683 122 : if (!input.empty() && input[0] == '@')
684 : {
685 : auto f =
686 2 : VSIVirtualHandleUniquePtr(VSIFOpenL(input.c_str() + 1, "r"));
687 2 : if (!f)
688 : {
689 0 : CPLError(CE_Failure, CPLE_FileIO, "Cannot open %s",
690 0 : input.c_str() + 1);
691 0 : return false;
692 : }
693 6 : while (const char *filename = CPLReadLineL(f.get()))
694 : {
695 4 : inputFilenames.push_back(filename);
696 4 : }
697 : }
698 : else
699 : {
700 120 : inputFilenames.push_back(input);
701 : }
702 : }
703 :
704 90 : return true;
705 : }
706 :
707 : /** Creates a VRT dataset with one or more derived raster bands containing
708 : * results of an expression.
709 : *
710 : * To make this work with muparser (which does not support vector types), we
711 : * do a simple parsing of the expression internally, transforming it into
712 : * multiple expressions with explicit band indices. For example, for a two-band
713 : * raster "X", the expression "X + 3" will be transformed into "X[1] + 3" and
714 : * "X[2] + 3". The use of brackets is for readability only; as far as the
715 : * expression engine is concerned, the variables "X[1]" and "X[2]" have nothing
716 : * to do with each other.
717 : *
718 : * @param inputs Either:
719 : * - a list of sources, expressed as NAME=DSN
720 : * - pointer to a single opened dataset
721 : * @param expressions A list of expressions to be evaluated
722 : * @param dialect Expression dialect
723 : * @param flatten Generate a single band output raster per expression, even if
724 : * input datasets are multiband.
725 : * @param noData NoData values to use for output bands, or "none", or ""
726 : * @param pixelFunctionArguments Pixel function arguments.
727 : * @param options flags controlling which checks should be performed on the inputs
728 : * @param[out] maxSourceBands Maximum number of bands in source dataset(s)
729 : * @param fakeSourceFilename If not empty, used instead of real input filenames.
730 : *
731 : * @return a newly created VRTDataset, or nullptr on error
732 : */
733 126 : static std::unique_ptr<GDALDataset> GDALCalcCreateVRTDerived(
734 : std::variant<GDALDataset *, const std::vector<std::string> *> inputs,
735 : const std::vector<std::string> &expressions, const std::string &dialect,
736 : bool flatten, const std::string &noData,
737 : const std::vector<std::vector<std::string>> &pixelFunctionArguments,
738 : const GDALCalcOptions &options, int &maxSourceBands,
739 : const std::string &fakeSourceFilename = std::string())
740 : {
741 252 : std::map<std::string, std::string> sources;
742 252 : std::map<std::string, SourceProperties> sourceProps;
743 126 : GDALDataset *pipelineInputDS = std::holds_alternative<GDALDataset *>(inputs)
744 126 : ? std::get<GDALDataset *>(inputs)
745 126 : : nullptr;
746 :
747 126 : maxSourceBands = 0;
748 :
749 : // Read properties from the first source
750 252 : SourceProperties out;
751 : {
752 0 : std::unique_ptr<GDALDataset> poTmpDS;
753 : const GDALDataset *poTemplateDS;
754 :
755 126 : if (pipelineInputDS)
756 : {
757 2 : poTemplateDS = pipelineInputDS;
758 : }
759 : else
760 : {
761 : const std::vector<std::string> &sourceDescriptors =
762 124 : *std::get<const std::vector<std::string> *>(inputs);
763 :
764 124 : if (sourceDescriptors.empty())
765 : {
766 7 : return nullptr;
767 : }
768 :
769 122 : const bool requireSourceNames = dialect != "builtin";
770 :
771 122 : std::string firstSource;
772 122 : if (!ParseSourceDescriptors(sourceDescriptors, sources, firstSource,
773 : requireSourceNames))
774 : {
775 5 : return nullptr;
776 : }
777 :
778 : // Use the first source provided to determine properties of the output
779 117 : const char *firstDSN = sources[firstSource].c_str();
780 :
781 117 : poTmpDS.reset(GDALDataset::Open(firstDSN, GDAL_OF_RASTER));
782 117 : if (!poTmpDS)
783 : {
784 0 : CPLError(CE_Failure, CPLE_AppDefined, "Failed to open %s",
785 : firstDSN);
786 0 : return nullptr;
787 : }
788 117 : poTemplateDS = poTmpDS.get();
789 : }
790 :
791 119 : out.nX = poTemplateDS->GetRasterXSize();
792 119 : out.nY = poTemplateDS->GetRasterYSize();
793 119 : out.nBands = 1;
794 238 : out.srs = OGRSpatialReferenceRefCountedPtr::makeClone(
795 238 : poTemplateDS->GetSpatialRef());
796 119 : out.hasGT = poTemplateDS->GetGeoTransform(out.gt) == CE_None;
797 :
798 119 : maxSourceBands = 0;
799 :
800 119 : if (pipelineInputDS)
801 : {
802 2 : sources[DEFAULT_SOURCE_NAME] = PIPELINE_INPUT_DSN;
803 2 : if (auto props =
804 2 : UpdateSourceProperties(out, pipelineInputDS, options))
805 : {
806 2 : sourceProps[DEFAULT_SOURCE_NAME] = props.value();
807 2 : maxSourceBands = props.value().nBands;
808 : }
809 : else
810 : {
811 0 : return nullptr; // error message emitted from UpdateSourceProperties
812 : }
813 : }
814 : else
815 : {
816 : // Collect properties of the different sources, and verify them for
817 : // consistency.
818 269 : for (const auto &[source_name, dsn] : sources)
819 : {
820 : // TODO avoid opening the first source twice.
821 : std::unique_ptr<GDALDataset> ds(
822 156 : GDALDataset::Open(dsn.c_str(), GDAL_OF_RASTER));
823 :
824 156 : if (!ds)
825 : {
826 0 : CPLError(CE_Failure, CPLE_AppDefined, "Failed to open %s",
827 : dsn.c_str());
828 0 : return nullptr;
829 : }
830 :
831 156 : auto props = UpdateSourceProperties(out, ds.get(), options);
832 156 : if (props.has_value())
833 : {
834 152 : maxSourceBands = std::max(maxSourceBands, props->nBands);
835 152 : sourceProps[source_name] = std::move(props.value());
836 : }
837 : else
838 : {
839 4 : return nullptr; // error message emitted from UpdateSourceProperties
840 : }
841 : }
842 : }
843 : }
844 :
845 115 : size_t iExpr = 0;
846 :
847 : auto poDS = VRTDataset::CreateVRTDataset("", out.nX, out.nY, 0,
848 230 : options.dstType, nullptr);
849 :
850 229 : for (const auto &origExpression : expressions)
851 : {
852 118 : GDALDataType bandType = options.dstType;
853 :
854 : // If output band type has not been specified, set it equal to the
855 : // input band type for certain pixel functions, if the inputs have
856 : // a consistent band type.
857 170 : if (bandType == GDT_Unknown && dialect == "builtin" &&
858 75 : (origExpression == "min" || origExpression == "max" ||
859 23 : origExpression == "mode"))
860 : {
861 12 : for (const auto &[_, props] : sourceProps)
862 : {
863 6 : if (bandType == GDT_Unknown)
864 : {
865 6 : bandType = props.eDT;
866 : }
867 0 : else if (props.eDT == GDT_Unknown || props.eDT != bandType)
868 : {
869 0 : bandType = GDT_Unknown;
870 0 : break;
871 : }
872 : }
873 : }
874 :
875 118 : if (!CreateVRTDerivedBand(
876 : poDS.get(), bandType, origExpression, dialect, flatten, noData,
877 118 : pixelFunctionArguments[iExpr], sources, sourceProps,
878 : fakeSourceFilename, pipelineInputDS))
879 : {
880 4 : return nullptr;
881 : }
882 114 : ++iExpr;
883 : }
884 :
885 111 : if (out.hasGT)
886 : {
887 59 : poDS->SetGeoTransform(out.gt);
888 : }
889 111 : if (out.srs)
890 : {
891 57 : poDS->SetSpatialRef(out.srs.get());
892 : }
893 :
894 111 : return poDS;
895 : }
896 :
897 : /************************************************************************/
898 : /* GDALRasterCalcAlgorithm::GDALRasterCalcAlgorithm() */
899 : /************************************************************************/
900 :
901 180 : GDALRasterCalcAlgorithm::GDALRasterCalcAlgorithm(bool standaloneStep) noexcept
902 : : GDALRasterPipelineStepAlgorithm(NAME, DESCRIPTION, HELP_URL,
903 540 : ConstructorOptions()
904 180 : .SetStandaloneStep(standaloneStep)
905 180 : .SetAddDefaultArguments(false)
906 180 : .SetAutoOpenInputDatasets(false)
907 360 : .SetInputDatasetMetaVar("INPUTS")
908 540 : .SetInputDatasetMaxCount(INT_MAX))
909 : {
910 180 : AddRasterInputArgs(false, false);
911 180 : if (standaloneStep)
912 : {
913 138 : AddProgressArg();
914 138 : AddRasterOutputArgs(false);
915 : }
916 :
917 180 : AddOutputDataTypeArg(&m_type);
918 :
919 : AddArg("no-check-crs", 0,
920 : _("Do not check consistency of input coordinate reference systems"),
921 360 : &m_noCheckCRS)
922 180 : .AddHiddenAlias("no-check-srs");
923 : AddArg("no-check-extent", 0, _("Do not check consistency of input extents"),
924 180 : &m_noCheckExtent);
925 :
926 : AddArg("propagate-nodata", 0,
927 : _("Whether to set pixels to the output NoData value if any of the "
928 : "input pixels is NoData"),
929 180 : &m_propagateNoData);
930 :
931 360 : AddArg("calc", 0, _("Expression(s) to evaluate"), &m_expr)
932 180 : .SetRequired()
933 180 : .SetPackedValuesAllowed(false)
934 180 : .SetMinCount(1)
935 : .SetAutoCompleteFunction(
936 4 : [this](const std::string ¤tValue)
937 : {
938 4 : std::vector<std::string> ret;
939 2 : if (m_dialect == "builtin")
940 : {
941 1 : if (currentValue.find('(') == std::string::npos)
942 1 : return VRTDerivedRasterBand::GetPixelFunctionNames();
943 : }
944 1 : return ret;
945 180 : });
946 :
947 360 : AddArg("dialect", 0, _("Expression dialect"), &m_dialect)
948 180 : .SetDefault(m_dialect)
949 180 : .SetChoices("muparser", "builtin");
950 :
951 : AddArg("flatten", 0,
952 : _("Generate a single band output raster per expression, even if "
953 : "input datasets are multiband"),
954 180 : &m_flatten);
955 :
956 180 : AddNodataArg(&m_nodata, true);
957 :
958 : // This is a hidden option only used by test_gdalalg_raster_calc_expression_rewriting()
959 : // for now
960 : AddArg("no-check-expression", 0,
961 : _("Whether to skip expression validity checks for virtual format "
962 : "output"),
963 360 : &m_noCheckExpression)
964 180 : .SetHidden();
965 :
966 180 : AddValidationAction(
967 188 : [this]()
968 : {
969 99 : GDALPipelineStepRunContext ctxt;
970 99 : return m_noCheckExpression || !IsGDALGOutput() || RunStep(ctxt);
971 : });
972 180 : }
973 :
974 : /************************************************************************/
975 : /* GDALRasterCalcAlgorithm::RunImpl() */
976 : /************************************************************************/
977 :
978 85 : bool GDALRasterCalcAlgorithm::RunImpl(GDALProgressFunc pfnProgress,
979 : void *pProgressData)
980 : {
981 85 : GDALPipelineStepRunContext stepCtxt;
982 85 : stepCtxt.m_pfnProgress = pfnProgress;
983 85 : stepCtxt.m_pProgressData = pProgressData;
984 85 : return RunPreStepPipelineValidations() && RunStep(stepCtxt);
985 : }
986 :
987 : /************************************************************************/
988 : /* GDALRasterCalcAlgorithm::RunStep() */
989 : /************************************************************************/
990 :
991 92 : bool GDALRasterCalcAlgorithm::RunStep(GDALPipelineStepRunContext &ctxt)
992 : {
993 92 : CPLAssert(!m_outputDataset.GetDatasetRef());
994 :
995 92 : GDALCalcOptions options;
996 92 : options.checkExtent = !m_noCheckExtent;
997 92 : options.checkCRS = !m_noCheckCRS;
998 92 : if (!m_type.empty())
999 : {
1000 5 : options.dstType = GDALGetDataTypeByName(m_type.c_str());
1001 : }
1002 :
1003 92 : GDALDataset *poPipelineInput = nullptr;
1004 184 : std::vector<std::string> inputFilenames;
1005 92 : if (m_inputDataset.size() == 1 && m_inputDataset[0].GetDatasetRef())
1006 : {
1007 2 : poPipelineInput = m_inputDataset[0].GetDatasetRef();
1008 : }
1009 : else
1010 : {
1011 90 : if (!ReadFileLists(m_inputDataset, inputFilenames))
1012 : {
1013 0 : return false;
1014 : }
1015 : }
1016 :
1017 184 : std::vector<std::vector<std::string>> pixelFunctionArgs;
1018 92 : if (m_dialect == "builtin")
1019 : {
1020 29 : for (std::string &expr : m_expr)
1021 : {
1022 : const CPLStringList aosTokens(
1023 : CSLTokenizeString2(expr.c_str(), "()",
1024 15 : CSLT_STRIPLEADSPACES | CSLT_STRIPENDSPACES));
1025 15 : const char *pszFunction = aosTokens[0];
1026 : const auto *pair =
1027 15 : VRTDerivedRasterBand::GetPixelFunction(pszFunction);
1028 15 : if (!pair)
1029 : {
1030 0 : ReportError(CE_Failure, CPLE_NotSupported,
1031 : "'%s' is a unknown builtin function", pszFunction);
1032 0 : return false;
1033 : }
1034 15 : if (aosTokens.size() == 2)
1035 : {
1036 2 : std::vector<std::string> validArguments;
1037 2 : AddOptionsSuggestions(pair->second.c_str(), 0, std::string(),
1038 : validArguments);
1039 6 : for (std::string &s : validArguments)
1040 : {
1041 4 : if (!s.empty() && s.back() == '=')
1042 4 : s.pop_back();
1043 : }
1044 :
1045 : const CPLStringList aosTokensArgs(CSLTokenizeString2(
1046 : aosTokens[1], ",",
1047 2 : CSLT_STRIPLEADSPACES | CSLT_STRIPENDSPACES));
1048 4 : for (const auto &[key, value] :
1049 6 : cpl::IterateNameValue(aosTokensArgs))
1050 : {
1051 2 : if (std::find(validArguments.begin(), validArguments.end(),
1052 2 : key) == validArguments.end())
1053 : {
1054 0 : if (validArguments.empty())
1055 : {
1056 0 : ReportError(
1057 : CE_Failure, CPLE_IllegalArg,
1058 : "'%s' is a unrecognized argument for builtin "
1059 : "function '%s'. It does not accept any "
1060 : "argument",
1061 : key, pszFunction);
1062 : }
1063 : else
1064 : {
1065 0 : std::string validArgumentsStr;
1066 0 : for (const std::string &s : validArguments)
1067 : {
1068 0 : if (!validArgumentsStr.empty())
1069 0 : validArgumentsStr += ", ";
1070 0 : validArgumentsStr += '\'';
1071 0 : validArgumentsStr += s;
1072 0 : validArgumentsStr += '\'';
1073 : }
1074 0 : ReportError(
1075 : CE_Failure, CPLE_IllegalArg,
1076 : "'%s' is a unrecognized argument for builtin "
1077 : "function '%s'. Only %s %s supported",
1078 : key, pszFunction,
1079 0 : validArguments.size() == 1 ? "is" : "are",
1080 : validArgumentsStr.c_str());
1081 : }
1082 0 : return false;
1083 : }
1084 2 : CPL_IGNORE_RET_VAL(value);
1085 : }
1086 2 : pixelFunctionArgs.emplace_back(aosTokensArgs);
1087 : }
1088 : else
1089 : {
1090 13 : pixelFunctionArgs.push_back(std::vector<std::string>());
1091 : }
1092 15 : expr = pszFunction;
1093 : }
1094 : }
1095 : else
1096 : {
1097 78 : pixelFunctionArgs.resize(m_expr.size());
1098 : }
1099 :
1100 92 : if (m_propagateNoData)
1101 : {
1102 2 : if (m_nodata == "none")
1103 : {
1104 0 : ReportError(CE_Failure, CPLE_AppDefined,
1105 : "Output NoData value must be specified to use "
1106 : "--propagate-nodata");
1107 0 : return false;
1108 : }
1109 4 : for (auto &args : pixelFunctionArgs)
1110 : {
1111 2 : args.push_back("propagateNoData=1");
1112 : }
1113 : }
1114 :
1115 92 : int maxSourceBands = 0;
1116 : const bool bIsVRT =
1117 236 : m_format == "VRT" ||
1118 91 : (m_format.empty() &&
1119 106 : EQUAL(CPLGetExtensionSafe(m_outputDataset.GetName().c_str()).c_str(),
1120 92 : "VRT"));
1121 :
1122 92 : std::variant<GDALDataset *, const std::vector<std::string> *> inputs;
1123 92 : if (poPipelineInput)
1124 : {
1125 2 : inputs = poPipelineInput;
1126 : }
1127 : else
1128 : {
1129 90 : inputs = &inputFilenames;
1130 : }
1131 :
1132 : auto vrt =
1133 92 : GDALCalcCreateVRTDerived(inputs, m_expr, m_dialect, m_flatten, m_nodata,
1134 184 : pixelFunctionArgs, options, maxSourceBands);
1135 92 : if (vrt == nullptr)
1136 : {
1137 13 : return false;
1138 : }
1139 :
1140 79 : if (!m_noCheckExpression)
1141 : {
1142 : const bool bIsGDALG =
1143 163 : m_format == "GDALG" ||
1144 65 : (m_format.empty() &&
1145 32 : cpl::ends_with(m_outputDataset.GetName(), ".gdalg.json"));
1146 66 : if (!m_standaloneStep || m_format == "stream" || bIsVRT || bIsGDALG)
1147 : {
1148 : // Try reading a single pixel to check formulas are valid.
1149 34 : std::vector<GByte> dummyData(vrt->GetRasterCount());
1150 :
1151 34 : auto poGTIFFDrv = GetGDALDriverManager()->GetDriverByName("GTiff");
1152 34 : std::string osTmpFilename;
1153 34 : if (poGTIFFDrv)
1154 : {
1155 : std::string osFilename =
1156 68 : VSIMemGenerateHiddenFilename("tmp.tif");
1157 : auto poDS = std::unique_ptr<GDALDataset>(
1158 : poGTIFFDrv->Create(osFilename.c_str(), 1, 1, maxSourceBands,
1159 68 : GDT_UInt8, nullptr));
1160 34 : if (poDS)
1161 34 : osTmpFilename = std::move(osFilename);
1162 : }
1163 34 : if (!osTmpFilename.empty())
1164 : {
1165 : auto fakeVRT = GDALCalcCreateVRTDerived(
1166 0 : &inputFilenames, m_expr, m_dialect, m_flatten, m_nodata,
1167 34 : pixelFunctionArgs, options, maxSourceBands, osTmpFilename);
1168 66 : if (fakeVRT &&
1169 32 : fakeVRT->RasterIO(GF_Read, 0, 0, 1, 1, dummyData.data(), 1,
1170 : 1, GDT_UInt8, vrt->GetRasterCount(),
1171 34 : nullptr, 0, 0, 0, nullptr) != CE_None)
1172 : {
1173 5 : return false;
1174 : }
1175 : }
1176 29 : if (bIsGDALG)
1177 : {
1178 1 : return true;
1179 : }
1180 : }
1181 : }
1182 :
1183 73 : if (m_format == "stream" || !m_standaloneStep)
1184 : {
1185 26 : m_outputDataset.Set(std::move(vrt));
1186 26 : return true;
1187 : }
1188 :
1189 47 : CPLStringList translateArgs;
1190 47 : if (!m_format.empty())
1191 : {
1192 9 : translateArgs.AddString("-of");
1193 9 : translateArgs.AddString(m_format.c_str());
1194 : }
1195 48 : for (const auto &co : m_creationOptions)
1196 : {
1197 1 : translateArgs.AddString("-co");
1198 1 : translateArgs.AddString(co.c_str());
1199 : }
1200 :
1201 47 : bool bOK = false;
1202 : GDALTranslateOptions *translateOptions =
1203 47 : GDALTranslateOptionsNew(translateArgs.List(), nullptr);
1204 47 : if (translateOptions)
1205 : {
1206 47 : GDALTranslateOptionsSetProgress(translateOptions, ctxt.m_pfnProgress,
1207 : ctxt.m_pProgressData);
1208 :
1209 : auto poOutDS =
1210 : std::unique_ptr<GDALDataset>(GDALDataset::FromHandle(GDALTranslate(
1211 47 : m_outputDataset.GetName().c_str(),
1212 94 : GDALDataset::ToHandle(vrt.get()), translateOptions, nullptr)));
1213 47 : GDALTranslateOptionsFree(translateOptions);
1214 :
1215 47 : bOK = poOutDS != nullptr;
1216 47 : m_outputDataset.Set(std::move(poOutDS));
1217 : }
1218 :
1219 47 : return bOK;
1220 : }
1221 :
1222 : GDALRasterCalcAlgorithmStandalone::~GDALRasterCalcAlgorithmStandalone() =
1223 : default;
1224 :
1225 : //! @endcond
|
