KJ135
/
microproduct
mirror of http://172.16.0.12:5000/WBJY/microproduct.git
4
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity
microproduct/dem-sentiral/ISCEApp/site-packages/osgeo_utils/samples/gdalinfo.py

611 lines
22 KiB
Python
Raw Permalink Blame History

#!/usr/bin/env python3
# *****************************************************************************
# $Id: gdalinfo.py e4fe7cc06270e5f38dfe78e6785a6bcca4e39e29 2021-04-01 21:02:04 +0300 Idan Miara $
#
# 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
from osgeo import osr
# **********************************************************************
# Usage()
# **********************************************************************
def Usage():
print("Usage: gdalinfo [--help-general] [-mm] [-stats] [-hist] [-nogcp] [-nomd]\n" +
" [-norat] [-noct] [-nofl] [-checksum] [-mdd domain]* datasetname")
return 1
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))
Reference in New Issue View Git Blame Copy Permalink