1 /******************************************************************************
2 *
3 * Project: GDAL
4 * Purpose: Command line raster query tool.
5 * Author: Frank Warmerdam <warmerdam@pobox.com>
6 *
7 ******************************************************************************
8 * Copyright (c) 2010, Frank Warmerdam <warmerdam@pobox.com>
9 * Copyright (c) 2010-2013, Even Rouault <even dot rouault at spatialys.com>
10 *
11 * Permission is hereby granted, free of charge, to any person obtaining a
12 * copy of this software and associated documentation files (the "Software"),
13 * to deal in the Software without restriction, including without limitation
14 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
15 * and/or sell copies of the Software, and to permit persons to whom the
16 * Software is furnished to do so, subject to the following conditions:
17 *
18 * The above copyright notice and this permission notice shall be included
19 * in all copies or substantial portions of the Software.
20 *
21 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
22 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
23 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
24 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
25 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
26 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
27 * DEALINGS IN THE SOFTWARE.
28 ****************************************************************************/
29
30 #include "cpl_string.h"
31 #include "cpl_minixml.h"
32 #include "gdal_version.h"
33 #include "gdal.h"
34 #include "commonutils.h"
35 #include "ogr_spatialref.h"
36 #include <vector>
37
38 #ifdef _WIN32
39 #include <io.h>
40 #else
41 #include <unistd.h>
42 #endif
43
44 CPL_CVSID("$Id: gdallocationinfo.cpp 8c3e4ef55212f20eec95aa7e12ba5d48dacfdc47 2020-10-01 21:20:51 +0200 Even Rouault $")
45
46 /************************************************************************/
47 /* Usage() */
48 /************************************************************************/
49
Usage()50 static void Usage()
51
52 {
53 printf( "Usage: gdallocationinfo [--help-general] [-xml] [-lifonly] [-valonly]\n"
54 " [-b band]* [-overview overview_level]\n"
55 " [-l_srs srs_def] [-geoloc] [-wgs84]\n"
56 " [-oo NAME=VALUE]* srcfile x y\n"
57 "\n" );
58 exit( 1 );
59 }
60
61 /************************************************************************/
62 /* SanitizeSRS */
63 /************************************************************************/
64
SanitizeSRS(const char * pszUserInput)65 static char *SanitizeSRS( const char *pszUserInput )
66
67 {
68 CPLErrorReset();
69
70 OGRSpatialReferenceH hSRS = OSRNewSpatialReference( nullptr );
71
72 char *pszResult = nullptr;
73 if( OSRSetFromUserInput( hSRS, pszUserInput ) == OGRERR_NONE )
74 OSRExportToWkt( hSRS, &pszResult );
75 else
76 {
77 CPLError( CE_Failure, CPLE_AppDefined,
78 "Translating source or target SRS failed:\n%s",
79 pszUserInput );
80 exit( 1 );
81 }
82
83 OSRDestroySpatialReference( hSRS );
84
85 return pszResult;
86 }
87
88 /************************************************************************/
89 /* main() */
90 /************************************************************************/
91
MAIN_START(argc,argv)92 MAIN_START(argc, argv)
93
94 {
95 const char *pszLocX = nullptr, *pszLocY = nullptr;
96 const char *pszSrcFilename = nullptr;
97 char *pszSourceSRS = nullptr;
98 std::vector<int> anBandList;
99 bool bAsXML = false, bLIFOnly = false;
100 bool bQuiet = false, bValOnly = false;
101 int nOverview = -1;
102 char **papszOpenOptions = nullptr;
103
104 GDALAllRegister();
105 argc = GDALGeneralCmdLineProcessor( argc, &argv, 0 );
106 if( argc < 1 )
107 exit( -argc );
108
109 /* -------------------------------------------------------------------- */
110 /* Parse arguments. */
111 /* -------------------------------------------------------------------- */
112 for( int i = 1; i < argc; i++ )
113 {
114 if( EQUAL(argv[i], "--utility_version") )
115 {
116 printf("%s was compiled against GDAL %s and is running against GDAL %s\n",
117 argv[0], GDAL_RELEASE_NAME, GDALVersionInfo("RELEASE_NAME"));
118 GDALDestroyDriverManager();
119 CSLDestroy(argv);
120 return 0;
121 }
122 else if( i < argc-1 && EQUAL(argv[i],"-b") )
123 {
124 anBandList.push_back( atoi(argv[++i]) );
125 }
126 else if( i < argc-1 && EQUAL(argv[i],"-overview") )
127 {
128 nOverview = atoi(argv[++i]) - 1;
129 }
130 else if( i < argc-1 && EQUAL(argv[i],"-l_srs") )
131 {
132 CPLFree(pszSourceSRS);
133 // coverity[tainted_data]
134 pszSourceSRS = SanitizeSRS(argv[++i]);
135 }
136 else if( EQUAL(argv[i],"-geoloc") )
137 {
138 CPLFree(pszSourceSRS);
139 pszSourceSRS = CPLStrdup("-geoloc");
140 }
141 else if( EQUAL(argv[i],"-wgs84") )
142 {
143 CPLFree(pszSourceSRS);
144 pszSourceSRS = SanitizeSRS("WGS84");
145 }
146 else if( EQUAL(argv[i],"-xml") )
147 {
148 bAsXML = true;
149 }
150 else if( EQUAL(argv[i],"-lifonly") )
151 {
152 bLIFOnly = true;
153 bQuiet = true;
154 }
155 else if( EQUAL(argv[i],"-valonly") )
156 {
157 bValOnly = true;
158 bQuiet = true;
159 }
160 else if( i < argc-1 && EQUAL(argv[i], "-oo") )
161 {
162 papszOpenOptions = CSLAddString( papszOpenOptions,
163 argv[++i] );
164 }
165 else if( argv[i][0] == '-' && !isdigit(argv[i][1]) )
166 Usage();
167
168 else if( pszSrcFilename == nullptr )
169 pszSrcFilename = argv[i];
170
171 else if( pszLocX == nullptr )
172 pszLocX = argv[i];
173
174 else if( pszLocY == nullptr )
175 pszLocY = argv[i];
176
177 else
178 Usage();
179 }
180
181 if( pszSrcFilename == nullptr || (pszLocX != nullptr && pszLocY == nullptr) )
182 Usage();
183
184 /* -------------------------------------------------------------------- */
185 /* Open source file. */
186 /* -------------------------------------------------------------------- */
187 GDALDatasetH hSrcDS
188 = GDALOpenEx( pszSrcFilename, GDAL_OF_RASTER | GDAL_OF_VERBOSE_ERROR,
189 nullptr,
190 papszOpenOptions, nullptr );
191 if( hSrcDS == nullptr )
192 exit( 1 );
193
194 /* -------------------------------------------------------------------- */
195 /* Setup coordinate transformation, if required */
196 /* -------------------------------------------------------------------- */
197 OGRSpatialReferenceH hSrcSRS = nullptr;
198 OGRCoordinateTransformationH hCT = nullptr;
199 if( pszSourceSRS != nullptr && !EQUAL(pszSourceSRS,"-geoloc") )
200 {
201
202 hSrcSRS = OSRNewSpatialReference( pszSourceSRS );
203 OSRSetAxisMappingStrategy(hSrcSRS, OAMS_TRADITIONAL_GIS_ORDER);
204 auto hTrgSRS = GDALGetSpatialRef( hSrcDS );
205 if( !hTrgSRS )
206 exit(1);
207
208 hCT = OCTNewCoordinateTransformation( hSrcSRS, hTrgSRS );
209 if( hCT == nullptr )
210 exit( 1 );
211 }
212
213 /* -------------------------------------------------------------------- */
214 /* If no bands were requested, we will query them all. */
215 /* -------------------------------------------------------------------- */
216 if( anBandList.empty() )
217 {
218 for( int i = 0; i < GDALGetRasterCount( hSrcDS ); i++ )
219 anBandList.push_back( i+1 );
220 }
221
222 /* -------------------------------------------------------------------- */
223 /* Turn the location into a pixel and line location. */
224 /* -------------------------------------------------------------------- */
225 int inputAvailable = 1;
226 double dfGeoX;
227 double dfGeoY;
228 CPLString osXML;
229
230 if( pszLocX == nullptr && pszLocY == nullptr )
231 {
232 // Is it an interactive terminal ?
233 if( isatty(static_cast<int>(fileno(stdin))) )
234 {
235 if( pszSourceSRS != nullptr )
236 {
237 fprintf(stderr, "Enter X Y values separated by space, and press Return.\n");
238 }
239 else
240 {
241 fprintf(stderr, "Enter pixel line values separated by space, and press Return.\n");
242 }
243 }
244
245 if (fscanf(stdin, "%lf %lf", &dfGeoX, &dfGeoY) != 2)
246 {
247 inputAvailable = 0;
248 }
249 }
250 else
251 {
252 dfGeoX = CPLAtof(pszLocX);
253 dfGeoY = CPLAtof(pszLocY);
254 }
255
256 while (inputAvailable)
257 {
258 int iPixel, iLine;
259
260 if (hCT)
261 {
262 if( !OCTTransform( hCT, 1, &dfGeoX, &dfGeoY, nullptr ) )
263 exit( 1 );
264 }
265
266 if( pszSourceSRS != nullptr )
267 {
268 double adfGeoTransform[6] = {};
269 if( GDALGetGeoTransform( hSrcDS, adfGeoTransform ) != CE_None )
270 {
271 CPLError(CE_Failure, CPLE_AppDefined, "Cannot get geotransform");
272 exit( 1 );
273 }
274
275 double adfInvGeoTransform[6] = {};
276 if( !GDALInvGeoTransform( adfGeoTransform, adfInvGeoTransform ) )
277 {
278 CPLError(CE_Failure, CPLE_AppDefined, "Cannot invert geotransform");
279 exit( 1 );
280 }
281
282 iPixel = static_cast<int>(floor(
283 adfInvGeoTransform[0]
284 + adfInvGeoTransform[1] * dfGeoX
285 + adfInvGeoTransform[2] * dfGeoY));
286 iLine = static_cast<int>(floor(
287 adfInvGeoTransform[3]
288 + adfInvGeoTransform[4] * dfGeoX
289 + adfInvGeoTransform[5] * dfGeoY));
290 }
291 else
292 {
293 iPixel = static_cast<int>(floor(dfGeoX));
294 iLine = static_cast<int>(floor(dfGeoY));
295 }
296
297 /* -------------------------------------------------------------------- */
298 /* Prepare report. */
299 /* -------------------------------------------------------------------- */
300 CPLString osLine;
301
302 if( bAsXML )
303 {
304 osLine.Printf( "<Report pixel=\"%d\" line=\"%d\">",
305 iPixel, iLine );
306 osXML += osLine;
307 }
308 else if( !bQuiet )
309 {
310 printf( "Report:\n" );
311 printf( " Location: (%dP,%dL)\n", iPixel, iLine );
312 }
313
314 bool bPixelReport = true;
315
316 if( iPixel < 0 || iLine < 0
317 || iPixel >= GDALGetRasterXSize( hSrcDS )
318 || iLine >= GDALGetRasterYSize( hSrcDS ) )
319 {
320 if( bAsXML )
321 osXML += "<Alert>Location is off this file! No further details to report.</Alert>";
322 else if( bValOnly )
323 printf("\n");
324 else if( !bQuiet )
325 printf( "\nLocation is off this file! No further details to report.\n");
326 bPixelReport = false;
327 }
328
329 /* -------------------------------------------------------------------- */
330 /* Process each band. */
331 /* -------------------------------------------------------------------- */
332 for( int i = 0; bPixelReport && i < static_cast<int>(anBandList.size());
333 i++ )
334 {
335 GDALRasterBandH hBand = GDALGetRasterBand( hSrcDS, anBandList[i] );
336
337 int iPixelToQuery = iPixel;
338 int iLineToQuery = iLine;
339
340 if (nOverview >= 0 && hBand != nullptr)
341 {
342 GDALRasterBandH hOvrBand = GDALGetOverview(hBand, nOverview);
343 if (hOvrBand != nullptr)
344 {
345 int nOvrXSize = GDALGetRasterBandXSize(hOvrBand);
346 int nOvrYSize = GDALGetRasterBandYSize(hOvrBand);
347 iPixelToQuery = static_cast<int>(
348 0.5 +
349 1.0 * iPixel / GDALGetRasterXSize(hSrcDS) * nOvrXSize);
350 iLineToQuery = static_cast<int>(
351 0.5 +
352 1.0 * iLine / GDALGetRasterYSize(hSrcDS) * nOvrYSize);
353 if (iPixelToQuery >= nOvrXSize)
354 iPixelToQuery = nOvrXSize - 1;
355 if (iLineToQuery >= nOvrYSize)
356 iLineToQuery = nOvrYSize - 1;
357 }
358 else
359 {
360 CPLError(CE_Failure, CPLE_AppDefined,
361 "Cannot get overview %d of band %d",
362 nOverview + 1, anBandList[i] );
363 }
364 hBand = hOvrBand;
365 }
366
367 if (hBand == nullptr)
368 continue;
369
370 if( bAsXML )
371 {
372 osLine.Printf( "<BandReport band=\"%d\">", anBandList[i] );
373 osXML += osLine;
374 }
375 else if( !bQuiet )
376 {
377 printf( " Band %d:\n", anBandList[i] );
378 }
379
380 /* -------------------------------------------------------------------- */
381 /* Request location info for this location. It is possible */
382 /* only the VRT driver actually supports this. */
383 /* -------------------------------------------------------------------- */
384 CPLString osItem;
385
386 osItem.Printf( "Pixel_%d_%d", iPixelToQuery, iLineToQuery );
387
388 const char *pszLI = GDALGetMetadataItem( hBand, osItem, "LocationInfo");
389
390 if( pszLI != nullptr )
391 {
392 if( bAsXML )
393 osXML += pszLI;
394 else if( !bQuiet )
395 printf( " %s\n", pszLI );
396 else if( bLIFOnly )
397 {
398 /* Extract all files, if any. */
399
400 CPLXMLNode *psRoot = CPLParseXMLString( pszLI );
401
402 if( psRoot != nullptr
403 && psRoot->psChild != nullptr
404 && psRoot->eType == CXT_Element
405 && EQUAL(psRoot->pszValue,"LocationInfo") )
406 {
407 for( CPLXMLNode *psNode = psRoot->psChild;
408 psNode != nullptr;
409 psNode = psNode->psNext )
410 {
411 if( psNode->eType == CXT_Element
412 && EQUAL(psNode->pszValue,"File")
413 && psNode->psChild != nullptr )
414 {
415 char* pszUnescaped = CPLUnescapeString(
416 psNode->psChild->pszValue, nullptr, CPLES_XML);
417 printf( "%s\n", pszUnescaped );
418 CPLFree(pszUnescaped);
419 }
420 }
421 }
422 CPLDestroyXMLNode( psRoot );
423 }
424 }
425
426 /* -------------------------------------------------------------------- */
427 /* Report the pixel value of this band. */
428 /* -------------------------------------------------------------------- */
429 double adfPixel[2];
430
431 if( GDALRasterIO( hBand, GF_Read, iPixelToQuery, iLineToQuery, 1, 1,
432 adfPixel, 1, 1, GDT_CFloat64, 0, 0) == CE_None )
433 {
434 CPLString osValue;
435
436 if( GDALDataTypeIsComplex( GDALGetRasterDataType( hBand ) ) )
437 osValue.Printf( "%.15g+%.15gi", adfPixel[0], adfPixel[1] );
438 else
439 osValue.Printf( "%.15g", adfPixel[0] );
440
441 if( bAsXML )
442 {
443 osXML += "<Value>";
444 osXML += osValue;
445 osXML += "</Value>";
446 }
447 else if( !bQuiet )
448 printf( " Value: %s\n", osValue.c_str() );
449 else if( bValOnly )
450 printf( "%s\n", osValue.c_str() );
451
452 // Report unscaled if we have scale/offset values.
453 int bSuccess;
454
455 double dfOffset = GDALGetRasterOffset( hBand, &bSuccess );
456 // TODO: Should we turn on checking of bSuccess?
457 // Alternatively, delete these checks and put a comment as to
458 // why checking bSuccess does not matter.
459 #if 0
460 if (bSuccess == FALSE)
461 {
462 CPLError( CE_Debug, CPLE_AppDefined,
463 "Unable to get raster offset." );
464 }
465 #endif
466 double dfScale = GDALGetRasterScale( hBand, &bSuccess );
467 #if 0
468 if (bSuccess == FALSE)
469 {
470 CPLError( CE_Debug, CPLE_AppDefined,
471 "Unable to get raster scale." );
472 }
473 #endif
474 if( dfOffset != 0.0 || dfScale != 1.0 )
475 {
476 adfPixel[0] = adfPixel[0] * dfScale + dfOffset;
477 adfPixel[1] = adfPixel[1] * dfScale + dfOffset;
478
479 if( GDALDataTypeIsComplex( GDALGetRasterDataType( hBand ) ) )
480 osValue.Printf( "%.15g+%.15gi", adfPixel[0], adfPixel[1] );
481 else
482 osValue.Printf( "%.15g", adfPixel[0] );
483
484 if( bAsXML )
485 {
486 osXML += "<DescaledValue>";
487 osXML += osValue;
488 osXML += "</DescaledValue>";
489 }
490 else if( !bQuiet )
491 printf( " Descaled Value: %s\n", osValue.c_str() );
492 }
493 }
494
495 if( bAsXML )
496 osXML += "</BandReport>";
497 }
498
499 osXML += "</Report>";
500
501 if( (pszLocX != nullptr && pszLocY != nullptr) ||
502 (fscanf(stdin, "%lf %lf", &dfGeoX, &dfGeoY) != 2) )
503 {
504 inputAvailable = 0;
505 }
506 }
507
508 /* -------------------------------------------------------------------- */
509 /* Finalize xml report and print. */
510 /* -------------------------------------------------------------------- */
511 if( bAsXML )
512 {
513 CPLXMLNode *psRoot = CPLParseXMLString( osXML );
514 char *pszFormattedXML = CPLSerializeXMLTree( psRoot );
515 CPLDestroyXMLNode( psRoot );
516
517 printf( "%s", pszFormattedXML );
518 CPLFree( pszFormattedXML );
519 }
520
521 /* -------------------------------------------------------------------- */
522 /* Cleanup */
523 /* -------------------------------------------------------------------- */
524 if (hCT) {
525 OSRDestroySpatialReference( hSrcSRS );
526 OCTDestroyCoordinateTransformation( hCT );
527 }
528
529 GDALClose(hSrcDS);
530
531 GDALDumpOpenDatasets( stderr );
532 GDALDestroyDriverManager();
533 CPLFree(pszSourceSRS);
534 CSLDestroy(papszOpenOptions);
535
536 CSLDestroy( argv );
537
538 return 0;
539 }
540 MAIN_END
541