Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: GDAL
4 : * Purpose: gdal "raster pixelinfo" subcommand
5 : * Author: Even Rouault <even dot rouault at spatialys.com>
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2025, Even Rouault <even dot rouault at spatialys.com>
9 : *
10 : * SPDX-License-Identifier: MIT
11 : ****************************************************************************/
12 :
13 : #include "gdalalg_raster_pixel_info.h"
14 :
15 : #include "cpl_conv.h"
16 : #include "cpl_json.h"
17 : #include "cpl_minixml.h"
18 : #include "gdal_priv.h"
19 : #include "ogrsf_frmts.h"
20 : #include "ogr_spatialref.h"
21 :
22 : #include <algorithm>
23 : #include <cmath>
24 : #include <limits>
25 : #include <set>
26 : #include <vector>
27 :
28 : //! @cond Doxygen_Suppress
29 :
30 : #ifndef _
31 : #define _(x) (x)
32 : #endif
33 :
34 : /************************************************************************/
35 : /* GDALRasterPixelInfoAlgorithm::GDALRasterPixelInfoAlgorithm() */
36 : /************************************************************************/
37 :
38 122 : GDALRasterPixelInfoAlgorithm::GDALRasterPixelInfoAlgorithm(bool standaloneStep)
39 : : GDALPipelineStepAlgorithm(
40 : NAME, DESCRIPTION, HELP_URL,
41 0 : ConstructorOptions()
42 122 : .SetStandaloneStep(standaloneStep)
43 122 : .SetAddAppendLayerArgument(false)
44 122 : .SetAddOverwriteLayerArgument(false)
45 122 : .SetAddUpdateArgument(false)
46 122 : .SetAddUpsertArgument(false)
47 122 : .SetAddSkipErrorsArgument(false)
48 366 : .SetOutputFormatCreateCapability(GDAL_DCAP_CREATE))
49 : {
50 122 : if (standaloneStep)
51 : {
52 92 : AddProgressArg(/* hidden = */ true);
53 : AddOutputFormatArg(&m_format, false, false,
54 : _("Output format (default is 'GeoJSON' if "
55 92 : "'position-dataset' not specified)"))
56 : .AddMetadataItem(GAAMDI_REQUIRED_CAPABILITIES,
57 368 : {GDAL_DCAP_VECTOR, GDAL_DCAP_CREATE})
58 92 : .SetAvailableInPipelineStep(false);
59 92 : AddOpenOptionsArg(&m_openOptions).SetAvailableInPipelineStep(false);
60 92 : AddInputFormatsArg(&m_inputFormats)
61 276 : .AddMetadataItem(GAAMDI_REQUIRED_CAPABILITIES, {GDAL_DCAP_RASTER})
62 92 : .SetAvailableInPipelineStep(false);
63 : }
64 :
65 122 : AddInputDatasetArg(&m_inputDataset, GDAL_OF_RASTER)
66 244 : .AddAlias("dataset")
67 122 : .SetMinCount(1)
68 122 : .SetMaxCount(1);
69 :
70 : {
71 : auto &coordinateDatasetArg =
72 : AddArg("position-dataset", '\0',
73 : _("Vector dataset with coordinates"), &m_vectorDataset,
74 244 : GDAL_OF_VECTOR)
75 122 : .SetMutualExclusionGroup("position-dataset-pos");
76 122 : if (!standaloneStep)
77 30 : coordinateDatasetArg.SetPositional();
78 :
79 122 : SetAutoCompleteFunctionForFilename(coordinateDatasetArg,
80 : GDAL_OF_VECTOR);
81 :
82 : auto &layerArg = AddArg(GDAL_ARG_NAME_INPUT_LAYER, 'l',
83 244 : _("Input layer name"), &m_inputLayerNames)
84 122 : .SetMaxCount(1)
85 122 : .AddAlias("layer");
86 122 : SetAutoCompleteFunctionForLayerName(layerArg, coordinateDatasetArg);
87 : }
88 :
89 : AddOutputDatasetArg(&m_outputDataset, GDAL_OF_VECTOR,
90 122 : /* positionalAndRequired = */ false)
91 122 : .SetHiddenForCLI(!standaloneStep);
92 122 : if (standaloneStep)
93 : {
94 92 : AddCreationOptionsArg(&m_creationOptions)
95 92 : .SetAvailableInPipelineStep(false);
96 92 : AddLayerCreationOptionsArg(&m_layerCreationOptions)
97 92 : .SetAvailableInPipelineStep(false);
98 92 : AddOverwriteArg(&m_overwrite).SetAvailableInPipelineStep(false);
99 92 : AddOutputStringArg(&m_output)
100 92 : .SetHiddenForCLI()
101 92 : .SetAvailableInPipelineStep(false);
102 : }
103 :
104 122 : AddBandArg(&m_band);
105 : AddArg("overview", 0, _("Which overview level of source file must be used"),
106 244 : &m_overview)
107 122 : .SetMinValueIncluded(0);
108 :
109 : auto &positionArg =
110 244 : AddArg("position", 'p', _("Pixel position"), &m_pos)
111 244 : .AddAlias("pos")
112 244 : .SetMetaVar("<column,line> or <X,Y>")
113 244 : .SetMutualExclusionGroup("position-dataset-pos")
114 : .AddValidationAction(
115 51 : [this]
116 : {
117 50 : if ((m_pos.size() % 2) != 0)
118 : {
119 1 : ReportError(
120 : CE_Failure, CPLE_IllegalArg,
121 : "An even number of values must be specified "
122 : "for 'position' argument");
123 1 : return false;
124 : }
125 49 : return true;
126 122 : });
127 122 : if (standaloneStep)
128 92 : positionArg.SetPositional();
129 :
130 : AddArg("position-crs", 0,
131 : _("CRS of position (default is 'pixel' if 'position-dataset' not "
132 : "specified)"),
133 244 : &m_posCrs)
134 488 : .SetIsCRSArg(false, {"pixel", "dataset"})
135 122 : .AddHiddenAlias("l_srs");
136 :
137 : AddArg("resampling", 'r', _("Resampling algorithm for interpolation"),
138 244 : &m_resampling)
139 122 : .SetDefault(m_resampling)
140 122 : .SetChoices("nearest", "bilinear", "cubic", "cubicspline")
141 122 : .SetHiddenChoices("near");
142 :
143 : AddArg(
144 : "promote-pixel-value-to-z", 0,
145 : _("Whether to set the pixel value as Z component of GeoJSON geometry"),
146 122 : &m_promotePixelValueToZ);
147 :
148 : AddArg("include-field", 0,
149 : _("Fields from coordinate dataset to include in output (special "
150 : "values: ALL and NONE)"),
151 244 : &m_includeFields)
152 122 : .SetDefault("ALL");
153 :
154 122 : AddValidationAction(
155 184 : [this]
156 : {
157 63 : if (m_inputDataset.size() == 1)
158 : {
159 59 : if (auto poSrcDS = m_inputDataset[0].GetDatasetRef())
160 : {
161 57 : if (poSrcDS->GetRasterCount() > 0)
162 : {
163 : const int nOvrCount =
164 56 : poSrcDS->GetRasterBand(1)->GetOverviewCount();
165 59 : if (m_overview >= 0 && poSrcDS->GetRasterCount() > 0 &&
166 3 : m_overview >= nOvrCount)
167 : {
168 2 : if (nOvrCount == 0)
169 : {
170 1 : ReportError(CE_Failure, CPLE_IllegalArg,
171 : "Source dataset has no overviews. "
172 : "Argument 'overview' must not be "
173 : "specified.");
174 : }
175 : else
176 : {
177 1 : ReportError(
178 : CE_Failure, CPLE_IllegalArg,
179 : "Source dataset has only %d overview "
180 : "level%s. "
181 : "'overview' "
182 : "value must be strictly lower than this "
183 : "number.",
184 : nOvrCount, nOvrCount > 1 ? "s" : "");
185 : }
186 2 : return false;
187 : }
188 : }
189 : else
190 : {
191 1 : ReportError(CE_Failure, CPLE_IllegalArg,
192 : "Source dataset has no raster band.");
193 1 : return false;
194 : }
195 : }
196 : }
197 60 : return true;
198 : });
199 122 : }
200 :
201 : /************************************************************************/
202 : /* GDALRasterPixelInfoAlgorithm::~GDALRasterPixelInfoAlgorithm() */
203 : /************************************************************************/
204 :
205 152 : GDALRasterPixelInfoAlgorithm::~GDALRasterPixelInfoAlgorithm()
206 : {
207 122 : if (!m_osTmpFilename.empty())
208 : {
209 6 : VSIRmdir(CPLGetPathSafe(m_osTmpFilename.c_str()).c_str());
210 6 : VSIUnlink(m_osTmpFilename.c_str());
211 : }
212 152 : }
213 :
214 : /************************************************************************/
215 : /* GDALRasterPixelInfoAlgorithm::RunImpl() */
216 : /************************************************************************/
217 :
218 43 : bool GDALRasterPixelInfoAlgorithm::RunImpl(GDALProgressFunc pfnProgress,
219 : void *pProgressData)
220 : {
221 43 : GDALPipelineStepRunContext stepCtxt;
222 43 : stepCtxt.m_pfnProgress = pfnProgress;
223 43 : stepCtxt.m_pProgressData = pProgressData;
224 43 : return RunPreStepPipelineValidations() && RunStep(stepCtxt);
225 : }
226 :
227 : /************************************************************************/
228 : /* GDALRasterPixelInfoAlgorithm::RunStep() */
229 : /************************************************************************/
230 :
231 49 : bool GDALRasterPixelInfoAlgorithm::RunStep(GDALPipelineStepRunContext &)
232 : {
233 49 : auto poSrcDS = m_inputDataset[0].GetDatasetRef();
234 49 : CPLAssert(poSrcDS);
235 :
236 49 : auto poVectorSrcDS = m_vectorDataset.GetDatasetRef();
237 :
238 49 : if (m_pos.empty() && !poVectorSrcDS && !IsCalledFromCommandLine())
239 : {
240 2 : ReportError(
241 : CE_Failure, CPLE_AppDefined,
242 : "Argument 'position' or 'position-dataset' must be specified.");
243 2 : return false;
244 : }
245 :
246 47 : if (!m_standaloneStep)
247 : {
248 5 : m_format = "MEM";
249 : }
250 42 : else if (m_outputDataset.GetName().empty() && m_format != "MEM")
251 : {
252 25 : if (m_format.empty())
253 : {
254 16 : m_format = "GeoJSON";
255 : }
256 9 : else if (!EQUAL(m_format.c_str(), "CSV") &&
257 0 : !EQUAL(m_format.c_str(), "GeoJSON"))
258 : {
259 0 : ReportError(CE_Failure, CPLE_AppDefined,
260 : "Only CSV or GeoJSON output format is allowed when "
261 : "'output' dataset is not specified.");
262 0 : return false;
263 : }
264 : }
265 17 : else if (m_format.empty())
266 : {
267 : const std::string osExt =
268 26 : CPLGetExtensionSafe(m_outputDataset.GetName().c_str());
269 13 : if (EQUAL(osExt.c_str(), "csv"))
270 1 : m_format = "CSV";
271 12 : else if (EQUAL(osExt.c_str(), "json"))
272 1 : m_format = "GeoJSON";
273 : }
274 :
275 47 : const bool bIsCSV = EQUAL(m_format.c_str(), "CSV");
276 47 : const bool bIsGeoJSON = EQUAL(m_format.c_str(), "GeoJSON");
277 :
278 47 : OGRLayer *poSrcLayer = nullptr;
279 94 : std::vector<int> anSrcFieldIndicesToInclude;
280 94 : std::vector<int> anMapSrcToDstFields;
281 :
282 47 : if (poVectorSrcDS)
283 : {
284 20 : if (m_inputLayerNames.empty())
285 : {
286 18 : const int nLayerCount = poVectorSrcDS->GetLayerCount();
287 18 : if (nLayerCount == 0)
288 : {
289 2 : ReportError(CE_Failure, CPLE_AppDefined,
290 : "Dataset '%s' has no vector layer",
291 2 : poVectorSrcDS->GetDescription());
292 2 : return false;
293 : }
294 16 : else if (nLayerCount > 1)
295 : {
296 1 : ReportError(CE_Failure, CPLE_AppDefined,
297 : "Dataset '%s' has more than one vector layer. "
298 : "Please specify the 'input-layer' argument",
299 1 : poVectorSrcDS->GetDescription());
300 1 : return false;
301 : }
302 15 : poSrcLayer = poVectorSrcDS->GetLayer(0);
303 : }
304 : else
305 : {
306 : poSrcLayer =
307 2 : poVectorSrcDS->GetLayerByName(m_inputLayerNames[0].c_str());
308 : }
309 17 : if (!poSrcLayer)
310 : {
311 2 : ReportError(
312 : CE_Failure, CPLE_AppDefined, "Cannot find layer '%s' in '%s'",
313 2 : m_inputLayerNames.empty() ? "" : m_inputLayerNames[0].c_str(),
314 1 : poVectorSrcDS->GetDescription());
315 1 : return false;
316 : }
317 16 : if (poSrcLayer->GetGeomType() == wkbNone)
318 : {
319 2 : ReportError(CE_Failure, CPLE_AppDefined,
320 : "Layer '%s' of '%s' has no geometry column",
321 2 : poSrcLayer->GetName(), poVectorSrcDS->GetDescription());
322 2 : return false;
323 : }
324 :
325 14 : if (!GetFieldIndices(m_includeFields, OGRLayer::ToHandle(poSrcLayer),
326 : anSrcFieldIndicesToInclude))
327 : {
328 1 : return false;
329 : }
330 :
331 13 : if (m_posCrs.empty())
332 : {
333 10 : const auto poVectorLayerSRS = poSrcLayer->GetSpatialRef();
334 10 : if (poVectorLayerSRS)
335 : {
336 4 : const char *const apszOptions[] = {"FORMAT=WKT2", nullptr};
337 4 : m_posCrs = poVectorLayerSRS->exportToWkt(apszOptions);
338 : }
339 : }
340 : }
341 :
342 40 : if (m_posCrs.empty())
343 28 : m_posCrs = "pixel";
344 :
345 : const GDALRIOResampleAlg eInterpolation =
346 40 : GDALRasterIOGetResampleAlg(m_resampling.c_str());
347 :
348 40 : const auto poSrcCRS = poSrcDS->GetSpatialRef();
349 40 : GDALGeoTransform gt;
350 40 : const bool bHasGT = poSrcDS->GetGeoTransform(gt) == CE_None;
351 40 : GDALGeoTransform invGT;
352 :
353 40 : if (m_posCrs != "pixel")
354 : {
355 10 : if (!poSrcCRS)
356 : {
357 1 : ReportError(CE_Failure, CPLE_AppDefined,
358 : "Dataset has no CRS. Only 'position-crs' = 'pixel' is "
359 : "supported.");
360 1 : return false;
361 : }
362 :
363 9 : if (!bHasGT)
364 : {
365 1 : ReportError(CE_Failure, CPLE_AppDefined, "Cannot get geotransform");
366 1 : return false;
367 : }
368 :
369 8 : if (!gt.GetInverse(invGT))
370 : {
371 1 : ReportError(CE_Failure, CPLE_AppDefined,
372 : "Cannot invert geotransform");
373 1 : return false;
374 : }
375 : }
376 :
377 37 : std::unique_ptr<OGRCoordinateTransformation> poCT;
378 74 : OGRSpatialReference oUserCRS;
379 37 : if (m_posCrs != "pixel" && m_posCrs != "dataset")
380 : {
381 6 : oUserCRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
382 : // Already validated due SetIsCRSArg()
383 6 : CPL_IGNORE_RET_VAL(oUserCRS.SetFromUserInput(m_posCrs.c_str()));
384 6 : poCT.reset(OGRCreateCoordinateTransformation(&oUserCRS, poSrcCRS));
385 6 : if (!poCT)
386 0 : return false;
387 : }
388 :
389 37 : if (m_band.empty())
390 : {
391 74 : for (int i = 1; i <= poSrcDS->GetRasterCount(); ++i)
392 37 : m_band.push_back(i);
393 : }
394 :
395 37 : std::unique_ptr<GDALDataset> poOutDS;
396 37 : OGRLayer *poOutLayer = nullptr;
397 37 : std::unique_ptr<OGRCoordinateTransformation> poCTSrcCRSToOutCRS;
398 :
399 : const bool isInteractive =
400 61 : m_pos.empty() && m_outputDataset.GetName().empty() &&
401 61 : IsCalledFromCommandLine() && CPLIsInteractive(stdin);
402 :
403 74 : std::string osOutFilename = m_outputDataset.GetName();
404 37 : if (osOutFilename.empty() && bIsCSV)
405 : {
406 6 : if (isInteractive)
407 : {
408 0 : osOutFilename = "/vsistdout/";
409 : }
410 : else
411 : {
412 6 : osOutFilename = VSIMemGenerateHiddenFilename("subdir");
413 6 : osOutFilename += "/out.csv";
414 6 : VSIMkdir(CPLGetPathSafe(osOutFilename.c_str()).c_str(), 0755);
415 6 : m_osTmpFilename = osOutFilename;
416 : }
417 : }
418 :
419 37 : if (!osOutFilename.empty() || m_format == "MEM" || !m_standaloneStep)
420 : {
421 25 : if (bIsGeoJSON)
422 : {
423 1 : m_outputFile.reset(VSIFOpenL(osOutFilename.c_str(), "wb"));
424 1 : if (!m_outputFile)
425 : {
426 0 : ReportError(CE_Failure, CPLE_FileIO, "Cannot create '%s'",
427 : osOutFilename.c_str());
428 0 : return false;
429 : }
430 : }
431 : else
432 : {
433 24 : if (m_format.empty())
434 : {
435 : const auto aosFormats =
436 : CPLStringList(GDALGetOutputDriversForDatasetName(
437 : osOutFilename.c_str(), GDAL_OF_VECTOR,
438 : /* bSingleMatch = */ true,
439 10 : /* bWarn = */ true));
440 10 : if (aosFormats.size() != 1)
441 : {
442 1 : ReportError(CE_Failure, CPLE_AppDefined,
443 : "Cannot guess driver for %s",
444 : osOutFilename.c_str());
445 1 : return false;
446 : }
447 9 : m_format = aosFormats[0];
448 : }
449 :
450 : auto poOutDrv =
451 23 : GetGDALDriverManager()->GetDriverByName(m_format.c_str());
452 23 : if (!poOutDrv)
453 : {
454 0 : ReportError(CE_Failure, CPLE_AppDefined,
455 : "Cannot find driver %s", m_format.c_str());
456 0 : return false;
457 : }
458 23 : poOutDS.reset(
459 : poOutDrv->Create(osOutFilename.c_str(), 0, 0, 0, GDT_Unknown,
460 46 : CPLStringList(m_creationOptions).List()));
461 23 : if (!poOutDS)
462 1 : return false;
463 :
464 22 : std::string osOutLayerName;
465 22 : if (EQUAL(m_format.c_str(), "ESRI Shapefile") || bIsCSV ||
466 : !poSrcLayer)
467 : {
468 10 : osOutLayerName = CPLGetBasenameSafe(osOutFilename.c_str());
469 : }
470 : else
471 12 : osOutLayerName = poSrcLayer->GetName();
472 :
473 22 : const OGRSpatialReference *poOutCRS = nullptr;
474 22 : if (!oUserCRS.IsEmpty())
475 5 : poOutCRS = &oUserCRS;
476 17 : else if (poSrcLayer)
477 : {
478 7 : poOutCRS = poSrcLayer->GetSpatialRef();
479 7 : if (!poOutCRS)
480 7 : poOutCRS = poSrcCRS;
481 : }
482 44 : poOutLayer = poOutDS->CreateLayer(
483 : osOutLayerName.c_str(), poOutCRS, wkbPoint,
484 44 : CPLStringList(m_layerCreationOptions).List());
485 22 : if (!poOutLayer)
486 : {
487 1 : ReportError(CE_Failure, CPLE_AppDefined,
488 : "Cannot create layer '%s' in '%s'",
489 : osOutLayerName.c_str(), osOutFilename.c_str());
490 1 : return false;
491 : }
492 21 : if (poSrcCRS && poOutCRS)
493 : {
494 12 : poCTSrcCRSToOutCRS.reset(
495 : OGRCreateCoordinateTransformation(poSrcCRS, poOutCRS));
496 12 : if (!poCTSrcCRSToOutCRS)
497 0 : return false;
498 : }
499 : }
500 : }
501 :
502 34 : if (poOutLayer)
503 : {
504 21 : bool bOK = true;
505 :
506 21 : if (bIsCSV)
507 : {
508 16 : OGRFieldDefn oFieldLine("geom_x", OFTReal);
509 8 : bOK &= poOutLayer->CreateField(&oFieldLine) == OGRERR_NONE;
510 :
511 8 : OGRFieldDefn oFieldColumn("geom_y", OFTReal);
512 8 : bOK &= poOutLayer->CreateField(&oFieldColumn) == OGRERR_NONE;
513 : }
514 :
515 21 : if (poSrcLayer)
516 : {
517 24 : CPLErrorStateBackuper oBackuper(CPLQuietErrorHandler);
518 12 : const auto poSrcLayerDefn = poSrcLayer->GetLayerDefn();
519 12 : auto poOutLayerDefn = poOutLayer->GetLayerDefn();
520 :
521 12 : anMapSrcToDstFields.resize(poSrcLayerDefn->GetFieldCount(), -1);
522 :
523 19 : for (int nSrcIdx : anSrcFieldIndicesToInclude)
524 : {
525 : const auto poSrcFieldDefn =
526 7 : poSrcLayerDefn->GetFieldDefn(nSrcIdx);
527 7 : if (poOutLayer->CreateField(poSrcFieldDefn) == OGRERR_NONE)
528 : {
529 7 : const int nDstIdx = poOutLayerDefn->GetFieldIndex(
530 7 : poSrcFieldDefn->GetNameRef());
531 7 : if (nDstIdx >= 0)
532 7 : anMapSrcToDstFields[nSrcIdx] = nDstIdx;
533 : }
534 : }
535 : }
536 : else
537 : {
538 9 : OGRFieldDefn oFieldExtraInput("extra_content", OFTString);
539 9 : bOK &= poOutLayer->CreateField(&oFieldExtraInput) == OGRERR_NONE;
540 : }
541 :
542 21 : OGRFieldDefn oFieldColumn("column", OFTReal);
543 21 : bOK &= poOutLayer->CreateField(&oFieldColumn) == OGRERR_NONE;
544 :
545 21 : OGRFieldDefn oFieldLine("line", OFTReal);
546 21 : bOK &= poOutLayer->CreateField(&oFieldLine) == OGRERR_NONE;
547 :
548 42 : for (int nBand : m_band)
549 : {
550 : auto hBand =
551 21 : GDALRasterBand::ToHandle(poSrcDS->GetRasterBand(nBand));
552 21 : const bool bIsComplex = CPL_TO_BOOL(
553 : GDALDataTypeIsComplex(GDALGetRasterDataType(hBand)));
554 21 : if (bIsComplex)
555 : {
556 : OGRFieldDefn oFieldReal(CPLSPrintf("band_%d_real_value", nBand),
557 6 : OFTReal);
558 3 : bOK &= poOutLayer->CreateField(&oFieldReal) == OGRERR_NONE;
559 :
560 : OGRFieldDefn oFieldImag(
561 3 : CPLSPrintf("band_%d_imaginary_value", nBand), OFTReal);
562 3 : bOK &= poOutLayer->CreateField(&oFieldImag) == OGRERR_NONE;
563 : }
564 : else
565 : {
566 : OGRFieldDefn oFieldRaw(CPLSPrintf("band_%d_raw_value", nBand),
567 36 : OFTReal);
568 18 : bOK &= poOutLayer->CreateField(&oFieldRaw) == OGRERR_NONE;
569 :
570 : OGRFieldDefn oFieldUnscaled(
571 18 : CPLSPrintf("band_%d_unscaled_value", nBand), OFTReal);
572 18 : bOK &= poOutLayer->CreateField(&oFieldUnscaled) == OGRERR_NONE;
573 : }
574 : }
575 :
576 21 : if (!bOK)
577 : {
578 0 : ReportError(CE_Failure, CPLE_AppDefined,
579 : "Cannot create fields in output layer");
580 0 : return false;
581 : }
582 : }
583 :
584 68 : CPLJSONObject oCollection;
585 34 : oCollection.Add("type", "FeatureCollection");
586 34 : std::unique_ptr<OGRCoordinateTransformation> poCTToWGS84;
587 34 : bool canOutputGeoJSONGeom = false;
588 34 : if (poSrcCRS && bHasGT)
589 : {
590 32 : const char *pszAuthName = poSrcCRS->GetAuthorityName();
591 32 : const char *pszAuthCode = poSrcCRS->GetAuthorityCode();
592 32 : if (pszAuthName && pszAuthCode && EQUAL(pszAuthName, "EPSG"))
593 : {
594 31 : canOutputGeoJSONGeom = true;
595 62 : CPLJSONObject jCRS;
596 31 : CPLJSONObject oProperties;
597 31 : if (EQUAL(pszAuthCode, "4326"))
598 2 : oProperties.Add("name", "urn:ogc:def:crs:OGC:1.3:CRS84");
599 : else
600 29 : oProperties.Add(
601 : "name",
602 58 : std::string("urn:ogc:def:crs:EPSG::").append(pszAuthCode));
603 31 : jCRS.Add("type", "name");
604 31 : jCRS.Add("properties", oProperties);
605 62 : oCollection.Add("crs", jCRS);
606 : }
607 : else
608 : {
609 2 : OGRSpatialReference oCRS;
610 1 : oCRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
611 1 : oCRS.importFromEPSG(4326);
612 1 : poCTToWGS84.reset(
613 : OGRCreateCoordinateTransformation(poSrcCRS, &oCRS));
614 1 : if (poCTToWGS84)
615 : {
616 1 : canOutputGeoJSONGeom = true;
617 2 : CPLJSONObject jCRS;
618 1 : CPLJSONObject oProperties;
619 1 : oProperties.Add("name", "urn:ogc:def:crs:OGC:1.3:CRS84");
620 1 : jCRS.Add("type", "name");
621 1 : jCRS.Add("properties", oProperties);
622 1 : oCollection.Add("crs", jCRS);
623 : }
624 : }
625 : }
626 68 : CPLJSONArray oFeatures;
627 34 : oCollection.Add("features", oFeatures);
628 :
629 : char szLine[1024];
630 34 : int nLine = 0;
631 34 : size_t iVal = 0;
632 16 : do
633 : {
634 50 : double x = 0, y = 0;
635 50 : std::string osExtraContent;
636 0 : std::unique_ptr<OGRFeature> poSrcFeature;
637 50 : if (poSrcLayer)
638 : {
639 26 : poSrcFeature.reset(poSrcLayer->GetNextFeature());
640 26 : if (!poSrcFeature)
641 13 : break;
642 13 : const auto poGeom = poSrcFeature->GetGeometryRef();
643 13 : if (!poGeom || poGeom->IsEmpty())
644 0 : continue;
645 13 : if (wkbFlatten(poGeom->getGeometryType()) == wkbPoint)
646 : {
647 13 : x = poGeom->toPoint()->getX();
648 13 : y = poGeom->toPoint()->getY();
649 : }
650 : else
651 : {
652 0 : OGRPoint oPoint;
653 0 : if (poGeom->Centroid(&oPoint) != OGRERR_NONE)
654 0 : return false;
655 0 : x = oPoint.getX();
656 0 : y = oPoint.getY();
657 : }
658 : }
659 24 : else if (iVal + 1 < m_pos.size())
660 : {
661 18 : x = m_pos[iVal++];
662 18 : y = m_pos[iVal++];
663 : }
664 : else
665 : {
666 6 : if (CPLIsInteractive(stdin))
667 : {
668 0 : if (m_posCrs != "pixel")
669 : {
670 0 : fprintf(stderr, "Enter X Y values separated by space, and "
671 : "press Return.\n");
672 : }
673 : else
674 : {
675 0 : fprintf(stderr,
676 : "Enter pixel line values separated by space, "
677 : "and press Return.\n");
678 : }
679 : }
680 :
681 6 : if (fgets(szLine, sizeof(szLine) - 1, stdin) && szLine[0] != '\n')
682 : {
683 3 : const CPLStringList aosTokens(CSLTokenizeString(szLine));
684 3 : const int nCount = aosTokens.size();
685 :
686 3 : ++nLine;
687 3 : if (nCount < 2)
688 : {
689 0 : fprintf(stderr, "Not enough values at line %d\n", nLine);
690 0 : return false;
691 : }
692 : else
693 : {
694 3 : x = CPLAtof(aosTokens[0]);
695 3 : y = CPLAtof(aosTokens[1]);
696 :
697 9 : for (int i = 2; i < nCount; ++i)
698 : {
699 6 : if (!osExtraContent.empty())
700 3 : osExtraContent += ' ';
701 6 : osExtraContent += aosTokens[i];
702 : }
703 6 : while (!osExtraContent.empty() &&
704 3 : isspace(static_cast<int>(osExtraContent.back())))
705 : {
706 0 : osExtraContent.pop_back();
707 : }
708 : }
709 : }
710 : else
711 : {
712 3 : break;
713 : }
714 : }
715 :
716 34 : const double xOri = x;
717 34 : const double yOri = y;
718 34 : double dfPixel{0}, dfLine{0};
719 :
720 34 : if (poCT)
721 : {
722 6 : if (!poCT->Transform(1, &x, &y, nullptr))
723 0 : return false;
724 : }
725 :
726 34 : if (m_posCrs != "pixel")
727 : {
728 7 : invGT.Apply(x, y, &dfPixel, &dfLine);
729 : }
730 : else
731 : {
732 27 : dfPixel = x;
733 27 : dfLine = y;
734 : }
735 : const int iPixel = static_cast<int>(
736 68 : std::clamp(std::floor(dfPixel), static_cast<double>(INT_MIN),
737 34 : static_cast<double>(INT_MAX)));
738 : const int iLine = static_cast<int>(
739 68 : std::clamp(std::floor(dfLine), static_cast<double>(INT_MIN),
740 34 : static_cast<double>(INT_MAX)));
741 :
742 34 : CPLJSONObject oFeature;
743 34 : CPLJSONObject oProperties;
744 0 : std::unique_ptr<OGRFeature> poFeature;
745 34 : if (bIsGeoJSON)
746 : {
747 13 : oFeature.Add("type", "Feature");
748 13 : oFeature.Add("properties", oProperties);
749 : {
750 13 : CPLJSONArray oArray;
751 13 : oArray.Add(xOri);
752 13 : oArray.Add(yOri);
753 13 : oProperties.Add("input_coordinate", oArray);
754 : }
755 13 : if (!osExtraContent.empty())
756 1 : oProperties.Add("extra_content", osExtraContent);
757 13 : oProperties.Add("column", dfPixel);
758 13 : oProperties.Add("line", dfLine);
759 :
760 13 : if (poSrcFeature)
761 : {
762 8 : for (int i : anSrcFieldIndicesToInclude)
763 : {
764 7 : const auto *poFieldDefn = poSrcFeature->GetFieldDefnRef(i);
765 7 : const char *pszName = poFieldDefn->GetNameRef();
766 7 : const auto eType = poFieldDefn->GetType();
767 7 : switch (eType)
768 : {
769 3 : case OFTInteger:
770 : case OFTInteger64:
771 : {
772 3 : if (poFieldDefn->GetSubType() == OFSTBoolean)
773 : {
774 1 : oProperties.Add(
775 : pszName,
776 1 : poSrcFeature->GetFieldAsInteger(i) != 0);
777 : }
778 : else
779 : {
780 2 : oProperties.Add(
781 : pszName,
782 : poSrcFeature->GetFieldAsInteger64(i));
783 : }
784 3 : break;
785 : }
786 :
787 1 : case OFTReal:
788 : {
789 1 : oProperties.Add(pszName,
790 : poSrcFeature->GetFieldAsDouble(i));
791 1 : break;
792 : }
793 :
794 2 : case OFTString:
795 : {
796 2 : if (poFieldDefn->GetSubType() != OFSTJSON)
797 : {
798 1 : oProperties.Add(
799 : pszName, poSrcFeature->GetFieldAsString(i));
800 1 : break;
801 : }
802 : else
803 : {
804 : [[fallthrough]];
805 : }
806 : }
807 :
808 : default:
809 : {
810 : char *pszJSON =
811 2 : poSrcFeature->GetFieldAsSerializedJSon(i);
812 4 : CPLJSONDocument oDoc;
813 2 : if (oDoc.LoadMemory(pszJSON))
814 2 : oProperties.Add(pszName, oDoc.GetRoot());
815 2 : CPLFree(pszJSON);
816 2 : break;
817 : }
818 : }
819 : }
820 : }
821 : }
822 : else
823 : {
824 21 : CPLAssert(poOutLayer);
825 : poFeature =
826 21 : std::make_unique<OGRFeature>(poOutLayer->GetLayerDefn());
827 :
828 21 : if (poSrcFeature)
829 : {
830 24 : CPLErrorStateBackuper oBackuper(CPLQuietErrorHandler);
831 24 : poFeature->SetFrom(poSrcFeature.get(),
832 12 : anMapSrcToDstFields.data());
833 : }
834 9 : else if (!osExtraContent.empty())
835 : {
836 2 : poFeature->SetField("extra_content", osExtraContent.c_str());
837 : }
838 :
839 21 : if (poOutLayer->GetSpatialRef() == nullptr)
840 : {
841 9 : if (bIsCSV)
842 : {
843 8 : poFeature->SetField("geom_x", xOri);
844 8 : poFeature->SetField("geom_y", yOri);
845 : }
846 : else
847 : {
848 2 : poFeature->SetGeometry(
849 2 : std::make_unique<OGRPoint>(xOri, yOri));
850 : }
851 : }
852 12 : else if (bHasGT && poCTSrcCRSToOutCRS)
853 : {
854 12 : gt.Apply(dfPixel, dfLine, &x, &y);
855 12 : if (poCTSrcCRSToOutCRS->Transform(1, &x, &y))
856 : {
857 12 : if (bIsCSV)
858 : {
859 0 : poFeature->SetField("geom_x", x);
860 0 : poFeature->SetField("geom_y", y);
861 : }
862 : else
863 : {
864 24 : poFeature->SetGeometry(
865 24 : std::make_unique<OGRPoint>(x, y));
866 : }
867 : }
868 : }
869 :
870 21 : poFeature->SetField("column", dfPixel);
871 21 : poFeature->SetField("line", dfLine);
872 : }
873 :
874 34 : CPLJSONArray oBands;
875 :
876 34 : double zValue = std::numeric_limits<double>::quiet_NaN();
877 68 : for (int nBand : m_band)
878 : {
879 34 : CPLJSONObject oBand;
880 34 : oBand.Add("band_number", nBand);
881 :
882 : auto hBand =
883 34 : GDALRasterBand::ToHandle(poSrcDS->GetRasterBand(nBand));
884 :
885 34 : int iPixelToQuery = iPixel;
886 34 : int iLineToQuery = iLine;
887 :
888 34 : double dfPixelToQuery = dfPixel;
889 34 : double dfLineToQuery = dfLine;
890 :
891 34 : if (m_overview >= 0 && hBand != nullptr)
892 : {
893 1 : GDALRasterBandH hOvrBand = GDALGetOverview(hBand, m_overview);
894 1 : if (hOvrBand != nullptr)
895 : {
896 1 : int nOvrXSize = GDALGetRasterBandXSize(hOvrBand);
897 1 : int nOvrYSize = GDALGetRasterBandYSize(hOvrBand);
898 1 : iPixelToQuery = static_cast<int>(
899 1 : 0.5 +
900 1 : 1.0 * iPixel / poSrcDS->GetRasterXSize() * nOvrXSize);
901 1 : iLineToQuery = static_cast<int>(
902 1 : 0.5 +
903 1 : 1.0 * iLine / poSrcDS->GetRasterYSize() * nOvrYSize);
904 1 : if (iPixelToQuery >= nOvrXSize)
905 0 : iPixelToQuery = nOvrXSize - 1;
906 1 : if (iLineToQuery >= nOvrYSize)
907 0 : iLineToQuery = nOvrYSize - 1;
908 1 : dfPixelToQuery =
909 1 : dfPixel / poSrcDS->GetRasterXSize() * nOvrXSize;
910 1 : dfLineToQuery =
911 1 : dfLine / poSrcDS->GetRasterYSize() * nOvrYSize;
912 : }
913 : else
914 : {
915 0 : ReportError(CE_Failure, CPLE_AppDefined,
916 : "Cannot get overview %d of band %d", m_overview,
917 : nBand);
918 0 : return false;
919 : }
920 1 : hBand = hOvrBand;
921 : }
922 :
923 34 : double adfPixel[2] = {0, 0};
924 34 : const bool bIsComplex = CPL_TO_BOOL(
925 : GDALDataTypeIsComplex(GDALGetRasterDataType(hBand)));
926 : int bIgnored;
927 34 : const double dfOffset = GDALGetRasterOffset(hBand, &bIgnored);
928 34 : const double dfScale = GDALGetRasterScale(hBand, &bIgnored);
929 34 : if (GDALRasterInterpolateAtPoint(
930 : hBand, dfPixelToQuery, dfLineToQuery, eInterpolation,
931 34 : &adfPixel[0], &adfPixel[1]) == CE_None)
932 : {
933 32 : if (!bIsComplex)
934 : {
935 29 : const double dfUnscaledVal =
936 29 : adfPixel[0] * dfScale + dfOffset;
937 29 : if (m_band.size() == 1 && m_promotePixelValueToZ)
938 1 : zValue = dfUnscaledVal;
939 29 : if (bIsGeoJSON)
940 : {
941 12 : if (GDALDataTypeIsInteger(GDALGetRasterDataType(hBand)))
942 : {
943 11 : oBand.Add("raw_value",
944 11 : static_cast<GInt64>(adfPixel[0]));
945 : }
946 : else
947 : {
948 1 : oBand.Add("raw_value", adfPixel[0]);
949 : }
950 :
951 12 : oBand.Add("unscaled_value", dfUnscaledVal);
952 : }
953 : else
954 : {
955 17 : poFeature->SetField(
956 : CPLSPrintf("band_%d_raw_value", nBand),
957 : adfPixel[0]);
958 17 : poFeature->SetField(
959 : CPLSPrintf("band_%d_unscaled_value", nBand),
960 : dfUnscaledVal);
961 : }
962 : }
963 : else
964 : {
965 3 : if (bIsGeoJSON)
966 : {
967 1 : CPLJSONObject oValue;
968 1 : oValue.Add("real", adfPixel[0]);
969 1 : oValue.Add("imaginary", adfPixel[1]);
970 1 : oBand.Add("value", oValue);
971 : }
972 : else
973 : {
974 2 : poFeature->SetField(
975 : CPLSPrintf("band_%d_real_value", nBand),
976 : adfPixel[0]);
977 2 : poFeature->SetField(
978 : CPLSPrintf("band_%d_imaginary_value", nBand),
979 : adfPixel[1]);
980 : }
981 : }
982 : }
983 :
984 : // Request location info for this location (just a few drivers,
985 : // like the VRT driver actually supports this).
986 68 : CPLString osItem;
987 34 : osItem.Printf("Pixel_%d_%d", iPixelToQuery, iLineToQuery);
988 :
989 34 : if (const char *pszLI =
990 34 : GDALGetMetadataItem(hBand, osItem, "LocationInfo"))
991 : {
992 4 : CPLXMLTreeCloser oTree(CPLParseXMLString(pszLI));
993 :
994 4 : if (oTree && oTree->psChild != nullptr &&
995 6 : oTree->eType == CXT_Element &&
996 2 : EQUAL(oTree->pszValue, "LocationInfo"))
997 : {
998 2 : CPLJSONArray oFiles;
999 :
1000 2 : for (const CPLXMLNode *psNode = oTree->psChild;
1001 4 : psNode != nullptr; psNode = psNode->psNext)
1002 : {
1003 2 : if (psNode->eType == CXT_Element &&
1004 2 : EQUAL(psNode->pszValue, "File") &&
1005 2 : psNode->psChild != nullptr)
1006 : {
1007 4 : char *pszUnescaped = CPLUnescapeString(
1008 2 : psNode->psChild->pszValue, nullptr, CPLES_XML);
1009 2 : oFiles.Add(pszUnescaped);
1010 2 : CPLFree(pszUnescaped);
1011 : }
1012 : }
1013 :
1014 2 : oBand.Add("files", oFiles);
1015 : }
1016 : else
1017 : {
1018 0 : oBand.Add("location_info", pszLI);
1019 : }
1020 : }
1021 :
1022 34 : oBands.Add(oBand);
1023 : }
1024 :
1025 34 : if (bIsGeoJSON)
1026 : {
1027 13 : oProperties.Add("bands", oBands);
1028 :
1029 13 : if (canOutputGeoJSONGeom)
1030 : {
1031 12 : x = dfPixel;
1032 12 : y = dfLine;
1033 :
1034 12 : gt.Apply(x, y, &x, &y);
1035 :
1036 12 : if (poCTToWGS84)
1037 1 : poCTToWGS84->Transform(1, &x, &y);
1038 :
1039 24 : CPLJSONObject oGeometry;
1040 12 : oFeature.Add("geometry", oGeometry);
1041 12 : oGeometry.Add("type", "Point");
1042 12 : CPLJSONArray oCoordinates;
1043 12 : oCoordinates.Add(x);
1044 12 : oCoordinates.Add(y);
1045 12 : if (!std::isnan(zValue))
1046 1 : oCoordinates.Add(zValue);
1047 12 : oGeometry.Add("coordinates", oCoordinates);
1048 : }
1049 : else
1050 : {
1051 1 : oFeature.AddNull("geometry");
1052 : }
1053 :
1054 13 : if (isInteractive)
1055 : {
1056 0 : CPLJSONDocument oDoc;
1057 0 : oDoc.SetRoot(oFeature);
1058 0 : printf("%s\n", oDoc.SaveAsString().c_str());
1059 : }
1060 : else
1061 : {
1062 13 : oFeatures.Add(oFeature);
1063 : }
1064 : }
1065 : else
1066 : {
1067 21 : if (poOutLayer->CreateFeature(std::move(poFeature)) != OGRERR_NONE)
1068 : {
1069 0 : return false;
1070 : }
1071 : }
1072 :
1073 34 : } while (m_pos.empty() || iVal + 1 < m_pos.size());
1074 :
1075 34 : if (bIsGeoJSON && !isInteractive)
1076 : {
1077 26 : CPLJSONDocument oDoc;
1078 13 : oDoc.SetRoot(oCollection);
1079 26 : std::string osRet = oDoc.SaveAsString();
1080 13 : if (m_outputFile)
1081 1 : m_outputFile->Write(osRet.data(), osRet.size());
1082 : else
1083 25 : m_output = std::move(osRet);
1084 : }
1085 21 : else if (poOutDS)
1086 : {
1087 27 : if (m_format != "MEM" && m_outputDataset.GetName().empty() &&
1088 6 : m_standaloneStep)
1089 : {
1090 6 : poOutDS.reset();
1091 6 : if (!isInteractive)
1092 : {
1093 6 : CPLAssert(!m_osTmpFilename.empty());
1094 : const GByte *pabyData =
1095 6 : VSIGetMemFileBuffer(m_osTmpFilename.c_str(), nullptr,
1096 : /* bUnlinkAndSeize = */ false);
1097 6 : m_output = reinterpret_cast<const char *>(pabyData);
1098 : }
1099 : }
1100 : else
1101 : {
1102 15 : m_outputDataset.Set(std::move(poOutDS));
1103 : }
1104 : }
1105 :
1106 34 : bool bRet = true;
1107 34 : if (m_outputFile)
1108 : {
1109 1 : bRet = m_outputFile->Close() == 0;
1110 1 : if (bRet)
1111 : {
1112 1 : poOutDS.reset(
1113 1 : GDALDataset::Open(m_outputDataset.GetName().c_str(),
1114 : GDAL_OF_VECTOR | GDAL_OF_VERBOSE_ERROR));
1115 1 : bRet = poOutDS != nullptr;
1116 1 : m_outputDataset.Set(std::move(poOutDS));
1117 : }
1118 : }
1119 :
1120 34 : return bRet;
1121 : }
1122 :
1123 : GDALRasterPixelInfoAlgorithmStandalone::
1124 : ~GDALRasterPixelInfoAlgorithmStandalone() = default;
1125 :
1126 : //! @endcond
|
