1#!/usr/bin/env python3
2# *****************************************************************************
3# $Id: gdalinfo.py e4fe7cc06270e5f38dfe78e6785a6bcca4e39e29 2021-04-01 21:02:04 +0300 Idan Miara $
4#
5# Project:  GDAL Utilities
6# Purpose:  Python port of Commandline application to list info about a file.
7# Author:   Even Rouault, <even dot rouault at spatialys.com>
8#
9# Port from gdalinfo.c whose author is Frank Warmerdam
10#
11# *****************************************************************************
12# Copyright (c) 2010-2011, Even Rouault <even dot rouault at spatialys.com>
13# Copyright (c) 1998, Frank Warmerdam
14#
15# Permission is hereby granted, free of charge, to any person obtaining a
16# copy of this software and associated documentation files (the "Software"),
17# to deal in the Software without restriction, including without limitation
18# the rights to use, copy, modify, merge, publish, distribute, sublicense,
19# and/or sell copies of the Software, and to permit persons to whom the
20# Software is furnished to do so, subject to the following conditions:
21#
22# The above copyright notice and this permission notice shall be included
23# in all copies or substantial portions of the Software.
24#
25# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
26# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
27# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
28# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
29# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
30# FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
31# DEALINGS IN THE SOFTWARE.
32# ***************************************************************************/
33
34import sys
35
36from osgeo import gdal
37from osgeo import osr
38
39# **********************************************************************
40#                               Usage()
41# **********************************************************************
42
43
44def Usage():
45    print("Usage: gdalinfo [--help-general] [-mm] [-stats] [-hist] [-nogcp] [-nomd]\n" +
46          "                [-norat] [-noct] [-nofl] [-checksum] [-mdd domain]* datasetname")
47    return 1
48
49
50def EQUAL(a, b):
51    return a.lower() == b.lower()
52
53
54def main(argv=None):
55    version_num = int(gdal.VersionInfo('VERSION_NUM'))
56    if version_num < 1800:  # because of GetGeoTransform(can_return_null)
57        print('ERROR: Python bindings of GDAL 1.8.0 or later required')
58        return 1
59
60    bComputeMinMax = False
61    bShowGCPs = True
62    bShowMetadata = True
63    bShowRAT = True
64    bStats = False
65    bApproxStats = True
66    bShowColorTable = True
67    bComputeChecksum = False
68    bReportHistograms = False
69    pszFilename = None
70    papszExtraMDDomains = []
71    pszProjection = None
72    hTransform = None
73    bShowFileList = True
74
75    # Must process GDAL_SKIP before GDALAllRegister(), but we can't call
76    # GDALGeneralCmdLineProcessor before it needs the drivers to be registered
77    # for the --format or --formats options
78    # for( i = 1; i < argc; i++ )
79    # {
80    #    if EQUAL(argv[i],"--config") and i + 2 < argc and EQUAL(argv[i + 1], "GDAL_SKIP"):
81    #    {
82    #        CPLSetConfigOption( argv[i+1], argv[i+2] );
83    #
84    #        i += 2;
85    #    }
86    # }
87    #
88    # GDALAllRegister();
89
90    if argv is None:
91        argv = sys.argv
92
93    argv = gdal.GeneralCmdLineProcessor(argv)
94
95    if argv is None:
96        return 1
97
98    nArgc = len(argv)
99# --------------------------------------------------------------------
100#      Parse arguments.
101# --------------------------------------------------------------------
102    i = 1
103    while i < nArgc:
104
105        if EQUAL(argv[i], "--utility_version"):
106            print("%s is running against GDAL %s" %
107                  (argv[0], gdal.VersionInfo("RELEASE_NAME")))
108            return 0
109        elif EQUAL(argv[i], "-mm"):
110            bComputeMinMax = True
111        elif EQUAL(argv[i], "-hist"):
112            bReportHistograms = True
113        elif EQUAL(argv[i], "-stats"):
114            bStats = True
115            bApproxStats = False
116        elif EQUAL(argv[i], "-approx_stats"):
117            bStats = True
118            bApproxStats = True
119        elif EQUAL(argv[i], "-checksum"):
120            bComputeChecksum = True
121        elif EQUAL(argv[i], "-nogcp"):
122            bShowGCPs = False
123        elif EQUAL(argv[i], "-nomd"):
124            bShowMetadata = False
125        elif EQUAL(argv[i], "-norat"):
126            bShowRAT = False
127        elif EQUAL(argv[i], "-noct"):
128            bShowColorTable = False
129        elif EQUAL(argv[i], "-mdd") and i < nArgc - 1:
130            i = i + 1
131            papszExtraMDDomains.append(argv[i])
132        elif EQUAL(argv[i], "-nofl"):
133            bShowFileList = False
134        elif argv[i][0] == '-':
135            return Usage()
136        elif pszFilename is None:
137            pszFilename = argv[i]
138        else:
139            return Usage()
140
141        i = i + 1
142
143    if pszFilename is None:
144        return Usage()
145
146# --------------------------------------------------------------------
147#      Open dataset.
148# --------------------------------------------------------------------
149    hDataset = gdal.Open(pszFilename, gdal.GA_ReadOnly)
150
151    if hDataset is None:
152
153        print("gdalinfo failed - unable to open '%s'." % pszFilename)
154
155        return 1
156
157# --------------------------------------------------------------------
158#      Report general info.
159# --------------------------------------------------------------------
160    hDriver = hDataset.GetDriver()
161    print("Driver: %s/%s" % (
162        hDriver.ShortName,
163        hDriver.LongName))
164
165    papszFileList = hDataset.GetFileList()
166    if papszFileList is None or not papszFileList:
167        print("Files: none associated")
168    else:
169        print("Files: %s" % papszFileList[0])
170        if bShowFileList:
171            for i in range(1, len(papszFileList)):
172                print("       %s" % papszFileList[i])
173
174    print("Size is %d, %d" % (hDataset.RasterXSize, hDataset.RasterYSize))
175
176# --------------------------------------------------------------------
177#      Report projection.
178# --------------------------------------------------------------------
179    pszProjection = hDataset.GetProjectionRef()
180    if pszProjection is not None:
181
182        hSRS = osr.SpatialReference()
183        if hSRS.ImportFromWkt(pszProjection) == gdal.CE_None:
184            pszPrettyWkt = hSRS.ExportToPrettyWkt(False)
185
186            print("Coordinate System is:\n%s" % pszPrettyWkt)
187        else:
188            print("Coordinate System is `%s'" % pszProjection)
189
190# --------------------------------------------------------------------
191#      Report Geotransform.
192# --------------------------------------------------------------------
193    adfGeoTransform = hDataset.GetGeoTransform(can_return_null=True)
194    if adfGeoTransform is not None:
195
196        if adfGeoTransform[2] == 0.0 and adfGeoTransform[4] == 0.0:
197            print("Origin = (%.15f,%.15f)" % (
198                adfGeoTransform[0], adfGeoTransform[3]))
199
200            print("Pixel Size = (%.15f,%.15f)" % (
201                adfGeoTransform[1], adfGeoTransform[5]))
202
203        else:
204            print("GeoTransform =\n"
205                  "  %.16g, %.16g, %.16g\n"
206                  "  %.16g, %.16g, %.16g" % (
207                      adfGeoTransform[0],
208                      adfGeoTransform[1],
209                      adfGeoTransform[2],
210                      adfGeoTransform[3],
211                      adfGeoTransform[4],
212                      adfGeoTransform[5]))
213
214# --------------------------------------------------------------------
215#      Report GCPs.
216# --------------------------------------------------------------------
217    if bShowGCPs and hDataset.GetGCPCount() > 0:
218
219        pszProjection = hDataset.GetGCPProjection()
220        if pszProjection is not None:
221
222            hSRS = osr.SpatialReference()
223            if hSRS.ImportFromWkt(pszProjection) == gdal.CE_None:
224                pszPrettyWkt = hSRS.ExportToPrettyWkt(False)
225                print("GCP Projection = \n%s" % pszPrettyWkt)
226
227            else:
228                print("GCP Projection = %s" %
229                      pszProjection)
230
231        gcps = hDataset.GetGCPs()
232        i = 0
233        for gcp in gcps:
234
235            print("GCP[%3d]: Id=%s, Info=%s\n"
236                  "          (%.15g,%.15g) -> (%.15g,%.15g,%.15g)" % (
237                      i, gcp.Id, gcp.Info,
238                      gcp.GCPPixel, gcp.GCPLine,
239                      gcp.GCPX, gcp.GCPY, gcp.GCPZ))
240            i = i + 1
241
242# --------------------------------------------------------------------
243#      Report metadata.
244# --------------------------------------------------------------------
245    if bShowMetadata:
246        papszMetadata = hDataset.GetMetadata_List()
247    else:
248        papszMetadata = None
249    if bShowMetadata and papszMetadata:
250        print("Metadata:")
251        for metadata in papszMetadata:
252            print("  %s" % metadata)
253
254    if bShowMetadata:
255        for extra_domain in papszExtraMDDomains:
256            papszMetadata = hDataset.GetMetadata_List(extra_domain)
257            if papszMetadata:
258                print("Metadata (%s):" % extra_domain)
259                for metadata in papszMetadata:
260                    print("  %s" % metadata)
261
262# --------------------------------------------------------------------
263#      Report "IMAGE_STRUCTURE" metadata.
264# --------------------------------------------------------------------
265    if bShowMetadata:
266        papszMetadata = hDataset.GetMetadata_List("IMAGE_STRUCTURE")
267    else:
268        papszMetadata = None
269    if bShowMetadata and papszMetadata:
270        print("Image Structure Metadata:")
271        for metadata in papszMetadata:
272            print("  %s" % metadata)
273
274# --------------------------------------------------------------------
275#      Report subdatasets.
276# --------------------------------------------------------------------
277    papszMetadata = hDataset.GetMetadata_List("SUBDATASETS")
278    if papszMetadata:
279        print("Subdatasets:")
280        for metadata in papszMetadata:
281            print("  %s" % metadata)
282
283# --------------------------------------------------------------------
284#      Report geolocation.
285# --------------------------------------------------------------------
286    if bShowMetadata:
287        papszMetadata = hDataset.GetMetadata_List("GEOLOCATION")
288    else:
289        papszMetadata = None
290    if bShowMetadata and papszMetadata:
291        print("Geolocation:")
292        for metadata in papszMetadata:
293            print("  %s" % metadata)
294
295# --------------------------------------------------------------------
296#      Report RPCs
297# --------------------------------------------------------------------
298    if bShowMetadata:
299        papszMetadata = hDataset.GetMetadata_List("RPC")
300    else:
301        papszMetadata = None
302    if bShowMetadata and papszMetadata:
303        print("RPC Metadata:")
304        for metadata in papszMetadata:
305            print("  %s" % metadata)
306
307# --------------------------------------------------------------------
308#      Setup projected to lat/long transform if appropriate.
309# --------------------------------------------------------------------
310    if pszProjection:
311        hProj = osr.SpatialReference(pszProjection)
312        if hProj is not None:
313            hLatLong = hProj.CloneGeogCS()
314
315        if hLatLong is not None:
316            gdal.PushErrorHandler('CPLQuietErrorHandler')
317            hTransform = osr.CoordinateTransformation(hProj, hLatLong)
318            gdal.PopErrorHandler()
319            if gdal.GetLastErrorMsg().find('Unable to load PROJ.4 library') != -1:
320                hTransform = None
321
322# --------------------------------------------------------------------
323#      Report corners.
324# --------------------------------------------------------------------
325    print("Corner Coordinates:")
326    GDALInfoReportCorner(hDataset, hTransform, "Upper Left",
327                         0.0, 0.0)
328    GDALInfoReportCorner(hDataset, hTransform, "Lower Left",
329                         0.0, hDataset.RasterYSize)
330    GDALInfoReportCorner(hDataset, hTransform, "Upper Right",
331                         hDataset.RasterXSize, 0.0)
332    GDALInfoReportCorner(hDataset, hTransform, "Lower Right",
333                         hDataset.RasterXSize,
334                         hDataset.RasterYSize)
335    GDALInfoReportCorner(hDataset, hTransform, "Center",
336                         hDataset.RasterXSize / 2.0,
337                         hDataset.RasterYSize / 2.0)
338
339# ====================================================================
340#      Loop over bands.
341# ====================================================================
342    for iBand in range(hDataset.RasterCount):
343
344        hBand = hDataset.GetRasterBand(iBand + 1)
345
346        # if( bSample )
347        # {
348        #    float afSample[10000];
349        #    int   nCount;
350        #
351        #    nCount = GDALGetRandomRasterSample( hBand, 10000, afSample );
352        #    print( "Got %d samples.\n", nCount );
353        # }
354
355        (nBlockXSize, nBlockYSize) = hBand.GetBlockSize()
356        print("Band %d Block=%dx%d Type=%s, ColorInterp=%s" % (iBand + 1,
357                                                               nBlockXSize, nBlockYSize,
358                                                               gdal.GetDataTypeName(hBand.DataType),
359                                                               gdal.GetColorInterpretationName(
360                                                                   hBand.GetRasterColorInterpretation())))
361
362        if hBand.GetDescription() is not None \
363                and hBand.GetDescription():
364            print("  Description = %s" % hBand.GetDescription())
365
366        dfMin = hBand.GetMinimum()
367        dfMax = hBand.GetMaximum()
368        if dfMin is not None or dfMax is not None or bComputeMinMax:
369
370            line = "  "
371            if dfMin is not None:
372                line = line + ("Min=%.3f " % dfMin)
373            if dfMax is not None:
374                line = line + ("Max=%.3f " % dfMax)
375
376            if bComputeMinMax:
377                gdal.ErrorReset()
378                adfCMinMax = hBand.ComputeRasterMinMax(False)
379                if gdal.GetLastErrorType() == gdal.CE_None:
380                    line = line + ("  Computed Min/Max=%.3f,%.3f" % (
381                        adfCMinMax[0], adfCMinMax[1]))
382
383            print(line)
384
385        stats = hBand.GetStatistics(bApproxStats, bStats)
386        # Dirty hack to recognize if stats are valid. If invalid, the returned
387        # stddev is negative
388        if stats[3] >= 0.0:
389            print("  Minimum=%.3f, Maximum=%.3f, Mean=%.3f, StdDev=%.3f" % (
390                stats[0], stats[1], stats[2], stats[3]))
391
392        if bReportHistograms:
393
394            hist = hBand.GetDefaultHistogram(force=True, callback=gdal.TermProgress_nocb)
395            if hist is not None:
396                dfMin = hist[0]
397                dfMax = hist[1]
398                nBucketCount = hist[2]
399                panHistogram = hist[3]
400
401                print("  %d buckets from %g to %g:" % (
402                    nBucketCount, dfMin, dfMax))
403                line = '  '
404                for bucket in panHistogram:
405                    line = line + ("%d " % bucket)
406
407                print(line)
408
409        if bComputeChecksum:
410            print("  Checksum=%d" % hBand.Checksum())
411
412        dfNoData = hBand.GetNoDataValue()
413        if dfNoData is not None:
414            if dfNoData != dfNoData:
415                print("  NoData Value=nan")
416            else:
417                print("  NoData Value=%.18g" % dfNoData)
418
419        if hBand.GetOverviewCount() > 0:
420
421            line = "  Overviews: "
422            for iOverview in range(hBand.GetOverviewCount()):
423
424                if iOverview != 0:
425                    line = line + ", "
426
427                hOverview = hBand.GetOverview(iOverview)
428                if hOverview is not None:
429
430                    line = line + ("%dx%d" % (hOverview.XSize, hOverview.YSize))
431
432                    pszResampling = \
433                        hOverview.GetMetadataItem("RESAMPLING", "")
434
435                    if pszResampling is not None \
436                       and len(pszResampling) >= 12 \
437                       and EQUAL(pszResampling[0:12], "AVERAGE_BIT2"):
438                        line = line + "*"
439
440                else:
441                    line = line + "(null)"
442
443            print(line)
444
445            if bComputeChecksum:
446
447                line = "  Overviews checksum: "
448                for iOverview in range(hBand.GetOverviewCount()):
449
450                    if iOverview != 0:
451                        line = line + ", "
452
453                    hOverview = hBand.GetOverview(iOverview)
454                    if hOverview is not None:
455                        line = line + ("%d" % hOverview.Checksum())
456                    else:
457                        line = line + "(null)"
458                print(line)
459
460        if hBand.HasArbitraryOverviews():
461            print("  Overviews: arbitrary")
462
463        nMaskFlags = hBand.GetMaskFlags()
464        if (nMaskFlags & (gdal.GMF_NODATA | gdal.GMF_ALL_VALID)) == 0:
465
466            hMaskBand = hBand.GetMaskBand()
467
468            line = "  Mask Flags: "
469            if (nMaskFlags & gdal.GMF_PER_DATASET) != 0:
470                line = line + "PER_DATASET "
471            if (nMaskFlags & gdal.GMF_ALPHA) != 0:
472                line = line + "ALPHA "
473            if (nMaskFlags & gdal.GMF_NODATA) != 0:
474                line = line + "NODATA "
475            if (nMaskFlags & gdal.GMF_ALL_VALID) != 0:
476                line = line + "ALL_VALID "
477            print(line)
478
479            if hMaskBand is not None and \
480                    hMaskBand.GetOverviewCount() > 0:
481
482                line = "  Overviews of mask band: "
483                for iOverview in range(hMaskBand.GetOverviewCount()):
484
485                    if iOverview != 0:
486                        line = line + ", "
487
488                    hOverview = hMaskBand.GetOverview(iOverview)
489                    if hOverview is not None:
490                        line = line + ("%d" % hOverview.Checksum())
491                    else:
492                        line = line + "(null)"
493
494        if hBand.GetUnitType():
495            print("  Unit Type: %s" % hBand.GetUnitType())
496
497        papszCategories = hBand.GetRasterCategoryNames()
498        if papszCategories is not None:
499
500            print("  Categories:")
501            i = 0
502            for category in papszCategories:
503                print("    %3d: %s" % (i, category))
504                i = i + 1
505
506        scale = hBand.GetScale()
507        if not scale:
508            scale = 1.0
509        offset = hBand.GetOffset()
510        if not offset:
511            offset = 0.0
512        if scale != 1.0 or offset != 0.0:
513            print("  Offset: %.15g,   Scale:%.15g" %
514                  (offset, scale))
515
516        if bShowMetadata:
517            papszMetadata = hBand.GetMetadata_List()
518        else:
519            papszMetadata = None
520        if bShowMetadata and papszMetadata:
521            print("  Metadata:")
522            for metadata in papszMetadata:
523                print("    %s" % metadata)
524
525        if bShowMetadata:
526            papszMetadata = hBand.GetMetadata_List("IMAGE_STRUCTURE")
527        else:
528            papszMetadata = None
529        if bShowMetadata and papszMetadata:
530            print("  Image Structure Metadata:")
531            for metadata in papszMetadata:
532                print("    %s" % metadata)
533
534        hTable = hBand.GetRasterColorTable()
535        if hBand.GetRasterColorInterpretation() == gdal.GCI_PaletteIndex  \
536                and hTable is not None:
537
538            print("  Color Table (%s with %d entries)" % (
539                gdal.GetPaletteInterpretationName(
540                    hTable.GetPaletteInterpretation()),
541                hTable.GetCount()))
542
543            if bShowColorTable:
544
545                for i in range(hTable.GetCount()):
546                    sEntry = hTable.GetColorEntry(i)
547                    print("  %3d: %d,%d,%d,%d" % (
548                        i,
549                        sEntry[0],
550                        sEntry[1],
551                        sEntry[2],
552                        sEntry[3]))
553
554        if bShowRAT:
555            pass
556            # hRAT = hBand.GetDefaultRAT()
557
558            # GDALRATDumpReadable( hRAT, None );
559
560    return 0
561
562# **********************************************************************
563#                        GDALInfoReportCorner()
564# **********************************************************************
565
566
567def GDALInfoReportCorner(hDataset, hTransform, corner_name, x, y):
568
569    line = "%-11s " % corner_name
570
571# --------------------------------------------------------------------
572#      Transform the point into georeferenced coordinates.
573# --------------------------------------------------------------------
574    adfGeoTransform = hDataset.GetGeoTransform(can_return_null=True)
575    if adfGeoTransform is not None:
576        dfGeoX = adfGeoTransform[0] + adfGeoTransform[1] * x \
577            + adfGeoTransform[2] * y
578        dfGeoY = adfGeoTransform[3] + adfGeoTransform[4] * x \
579            + adfGeoTransform[5] * y
580
581    else:
582        line = line + ("(%7.1f,%7.1f)" % (x, y))
583        print(line)
584        return False
585
586# --------------------------------------------------------------------
587#      Report the georeferenced coordinates.
588# --------------------------------------------------------------------
589    if abs(dfGeoX) < 181 and abs(dfGeoY) < 91:
590        line = line + ("(%12.7f,%12.7f) " % (dfGeoX, dfGeoY))
591
592    else:
593        line = line + ("(%12.3f,%12.3f) " % (dfGeoX, dfGeoY))
594
595# --------------------------------------------------------------------
596#      Transform to latlong and report.
597# --------------------------------------------------------------------
598    if hTransform is not None:
599        pnt = hTransform.TransformPoint(dfGeoX, dfGeoY, 0)
600        if pnt is not None:
601            line = line + ("(%s," % gdal.DecToDMS(pnt[0], "Long", 2))
602            line = line + ("%s)" % gdal.DecToDMS(pnt[1], "Lat", 2))
603
604    print(line)
605
606    return True
607
608
609if __name__ == '__main__':
610    sys.exit(main(sys.argv))
611