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 : //! @cond Doxygen_Suppress
20 :
21 : #ifndef _
22 : #define _(x) (x)
23 : #endif
24 :
25 : /************************************************************************/
26 : /* GetGCPFilename() */
27 : /************************************************************************/
28 :
29 20 : static std::string GetGCPFilename(const std::vector<std::string> &gcps)
30 : {
31 20 : if (gcps.size() == 1 && !gcps[0].empty() && gcps[0][0] == '@')
32 : {
33 8 : return gcps[0].substr(1);
34 : }
35 12 : return std::string();
36 : }
37 :
38 : /************************************************************************/
39 : /* GDALRasterEditAlgorithm::GDALRasterEditAlgorithm() */
40 : /************************************************************************/
41 :
42 64 : GDALRasterEditAlgorithm::GDALRasterEditAlgorithm(bool standaloneStep)
43 : : GDALRasterPipelineStepAlgorithm(
44 : NAME, DESCRIPTION, HELP_URL,
45 64 : ConstructorOptions().SetAddDefaultArguments(false))
46 : {
47 64 : if (standaloneStep)
48 : {
49 31 : AddProgressArg();
50 :
51 : AddArg("dataset", 0,
52 : _("Dataset (to be updated in-place, unless --auxiliary)"),
53 62 : &m_dataset, GDAL_OF_RASTER | GDAL_OF_UPDATE)
54 31 : .SetPositional()
55 31 : .SetRequired();
56 : AddArg("auxiliary", 0,
57 : _("Ask for an auxiliary .aux.xml file to be edited"),
58 62 : &m_readOnly)
59 62 : .AddHiddenAlias("ro")
60 31 : .AddHiddenAlias(GDAL_ARG_NAME_READ_ONLY);
61 : }
62 :
63 128 : AddArg("crs", 0, _("Override CRS (without reprojection)"), &m_overrideCrs)
64 128 : .AddHiddenAlias("a_srs")
65 128 : .AddHiddenAlias("srs")
66 64 : .SetIsCRSArg(/*noneAllowed=*/true);
67 :
68 64 : AddBBOXArg(&m_bbox);
69 :
70 64 : AddNodataArg(&m_nodata, /* noneAllowed = */ true);
71 :
72 : {
73 : auto &arg = AddArg("metadata", 0, _("Add/update dataset metadata item"),
74 128 : &m_metadata)
75 128 : .SetMetaVar("<KEY>=<VALUE>")
76 64 : .SetPackedValuesAllowed(false);
77 3 : arg.AddValidationAction([this, &arg]()
78 67 : { return ParseAndValidateKeyValue(arg); });
79 64 : arg.AddHiddenAlias("mo");
80 : }
81 :
82 : AddArg("unset-metadata", 0, _("Remove dataset metadata item"),
83 128 : &m_unsetMetadata)
84 64 : .SetMetaVar("<KEY>");
85 :
86 : AddArg("gcp", 0,
87 : _("Add ground control point, formatted as "
88 : "pixel,line,easting,northing[,elevation], or @filename"),
89 128 : &m_gcps)
90 64 : .SetPackedValuesAllowed(false)
91 : .AddValidationAction(
92 22 : [this]()
93 : {
94 12 : if (GetGCPFilename(m_gcps).empty())
95 : {
96 16 : for (const std::string &gcp : m_gcps)
97 : {
98 : const CPLStringList aosTokens(
99 10 : CSLTokenizeString2(gcp.c_str(), ",", 0));
100 10 : if (aosTokens.size() != 4 && aosTokens.size() != 5)
101 : {
102 1 : ReportError(CE_Failure, CPLE_IllegalArg,
103 : "Bad format for %s", gcp.c_str());
104 1 : return false;
105 : }
106 47 : for (int i = 0; i < aosTokens.size(); ++i)
107 : {
108 39 : if (CPLGetValueType(aosTokens[i]) ==
109 : CPL_VALUE_STRING)
110 : {
111 1 : ReportError(CE_Failure, CPLE_IllegalArg,
112 : "Bad format for %s", gcp.c_str());
113 1 : return false;
114 : }
115 : }
116 : }
117 : }
118 10 : return true;
119 64 : });
120 :
121 64 : if (standaloneStep)
122 : {
123 62 : AddArg("stats", 0, _("Compute statistics, using all pixels"), &m_stats)
124 31 : .SetMutualExclusionGroup("stats");
125 : AddArg("approx-stats", 0,
126 : _("Compute statistics, using a subset of pixels"),
127 62 : &m_approxStats)
128 31 : .SetMutualExclusionGroup("stats");
129 31 : AddArg("hist", 0, _("Compute histogram"), &m_hist);
130 : }
131 64 : }
132 :
133 : /************************************************************************/
134 : /* GDALRasterEditAlgorithm::~GDALRasterEditAlgorithm() */
135 : /************************************************************************/
136 :
137 : GDALRasterEditAlgorithm::~GDALRasterEditAlgorithm() = default;
138 :
139 : /************************************************************************/
140 : /* ParseGCPs() */
141 : /************************************************************************/
142 :
143 8 : std::vector<gdal::GCP> GDALRasterEditAlgorithm::ParseGCPs() const
144 : {
145 8 : std::vector<gdal::GCP> ret;
146 16 : const std::string osGCPFilename = GetGCPFilename(m_gcps);
147 8 : if (!osGCPFilename.empty())
148 : {
149 : auto poDS = std::unique_ptr<GDALDataset>(GDALDataset::Open(
150 4 : osGCPFilename.c_str(), GDAL_OF_VECTOR | GDAL_OF_VERBOSE_ERROR));
151 4 : if (!poDS)
152 1 : return ret;
153 3 : if (poDS->GetLayerCount() != 1)
154 : {
155 1 : ReportError(CE_Failure, CPLE_AppDefined,
156 : "GCPs can only be specified for single-layer datasets");
157 1 : return ret;
158 : }
159 2 : auto poLayer = poDS->GetLayer(0);
160 2 : const auto poLayerDefn = poLayer->GetLayerDefn();
161 2 : int nIdIdx = -1, nInfoIdx = -1, nColIdx = -1, nLineIdx = -1, nXIdx = -1,
162 2 : nYIdx = -1, nZIdx = -1;
163 :
164 : const struct
165 : {
166 : int &idx;
167 : const char *name;
168 : bool required;
169 2 : } aFields[] = {
170 : {nIdIdx, "id", false}, {nInfoIdx, "info", false},
171 : {nColIdx, "column", true}, {nLineIdx, "line", true},
172 : {nXIdx, "x", true}, {nYIdx, "y", true},
173 : {nZIdx, "z", false},
174 2 : };
175 :
176 11 : for (auto &field : aFields)
177 : {
178 10 : field.idx = poLayerDefn->GetFieldIndex(field.name);
179 10 : if (field.idx < 0 && field.required)
180 : {
181 3 : ReportError(CE_Failure, CPLE_AppDefined,
182 1 : "Field '%s' cannot be found in '%s'", field.name,
183 1 : poDS->GetDescription());
184 1 : return ret;
185 : }
186 : }
187 2 : for (auto &&poFeature : poLayer)
188 : {
189 2 : gdal::GCP gcp;
190 1 : if (nIdIdx >= 0)
191 1 : gcp.SetId(poFeature->GetFieldAsString(nIdIdx));
192 1 : if (nInfoIdx >= 0)
193 1 : gcp.SetInfo(poFeature->GetFieldAsString(nInfoIdx));
194 1 : gcp.Pixel() = poFeature->GetFieldAsDouble(nColIdx);
195 1 : gcp.Line() = poFeature->GetFieldAsDouble(nLineIdx);
196 1 : gcp.X() = poFeature->GetFieldAsDouble(nXIdx);
197 1 : gcp.Y() = poFeature->GetFieldAsDouble(nYIdx);
198 1 : if (nZIdx >= 0 && poFeature->IsFieldSetAndNotNull(nZIdx))
199 1 : gcp.Z() = poFeature->GetFieldAsDouble(nZIdx);
200 1 : ret.push_back(std::move(gcp));
201 : }
202 : }
203 : else
204 : {
205 10 : for (const std::string &gcpStr : m_gcps)
206 : {
207 : const CPLStringList aosTokens(
208 12 : CSLTokenizeString2(gcpStr.c_str(), ",", 0));
209 : // Verified by validation action
210 6 : CPLAssert(aosTokens.size() == 4 || aosTokens.size() == 5);
211 12 : gdal::GCP gcp;
212 6 : gcp.Pixel() = CPLAtof(aosTokens[0]);
213 6 : gcp.Line() = CPLAtof(aosTokens[1]);
214 6 : gcp.X() = CPLAtof(aosTokens[2]);
215 6 : gcp.Y() = CPLAtof(aosTokens[3]);
216 6 : if (aosTokens.size() == 5)
217 3 : gcp.Z() = CPLAtof(aosTokens[4]);
218 6 : ret.push_back(std::move(gcp));
219 : }
220 : }
221 5 : return ret;
222 : }
223 :
224 : /************************************************************************/
225 : /* GDALRasterEditAlgorithm::RunStep() */
226 : /************************************************************************/
227 :
228 33 : bool GDALRasterEditAlgorithm::RunStep(GDALPipelineStepRunContext &ctxt)
229 : {
230 33 : GDALDataset *poDS = m_dataset.GetDatasetRef();
231 33 : if (poDS)
232 : {
233 24 : if (poDS->GetAccess() != GA_Update && !m_readOnly)
234 : {
235 1 : ReportError(CE_Failure, CPLE_AppDefined,
236 : "Dataset should be opened in update mode unless "
237 : "--auxiliary is set");
238 1 : return false;
239 : }
240 : }
241 : else
242 : {
243 9 : const auto poSrcDS = m_inputDataset[0].GetDatasetRef();
244 9 : CPLAssert(poSrcDS);
245 9 : CPLAssert(m_outputDataset.GetName().empty());
246 9 : CPLAssert(!m_outputDataset.GetDatasetRef());
247 :
248 18 : CPLStringList aosOptions;
249 9 : aosOptions.push_back("-of");
250 9 : aosOptions.push_back("VRT");
251 : GDALTranslateOptions *psOptions =
252 9 : GDALTranslateOptionsNew(aosOptions.List(), nullptr);
253 9 : GDALDatasetH hSrcDS = GDALDataset::ToHandle(poSrcDS);
254 9 : auto poRetDS = GDALDataset::FromHandle(
255 : GDALTranslate("", hSrcDS, psOptions, nullptr));
256 9 : GDALTranslateOptionsFree(psOptions);
257 9 : m_outputDataset.Set(std::unique_ptr<GDALDataset>(poRetDS));
258 9 : poDS = m_outputDataset.GetDatasetRef();
259 : }
260 :
261 32 : bool ret = poDS != nullptr;
262 :
263 32 : if (poDS)
264 : {
265 32 : if (m_overrideCrs == "null" || m_overrideCrs == "none")
266 : {
267 3 : if (poDS->SetSpatialRef(nullptr) != CE_None)
268 : {
269 1 : ReportError(CE_Failure, CPLE_AppDefined,
270 : "SetSpatialRef(%s) failed", m_overrideCrs.c_str());
271 10 : return false;
272 : }
273 : }
274 29 : else if (!m_overrideCrs.empty() && m_gcps.empty())
275 : {
276 3 : OGRSpatialReference oSRS;
277 3 : oSRS.SetFromUserInput(m_overrideCrs.c_str());
278 3 : oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
279 3 : if (poDS->SetSpatialRef(&oSRS) != CE_None)
280 : {
281 1 : ReportError(CE_Failure, CPLE_AppDefined,
282 : "SetSpatialRef(%s) failed", m_overrideCrs.c_str());
283 1 : return false;
284 : }
285 : }
286 :
287 30 : if (!m_bbox.empty())
288 : {
289 4 : if (poDS->GetRasterXSize() == 0 || poDS->GetRasterYSize() == 0)
290 : {
291 1 : ReportError(CE_Failure, CPLE_AppDefined,
292 : "Cannot set extent because one of dataset height "
293 : "or width is null");
294 2 : return false;
295 : }
296 : double adfGT[6];
297 3 : adfGT[0] = m_bbox[0];
298 3 : adfGT[1] = (m_bbox[2] - m_bbox[0]) / poDS->GetRasterXSize();
299 3 : adfGT[2] = 0;
300 3 : adfGT[3] = m_bbox[3];
301 3 : adfGT[4] = 0;
302 3 : adfGT[5] = -(m_bbox[3] - m_bbox[1]) / poDS->GetRasterYSize();
303 3 : if (poDS->SetGeoTransform(adfGT) != CE_None)
304 : {
305 1 : ReportError(CE_Failure, CPLE_AppDefined,
306 : "Setting extent failed");
307 1 : return false;
308 : }
309 : }
310 :
311 28 : if (!m_nodata.empty())
312 : {
313 8 : for (int i = 0; i < poDS->GetRasterCount(); ++i)
314 : {
315 4 : if (EQUAL(m_nodata.c_str(), "none"))
316 2 : poDS->GetRasterBand(i + 1)->DeleteNoDataValue();
317 : else
318 4 : poDS->GetRasterBand(i + 1)->SetNoDataValue(
319 2 : CPLAtof(m_nodata.c_str()));
320 : }
321 : }
322 :
323 28 : const CPLStringList aosMD(m_metadata);
324 32 : for (const auto &[key, value] : cpl::IterateNameValue(aosMD))
325 : {
326 5 : if (poDS->SetMetadataItem(key, value) != CE_None)
327 : {
328 1 : ReportError(CE_Failure, CPLE_AppDefined,
329 : "SetMetadataItem('%s', '%s') failed", key, value);
330 1 : return false;
331 : }
332 : }
333 :
334 29 : for (const std::string &key : m_unsetMetadata)
335 : {
336 3 : if (poDS->SetMetadataItem(key.c_str(), nullptr) != CE_None)
337 : {
338 1 : ReportError(CE_Failure, CPLE_AppDefined,
339 : "SetMetadataItem('%s', NULL) failed", key.c_str());
340 1 : return false;
341 : }
342 : }
343 :
344 26 : if (!m_gcps.empty())
345 : {
346 8 : const auto gcps = ParseGCPs();
347 8 : if (gcps.empty())
348 3 : return false; // error already emitted by ParseGCPs()
349 :
350 5 : OGRSpatialReference oSRS;
351 5 : if (!m_overrideCrs.empty())
352 : {
353 1 : oSRS.SetFromUserInput(m_overrideCrs.c_str());
354 1 : oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
355 : }
356 :
357 5 : if (poDS->SetGCPs(static_cast<int>(gcps.size()), gcps[0].c_ptr(),
358 10 : oSRS.IsEmpty() ? nullptr : &oSRS) != CE_None)
359 : {
360 1 : ReportError(CE_Failure, CPLE_AppDefined, "Setting GCPs failed");
361 1 : return false;
362 : }
363 : }
364 :
365 22 : const int nBands = poDS->GetRasterCount();
366 22 : int nCurProgress = 0;
367 22 : const double dfTotalProgress =
368 22 : ((m_stats || m_approxStats) ? nBands : 0) + (m_hist ? nBands : 0);
369 22 : if (m_stats || m_approxStats)
370 : {
371 4 : for (int i = 0; (i < nBands) && ret; ++i)
372 : {
373 4 : void *pScaledProgress = GDALCreateScaledProgress(
374 : nCurProgress / dfTotalProgress,
375 2 : (nCurProgress + 1) / dfTotalProgress, ctxt.m_pfnProgress,
376 : ctxt.m_pProgressData);
377 2 : ++nCurProgress;
378 2 : double dfMin = 0.0;
379 2 : double dfMax = 0.0;
380 2 : double dfMean = 0.0;
381 2 : double dfStdDev = 0.0;
382 4 : ret = poDS->GetRasterBand(i + 1)->ComputeStatistics(
383 2 : m_approxStats, &dfMin, &dfMax, &dfMean, &dfStdDev,
384 : pScaledProgress ? GDALScaledProgress : nullptr,
385 2 : pScaledProgress) == CE_None;
386 2 : GDALDestroyScaledProgress(pScaledProgress);
387 : }
388 : }
389 :
390 22 : if (m_hist)
391 : {
392 2 : for (int i = 0; (i < nBands) && ret; ++i)
393 : {
394 2 : void *pScaledProgress = GDALCreateScaledProgress(
395 : nCurProgress / dfTotalProgress,
396 1 : (nCurProgress + 1) / dfTotalProgress, ctxt.m_pfnProgress,
397 : ctxt.m_pProgressData);
398 1 : ++nCurProgress;
399 1 : double dfMin = 0.0;
400 1 : double dfMax = 0.0;
401 1 : int nBucketCount = 0;
402 1 : GUIntBig *panHistogram = nullptr;
403 2 : ret = poDS->GetRasterBand(i + 1)->GetDefaultHistogram(
404 : &dfMin, &dfMax, &nBucketCount, &panHistogram, TRUE,
405 : pScaledProgress ? GDALScaledProgress : nullptr,
406 1 : pScaledProgress) == CE_None;
407 1 : if (ret)
408 : {
409 2 : ret = poDS->GetRasterBand(i + 1)->SetDefaultHistogram(
410 1 : dfMin, dfMax, nBucketCount, panHistogram) ==
411 : CE_None;
412 : }
413 1 : CPLFree(panHistogram);
414 1 : GDALDestroyScaledProgress(pScaledProgress);
415 : }
416 : }
417 : }
418 :
419 22 : return poDS != nullptr;
420 : }
421 :
422 : GDALRasterEditAlgorithmStandalone::~GDALRasterEditAlgorithmStandalone() =
423 : default;
424 :
425 : //! @endcond
|
