JezK
Edit File: gdalinfo.py
#!/usr/bin/env python3 # ***************************************************************************** # $Id$ # # Project: GDAL Utilities # Purpose: Python port of Commandline application to list info about a file. # Author: Even Rouault, <even dot rouault at spatialys.com> # # Port from gdalinfo.c whose author is Frank Warmerdam # # ***************************************************************************** # Copyright (c) 2010-2011, Even Rouault <even dot rouault at spatialys.com> # Copyright (c) 1998, Frank Warmerdam # # Permission is hereby granted, free of charge, to any person obtaining a # copy of this software and associated documentation files (the "Software"), # to deal in the Software without restriction, including without limitation # the rights to use, copy, modify, merge, publish, distribute, sublicense, # and/or sell copies of the Software, and to permit persons to whom the # Software is furnished to do so, subject to the following conditions: # # The above copyright notice and this permission notice shall be included # in all copies or substantial portions of the Software. # # THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS # OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, # FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL # THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER # LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING # FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER # DEALINGS IN THE SOFTWARE. # ***************************************************************************/ import sys from osgeo import gdal, osr # ********************************************************************** # Usage() # ********************************************************************** def Usage(): print( "Usage: gdalinfo [--help-general] [-mm] [-stats] [-hist] [-nogcp] [-nomd]\n" + " [-norat] [-noct] [-nofl] [-checksum] [-mdd domain]* datasetname" ) return 2 def EQUAL(a, b): return a.lower() == b.lower() def main(argv=None): version_num = int(gdal.VersionInfo("VERSION_NUM")) if version_num < 1800: # because of GetGeoTransform(can_return_null) print("ERROR: Python bindings of GDAL 1.8.0 or later required") return 1 bComputeMinMax = False bShowGCPs = True bShowMetadata = True bShowRAT = True bStats = False bApproxStats = True bShowColorTable = True bComputeChecksum = False bReportHistograms = False pszFilename = None papszExtraMDDomains = [] pszProjection = None hTransform = None bShowFileList = True # Must process GDAL_SKIP before GDALAllRegister(), but we can't call # GDALGeneralCmdLineProcessor before it needs the drivers to be registered # for the --format or --formats options # for( i = 1; i < argc; i++ ) # { # if EQUAL(argv[i],"--config") and i + 2 < argc and EQUAL(argv[i + 1], "GDAL_SKIP"): # { # CPLSetConfigOption( argv[i+1], argv[i+2] ); # # i += 2; # } # } # # GDALAllRegister(); if argv is None: argv = sys.argv argv = gdal.GeneralCmdLineProcessor(argv) if argv is None: return 1 nArgc = len(argv) # -------------------------------------------------------------------- # Parse arguments. # -------------------------------------------------------------------- i = 1 while i < nArgc: if EQUAL(argv[i], "--utility_version"): print( "%s is running against GDAL %s" % (argv[0], gdal.VersionInfo("RELEASE_NAME")) ) return 0 elif EQUAL(argv[i], "-mm"): bComputeMinMax = True elif EQUAL(argv[i], "-hist"): bReportHistograms = True elif EQUAL(argv[i], "-stats"): bStats = True bApproxStats = False elif EQUAL(argv[i], "-approx_stats"): bStats = True bApproxStats = True elif EQUAL(argv[i], "-checksum"): bComputeChecksum = True elif EQUAL(argv[i], "-nogcp"): bShowGCPs = False elif EQUAL(argv[i], "-nomd"): bShowMetadata = False elif EQUAL(argv[i], "-norat"): bShowRAT = False elif EQUAL(argv[i], "-noct"): bShowColorTable = False elif EQUAL(argv[i], "-mdd") and i < nArgc - 1: i = i + 1 papszExtraMDDomains.append(argv[i]) elif EQUAL(argv[i], "-nofl"): bShowFileList = False elif argv[i][0] == "-": return Usage() elif pszFilename is None: pszFilename = argv[i] else: return Usage() i = i + 1 if pszFilename is None: return Usage() # -------------------------------------------------------------------- # Open dataset. # -------------------------------------------------------------------- hDataset = gdal.Open(pszFilename, gdal.GA_ReadOnly) if hDataset is None: print("gdalinfo failed - unable to open '%s'." % pszFilename) return 1 # -------------------------------------------------------------------- # Report general info. # -------------------------------------------------------------------- hDriver = hDataset.GetDriver() print("Driver: %s/%s" % (hDriver.ShortName, hDriver.LongName)) papszFileList = hDataset.GetFileList() if papszFileList is None or not papszFileList: print("Files: none associated") else: print("Files: %s" % papszFileList[0]) if bShowFileList: for i in range(1, len(papszFileList)): print(" %s" % papszFileList[i]) print("Size is %d, %d" % (hDataset.RasterXSize, hDataset.RasterYSize)) # -------------------------------------------------------------------- # Report projection. # -------------------------------------------------------------------- pszProjection = hDataset.GetProjectionRef() if pszProjection is not None: hSRS = osr.SpatialReference() if hSRS.ImportFromWkt(pszProjection) == gdal.CE_None: pszPrettyWkt = hSRS.ExportToPrettyWkt(False) print("Coordinate System is:\n%s" % pszPrettyWkt) else: print("Coordinate System is `%s'" % pszProjection) # -------------------------------------------------------------------- # Report Geotransform. # -------------------------------------------------------------------- adfGeoTransform = hDataset.GetGeoTransform(can_return_null=True) if adfGeoTransform is not None: if adfGeoTransform[2] == 0.0 and adfGeoTransform[4] == 0.0: print("Origin = (%.15f,%.15f)" % (adfGeoTransform[0], adfGeoTransform[3])) print( "Pixel Size = (%.15f,%.15f)" % (adfGeoTransform[1], adfGeoTransform[5]) ) else: print( "GeoTransform =\n" " %.16g, %.16g, %.16g\n" " %.16g, %.16g, %.16g" % ( adfGeoTransform[0], adfGeoTransform[1], adfGeoTransform[2], adfGeoTransform[3], adfGeoTransform[4], adfGeoTransform[5], ) ) # -------------------------------------------------------------------- # Report GCPs. # -------------------------------------------------------------------- if bShowGCPs and hDataset.GetGCPCount() > 0: pszProjection = hDataset.GetGCPProjection() if pszProjection is not None: hSRS = osr.SpatialReference() if hSRS.ImportFromWkt(pszProjection) == gdal.CE_None: pszPrettyWkt = hSRS.ExportToPrettyWkt(False) print("GCP Projection = \n%s" % pszPrettyWkt) else: print("GCP Projection = %s" % pszProjection) gcps = hDataset.GetGCPs() i = 0 for gcp in gcps: print( "GCP[%3d]: Id=%s, Info=%s\n" " (%.15g,%.15g) -> (%.15g,%.15g,%.15g)" % ( i, gcp.Id, gcp.Info, gcp.GCPPixel, gcp.GCPLine, gcp.GCPX, gcp.GCPY, gcp.GCPZ, ) ) i = i + 1 # -------------------------------------------------------------------- # Report metadata. # -------------------------------------------------------------------- if bShowMetadata: papszMetadata = hDataset.GetMetadata_List() else: papszMetadata = None if bShowMetadata and papszMetadata: print("Metadata:") for metadata in papszMetadata: print(" %s" % metadata) if bShowMetadata: for extra_domain in papszExtraMDDomains: papszMetadata = hDataset.GetMetadata_List(extra_domain) if papszMetadata: print("Metadata (%s):" % extra_domain) for metadata in papszMetadata: print(" %s" % metadata) # -------------------------------------------------------------------- # Report "IMAGE_STRUCTURE" metadata. # -------------------------------------------------------------------- if bShowMetadata: papszMetadata = hDataset.GetMetadata_List("IMAGE_STRUCTURE") else: papszMetadata = None if bShowMetadata and papszMetadata: print("Image Structure Metadata:") for metadata in papszMetadata: print(" %s" % metadata) # -------------------------------------------------------------------- # Report subdatasets. # -------------------------------------------------------------------- papszMetadata = hDataset.GetMetadata_List("SUBDATASETS") if papszMetadata: print("Subdatasets:") for metadata in papszMetadata: print(" %s" % metadata) # -------------------------------------------------------------------- # Report geolocation. # -------------------------------------------------------------------- if bShowMetadata: papszMetadata = hDataset.GetMetadata_List("GEOLOCATION") else: papszMetadata = None if bShowMetadata and papszMetadata: print("Geolocation:") for metadata in papszMetadata: print(" %s" % metadata) # -------------------------------------------------------------------- # Report RPCs # -------------------------------------------------------------------- if bShowMetadata: papszMetadata = hDataset.GetMetadata_List("RPC") else: papszMetadata = None if bShowMetadata and papszMetadata: print("RPC Metadata:") for metadata in papszMetadata: print(" %s" % metadata) # -------------------------------------------------------------------- # Setup projected to lat/long transform if appropriate. # -------------------------------------------------------------------- if pszProjection: hProj = osr.SpatialReference(pszProjection) if hProj is not None: hLatLong = hProj.CloneGeogCS() if hLatLong is not None: gdal.PushErrorHandler("CPLQuietErrorHandler") hTransform = osr.CoordinateTransformation(hProj, hLatLong) gdal.PopErrorHandler() if gdal.GetLastErrorMsg().find("Unable to load PROJ.4 library") != -1: hTransform = None # -------------------------------------------------------------------- # Report corners. # -------------------------------------------------------------------- print("Corner Coordinates:") GDALInfoReportCorner(hDataset, hTransform, "Upper Left", 0.0, 0.0) GDALInfoReportCorner(hDataset, hTransform, "Lower Left", 0.0, hDataset.RasterYSize) GDALInfoReportCorner(hDataset, hTransform, "Upper Right", hDataset.RasterXSize, 0.0) GDALInfoReportCorner( hDataset, hTransform, "Lower Right", hDataset.RasterXSize, hDataset.RasterYSize ) GDALInfoReportCorner( hDataset, hTransform, "Center", hDataset.RasterXSize / 2.0, hDataset.RasterYSize / 2.0, ) # ==================================================================== # Loop over bands. # ==================================================================== for iBand in range(hDataset.RasterCount): hBand = hDataset.GetRasterBand(iBand + 1) # if( bSample ) # { # float afSample[10000]; # int nCount; # # nCount = GDALGetRandomRasterSample( hBand, 10000, afSample ); # print( "Got %d samples.\n", nCount ); # } (nBlockXSize, nBlockYSize) = hBand.GetBlockSize() print( "Band %d Block=%dx%d Type=%s, ColorInterp=%s" % ( iBand + 1, nBlockXSize, nBlockYSize, gdal.GetDataTypeName(hBand.DataType), gdal.GetColorInterpretationName(hBand.GetRasterColorInterpretation()), ) ) if hBand.GetDescription() is not None and hBand.GetDescription(): print(" Description = %s" % hBand.GetDescription()) dfMin = hBand.GetMinimum() dfMax = hBand.GetMaximum() if dfMin is not None or dfMax is not None or bComputeMinMax: line = " " if dfMin is not None: line = line + ("Min=%.3f " % dfMin) if dfMax is not None: line = line + ("Max=%.3f " % dfMax) if bComputeMinMax: gdal.ErrorReset() adfCMinMax = hBand.ComputeRasterMinMax(False) if gdal.GetLastErrorType() == gdal.CE_None: line = line + ( " Computed Min/Max=%.3f,%.3f" % (adfCMinMax[0], adfCMinMax[1]) ) print(line) stats = hBand.GetStatistics(bApproxStats, bStats) # Dirty hack to recognize if stats are valid. If invalid, the returned # stddev is negative if stats[3] >= 0.0: print( " Minimum=%.3f, Maximum=%.3f, Mean=%.3f, StdDev=%.3f" % (stats[0], stats[1], stats[2], stats[3]) ) if bReportHistograms: hist = hBand.GetDefaultHistogram( force=True, callback=gdal.TermProgress_nocb ) if hist is not None: dfMin = hist[0] dfMax = hist[1] nBucketCount = hist[2] panHistogram = hist[3] print(" %d buckets from %g to %g:" % (nBucketCount, dfMin, dfMax)) line = " " for bucket in panHistogram: line = line + ("%d " % bucket) print(line) if bComputeChecksum: print(" Checksum=%d" % hBand.Checksum()) dfNoData = hBand.GetNoDataValue() if dfNoData is not None: if dfNoData != dfNoData: print(" NoData Value=nan") else: print(" NoData Value=%.18g" % dfNoData) if hBand.GetOverviewCount() > 0: line = " Overviews: " for iOverview in range(hBand.GetOverviewCount()): if iOverview != 0: line = line + ", " hOverview = hBand.GetOverview(iOverview) if hOverview is not None: line = line + ("%dx%d" % (hOverview.XSize, hOverview.YSize)) pszResampling = hOverview.GetMetadataItem("RESAMPLING", "") if ( pszResampling is not None and len(pszResampling) >= 12 and EQUAL(pszResampling[0:12], "AVERAGE_BIT2") ): line = line + "*" else: line = line + "(null)" print(line) if bComputeChecksum: line = " Overviews checksum: " for iOverview in range(hBand.GetOverviewCount()): if iOverview != 0: line = line + ", " hOverview = hBand.GetOverview(iOverview) if hOverview is not None: line = line + ("%d" % hOverview.Checksum()) else: line = line + "(null)" print(line) if hBand.HasArbitraryOverviews(): print(" Overviews: arbitrary") nMaskFlags = hBand.GetMaskFlags() if (nMaskFlags & (gdal.GMF_NODATA | gdal.GMF_ALL_VALID)) == 0: hMaskBand = hBand.GetMaskBand() line = " Mask Flags: " if (nMaskFlags & gdal.GMF_PER_DATASET) != 0: line = line + "PER_DATASET " if (nMaskFlags & gdal.GMF_ALPHA) != 0: line = line + "ALPHA " if (nMaskFlags & gdal.GMF_NODATA) != 0: line = line + "NODATA " if (nMaskFlags & gdal.GMF_ALL_VALID) != 0: line = line + "ALL_VALID " print(line) if hMaskBand is not None and hMaskBand.GetOverviewCount() > 0: line = " Overviews of mask band: " for iOverview in range(hMaskBand.GetOverviewCount()): if iOverview != 0: line = line + ", " hOverview = hMaskBand.GetOverview(iOverview) if hOverview is not None: line = line + ("%d" % hOverview.Checksum()) else: line = line + "(null)" if hBand.GetUnitType(): print(" Unit Type: %s" % hBand.GetUnitType()) papszCategories = hBand.GetRasterCategoryNames() if papszCategories is not None: print(" Categories:") i = 0 for category in papszCategories: print(" %3d: %s" % (i, category)) i = i + 1 scale = hBand.GetScale() if not scale: scale = 1.0 offset = hBand.GetOffset() if not offset: offset = 0.0 if scale != 1.0 or offset != 0.0: print(" Offset: %.15g, Scale:%.15g" % (offset, scale)) if bShowMetadata: papszMetadata = hBand.GetMetadata_List() else: papszMetadata = None if bShowMetadata and papszMetadata: print(" Metadata:") for metadata in papszMetadata: print(" %s" % metadata) if bShowMetadata: papszMetadata = hBand.GetMetadata_List("IMAGE_STRUCTURE") else: papszMetadata = None if bShowMetadata and papszMetadata: print(" Image Structure Metadata:") for metadata in papszMetadata: print(" %s" % metadata) hTable = hBand.GetRasterColorTable() if ( hBand.GetRasterColorInterpretation() == gdal.GCI_PaletteIndex and hTable is not None ): print( " Color Table (%s with %d entries)" % ( gdal.GetPaletteInterpretationName( hTable.GetPaletteInterpretation() ), hTable.GetCount(), ) ) if bShowColorTable: for i in range(hTable.GetCount()): sEntry = hTable.GetColorEntry(i) print( " %3d: %d,%d,%d,%d" % (i, sEntry[0], sEntry[1], sEntry[2], sEntry[3]) ) if bShowRAT: pass # hRAT = hBand.GetDefaultRAT() # GDALRATDumpReadable( hRAT, None ); return 0 # ********************************************************************** # GDALInfoReportCorner() # ********************************************************************** def GDALInfoReportCorner(hDataset, hTransform, corner_name, x, y): line = "%-11s " % corner_name # -------------------------------------------------------------------- # Transform the point into georeferenced coordinates. # -------------------------------------------------------------------- adfGeoTransform = hDataset.GetGeoTransform(can_return_null=True) if adfGeoTransform is not None: dfGeoX = adfGeoTransform[0] + adfGeoTransform[1] * x + adfGeoTransform[2] * y dfGeoY = adfGeoTransform[3] + adfGeoTransform[4] * x + adfGeoTransform[5] * y else: line = line + ("(%7.1f,%7.1f)" % (x, y)) print(line) return False # -------------------------------------------------------------------- # Report the georeferenced coordinates. # -------------------------------------------------------------------- if abs(dfGeoX) < 181 and abs(dfGeoY) < 91: line = line + ("(%12.7f,%12.7f) " % (dfGeoX, dfGeoY)) else: line = line + ("(%12.3f,%12.3f) " % (dfGeoX, dfGeoY)) # -------------------------------------------------------------------- # Transform to latlong and report. # -------------------------------------------------------------------- if hTransform is not None: pnt = hTransform.TransformPoint(dfGeoX, dfGeoY, 0) if pnt is not None: line = line + ("(%s," % gdal.DecToDMS(pnt[0], "Long", 2)) line = line + ("%s)" % gdal.DecToDMS(pnt[1], "Lat", 2)) print(line) return True if __name__ == "__main__": sys.exit(main(sys.argv))