Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: GDAL
4 : * Purpose: "edit" step of "raster pipeline"
5 : * Author: Even Rouault <even dot rouault at spatialys.com>
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2024, Even Rouault <even dot rouault at spatialys.com>
9 : *
10 : * SPDX-License-Identifier: MIT
11 : ****************************************************************************/
12 :
13 : #include "gdalalg_raster_edit.h"
14 :
15 : #include "gdal_priv.h"
16 : #include "gdal_utils.h"
17 : #include "ogrsf_frmts.h"
18 :
19 : #include <optional>
20 :
21 : //! @cond Doxygen_Suppress
22 :
23 : #ifndef _
24 : #define _(x) (x)
25 : #endif
26 :
27 : /************************************************************************/
28 : /* GetGCPFilename() */
29 : /************************************************************************/
30 :
31 29 : static std::string GetGCPFilename(const std::vector<std::string> &gcps)
32 : {
33 29 : if (gcps.size() == 1 && !gcps[0].empty() && gcps[0][0] == '@')
34 : {
35 12 : return gcps[0].substr(1);
36 : }
37 17 : return std::string();
38 : }
39 :
40 : /************************************************************************/
41 : /* GDALRasterEditAlgorithm::GDALRasterEditAlgorithm() */
42 : /************************************************************************/
43 :
44 185 : GDALRasterEditAlgorithm::GDALRasterEditAlgorithm(bool standaloneStep)
45 : : GDALRasterPipelineStepAlgorithm(
46 : NAME, DESCRIPTION, HELP_URL,
47 185 : ConstructorOptions().SetAddDefaultArguments(false))
48 : {
49 185 : if (standaloneStep)
50 : {
51 96 : AddProgressArg();
52 :
53 : AddArg("dataset", 0,
54 : _("Dataset (to be updated in-place, unless --auxiliary)"),
55 192 : &m_dataset, GDAL_OF_RASTER | GDAL_OF_UPDATE)
56 96 : .SetPositional()
57 96 : .SetRequired();
58 96 : AddOpenOptionsArg(&m_openOptions);
59 : AddArg("auxiliary", 0,
60 : _("Ask for an auxiliary .aux.xml file to be edited"),
61 192 : &m_readOnly)
62 192 : .AddHiddenAlias("ro")
63 96 : .AddHiddenAlias(GDAL_ARG_NAME_READ_ONLY);
64 : }
65 : else
66 : {
67 89 : AddRasterHiddenInputDatasetArg();
68 : }
69 :
70 370 : AddArg("crs", 0, _("Override CRS (without reprojection)"), &m_overrideCrs)
71 370 : .AddHiddenAlias("a_srs")
72 370 : .AddHiddenAlias("srs")
73 185 : .SetIsCRSArg(/*noneAllowed=*/true);
74 :
75 185 : AddBBOXArg(&m_bbox);
76 :
77 185 : AddNodataArg(&m_nodata, /* noneAllowed = */ true);
78 :
79 : AddArg("color-interpretation", 0, _("Set band color interpretation"),
80 370 : &m_colorInterpretation)
81 : .SetAutoCompleteFunction(
82 5 : [this](const std::string &s)
83 : {
84 3 : std::vector<std::string> ret;
85 3 : int nValues = 0;
86 3 : const auto paeVals = GDALGetColorInterpretationList(&nValues);
87 3 : if (s.find('=') == std::string::npos)
88 : {
89 2 : ret.push_back("all=");
90 2 : if (auto poDS = m_dataset.GetDatasetRef())
91 : {
92 2 : for (int i = 0; i < poDS->GetRasterCount(); ++i)
93 1 : ret.push_back(std::to_string(i + 1).append("="));
94 : }
95 70 : for (int i = 0; i < nValues; ++i)
96 136 : ret.push_back(
97 68 : GDALGetColorInterpretationName(paeVals[i]));
98 : }
99 : else
100 : {
101 35 : for (int i = 0; i < nValues; ++i)
102 68 : ret.push_back(
103 34 : GDALGetColorInterpretationName(paeVals[i]));
104 : }
105 6 : return ret;
106 185 : });
107 :
108 : {
109 : auto &arg = AddArg("metadata", 0, _("Add/update dataset metadata item"),
110 370 : &m_metadata)
111 370 : .SetMetaVar("<KEY>=<VALUE>")
112 185 : .SetPackedValuesAllowed(false);
113 43 : arg.AddValidationAction([this, &arg]()
114 228 : { return ParseAndValidateKeyValue(arg); });
115 185 : arg.AddHiddenAlias("mo");
116 : }
117 :
118 : AddArg("unset-metadata", 0, _("Remove dataset metadata item(s)"),
119 370 : &m_unsetMetadata)
120 185 : .SetMetaVar("<KEY>");
121 :
122 : AddArg("unset-metadata-domain", 0, _("Remove dataset metadata domain(s)"),
123 370 : &m_unsetMetadataDomain)
124 185 : .SetMetaVar("<DOMAIN>");
125 :
126 : AddArg("gcp", 0,
127 : _("Add ground control point, formatted as "
128 : "pixel,line,easting,northing[,elevation], or @filename"),
129 370 : &m_gcps)
130 185 : .SetPackedValuesAllowed(false)
131 : .AddValidationAction(
132 36 : [this]()
133 : {
134 21 : if (GetGCPFilename(m_gcps).empty())
135 : {
136 28 : for (const std::string &gcp : m_gcps)
137 : {
138 : const CPLStringList aosTokens(
139 17 : CSLTokenizeString2(gcp.c_str(), ",", 0));
140 17 : if (aosTokens.size() != 4 && aosTokens.size() != 5)
141 : {
142 1 : ReportError(CE_Failure, CPLE_IllegalArg,
143 : "Bad format for %s", gcp.c_str());
144 1 : return false;
145 : }
146 85 : for (int i = 0; i < aosTokens.size(); ++i)
147 : {
148 70 : if (CPLGetValueType(aosTokens[i]) ==
149 : CPL_VALUE_STRING)
150 : {
151 1 : ReportError(CE_Failure, CPLE_IllegalArg,
152 : "Bad format for %s", gcp.c_str());
153 1 : return false;
154 : }
155 : }
156 : }
157 : }
158 19 : return true;
159 185 : });
160 :
161 185 : if (standaloneStep)
162 : {
163 192 : AddArg("stats", 0, _("Compute statistics, using all pixels"), &m_stats)
164 96 : .SetMutualExclusionGroup("stats");
165 : AddArg("approx-stats", 0,
166 : _("Compute statistics, using a subset of pixels"),
167 192 : &m_approxStats)
168 96 : .SetMutualExclusionGroup("stats");
169 96 : AddArg("hist", 0, _("Compute histogram"), &m_hist);
170 : }
171 185 : }
172 :
173 : /************************************************************************/
174 : /* GDALRasterEditAlgorithm::~GDALRasterEditAlgorithm() */
175 : /************************************************************************/
176 :
177 : GDALRasterEditAlgorithm::~GDALRasterEditAlgorithm() = default;
178 :
179 : /************************************************************************/
180 : /* ParseGCPs() */
181 : /************************************************************************/
182 :
183 8 : std::vector<gdal::GCP> GDALRasterEditAlgorithm::ParseGCPs() const
184 : {
185 8 : std::vector<gdal::GCP> ret;
186 16 : const std::string osGCPFilename = GetGCPFilename(m_gcps);
187 8 : if (!osGCPFilename.empty())
188 : {
189 : auto poDS = std::unique_ptr<GDALDataset>(GDALDataset::Open(
190 4 : osGCPFilename.c_str(), GDAL_OF_VECTOR | GDAL_OF_VERBOSE_ERROR));
191 4 : if (!poDS)
192 1 : return ret;
193 3 : if (poDS->GetLayerCount() != 1)
194 : {
195 1 : ReportError(CE_Failure, CPLE_AppDefined,
196 : "GCPs can only be specified for single-layer datasets");
197 1 : return ret;
198 : }
199 2 : auto poLayer = poDS->GetLayer(0);
200 2 : const auto poLayerDefn = poLayer->GetLayerDefn();
201 2 : int nIdIdx = -1, nInfoIdx = -1, nColIdx = -1, nLineIdx = -1, nXIdx = -1,
202 2 : nYIdx = -1, nZIdx = -1;
203 :
204 : const struct
205 : {
206 : int &idx;
207 : const char *name;
208 : bool required;
209 2 : } aFields[] = {
210 : {nIdIdx, "id", false}, {nInfoIdx, "info", false},
211 : {nColIdx, "column", true}, {nLineIdx, "line", true},
212 : {nXIdx, "x", true}, {nYIdx, "y", true},
213 : {nZIdx, "z", false},
214 2 : };
215 :
216 11 : for (auto &field : aFields)
217 : {
218 10 : field.idx = poLayerDefn->GetFieldIndex(field.name);
219 10 : if (field.idx < 0 && field.required)
220 : {
221 3 : ReportError(CE_Failure, CPLE_AppDefined,
222 1 : "Field '%s' cannot be found in '%s'", field.name,
223 1 : poDS->GetDescription());
224 1 : return ret;
225 : }
226 : }
227 2 : for (auto &&poFeature : poLayer)
228 : {
229 2 : gdal::GCP gcp;
230 1 : if (nIdIdx >= 0)
231 1 : gcp.SetId(poFeature->GetFieldAsString(nIdIdx));
232 1 : if (nInfoIdx >= 0)
233 1 : gcp.SetInfo(poFeature->GetFieldAsString(nInfoIdx));
234 1 : gcp.Pixel() = poFeature->GetFieldAsDouble(nColIdx);
235 1 : gcp.Line() = poFeature->GetFieldAsDouble(nLineIdx);
236 1 : gcp.X() = poFeature->GetFieldAsDouble(nXIdx);
237 1 : gcp.Y() = poFeature->GetFieldAsDouble(nYIdx);
238 1 : if (nZIdx >= 0 && poFeature->IsFieldSetAndNotNull(nZIdx))
239 1 : gcp.Z() = poFeature->GetFieldAsDouble(nZIdx);
240 1 : ret.push_back(std::move(gcp));
241 : }
242 : }
243 : else
244 : {
245 10 : for (const std::string &gcpStr : m_gcps)
246 : {
247 : const CPLStringList aosTokens(
248 12 : CSLTokenizeString2(gcpStr.c_str(), ",", 0));
249 : // Verified by validation action
250 6 : CPLAssert(aosTokens.size() == 4 || aosTokens.size() == 5);
251 12 : gdal::GCP gcp;
252 6 : gcp.Pixel() = CPLAtof(aosTokens[0]);
253 6 : gcp.Line() = CPLAtof(aosTokens[1]);
254 6 : gcp.X() = CPLAtof(aosTokens[2]);
255 6 : gcp.Y() = CPLAtof(aosTokens[3]);
256 6 : if (aosTokens.size() == 5)
257 3 : gcp.Z() = CPLAtof(aosTokens[4]);
258 6 : ret.push_back(std::move(gcp));
259 : }
260 : }
261 5 : return ret;
262 : }
263 :
264 : /************************************************************************/
265 : /* GDALRasterEditAlgorithm::RunStep() */
266 : /************************************************************************/
267 :
268 62 : bool GDALRasterEditAlgorithm::RunStep(GDALPipelineStepRunContext &ctxt)
269 : {
270 62 : GDALDataset *poDS = m_dataset.GetDatasetRef();
271 62 : if (poDS)
272 : {
273 44 : if (poDS->GetAccess() != GA_Update && !m_readOnly)
274 : {
275 1 : ReportError(CE_Failure, CPLE_AppDefined,
276 : "Dataset should be opened in update mode unless "
277 : "--auxiliary is set");
278 1 : return false;
279 : }
280 : }
281 : else
282 : {
283 18 : const auto poSrcDS = m_inputDataset[0].GetDatasetRef();
284 18 : CPLAssert(poSrcDS);
285 18 : CPLAssert(m_outputDataset.GetName().empty());
286 18 : CPLAssert(!m_outputDataset.GetDatasetRef());
287 :
288 36 : CPLStringList aosOptions;
289 18 : aosOptions.push_back("-of");
290 18 : aosOptions.push_back("VRT");
291 : GDALTranslateOptions *psOptions =
292 18 : GDALTranslateOptionsNew(aosOptions.List(), nullptr);
293 18 : GDALDatasetH hSrcDS = GDALDataset::ToHandle(poSrcDS);
294 18 : auto poRetDS = GDALDataset::FromHandle(
295 : GDALTranslate("", hSrcDS, psOptions, nullptr));
296 18 : GDALTranslateOptionsFree(psOptions);
297 18 : m_outputDataset.Set(std::unique_ptr<GDALDataset>(poRetDS));
298 18 : poDS = m_outputDataset.GetDatasetRef();
299 : }
300 :
301 61 : bool ret = poDS != nullptr;
302 :
303 61 : if (poDS)
304 : {
305 61 : if (m_overrideCrs == "null" || m_overrideCrs == "none")
306 : {
307 3 : if (poDS->SetSpatialRef(nullptr) != CE_None)
308 : {
309 1 : ReportError(CE_Failure, CPLE_AppDefined,
310 : "SetSpatialRef(%s) failed", m_overrideCrs.c_str());
311 22 : return false;
312 : }
313 : }
314 58 : else if (!m_overrideCrs.empty() && m_gcps.empty())
315 : {
316 4 : OGRSpatialReference oSRS;
317 4 : oSRS.SetFromUserInput(m_overrideCrs.c_str());
318 4 : oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
319 4 : if (poDS->SetSpatialRef(&oSRS) != CE_None)
320 : {
321 1 : ReportError(CE_Failure, CPLE_AppDefined,
322 : "SetSpatialRef(%s) failed", m_overrideCrs.c_str());
323 1 : return false;
324 : }
325 : }
326 :
327 59 : if (!m_bbox.empty())
328 : {
329 4 : if (poDS->GetRasterXSize() == 0 || poDS->GetRasterYSize() == 0)
330 : {
331 1 : ReportError(CE_Failure, CPLE_AppDefined,
332 : "Cannot set extent because one of dataset height "
333 : "or width is null");
334 2 : return false;
335 : }
336 3 : GDALGeoTransform gt;
337 3 : gt.xorig = m_bbox[0];
338 3 : gt.xscale = (m_bbox[2] - m_bbox[0]) / poDS->GetRasterXSize();
339 3 : gt.xrot = 0;
340 3 : gt.yorig = m_bbox[3];
341 3 : gt.yrot = 0;
342 3 : gt.yscale = -(m_bbox[3] - m_bbox[1]) / poDS->GetRasterYSize();
343 3 : if (poDS->SetGeoTransform(gt) != CE_None)
344 : {
345 1 : ReportError(CE_Failure, CPLE_AppDefined,
346 : "Setting extent failed");
347 1 : return false;
348 : }
349 : }
350 :
351 57 : if (!m_nodata.empty())
352 : {
353 18 : for (int i = 0; i < poDS->GetRasterCount(); ++i)
354 : {
355 10 : if (EQUAL(m_nodata.c_str(), "none"))
356 2 : poDS->GetRasterBand(i + 1)->DeleteNoDataValue();
357 : else
358 16 : poDS->GetRasterBand(i + 1)->SetNoDataValue(
359 8 : CPLAtof(m_nodata.c_str()));
360 : }
361 : }
362 :
363 57 : if (!m_colorInterpretation.empty())
364 : {
365 : const auto GetColorInterp =
366 28 : [this](const char *pszStr) -> std::optional<GDALColorInterp>
367 : {
368 25 : if (EQUAL(pszStr, "undefined"))
369 0 : return GCI_Undefined;
370 : const GDALColorInterp eInterp =
371 25 : GDALGetColorInterpretationByName(pszStr);
372 25 : if (eInterp != GCI_Undefined)
373 22 : return eInterp;
374 3 : ReportError(CE_Failure, CPLE_NotSupported,
375 : "Unsupported color interpretation: %s", pszStr);
376 3 : return {};
377 18 : };
378 :
379 18 : if (m_colorInterpretation.size() == 1 &&
380 20 : poDS->GetRasterCount() > 1 &&
381 2 : !cpl::starts_with(m_colorInterpretation[0], "all="))
382 : {
383 1 : ReportError(
384 : CE_Failure, CPLE_NotSupported,
385 : "With several bands, specify as many color interpretation "
386 : "as bands, one or many values of the form "
387 : "<band_number>=<color> or a single value all=<color>");
388 12 : return false;
389 : }
390 : else
391 : {
392 17 : int nBandIter = 0;
393 17 : bool bSyntaxAll = false;
394 17 : bool bSyntaxExplicitBand = false;
395 17 : bool bSyntaxImplicitBand = false;
396 39 : for (const std::string &token : m_colorInterpretation)
397 : {
398 : const CPLStringList aosTokens(
399 29 : CSLTokenizeString2(token.c_str(), "=", 0));
400 29 : if (aosTokens.size() == 2 && EQUAL(aosTokens[0], "all"))
401 : {
402 4 : bSyntaxAll = true;
403 4 : const auto eColorInterp = GetColorInterp(aosTokens[1]);
404 4 : if (!eColorInterp)
405 : {
406 1 : return false;
407 : }
408 : else
409 : {
410 10 : for (int i = 0; i < poDS->GetRasterCount(); ++i)
411 : {
412 7 : if (poDS->GetRasterBand(i + 1)
413 7 : ->SetColorInterpretation(
414 14 : *eColorInterp) != CE_None)
415 0 : return false;
416 : }
417 : }
418 : }
419 25 : else if (aosTokens.size() == 2)
420 : {
421 9 : bSyntaxExplicitBand = true;
422 9 : const int nBand = atoi(aosTokens[0]);
423 9 : if (nBand <= 0 || nBand > poDS->GetRasterCount())
424 : {
425 2 : ReportError(CE_Failure, CPLE_NotSupported,
426 : "Invalid band number '%s' in '%s'",
427 : aosTokens[0], token.c_str());
428 3 : return false;
429 : }
430 7 : const auto eColorInterp = GetColorInterp(aosTokens[1]);
431 7 : if (!eColorInterp)
432 : {
433 1 : return false;
434 : }
435 6 : else if (poDS->GetRasterBand(nBand)
436 6 : ->SetColorInterpretation(*eColorInterp) !=
437 : CE_None)
438 : {
439 0 : return false;
440 : }
441 : }
442 : else
443 : {
444 16 : bSyntaxImplicitBand = true;
445 16 : ++nBandIter;
446 16 : if (nBandIter > poDS->GetRasterCount())
447 : {
448 2 : ReportError(CE_Failure, CPLE_IllegalArg,
449 : "More color interpretation values "
450 : "specified than bands in the dataset");
451 3 : return false;
452 : }
453 14 : const auto eColorInterp = GetColorInterp(token.c_str());
454 14 : if (!eColorInterp)
455 : {
456 1 : return false;
457 : }
458 13 : else if (poDS->GetRasterBand(nBandIter)
459 13 : ->SetColorInterpretation(*eColorInterp) !=
460 : CE_None)
461 : {
462 0 : return false;
463 : }
464 : }
465 : }
466 20 : if ((bSyntaxAll ? 1 : 0) + (bSyntaxExplicitBand ? 1 : 0) +
467 10 : (bSyntaxImplicitBand ? 1 : 0) !=
468 : 1)
469 : {
470 3 : ReportError(CE_Failure, CPLE_IllegalArg,
471 : "Mix of different syntaxes to specify color "
472 : "interpretation");
473 3 : return false;
474 : }
475 7 : if (bSyntaxImplicitBand && nBandIter != poDS->GetRasterCount())
476 : {
477 1 : ReportError(CE_Failure, CPLE_IllegalArg,
478 : "Less color interpretation values specified "
479 : "than bands in the dataset");
480 1 : return false;
481 : }
482 : }
483 : }
484 :
485 45 : const CPLStringList aosMD(m_metadata);
486 56 : for (const auto &[key, value] : cpl::IterateNameValue(aosMD))
487 : {
488 12 : if (poDS->SetMetadataItem(key, value) != CE_None)
489 : {
490 1 : ReportError(CE_Failure, CPLE_AppDefined,
491 : "SetMetadataItem('%s', '%s') failed", key, value);
492 1 : return false;
493 : }
494 : }
495 :
496 46 : for (const std::string &key : m_unsetMetadata)
497 : {
498 3 : if (poDS->SetMetadataItem(key.c_str(), nullptr) != CE_None)
499 : {
500 1 : ReportError(CE_Failure, CPLE_AppDefined,
501 : "SetMetadataItem('%s', NULL) failed", key.c_str());
502 1 : return false;
503 : }
504 : }
505 :
506 44 : for (const std::string &domain : m_unsetMetadataDomain)
507 : {
508 1 : if (poDS->SetMetadata(nullptr, domain.c_str()) != CE_None)
509 : {
510 0 : ReportError(CE_Failure, CPLE_AppDefined,
511 : "SetMetadata(NULL, '%s') failed", domain.c_str());
512 0 : return false;
513 : }
514 : }
515 :
516 43 : if (!m_gcps.empty())
517 : {
518 8 : const auto gcps = ParseGCPs();
519 8 : if (gcps.empty())
520 3 : return false; // error already emitted by ParseGCPs()
521 :
522 5 : OGRSpatialReference oSRS;
523 5 : if (!m_overrideCrs.empty())
524 : {
525 1 : oSRS.SetFromUserInput(m_overrideCrs.c_str());
526 1 : oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
527 : }
528 :
529 5 : if (poDS->SetGCPs(static_cast<int>(gcps.size()), gcps[0].c_ptr(),
530 10 : oSRS.IsEmpty() ? nullptr : &oSRS) != CE_None)
531 : {
532 1 : ReportError(CE_Failure, CPLE_AppDefined, "Setting GCPs failed");
533 1 : return false;
534 : }
535 : }
536 :
537 39 : const int nBands = poDS->GetRasterCount();
538 39 : int nCurProgress = 0;
539 39 : const double dfTotalProgress =
540 39 : ((m_stats || m_approxStats) ? nBands : 0) + (m_hist ? nBands : 0);
541 39 : if (m_stats || m_approxStats)
542 : {
543 4 : for (int i = 0; (i < nBands) && ret; ++i)
544 : {
545 4 : void *pScaledProgress = GDALCreateScaledProgress(
546 : nCurProgress / dfTotalProgress,
547 2 : (nCurProgress + 1) / dfTotalProgress, ctxt.m_pfnProgress,
548 : ctxt.m_pProgressData);
549 2 : ++nCurProgress;
550 2 : double dfMin = 0.0;
551 2 : double dfMax = 0.0;
552 2 : double dfMean = 0.0;
553 2 : double dfStdDev = 0.0;
554 4 : ret = poDS->GetRasterBand(i + 1)->ComputeStatistics(
555 2 : m_approxStats, &dfMin, &dfMax, &dfMean, &dfStdDev,
556 : pScaledProgress ? GDALScaledProgress : nullptr,
557 2 : pScaledProgress) == CE_None;
558 2 : GDALDestroyScaledProgress(pScaledProgress);
559 : }
560 : }
561 :
562 39 : if (m_hist)
563 : {
564 2 : for (int i = 0; (i < nBands) && ret; ++i)
565 : {
566 2 : void *pScaledProgress = GDALCreateScaledProgress(
567 : nCurProgress / dfTotalProgress,
568 1 : (nCurProgress + 1) / dfTotalProgress, ctxt.m_pfnProgress,
569 : ctxt.m_pProgressData);
570 1 : ++nCurProgress;
571 1 : double dfMin = 0.0;
572 1 : double dfMax = 0.0;
573 1 : int nBucketCount = 0;
574 1 : GUIntBig *panHistogram = nullptr;
575 2 : ret = poDS->GetRasterBand(i + 1)->GetDefaultHistogram(
576 : &dfMin, &dfMax, &nBucketCount, &panHistogram, TRUE,
577 : pScaledProgress ? GDALScaledProgress : nullptr,
578 1 : pScaledProgress) == CE_None;
579 1 : if (ret)
580 : {
581 2 : ret = poDS->GetRasterBand(i + 1)->SetDefaultHistogram(
582 1 : dfMin, dfMax, nBucketCount, panHistogram) ==
583 : CE_None;
584 : }
585 1 : CPLFree(panHistogram);
586 1 : GDALDestroyScaledProgress(pScaledProgress);
587 : }
588 : }
589 : }
590 :
591 39 : return poDS != nullptr;
592 : }
593 :
594 : GDALRasterEditAlgorithmStandalone::~GDALRasterEditAlgorithmStandalone() =
595 : default;
596 :
597 : //! @endcond
|
