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