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 :
23 : #include <algorithm>
24 : #include <array>
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_Float64};
36 : bool checkSRS{true};
37 : bool checkExtent{true};
38 : };
39 :
40 101 : static bool MatchIsCompleteVariableNameWithNoIndex(const std::string &str,
41 : size_t from, size_t to)
42 : {
43 101 : 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 130 : if (std::isalnum(str[to]) || str[to] == '_' || str[to] == '[' ||
50 46 : str[to] == '(')
51 : {
52 39 : return false;
53 : }
54 : }
55 62 : 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 39 : if (std::isalnum(str[from - 1]) || str[from - 1] == '_')
60 : {
61 3 : return false;
62 : }
63 : }
64 :
65 59 : 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 65 : static std::string SetBandIndices(const std::string &origExpression,
74 : const std::string &variable, int band,
75 : bool &expressionChanged)
76 : {
77 65 : std::string expression = origExpression;
78 65 : expressionChanged = false;
79 :
80 65 : std::string::size_type seekPos = 0;
81 65 : auto pos = expression.find(variable, seekPos);
82 166 : while (pos != std::string::npos)
83 : {
84 101 : auto end = pos + variable.size();
85 :
86 101 : if (MatchIsCompleteVariableNameWithNoIndex(expression, pos, end))
87 : {
88 : // No index specified for variable
89 118 : expression = expression.substr(0, pos + variable.size()) + '[' +
90 177 : std::to_string(band) + ']' + expression.substr(end);
91 59 : expressionChanged = true;
92 : }
93 :
94 101 : seekPos = end;
95 101 : pos = expression.find(variable, seekPos);
96 : }
97 :
98 65 : return expression;
99 : }
100 :
101 : struct SourceProperties
102 : {
103 : int nBands{0};
104 : int nX{0};
105 : int nY{0};
106 : std::array<double, 6> gt{};
107 : std::unique_ptr<OGRSpatialReference, OGRSpatialReferenceReleaser> srs{
108 : nullptr};
109 : };
110 :
111 : static std::optional<SourceProperties>
112 55 : UpdateSourceProperties(SourceProperties &out, const std::string &dsn,
113 : const GDALCalcOptions &options)
114 : {
115 110 : SourceProperties source;
116 55 : bool srsMismatch = false;
117 55 : bool extentMismatch = false;
118 55 : bool dimensionMismatch = false;
119 :
120 : {
121 : std::unique_ptr<GDALDataset> ds(
122 55 : GDALDataset::Open(dsn.c_str(), GDAL_OF_RASTER));
123 :
124 55 : if (!ds)
125 : {
126 0 : CPLError(CE_Failure, CPLE_AppDefined, "Failed to open %s",
127 : dsn.c_str());
128 0 : return std::nullopt;
129 : }
130 :
131 55 : source.nBands = ds->GetRasterCount();
132 55 : source.nX = ds->GetRasterXSize();
133 55 : source.nY = ds->GetRasterYSize();
134 :
135 55 : if (options.checkExtent)
136 : {
137 51 : ds->GetGeoTransform(source.gt.data());
138 : }
139 :
140 55 : if (options.checkSRS && out.srs)
141 : {
142 12 : const OGRSpatialReference *srs = ds->GetSpatialRef();
143 12 : srsMismatch = srs && !srs->IsSame(out.srs.get());
144 : }
145 : }
146 :
147 55 : if (source.nX != out.nX || source.nY != out.nY)
148 : {
149 3 : dimensionMismatch = true;
150 : }
151 :
152 110 : if (source.gt[0] != out.gt[0] || source.gt[2] != out.gt[2] ||
153 110 : source.gt[3] != out.gt[3] || source.gt[4] != out.gt[4])
154 : {
155 5 : extentMismatch = true;
156 : }
157 55 : if (source.gt[1] != out.gt[1] || source.gt[5] != out.gt[5])
158 : {
159 : // Resolutions are different. Are the extents the same?
160 8 : double xmaxOut = out.gt[0] + out.nX * out.gt[1] + out.nY * out.gt[2];
161 8 : double yminOut = out.gt[3] + out.nX * out.gt[4] + out.nY * out.gt[5];
162 :
163 : double xmax =
164 8 : source.gt[0] + source.nX * source.gt[1] + source.nY * source.gt[2];
165 : double ymin =
166 8 : source.gt[3] + source.nX * source.gt[4] + source.nY * source.gt[5];
167 :
168 : // Max allowable extent misalignment, expressed as fraction of a pixel
169 8 : constexpr double EXTENT_RTOL = 1e-3;
170 :
171 11 : if (std::abs(xmax - xmaxOut) > EXTENT_RTOL * std::abs(source.gt[1]) ||
172 3 : std::abs(ymin - yminOut) > EXTENT_RTOL * std::abs(source.gt[5]))
173 : {
174 5 : extentMismatch = true;
175 : }
176 : }
177 :
178 55 : if (options.checkExtent && extentMismatch)
179 : {
180 1 : CPLError(CE_Failure, CPLE_AppDefined,
181 : "Input extents are inconsistent.");
182 1 : return std::nullopt;
183 : }
184 :
185 54 : if (!options.checkExtent && dimensionMismatch)
186 : {
187 1 : CPLError(CE_Failure, CPLE_AppDefined,
188 : "Inputs do not have the same dimensions.");
189 1 : return std::nullopt;
190 : }
191 :
192 : // Find a common resolution
193 53 : if (source.nX > out.nX)
194 : {
195 1 : auto dx = CPLGreatestCommonDivisor(out.gt[1], source.gt[1]);
196 1 : if (dx == 0)
197 : {
198 0 : CPLError(CE_Failure, CPLE_AppDefined,
199 : "Failed to find common resolution for inputs.");
200 0 : return std::nullopt;
201 : }
202 1 : out.nX = static_cast<int>(
203 1 : std::round(static_cast<double>(out.nX) * out.gt[1] / dx));
204 1 : out.gt[1] = dx;
205 : }
206 53 : if (source.nY > out.nY)
207 : {
208 1 : auto dy = CPLGreatestCommonDivisor(out.gt[5], source.gt[5]);
209 1 : if (dy == 0)
210 : {
211 0 : CPLError(CE_Failure, CPLE_AppDefined,
212 : "Failed to find common resolution for inputs.");
213 0 : return std::nullopt;
214 : }
215 1 : out.nY = static_cast<int>(
216 1 : std::round(static_cast<double>(out.nY) * out.gt[5] / dy));
217 1 : out.gt[5] = dy;
218 : }
219 :
220 53 : if (srsMismatch)
221 : {
222 1 : CPLError(CE_Failure, CPLE_AppDefined,
223 : "Input spatial reference systems are inconsistent.");
224 1 : return std::nullopt;
225 : }
226 :
227 52 : return source;
228 : }
229 :
230 : /** Create XML nodes for one or more derived bands resulting from the evaluation
231 : * of a single expression
232 : *
233 : * @param root VRTDataset node to which the band nodes should be added
234 : * @param nXOut Number of columns in VRT dataset
235 : * @param nYOut Number of rows in VRT dataset
236 : * @param expression Expression for which band(s) should be added
237 : * @param sources Mapping of source names to DSNs
238 : * @param sourceProps Mapping of source names to properties
239 : * @return true if the band(s) were added, false otherwise
240 : */
241 : static bool
242 38 : CreateDerivedBandXML(CPLXMLNode *root, int nXOut, int nYOut,
243 : GDALDataType bandType, const std::string &expression,
244 : const std::map<std::string, std::string> &sources,
245 : const std::map<std::string, SourceProperties> &sourceProps)
246 : {
247 38 : int nOutBands = 1; // By default, each expression produces a single output
248 : // band. When processing the expression below, we may
249 : // discover that the expression produces multiple bands,
250 : // in which case this will be updated.
251 84 : for (int nOutBand = 1; nOutBand <= nOutBands; nOutBand++)
252 : {
253 : // Copy the expression for each output band, because we may modify it
254 : // when adding band indices (e.g., X -> X[1]) to the variables in the
255 : // expression.
256 48 : std::string bandExpression = expression;
257 :
258 48 : CPLXMLNode *band = CPLCreateXMLNode(root, CXT_Element, "VRTRasterBand");
259 48 : CPLAddXMLAttributeAndValue(band, "subClass", "VRTDerivedRasterBand");
260 48 : CPLAddXMLAttributeAndValue(band, "dataType",
261 : GDALGetDataTypeName(bandType));
262 :
263 : CPLXMLNode *sourceTransferType =
264 48 : CPLCreateXMLNode(band, CXT_Element, "SourceTransferType");
265 48 : CPLCreateXMLNode(sourceTransferType, CXT_Text,
266 : GDALGetDataTypeName(GDT_Float64));
267 :
268 : CPLXMLNode *pixelFunctionType =
269 48 : CPLCreateXMLNode(band, CXT_Element, "PixelFunctionType");
270 48 : CPLCreateXMLNode(pixelFunctionType, CXT_Text, "expression");
271 : CPLXMLNode *arguments =
272 48 : CPLCreateXMLNode(band, CXT_Element, "PixelFunctionArguments");
273 :
274 111 : for (const auto &[source_name, dsn] : sources)
275 : {
276 65 : auto it = sourceProps.find(source_name);
277 65 : CPLAssert(it != sourceProps.end());
278 65 : const auto &props = it->second;
279 :
280 : {
281 65 : const int nDefaultInBand = std::min(props.nBands, nOutBand);
282 :
283 65 : CPLString expressionBandVariable;
284 : expressionBandVariable.Printf("%s[%d]", source_name.c_str(),
285 65 : nDefaultInBand);
286 :
287 65 : bool expressionUsesAllBands = false;
288 : bandExpression =
289 130 : SetBandIndices(bandExpression, source_name, nDefaultInBand,
290 65 : expressionUsesAllBands);
291 :
292 65 : if (expressionUsesAllBands)
293 : {
294 51 : if (nOutBands <= 1)
295 : {
296 37 : nOutBands = props.nBands;
297 : }
298 14 : else if (props.nBands != 1 && props.nBands != nOutBands)
299 : {
300 2 : CPLError(CE_Failure, CPLE_AppDefined,
301 : "Expression cannot operate on all bands of "
302 : "rasters with incompatible numbers of bands "
303 : "(source %s has %d bands but expected to have "
304 : "1 or %d bands).",
305 2 : source_name.c_str(), props.nBands, nOutBands);
306 2 : return false;
307 : }
308 : }
309 : }
310 :
311 : // Create a <SimpleSource> for each input band that is used in
312 : // the expression.
313 170 : for (int nInBand = 1; nInBand <= props.nBands; nInBand++)
314 : {
315 107 : CPLString inBandVariable;
316 107 : inBandVariable.Printf("%s[%d]", source_name.c_str(), nInBand);
317 107 : if (bandExpression.find(inBandVariable) == std::string::npos)
318 : {
319 31 : continue;
320 : }
321 :
322 : CPLXMLNode *source =
323 76 : CPLCreateXMLNode(band, CXT_Element, "SimpleSource");
324 76 : CPLAddXMLAttributeAndValue(source, "name",
325 : inBandVariable.c_str());
326 :
327 : CPLXMLNode *sourceFilename =
328 76 : CPLCreateXMLNode(source, CXT_Element, "SourceFilename");
329 76 : CPLAddXMLAttributeAndValue(sourceFilename, "relativeToVRT",
330 : "0");
331 76 : CPLCreateXMLNode(sourceFilename, CXT_Text, dsn.c_str());
332 :
333 : CPLXMLNode *sourceBand =
334 76 : CPLCreateXMLNode(source, CXT_Element, "SourceBand");
335 76 : CPLCreateXMLNode(sourceBand, CXT_Text,
336 152 : std::to_string(nInBand).c_str());
337 :
338 : // TODO add <SourceProperties> ?
339 :
340 : CPLXMLNode *srcRect =
341 76 : CPLCreateXMLNode(source, CXT_Element, "SrcRect");
342 76 : CPLAddXMLAttributeAndValue(srcRect, "xOff", "0");
343 76 : CPLAddXMLAttributeAndValue(srcRect, "yOff", "0");
344 76 : CPLAddXMLAttributeAndValue(srcRect, "xSize",
345 152 : std::to_string(props.nX).c_str());
346 76 : CPLAddXMLAttributeAndValue(srcRect, "ySize",
347 152 : std::to_string(props.nY).c_str());
348 :
349 : CPLXMLNode *dstRect =
350 76 : CPLCreateXMLNode(source, CXT_Element, "DstRect");
351 76 : CPLAddXMLAttributeAndValue(dstRect, "xOff", "0");
352 76 : CPLAddXMLAttributeAndValue(dstRect, "yOff", "0");
353 76 : CPLAddXMLAttributeAndValue(dstRect, "xSize",
354 152 : std::to_string(nXOut).c_str());
355 76 : CPLAddXMLAttributeAndValue(dstRect, "ySize",
356 152 : std::to_string(nYOut).c_str());
357 : }
358 : }
359 :
360 : // Add the expression as a last step, because we may modify the
361 : // expression as we iterate through the bands.
362 46 : CPLAddXMLAttributeAndValue(arguments, "expression",
363 : bandExpression.c_str());
364 46 : CPLAddXMLAttributeAndValue(arguments, "dialect", "muparser");
365 : }
366 :
367 36 : return true;
368 : }
369 :
370 40 : static bool ParseSourceDescriptors(const std::vector<std::string> &inputs,
371 : std::map<std::string, std::string> &datasets,
372 : std::string &firstSourceName)
373 : {
374 40 : bool isFirst = true;
375 :
376 97 : for (const auto &input : inputs)
377 : {
378 57 : std::string name = "";
379 :
380 57 : auto pos = input.find('=');
381 57 : if (pos == std::string::npos)
382 : {
383 10 : if (inputs.size() > 1)
384 : {
385 0 : CPLError(CE_Failure, CPLE_AppDefined,
386 : "Inputs must be named when more than one input is "
387 : "provided.");
388 0 : return false;
389 : }
390 10 : name = "X";
391 : }
392 : else
393 : {
394 47 : name = input.substr(0, pos);
395 : }
396 :
397 : std::string dsn =
398 114 : (pos == std::string::npos) ? input : input.substr(pos + 1);
399 57 : datasets[name] = std::move(dsn);
400 :
401 57 : if (isFirst)
402 : {
403 40 : firstSourceName = name;
404 40 : isFirst = false;
405 : }
406 : }
407 :
408 40 : return true;
409 : }
410 :
411 40 : static bool ReadFileLists(std::vector<std::string> &inputs)
412 : {
413 97 : for (std::size_t i = 0; i < inputs.size(); i++)
414 : {
415 57 : const auto &input = inputs[i];
416 57 : if (input[0] == '@')
417 : {
418 : auto f =
419 1 : VSIVirtualHandleUniquePtr(VSIFOpenL(input.c_str() + 1, "r"));
420 1 : if (!f)
421 : {
422 0 : CPLError(CE_Failure, CPLE_FileIO, "Cannot open %s",
423 0 : input.c_str() + 1);
424 0 : return false;
425 : }
426 1 : std::vector<std::string> sources;
427 3 : while (const char *filename = CPLReadLineL(f.get()))
428 : {
429 2 : sources.push_back(filename);
430 2 : }
431 1 : inputs.erase(inputs.begin() + static_cast<int>(i));
432 1 : inputs.insert(inputs.end(), sources.begin(), sources.end());
433 : }
434 : }
435 :
436 40 : return true;
437 : }
438 :
439 : /** Creates a VRT datasource with one or more derived raster bands containing
440 : * results of an expression.
441 : *
442 : * To make this work with muparser (which does not support vector types), we
443 : * do a simple parsing of the expression internally, transforming it into
444 : * multiple expressions with explicit band indices. For example, for a two-band
445 : * raster "X", the expression "X + 3" will be transformed into "X[1] + 3" and
446 : * "X[2] + 3". The use of brackets is for readability only; as far as the
447 : * expression engine is concerned, the variables "X[1]" and "X[2]" have nothing
448 : * to do with each other.
449 : *
450 : * @param inputs A list of sources, expressed as NAME=DSN
451 : * @param expressions A list of expressions to be evaluated
452 : * @param options flags controlling which checks should be performed on the inputs
453 : *
454 : * @return a newly created VRTDataset, or nullptr on error
455 : */
456 : static std::unique_ptr<GDALDataset>
457 40 : GDALCalcCreateVRTDerived(const std::vector<std::string> &inputs,
458 : const std::vector<std::string> &expressions,
459 : const GDALCalcOptions &options)
460 : {
461 40 : if (inputs.empty())
462 : {
463 0 : return nullptr;
464 : }
465 :
466 80 : std::map<std::string, std::string> sources;
467 80 : std::string firstSource;
468 40 : if (!ParseSourceDescriptors(inputs, sources, firstSource))
469 : {
470 0 : return nullptr;
471 : }
472 :
473 : // Use the first source provided to determine properties of the output
474 40 : const char *firstDSN = sources[firstSource].c_str();
475 :
476 : // Read properties from the first source
477 80 : SourceProperties out;
478 : {
479 : std::unique_ptr<GDALDataset> ds(
480 40 : GDALDataset::Open(firstDSN, GDAL_OF_RASTER));
481 :
482 40 : if (!ds)
483 : {
484 0 : CPLError(CE_Failure, CPLE_AppDefined, "Failed to open %s",
485 : firstDSN);
486 0 : return nullptr;
487 : }
488 :
489 40 : out.nX = ds->GetRasterXSize();
490 40 : out.nY = ds->GetRasterYSize();
491 40 : out.nBands = 1;
492 40 : out.srs.reset(ds->GetSpatialRef() ? ds->GetSpatialRef()->Clone()
493 : : nullptr);
494 40 : ds->GetGeoTransform(out.gt.data());
495 : }
496 :
497 80 : CPLXMLTreeCloser root(CPLCreateXMLNode(nullptr, CXT_Element, "VRTDataset"));
498 :
499 : // Collect properties of the different sources, and verity them for
500 : // consistency.
501 80 : std::map<std::string, SourceProperties> sourceProps;
502 92 : for (const auto &[source_name, dsn] : sources)
503 : {
504 : // TODO avoid opening the first source twice.
505 55 : auto props = UpdateSourceProperties(out, dsn, options);
506 55 : if (props.has_value())
507 : {
508 52 : sourceProps[source_name] = std::move(props.value());
509 : }
510 : else
511 : {
512 3 : return nullptr;
513 : }
514 : }
515 :
516 73 : for (const auto &origExpression : expressions)
517 : {
518 38 : if (!CreateDerivedBandXML(root.get(), out.nX, out.nY, options.dstType,
519 : origExpression, sources, sourceProps))
520 : {
521 2 : return nullptr;
522 : }
523 : }
524 :
525 : //CPLDebug("VRT", "%s", CPLSerializeXMLTree(root.get()));
526 :
527 70 : auto ds = std::make_unique<VRTDataset>(out.nX, out.nY);
528 35 : if (ds->XMLInit(root.get(), "") != CE_None)
529 : {
530 0 : return nullptr;
531 : };
532 35 : ds->SetGeoTransform(out.gt.data());
533 35 : if (out.srs)
534 : {
535 12 : ds->SetSpatialRef(out.srs.get());
536 : }
537 :
538 35 : return ds;
539 : }
540 :
541 : /************************************************************************/
542 : /* GDALRasterCalcAlgorithm::GDALRasterCalcAlgorithm() */
543 : /************************************************************************/
544 :
545 38 : GDALRasterCalcAlgorithm::GDALRasterCalcAlgorithm() noexcept
546 38 : : GDALAlgorithm(NAME, DESCRIPTION, HELP_URL)
547 : {
548 38 : m_supportsStreamedOutput = true;
549 :
550 38 : AddProgressArg();
551 :
552 76 : AddArg(GDAL_ARG_NAME_INPUT, 'i', _("Input raster datasets"), &m_inputs)
553 38 : .SetPositional()
554 38 : .SetMinCount(1)
555 38 : .SetAutoOpenDataset(false)
556 38 : .SetMetaVar("INPUTS");
557 :
558 : AddOutputFormatArg(&m_format, /* bStreamAllowed = */ true,
559 38 : /* bGDALGAllowed = */ true);
560 38 : AddOutputDatasetArg(&m_outputDataset, GDAL_OF_RASTER);
561 38 : AddCreationOptionsArg(&m_creationOptions);
562 38 : AddOverwriteArg(&m_overwrite);
563 38 : AddOutputDataTypeArg(&m_type);
564 :
565 : AddArg("no-check-srs", 0,
566 : _("Do not check consistency of input spatial reference systems"),
567 38 : &m_NoCheckSRS);
568 : AddArg("no-check-extent", 0, _("Do not check consistency of input extents"),
569 38 : &m_NoCheckExtent);
570 :
571 76 : AddArg("calc", 0, _("Expression(s) to evaluate"), &m_expr)
572 38 : .SetRequired()
573 38 : .SetMinCount(1);
574 38 : }
575 :
576 : /************************************************************************/
577 : /* GDALRasterCalcAlgorithm::RunImpl() */
578 : /************************************************************************/
579 :
580 40 : bool GDALRasterCalcAlgorithm::RunImpl(GDALProgressFunc pfnProgress,
581 : void *pProgressData)
582 : {
583 40 : CPLAssert(!m_outputDataset.GetDatasetRef());
584 :
585 40 : GDALCalcOptions options;
586 40 : options.checkExtent = !m_NoCheckExtent;
587 40 : options.checkSRS = !m_NoCheckSRS;
588 40 : if (!m_type.empty())
589 : {
590 1 : options.dstType = GDALGetDataTypeByName(m_type.c_str());
591 : }
592 :
593 40 : if (!ReadFileLists(m_inputs))
594 : {
595 0 : return false;
596 : }
597 :
598 80 : auto vrt = GDALCalcCreateVRTDerived(m_inputs, m_expr, options);
599 40 : if (vrt == nullptr)
600 : {
601 5 : return false;
602 : }
603 :
604 35 : if (m_format == "stream")
605 : {
606 1 : m_outputDataset.Set(std::move(vrt));
607 1 : return true;
608 : }
609 :
610 68 : CPLStringList translateArgs;
611 34 : if (!m_format.empty())
612 : {
613 2 : translateArgs.AddString("-of");
614 2 : translateArgs.AddString(m_format.c_str());
615 : }
616 35 : for (const auto &co : m_creationOptions)
617 : {
618 1 : translateArgs.AddString("-co");
619 1 : translateArgs.AddString(co.c_str());
620 : }
621 :
622 : GDALTranslateOptions *translateOptions =
623 34 : GDALTranslateOptionsNew(translateArgs.List(), nullptr);
624 34 : GDALTranslateOptionsSetProgress(translateOptions, pfnProgress,
625 : pProgressData);
626 :
627 : auto poOutDS =
628 : std::unique_ptr<GDALDataset>(GDALDataset::FromHandle(GDALTranslate(
629 34 : m_outputDataset.GetName().c_str(), GDALDataset::ToHandle(vrt.get()),
630 68 : translateOptions, nullptr)));
631 34 : GDALTranslateOptionsFree(translateOptions);
632 :
633 34 : const bool bOK = poOutDS != nullptr;
634 34 : m_outputDataset.Set(std::move(poOutDS));
635 :
636 34 : return bOK;
637 : }
638 :
639 : //! @endcond
|
