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