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