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