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 "ogr_spatialref.h"
20 :
21 : #include <algorithm>
22 : #include <cmath>
23 : #include <limits>
24 :
25 : //! @cond Doxygen_Suppress
26 :
27 : #ifndef _
28 : #define _(x) (x)
29 : #endif
30 :
31 : /************************************************************************/
32 : /* GDALRasterPixelInfoAlgorithm::GDALRasterPixelInfoAlgorithm() */
33 : /************************************************************************/
34 :
35 32 : GDALRasterPixelInfoAlgorithm::GDALRasterPixelInfoAlgorithm()
36 32 : : GDALAlgorithm(NAME, DESCRIPTION, HELP_URL)
37 : {
38 32 : AddOutputFormatArg(&m_format)
39 32 : .SetDefault("geojson")
40 32 : .SetChoices("geojson", "csv")
41 32 : .SetHiddenChoices("json")
42 32 : .SetDefault(m_format);
43 32 : AddOpenOptionsArg(&m_openOptions);
44 32 : AddInputFormatsArg(&m_inputFormats)
45 64 : .AddMetadataItem(GAAMDI_REQUIRED_CAPABILITIES, {GDAL_DCAP_RASTER});
46 32 : AddInputDatasetArg(&m_dataset, GDAL_OF_RASTER).AddAlias("dataset");
47 32 : AddOutputStringArg(&m_output);
48 :
49 32 : AddBandArg(&m_band);
50 : AddArg("overview", 0, _("Which overview level of source file must be used"),
51 64 : &m_overview)
52 32 : .SetMinValueIncluded(0);
53 :
54 64 : AddArg("position", 'p', _("Pixel position"), &m_pos)
55 64 : .AddAlias("pos")
56 64 : .SetMetaVar("<column,line> or <X,Y>")
57 32 : .SetPositional()
58 : .AddValidationAction(
59 38 : [this]
60 : {
61 37 : if ((m_pos.size() % 2) != 0)
62 : {
63 1 : ReportError(CE_Failure, CPLE_IllegalArg,
64 : "An even number of values must be specified "
65 : "for 'position' argument");
66 1 : return false;
67 : }
68 36 : return true;
69 32 : });
70 64 : AddArg("position-crs", 0, _("CRS of position"), &m_posCrs)
71 128 : .SetIsCRSArg(false, {"pixel", "dataset"})
72 32 : .SetDefault("pixel")
73 32 : .AddHiddenAlias("l_srs");
74 :
75 : AddArg("resampling", 'r', _("Resampling algorithm for interpolation"),
76 64 : &m_resampling)
77 32 : .SetDefault(m_resampling)
78 32 : .SetChoices("nearest", "bilinear", "cubic", "cubicspline")
79 32 : .SetHiddenChoices("near");
80 :
81 : AddArg(
82 : "promote-pixel-value-to-z", 0,
83 : _("Whether to set the pixel value as Z component of GeoJSON geometry"),
84 32 : &m_promotePixelValueToZ);
85 :
86 32 : AddValidationAction(
87 55 : [this]
88 : {
89 25 : if (auto poSrcDS = m_dataset.GetDatasetRef())
90 : {
91 : const int nOvrCount =
92 25 : poSrcDS->GetRasterBand(1)->GetOverviewCount();
93 28 : if (m_overview >= 0 && poSrcDS->GetRasterCount() > 0 &&
94 3 : m_overview >= nOvrCount)
95 : {
96 2 : if (nOvrCount == 0)
97 : {
98 1 : ReportError(
99 : CE_Failure, CPLE_IllegalArg,
100 : "Source dataset has no overviews. "
101 : "Argument 'overview' must not be specified.");
102 : }
103 : else
104 : {
105 1 : ReportError(
106 : CE_Failure, CPLE_IllegalArg,
107 : "Source dataset has only %d overview level%s. "
108 : "'overview' "
109 : "value must be strictly lower than this number.",
110 : nOvrCount, nOvrCount > 1 ? "s" : "");
111 : }
112 2 : return false;
113 : }
114 : }
115 23 : return true;
116 : });
117 32 : }
118 :
119 : /************************************************************************/
120 : /* GDALRasterPixelInfoAlgorithm::PrintLine() */
121 : /************************************************************************/
122 :
123 12 : void GDALRasterPixelInfoAlgorithm::PrintLine(const std::string &str)
124 : {
125 12 : if (IsCalledFromCommandLine())
126 : {
127 2 : printf("%s\n", str.c_str());
128 : }
129 : else
130 : {
131 10 : m_output += str;
132 10 : m_output += '\n';
133 : }
134 12 : }
135 :
136 : /************************************************************************/
137 : /* GDALRasterPixelInfoAlgorithm::RunImpl() */
138 : /************************************************************************/
139 :
140 21 : bool GDALRasterPixelInfoAlgorithm::RunImpl(GDALProgressFunc, void *)
141 : {
142 21 : auto poSrcDS = m_dataset.GetDatasetRef();
143 21 : CPLAssert(poSrcDS);
144 :
145 21 : if (m_pos.empty() && !IsCalledFromCommandLine())
146 : {
147 1 : ReportError(CE_Failure, CPLE_AppDefined,
148 : "Argument 'position' must be specified.");
149 1 : return false;
150 : }
151 :
152 : const GDALRIOResampleAlg eInterpolation =
153 20 : GDALRasterIOGetResampleAlg(m_resampling.c_str());
154 :
155 20 : const auto poSrcCRS = poSrcDS->GetSpatialRef();
156 20 : GDALGeoTransform gt;
157 20 : const bool bHasGT = poSrcDS->GetGeoTransform(gt) == CE_None;
158 20 : GDALGeoTransform invGT;
159 :
160 20 : if (m_posCrs != "pixel")
161 : {
162 5 : if (!poSrcCRS)
163 : {
164 1 : ReportError(CE_Failure, CPLE_AppDefined,
165 : "Dataset has no CRS. Only 'position-crs' = 'pixel' is "
166 : "supported.");
167 1 : return false;
168 : }
169 :
170 4 : if (!bHasGT)
171 : {
172 1 : ReportError(CE_Failure, CPLE_AppDefined, "Cannot get geotransform");
173 1 : return false;
174 : }
175 :
176 3 : if (!gt.GetInverse(invGT))
177 : {
178 1 : ReportError(CE_Failure, CPLE_AppDefined,
179 : "Cannot invert geotransform");
180 1 : return false;
181 : }
182 : }
183 :
184 17 : std::unique_ptr<OGRCoordinateTransformation> poCT;
185 17 : if (m_posCrs != "pixel" && m_posCrs != "dataset")
186 : {
187 1 : OGRSpatialReference oUserCRS;
188 1 : oUserCRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
189 : // Already validated due SetIsCRSArg()
190 1 : CPL_IGNORE_RET_VAL(oUserCRS.SetFromUserInput(m_posCrs.c_str()));
191 1 : poCT.reset(OGRCreateCoordinateTransformation(&oUserCRS, poSrcCRS));
192 1 : if (!poCT)
193 0 : return false;
194 : }
195 :
196 17 : if (m_band.empty())
197 : {
198 34 : for (int i = 1; i <= poSrcDS->GetRasterCount(); ++i)
199 17 : m_band.push_back(i);
200 : }
201 :
202 17 : if (m_format == "csv")
203 : {
204 12 : std::string line = "input_x,input_y,extra_input,column,line";
205 12 : for (int nBand : m_band)
206 : {
207 : auto hBand =
208 6 : GDALRasterBand::ToHandle(poSrcDS->GetRasterBand(nBand));
209 6 : const bool bIsComplex = CPL_TO_BOOL(
210 : GDALDataTypeIsComplex(GDALGetRasterDataType(hBand)));
211 6 : if (bIsComplex)
212 : line +=
213 : CPLSPrintf(",band_%d_real_value,band_%d_imaginary_value",
214 2 : nBand, nBand);
215 : else
216 : {
217 : line += CPLSPrintf(",band_%d_raw_value,band_%d_unscaled_value",
218 4 : nBand, nBand);
219 : }
220 : }
221 6 : PrintLine(line);
222 : }
223 :
224 : const bool isInteractive =
225 17 : m_pos.empty() && IsCalledFromCommandLine() && CPLIsInteractive(stdin);
226 :
227 34 : CPLJSONObject oCollection;
228 17 : oCollection.Add("type", "FeatureCollection");
229 17 : std::unique_ptr<OGRCoordinateTransformation> poCTToWGS84;
230 17 : bool canOutputGeoJSONGeom = false;
231 17 : if (poSrcCRS && bHasGT)
232 : {
233 15 : const char *pszAuthName = poSrcCRS->GetAuthorityName(nullptr);
234 15 : const char *pszAuthCode = poSrcCRS->GetAuthorityCode(nullptr);
235 15 : if (pszAuthName && pszAuthCode && EQUAL(pszAuthName, "EPSG"))
236 : {
237 14 : canOutputGeoJSONGeom = true;
238 28 : CPLJSONObject jCRS;
239 14 : CPLJSONObject oProperties;
240 14 : if (EQUAL(pszAuthCode, "4326"))
241 2 : oProperties.Add("name", "urn:ogc:def:crs:OGC:1.3:CRS84");
242 : else
243 12 : oProperties.Add(
244 : "name",
245 24 : std::string("urn:ogc:def:crs:EPSG::").append(pszAuthCode));
246 14 : jCRS.Add("type", "name");
247 14 : jCRS.Add("properties", oProperties);
248 28 : oCollection.Add("crs", jCRS);
249 : }
250 : else
251 : {
252 2 : OGRSpatialReference oCRS;
253 1 : oCRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
254 1 : oCRS.importFromEPSG(4326);
255 1 : poCTToWGS84.reset(
256 : OGRCreateCoordinateTransformation(poSrcCRS, &oCRS));
257 1 : if (poCTToWGS84)
258 : {
259 1 : canOutputGeoJSONGeom = true;
260 2 : CPLJSONObject jCRS;
261 1 : CPLJSONObject oProperties;
262 1 : oProperties.Add("name", "urn:ogc:def:crs:OGC:1.3:CRS84");
263 1 : jCRS.Add("type", "name");
264 1 : jCRS.Add("properties", oProperties);
265 1 : oCollection.Add("crs", jCRS);
266 : }
267 : }
268 : }
269 34 : CPLJSONArray oFeatures;
270 17 : oCollection.Add("features", oFeatures);
271 :
272 : char szLine[1024];
273 17 : int nLine = 0;
274 17 : size_t iVal = 0;
275 2 : do
276 : {
277 19 : double x = 0, y = 0;
278 19 : std::string osExtraContent;
279 19 : if (iVal + 1 < m_pos.size())
280 : {
281 15 : x = m_pos[iVal++];
282 15 : y = m_pos[iVal++];
283 : }
284 : else
285 : {
286 4 : if (CPLIsInteractive(stdin))
287 : {
288 0 : if (m_posCrs != "pixel")
289 : {
290 0 : fprintf(stderr, "Enter X Y values separated by space, and "
291 : "press Return.\n");
292 : }
293 : else
294 : {
295 0 : fprintf(stderr,
296 : "Enter pixel line values separated by space, "
297 : "and press Return.\n");
298 : }
299 : }
300 :
301 4 : if (fgets(szLine, sizeof(szLine) - 1, stdin) && szLine[0] != '\n')
302 : {
303 2 : const CPLStringList aosTokens(CSLTokenizeString(szLine));
304 2 : const int nCount = aosTokens.size();
305 :
306 2 : ++nLine;
307 2 : if (nCount < 2)
308 : {
309 0 : fprintf(stderr, "Not enough values at line %d\n", nLine);
310 0 : return false;
311 : }
312 : else
313 : {
314 2 : x = CPLAtof(aosTokens[0]);
315 2 : y = CPLAtof(aosTokens[1]);
316 :
317 6 : for (int i = 2; i < nCount; ++i)
318 : {
319 4 : if (!osExtraContent.empty())
320 2 : osExtraContent += ' ';
321 4 : osExtraContent += aosTokens[i];
322 : }
323 4 : while (!osExtraContent.empty() &&
324 2 : isspace(static_cast<int>(osExtraContent.back())))
325 : {
326 0 : osExtraContent.pop_back();
327 : }
328 : }
329 : }
330 : else
331 : {
332 2 : break;
333 : }
334 : }
335 :
336 17 : const double xOri = x;
337 17 : const double yOri = y;
338 17 : double dfPixel{0}, dfLine{0};
339 :
340 17 : if (poCT)
341 : {
342 1 : if (!poCT->Transform(1, &x, &y, nullptr))
343 0 : return false;
344 : }
345 :
346 17 : if (m_posCrs != "pixel")
347 : {
348 2 : invGT.Apply(x, y, &dfPixel, &dfLine);
349 : }
350 : else
351 : {
352 15 : dfPixel = x;
353 15 : dfLine = y;
354 : }
355 : const int iPixel = static_cast<int>(
356 34 : std::clamp(std::floor(dfPixel), static_cast<double>(INT_MIN),
357 17 : static_cast<double>(INT_MAX)));
358 : const int iLine = static_cast<int>(
359 34 : std::clamp(std::floor(dfLine), static_cast<double>(INT_MIN),
360 17 : static_cast<double>(INT_MAX)));
361 :
362 17 : std::string line;
363 17 : CPLJSONObject oFeature;
364 17 : CPLJSONObject oProperties;
365 17 : if (m_format == "csv")
366 : {
367 6 : line = CPLSPrintf("%.17g,%.17g", xOri, yOri);
368 6 : line += ",\"";
369 6 : line += CPLString(osExtraContent).replaceAll('"', "\"\"");
370 6 : line += '"';
371 6 : line += CPLSPrintf(",%.17g,%.17g", dfPixel, dfLine);
372 : }
373 : else
374 : {
375 11 : oFeature.Add("type", "Feature");
376 11 : oFeature.Add("properties", oProperties);
377 : {
378 11 : CPLJSONArray oArray;
379 11 : oArray.Add(xOri);
380 11 : oArray.Add(yOri);
381 11 : oProperties.Add("input_coordinate", oArray);
382 : }
383 11 : if (!osExtraContent.empty())
384 1 : oProperties.Add("extra_content", osExtraContent);
385 11 : oProperties.Add("column", dfPixel);
386 11 : oProperties.Add("line", dfLine);
387 : }
388 :
389 17 : CPLJSONArray oBands;
390 :
391 17 : double zValue = std::numeric_limits<double>::quiet_NaN();
392 34 : for (int nBand : m_band)
393 : {
394 17 : CPLJSONObject oBand;
395 17 : oBand.Add("band_number", nBand);
396 :
397 : auto hBand =
398 17 : GDALRasterBand::ToHandle(poSrcDS->GetRasterBand(nBand));
399 :
400 17 : int iPixelToQuery = iPixel;
401 17 : int iLineToQuery = iLine;
402 :
403 17 : double dfPixelToQuery = dfPixel;
404 17 : double dfLineToQuery = dfLine;
405 :
406 17 : if (m_overview >= 0 && hBand != nullptr)
407 : {
408 1 : GDALRasterBandH hOvrBand = GDALGetOverview(hBand, m_overview);
409 1 : if (hOvrBand != nullptr)
410 : {
411 1 : int nOvrXSize = GDALGetRasterBandXSize(hOvrBand);
412 1 : int nOvrYSize = GDALGetRasterBandYSize(hOvrBand);
413 1 : iPixelToQuery = static_cast<int>(
414 1 : 0.5 +
415 1 : 1.0 * iPixel / poSrcDS->GetRasterXSize() * nOvrXSize);
416 1 : iLineToQuery = static_cast<int>(
417 1 : 0.5 +
418 1 : 1.0 * iLine / poSrcDS->GetRasterYSize() * nOvrYSize);
419 1 : if (iPixelToQuery >= nOvrXSize)
420 0 : iPixelToQuery = nOvrXSize - 1;
421 1 : if (iLineToQuery >= nOvrYSize)
422 0 : iLineToQuery = nOvrYSize - 1;
423 1 : dfPixelToQuery =
424 1 : dfPixel / poSrcDS->GetRasterXSize() * nOvrXSize;
425 1 : dfLineToQuery =
426 1 : dfLine / poSrcDS->GetRasterYSize() * nOvrYSize;
427 : }
428 : else
429 : {
430 0 : ReportError(CE_Failure, CPLE_AppDefined,
431 : "Cannot get overview %d of band %d", m_overview,
432 : nBand);
433 0 : return false;
434 : }
435 1 : hBand = hOvrBand;
436 : }
437 :
438 17 : double adfPixel[2] = {0, 0};
439 17 : const bool bIsComplex = CPL_TO_BOOL(
440 : GDALDataTypeIsComplex(GDALGetRasterDataType(hBand)));
441 : int bIgnored;
442 17 : const double dfOffset = GDALGetRasterOffset(hBand, &bIgnored);
443 17 : const double dfScale = GDALGetRasterScale(hBand, &bIgnored);
444 17 : if (GDALRasterInterpolateAtPoint(
445 : hBand, dfPixelToQuery, dfLineToQuery, eInterpolation,
446 17 : &adfPixel[0], &adfPixel[1]) == CE_None)
447 : {
448 15 : if (!bIsComplex)
449 : {
450 13 : const double dfUnscaledVal =
451 13 : adfPixel[0] * dfScale + dfOffset;
452 13 : if (m_band.size() == 1 && m_promotePixelValueToZ)
453 1 : zValue = dfUnscaledVal;
454 13 : if (m_format == "csv")
455 : {
456 3 : line += CPLSPrintf(",%.17g", adfPixel[0]);
457 3 : line += CPLSPrintf(",%.17g", dfUnscaledVal);
458 : }
459 : else
460 : {
461 10 : if (GDALDataTypeIsInteger(GDALGetRasterDataType(hBand)))
462 : {
463 9 : oBand.Add("raw_value",
464 9 : static_cast<GInt64>(adfPixel[0]));
465 : }
466 : else
467 : {
468 1 : oBand.Add("raw_value", adfPixel[0]);
469 : }
470 :
471 10 : oBand.Add("unscaled_value", dfUnscaledVal);
472 : }
473 : }
474 : else
475 : {
476 2 : if (m_format == "csv")
477 : {
478 : line += CPLSPrintf(",%.17g,%.17g", adfPixel[0],
479 1 : adfPixel[1]);
480 : }
481 : else
482 : {
483 1 : CPLJSONObject oValue;
484 1 : oValue.Add("real", adfPixel[0]);
485 1 : oValue.Add("imaginary", adfPixel[1]);
486 1 : oBand.Add("value", oValue);
487 : }
488 : }
489 : }
490 2 : else if (m_format == "csv")
491 : {
492 2 : line += ",,";
493 : }
494 :
495 : // Request location info for this location (just a few drivers,
496 : // like the VRT driver actually supports this).
497 34 : CPLString osItem;
498 17 : osItem.Printf("Pixel_%d_%d", iPixelToQuery, iLineToQuery);
499 :
500 17 : if (const char *pszLI =
501 17 : GDALGetMetadataItem(hBand, osItem, "LocationInfo"))
502 : {
503 4 : CPLXMLTreeCloser oTree(CPLParseXMLString(pszLI));
504 :
505 4 : if (oTree && oTree->psChild != nullptr &&
506 6 : oTree->eType == CXT_Element &&
507 2 : EQUAL(oTree->pszValue, "LocationInfo"))
508 : {
509 2 : CPLJSONArray oFiles;
510 :
511 2 : for (const CPLXMLNode *psNode = oTree->psChild;
512 4 : psNode != nullptr; psNode = psNode->psNext)
513 : {
514 2 : if (psNode->eType == CXT_Element &&
515 2 : EQUAL(psNode->pszValue, "File") &&
516 2 : psNode->psChild != nullptr)
517 : {
518 4 : char *pszUnescaped = CPLUnescapeString(
519 2 : psNode->psChild->pszValue, nullptr, CPLES_XML);
520 2 : oFiles.Add(pszUnescaped);
521 2 : CPLFree(pszUnescaped);
522 : }
523 : }
524 :
525 2 : oBand.Add("files", oFiles);
526 : }
527 : else
528 : {
529 0 : oBand.Add("location_info", pszLI);
530 : }
531 : }
532 :
533 17 : oBands.Add(oBand);
534 : }
535 :
536 17 : if (m_format == "csv")
537 : {
538 6 : PrintLine(line);
539 : }
540 : else
541 : {
542 11 : oProperties.Add("bands", oBands);
543 :
544 11 : if (canOutputGeoJSONGeom)
545 : {
546 10 : x = dfPixel;
547 10 : y = dfLine;
548 :
549 10 : gt.Apply(x, y, &x, &y);
550 :
551 10 : if (poCTToWGS84)
552 1 : poCTToWGS84->Transform(1, &x, &y);
553 :
554 20 : CPLJSONObject oGeometry;
555 10 : oFeature.Add("geometry", oGeometry);
556 10 : oGeometry.Add("type", "Point");
557 10 : CPLJSONArray oCoordinates;
558 10 : oCoordinates.Add(x);
559 10 : oCoordinates.Add(y);
560 10 : if (!std::isnan(zValue))
561 1 : oCoordinates.Add(zValue);
562 10 : oGeometry.Add("coordinates", oCoordinates);
563 : }
564 : else
565 : {
566 1 : oFeature.AddNull("geometry");
567 : }
568 :
569 11 : if (isInteractive)
570 : {
571 0 : CPLJSONDocument oDoc;
572 0 : oDoc.SetRoot(oFeature);
573 0 : printf("%s\n", oDoc.SaveAsString().c_str());
574 : }
575 : else
576 : {
577 11 : oFeatures.Add(oFeature);
578 : }
579 : }
580 :
581 17 : } while (m_pos.empty() || iVal + 1 < m_pos.size());
582 :
583 17 : if (m_format != "csv" && !isInteractive)
584 : {
585 11 : CPLJSONDocument oDoc;
586 11 : oDoc.SetRoot(oCollection);
587 11 : m_output = oDoc.SaveAsString();
588 : }
589 :
590 17 : return true;
591 : }
592 :
593 : //! @endcond
|
