1 : /******************************************************************************
2 : * $Id: gdallocationinfo.cpp 23726 2012-01-07 22:57:12Z rouault $
3 : *
4 : * Project: GDAL
5 : * Purpose: Commandline raster query tool.
6 : * Author: Frank Warmerdam <warmerdam@pobox.com>
7 : *
8 : ******************************************************************************
9 : * Copyright (c) 2010, Frank Warmerdam <warmerdam@pobox.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 "gdal.h"
31 : #include "cpl_string.h"
32 : #include "ogr_spatialref.h"
33 : #include "cpl_minixml.h"
34 : #include <vector>
35 :
36 : CPL_CVSID("$Id: gdallocationinfo.cpp 23726 2012-01-07 22:57:12Z rouault $");
37 :
38 : /******************************************************************************/
39 : /*! \page gdallocationinfo gdallocationinfo
40 :
41 : raster query tool
42 :
43 : \section gdallocationinfo_synopsis SYNOPSIS
44 :
45 : \htmlonly
46 : Usage:
47 : \endhtmlonly
48 :
49 : \verbatim
50 : Usage: gdallocationinfo [--help-general] [-xml] [-lifonly] [-valonly]
51 : [-b band]* [-overview overview_level]
52 : [-l_srs srs_def] [-geoloc] [-wgs84]
53 : srcfile [x y]
54 : \endverbatim
55 :
56 : \section gdallocationinfo_description DESCRIPTION
57 :
58 : <p>
59 : The gdallocationinfo utility provide a mechanism to query information about
60 : a pixel given it's location in one of a variety of coordinate systems. Several
61 : reporting options are provided.
62 :
63 : <p>
64 : <dl>
65 : <dt> <b>-xml</b>:</dt>
66 : <dd> The output report will be XML formatted for convenient post processing.</dd>
67 :
68 : <dt> <b>-lifonly</b>:</dt>
69 : <dd>The only output is filenames production from the LocationInfo request
70 : against the database (ie. for identifying impacted file from VRT).</dd>
71 :
72 : <dt> <b>-valonly</b>:</dt>
73 : <dd>The only output is the pixel values of the selected pixel on each of
74 : the selected bands.</dd>
75 :
76 : <dt> <b>-b</b> <em>band</em>:</dt>
77 : <dd>Selects a band to query. Multiple bands can be listed. By default all
78 : bands are queried.</dd>
79 :
80 : <dt> <b>-overview</b> <em>overview_level</em>:</dt>
81 : <dd>Query the (overview_level)th overview (overview_level=1 is the 1st overview),
82 : instead of the base band. Note that the x,y location (if the coordinate system is
83 : pixel/line) must still be given with respect to the base band.</dd>
84 :
85 : <dt> <b>-l_srs</b> <em>srs def</em>:</dt>
86 : <dd> The coordinate system of the input x, y location.</dd>
87 :
88 : <dt> <b>-geoloc</b>:</dt>
89 : <dd> Indicates input x,y points are in the georeferencing system of the image.</dd>
90 :
91 : <dt> <b>-wgs84</b>:</dt>
92 : <dd> Indicates input x,y points are WGS84 long, lat.</dd>
93 :
94 : <dt> <em>srcfile</em>:</dt><dd> The source GDAL raster datasource name.</dd>
95 :
96 : <dt> <em>x</em>:</dt><dd> X location of target pixel. By default the
97 : coordinate system is pixel/line unless -l_srs, -wgs84 or -geoloc supplied. </dd>
98 :
99 : <dt> <em>y</em>:</dt><dd> Y location of target pixel. By default the
100 : coordinate system is pixel/line unless -l_srs, -wgs84 or -geoloc supplied. </dd>
101 :
102 : </dl>
103 :
104 : This utility is intended to provide a variety of information about a
105 : pixel. Currently it reports three things:
106 :
107 : <ul>
108 : <li> The location of the pixel in pixel/line space.
109 : <li> The result of a LocationInfo metadata query against the datasource -
110 : currently this is only implemented for VRT files which will report the
111 : file(s) used to satisfy requests for that pixel.
112 : <li> The raster pixel value of that pixel for all or a subset of the bands.
113 : <li> The unscaled pixel value if a Scale and/or Offset apply to the band.
114 : </ul>
115 :
116 : The pixel selected is requested by x/y coordinate on the commandline, or read
117 : from stdin. More than one coordinate pair can be supplied when reading
118 : coordinatesis from stdin. By default pixel/line coordinates are expected.
119 : However with use of the -geoloc, -wgs84, or -l_srs switches it is possible
120 : to specify the location in other coordinate systems.
121 :
122 : The default report is in a human readable text format. It is possible to
123 : instead request xml output with the -xml switch.
124 :
125 : For scripting purposes, the -valonly and -lifonly switches are provided to
126 : restrict output to the actual pixel values, or the LocationInfo files
127 : identified for the pixel.
128 :
129 : It is anticipated that additional reporting capabilities will be added to
130 : gdallocationinfo in the future.
131 :
132 : <p>
133 : \section gdallocationinfo_example EXAMPLE
134 :
135 : Simple example reporting on pixel (256,256) on the file utm.tif.
136 :
137 : \verbatim
138 : $ gdallocationinfo utm.tif 256 256
139 : Report:
140 : Location: (256P,256L)
141 : Band 1:
142 : Value: 115
143 : \endverbatim
144 :
145 : Query a VRT file providing the location in WGS84, and getting the result in xml.
146 :
147 : \verbatim
148 : $ gdallocationinfo -xml -wgs84 utm.vrt -117.5 33.75
149 : <Report pixel="217" line="282">
150 : <BandReport band="1">
151 : <LocationInfo>
152 : <File>utm.tif</File>
153 : </LocationInfo>
154 : <Value>16</Value>
155 : </BandReport>
156 : </Report>
157 : \endverbatim
158 :
159 : \if man
160 : \section gdallocationinfo_author AUTHORS
161 : Frank Warmerdam <warmerdam@pobox.com>
162 : \endif
163 : */
164 :
165 : /************************************************************************/
166 : /* Usage() */
167 : /************************************************************************/
168 :
169 0 : static void Usage()
170 :
171 : {
172 : printf( "Usage: gdallocationinfo [--help-general] [-xml] [-lifonly] [-valonly]\n"
173 : " [-b band]* [-overview overview_level]\n"
174 : " [-l_srs srs_def] [-geoloc] [-wgs84]\n"
175 : " srcfile x y\n"
176 0 : "\n" );
177 0 : exit( 1 );
178 : }
179 :
180 : /************************************************************************/
181 : /* SanitizeSRS */
182 : /************************************************************************/
183 :
184 0 : char *SanitizeSRS( const char *pszUserInput )
185 :
186 : {
187 : OGRSpatialReferenceH hSRS;
188 0 : char *pszResult = NULL;
189 :
190 0 : CPLErrorReset();
191 :
192 0 : hSRS = OSRNewSpatialReference( NULL );
193 0 : if( OSRSetFromUserInput( hSRS, pszUserInput ) == OGRERR_NONE )
194 0 : OSRExportToWkt( hSRS, &pszResult );
195 : else
196 : {
197 : CPLError( CE_Failure, CPLE_AppDefined,
198 : "Translating source or target SRS failed:\n%s",
199 0 : pszUserInput );
200 0 : exit( 1 );
201 : }
202 :
203 0 : OSRDestroySpatialReference( hSRS );
204 :
205 0 : return pszResult;
206 : }
207 :
208 : /************************************************************************/
209 : /* main() */
210 : /************************************************************************/
211 :
212 7 : int main( int argc, char ** argv )
213 :
214 : {
215 7 : const char *pszLocX = NULL, *pszLocY = NULL;
216 7 : const char *pszSrcFilename = NULL;
217 7 : char *pszSourceSRS = NULL;
218 7 : std::vector<int> anBandList;
219 7 : bool bAsXML = false, bLIFOnly = false;
220 7 : bool bQuiet = false, bValOnly = false;
221 7 : int nOverview = -1;
222 :
223 7 : GDALAllRegister();
224 7 : argc = GDALGeneralCmdLineProcessor( argc, &argv, 0 );
225 7 : if( argc < 1 )
226 0 : exit( -argc );
227 :
228 : /* -------------------------------------------------------------------- */
229 : /* Parse arguments. */
230 : /* -------------------------------------------------------------------- */
231 : int i;
232 :
233 31 : for( i = 1; i < argc; i++ )
234 : {
235 25 : if( EQUAL(argv[i], "--utility_version") )
236 : {
237 : printf("%s was compiled against GDAL %s and is running against GDAL %s\n",
238 1 : argv[0], GDAL_RELEASE_NAME, GDALVersionInfo("RELEASE_NAME"));
239 1 : return 0;
240 : }
241 25 : else if( EQUAL(argv[i],"-b") && i < argc-1 )
242 : {
243 1 : anBandList.push_back( atoi(argv[++i]) );
244 : }
245 24 : else if( EQUAL(argv[i],"-overview") && i < argc-1 )
246 : {
247 1 : nOverview = atoi(argv[++i]) - 1;
248 : }
249 22 : else if( EQUAL(argv[i],"-l_srs") && i < argc-1 )
250 : {
251 0 : CPLFree(pszSourceSRS);
252 0 : pszSourceSRS = SanitizeSRS(argv[++i]);
253 : }
254 22 : else if( EQUAL(argv[i],"-geoloc") )
255 : {
256 1 : CPLFree(pszSourceSRS);
257 1 : pszSourceSRS = CPLStrdup("-geoloc");
258 : }
259 21 : else if( EQUAL(argv[i],"-wgs84") )
260 : {
261 0 : CPLFree(pszSourceSRS);
262 0 : pszSourceSRS = SanitizeSRS("WGS84");
263 : }
264 21 : else if( EQUAL(argv[i],"-xml") )
265 : {
266 1 : bAsXML = true;
267 : }
268 20 : else if( EQUAL(argv[i],"-lifonly") )
269 : {
270 1 : bLIFOnly = true;
271 1 : bQuiet = true;
272 : }
273 19 : else if( EQUAL(argv[i],"-valonly") )
274 : {
275 1 : bValOnly = true;
276 1 : bQuiet = true;
277 : }
278 18 : else if( argv[i][0] == '-' && !isdigit(argv[i][1]) )
279 0 : Usage();
280 :
281 18 : else if( pszSrcFilename == NULL )
282 6 : pszSrcFilename = argv[i];
283 :
284 12 : else if( pszLocX == NULL )
285 6 : pszLocX = argv[i];
286 :
287 6 : else if( pszLocY == NULL )
288 6 : pszLocY = argv[i];
289 :
290 : else
291 0 : Usage();
292 : }
293 :
294 6 : if( pszSrcFilename == NULL || (pszLocX != NULL && pszLocY == NULL) )
295 0 : Usage();
296 :
297 : /* -------------------------------------------------------------------- */
298 : /* Open source file. */
299 : /* -------------------------------------------------------------------- */
300 6 : GDALDatasetH hSrcDS = NULL;
301 :
302 6 : hSrcDS = GDALOpen( pszSrcFilename, GA_ReadOnly );
303 6 : if( hSrcDS == NULL )
304 0 : exit( 1 );
305 :
306 : /* -------------------------------------------------------------------- */
307 : /* Setup coordinate transformation, if required */
308 : /* -------------------------------------------------------------------- */
309 6 : OGRSpatialReferenceH hSrcSRS = NULL, hTrgSRS = NULL;
310 6 : OGRCoordinateTransformationH hCT = NULL;
311 6 : if( pszSourceSRS != NULL && !EQUAL(pszSourceSRS,"-geoloc") )
312 : {
313 :
314 0 : hSrcSRS = OSRNewSpatialReference( pszSourceSRS );
315 0 : hTrgSRS = OSRNewSpatialReference( GDALGetProjectionRef( hSrcDS ) );
316 :
317 0 : hCT = OCTNewCoordinateTransformation( hSrcSRS, hTrgSRS );
318 0 : if( hCT == NULL )
319 0 : exit( 1 );
320 : }
321 :
322 : /* -------------------------------------------------------------------- */
323 : /* If no bands were requested, we will query them all. */
324 : /* -------------------------------------------------------------------- */
325 6 : if( anBandList.size() == 0 )
326 : {
327 10 : for( i = 0; i < GDALGetRasterCount( hSrcDS ); i++ )
328 5 : anBandList.push_back( i+1 );
329 : }
330 :
331 : /* -------------------------------------------------------------------- */
332 : /* Turn the location into a pixel and line location. */
333 : /* -------------------------------------------------------------------- */
334 6 : int inputAvailable = 1;
335 : double dfGeoX;
336 : double dfGeoY;
337 6 : CPLString osXML;
338 :
339 6 : if( pszLocX == NULL && pszLocY == NULL )
340 : {
341 0 : if (fscanf(stdin, "%lf %lf", &dfGeoX, &dfGeoY) != 2)
342 : {
343 0 : inputAvailable = 0;
344 : }
345 : }
346 : else
347 : {
348 6 : dfGeoX = CPLAtof(pszLocX);
349 6 : dfGeoY = CPLAtof(pszLocY);
350 : }
351 :
352 12 : while (inputAvailable)
353 : {
354 : int iPixel, iLine;
355 :
356 6 : if (hCT)
357 : {
358 0 : if( !OCTTransform( hCT, 1, &dfGeoX, &dfGeoY, NULL ) )
359 0 : exit( 1 );
360 : }
361 :
362 6 : if( pszSourceSRS != NULL )
363 : {
364 : double adfGeoTransform[6], adfInvGeoTransform[6];
365 :
366 1 : if( GDALGetGeoTransform( hSrcDS, adfGeoTransform ) != CE_None )
367 0 : exit( 1 );
368 :
369 1 : GDALInvGeoTransform( adfGeoTransform, adfInvGeoTransform );
370 :
371 : iPixel = (int) floor(
372 1 : adfInvGeoTransform[0]
373 1 : + adfInvGeoTransform[1] * dfGeoX
374 1 : + adfInvGeoTransform[2] * dfGeoY );
375 : iLine = (int) floor(
376 1 : adfInvGeoTransform[3]
377 1 : + adfInvGeoTransform[4] * dfGeoX
378 1 : + adfInvGeoTransform[5] * dfGeoY );
379 : }
380 : else
381 : {
382 5 : iPixel = (int) floor(dfGeoX);
383 5 : iLine = (int) floor(dfGeoY);
384 : }
385 :
386 : /* -------------------------------------------------------------------- */
387 : /* Prepare report. */
388 : /* -------------------------------------------------------------------- */
389 6 : CPLString osLine;
390 :
391 6 : if( bAsXML )
392 : {
393 : osLine.Printf( "<Report pixel=\"%d\" line=\"%d\">",
394 1 : iPixel, iLine );
395 1 : osXML += osLine;
396 : }
397 5 : else if( !bQuiet )
398 : {
399 3 : printf( "Report:\n" );
400 3 : printf( " Location: (%dP,%dL)\n", iPixel, iLine );
401 : }
402 :
403 6 : int bPixelReport = TRUE;
404 :
405 6 : if( iPixel < 0 || iLine < 0
406 : || iPixel >= GDALGetRasterXSize( hSrcDS )
407 : || iLine >= GDALGetRasterYSize( hSrcDS ) )
408 : {
409 0 : if( bAsXML )
410 0 : osXML += "<Alert>Location is off this file! No further details to report.</Alert>";
411 0 : else if( bValOnly )
412 0 : printf("\n");
413 0 : else if( !bQuiet )
414 0 : printf( "\nLocation is off this file! No further details to report.\n");
415 0 : bPixelReport = FALSE;
416 : }
417 :
418 : /* -------------------------------------------------------------------- */
419 : /* Process each band. */
420 : /* -------------------------------------------------------------------- */
421 12 : for( i = 0; bPixelReport && i < (int) anBandList.size(); i++ )
422 : {
423 6 : GDALRasterBandH hBand = GDALGetRasterBand( hSrcDS, anBandList[i] );
424 :
425 6 : int iPixelToQuery = iPixel;
426 6 : int iLineToQuery = iLine;
427 :
428 6 : if (nOverview >= 0 && hBand != NULL)
429 : {
430 1 : GDALRasterBandH hOvrBand = GDALGetOverview(hBand, nOverview);
431 1 : if (hOvrBand != NULL)
432 : {
433 1 : int nOvrXSize = GDALGetRasterBandXSize(hOvrBand);
434 1 : int nOvrYSize = GDALGetRasterBandYSize(hOvrBand);
435 1 : iPixelToQuery = (int)(0.5 + 1.0 * iPixel / GDALGetRasterXSize( hSrcDS ) * nOvrXSize);
436 1 : iLineToQuery = (int)(0.5 + 1.0 * iLine / GDALGetRasterYSize( hSrcDS ) * nOvrYSize);
437 1 : if (iPixelToQuery >= nOvrXSize)
438 0 : iPixelToQuery = nOvrXSize - 1;
439 1 : if (iLineToQuery >= nOvrYSize)
440 0 : iLineToQuery = nOvrYSize - 1;
441 : }
442 : else
443 : {
444 : CPLError(CE_Failure, CPLE_AppDefined,
445 : "Cannot get overview %d of band %d",
446 0 : nOverview + 1, anBandList[i] );
447 : }
448 1 : hBand = hOvrBand;
449 : }
450 :
451 6 : if (hBand == NULL)
452 0 : continue;
453 :
454 6 : if( bAsXML )
455 : {
456 1 : osLine.Printf( "<BandReport band=\"%d\">", anBandList[i] );
457 1 : osXML += osLine;
458 : }
459 5 : else if( !bQuiet )
460 : {
461 3 : printf( " Band %d:\n", anBandList[i] );
462 : }
463 :
464 : /* -------------------------------------------------------------------- */
465 : /* Request location info for this location. It is possible */
466 : /* only the VRT driver actually supports this. */
467 : /* -------------------------------------------------------------------- */
468 6 : CPLString osItem;
469 :
470 6 : osItem.Printf( "Pixel_%d_%d", iPixelToQuery, iLineToQuery );
471 :
472 6 : const char *pszLI = GDALGetMetadataItem( hBand, osItem, "LocationInfo");
473 :
474 6 : if( pszLI != NULL )
475 : {
476 1 : if( bAsXML )
477 0 : osXML += pszLI;
478 1 : else if( !bQuiet )
479 0 : printf( " %s\n", pszLI );
480 1 : else if( bLIFOnly )
481 : {
482 : /* Extract all files, if any. */
483 :
484 1 : CPLXMLNode *psRoot = CPLParseXMLString( pszLI );
485 :
486 1 : if( psRoot != NULL
487 : && psRoot->psChild != NULL
488 : && psRoot->eType == CXT_Element
489 : && EQUAL(psRoot->pszValue,"LocationInfo") )
490 : {
491 : CPLXMLNode *psNode;
492 :
493 2 : for( psNode = psRoot->psChild;
494 : psNode != NULL;
495 : psNode = psNode->psNext )
496 : {
497 1 : if( psNode->eType == CXT_Element
498 : && EQUAL(psNode->pszValue,"File")
499 : && psNode->psChild != NULL )
500 : {
501 : char* pszUnescaped = CPLUnescapeString(
502 1 : psNode->psChild->pszValue, NULL, CPLES_XML);
503 1 : printf( "%s\n", pszUnescaped );
504 1 : CPLFree(pszUnescaped);
505 : }
506 : }
507 : }
508 1 : CPLDestroyXMLNode( psRoot );
509 : }
510 : }
511 :
512 : /* -------------------------------------------------------------------- */
513 : /* Report the pixel value of this band. */
514 : /* -------------------------------------------------------------------- */
515 : double adfPixel[2];
516 :
517 6 : if( GDALRasterIO( hBand, GF_Read, iPixelToQuery, iLineToQuery, 1, 1,
518 : adfPixel, 1, 1, GDT_CFloat64, 0, 0) == CE_None )
519 : {
520 6 : CPLString osValue;
521 :
522 6 : if( GDALDataTypeIsComplex( GDALGetRasterDataType( hBand ) ) )
523 0 : osValue.Printf( "%.15g+%.15gi", adfPixel[0], adfPixel[1] );
524 : else
525 6 : osValue.Printf( "%.15g", adfPixel[0] );
526 :
527 6 : if( bAsXML )
528 : {
529 1 : osXML += "<Value>";
530 1 : osXML += osValue;
531 1 : osXML += "</Value>";
532 : }
533 5 : else if( !bQuiet )
534 3 : printf( " Value: %s\n", osValue.c_str() );
535 2 : else if( bValOnly )
536 1 : printf( "%s\n", osValue.c_str() );
537 :
538 : // Report unscaled if we have scale/offset values.
539 : int bSuccess;
540 :
541 6 : double dfOffset = GDALGetRasterOffset( hBand, &bSuccess );
542 6 : double dfScale = GDALGetRasterScale( hBand, &bSuccess );
543 :
544 6 : if( dfOffset != 0.0 || dfScale != 1.0 )
545 : {
546 0 : adfPixel[0] = adfPixel[0] * dfScale + dfOffset;
547 0 : adfPixel[1] = adfPixel[1] * dfScale + dfOffset;
548 :
549 0 : if( GDALDataTypeIsComplex( GDALGetRasterDataType( hBand ) ) )
550 0 : osValue.Printf( "%.15g+%.15gi", adfPixel[0], adfPixel[1] );
551 : else
552 0 : osValue.Printf( "%.15g", adfPixel[0] );
553 :
554 0 : if( bAsXML )
555 : {
556 0 : osXML += "<DescaledValue>";
557 0 : osXML += osValue;
558 0 : osXML += "</DescaledValue>";
559 : }
560 0 : else if( !bQuiet )
561 0 : printf( " Descaled Value: %s\n", osValue.c_str() );
562 6 : }
563 : }
564 :
565 6 : if( bAsXML )
566 1 : osXML += "</BandReport>";
567 : }
568 :
569 6 : osXML += "</Report>";
570 :
571 6 : if( (pszLocX != NULL && pszLocY != NULL) ||
572 : (fscanf(stdin, "%lf %lf", &dfGeoX, &dfGeoY) != 2) )
573 : {
574 6 : inputAvailable = 0;
575 : }
576 :
577 : }
578 :
579 : /* -------------------------------------------------------------------- */
580 : /* Finalize xml report and print. */
581 : /* -------------------------------------------------------------------- */
582 6 : if( bAsXML )
583 : {
584 : CPLXMLNode *psRoot;
585 : char *pszFormattedXML;
586 :
587 :
588 1 : psRoot = CPLParseXMLString( osXML );
589 1 : pszFormattedXML = CPLSerializeXMLTree( psRoot );
590 1 : CPLDestroyXMLNode( psRoot );
591 :
592 1 : printf( "%s", pszFormattedXML );
593 1 : CPLFree( pszFormattedXML );
594 : }
595 :
596 : /* -------------------------------------------------------------------- */
597 : /* Cleanup */
598 : /* -------------------------------------------------------------------- */
599 6 : if (hCT) {
600 0 : OSRDestroySpatialReference( hSrcSRS );
601 0 : OSRDestroySpatialReference( hTrgSRS );
602 0 : OCTDestroyCoordinateTransformation( hCT );
603 : }
604 :
605 6 : if (hSrcDS)
606 6 : GDALClose(hSrcDS);
607 :
608 6 : GDALDumpOpenDatasets( stderr );
609 6 : GDALDestroyDriverManager();
610 6 : CPLFree(pszSourceSRS);
611 :
612 6 : CSLDestroy( argv );
613 :
614 6 : return 0;
615 : }
|
