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