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