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