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