用GDAL打开从USGS下载的img影像文件
来源于我的独立博客：用GDAL打开从USGS下载的img影像文件
缘由：我想找全球DEM数据，在这里发现了5种免费的数据库，其中一个是SRTM DEM由USGS提供，可以在这里下载。
下载的数据是GeoTIFF格式，如不了解，可以看wikipedia的解释：GeoTIFF。这个格式的数据matlab有相应的函数可以打开，名叫：geotiffread ，但是要高版本才行，我用的2014b不行。
当然最简单的方式是使用现有软件，将数据直接拖进软件窗口即可显示，常用软件有ArcGis，它虽然是全球最好的GIS软件，但是收费的，免费的有QGIS 。
但如果想用代码的方式来读数据、显示怎么办？matlab是个不错的选择，但我的2014不支持，并且matlab也是收费的，这里用GDAL。
接下来跟我一起实验一下吧。
首先在USGS下载一个GeoTIFF格式的影像，（我这里用的是DEM数据，它也是以影像的形式呈现的，只有一个波段）。
确认已经安装了Python，我这里用的Anaconda，python27. gdal我包直接在Anaconda中安装。
参考gdal的指导：
# 导入包
from osgeo import gdal
import matplotlib.pyplot as plt
# 使用gdal打开img文件
dataset = gdal.Open('ASTGTM_S34E141F.img')查看driver，可以理解为看它的格式:
print("Driver: {}/{}".format(dataset.GetDriver().ShortName,
dataset.GetDriver().LongName))输出：
Driver: HFA/Erdas Imagine Images (.img)这里附带提一下，我们知道ArcGrid格式是ESRI的二进制格式，GridFloat是USGS的DEM数据格式，IMG是Erdas的。
# 查看大小、波段
print("Size is {} x {} x {}".format(dataset.RasterXSize,
dataset.RasterYSize,
dataset.RasterCount))输出：
Size is 3115 x 3712 x 1
# 行数 列数 波段数查看投影方式：
print("Projection is {}".format(dataset.GetProjection()))地理坐标转换的参数:
geotransform = dataset.GetGeoTransform()
if geotransform:
print("Origin = ({}, {})".format(geotransform[0], geotransform[3]))
print("Pixel Size = ({}, {})".format(geotransform[1], geotransform[5]))输出：
Origin = (499987.025558, 6348728.45327)
Pixel Size = (30.0, -30.0)查看数据类型:
band = dataset.GetRasterBand(1)
print("Band Type={}".format(gdal.GetDataTypeName(band.DataType)))输出：
Band Type=Int16还有一些：
print '元数据：'
print dataset.GetMetadata()
print 'space: (Mb)'
print (lc.nbytes/(8*1024*1024.))
print '数组的统计数据：'
print (lc.min(), lc.max(), lc.mean(), lc.std())完整代码：
# -*- coding: utf-8 -*-
"""
Created on Thu Jan 04 19:55:49 2018
@author: Administrator
"""
# 导入包
from osgeo import gdal
import matplotlib.pyplot as plt
# 使用gdal打开img文件
dataset = gdal.Open('ASTGTM_S34E141F.img')
# 查看driver，可以理解为看它的格式
print("Driver: {}/{}".format(dataset.GetDriver().ShortName,
dataset.GetDriver().LongName))
# 查看大小、波段
print("Size is {} x {} x {}".format(dataset.RasterXSize,
dataset.RasterYSize,
dataset.RasterCount))
# 查看投影方式
#print("Projection is {}".format(dataset.GetProjection()))
# 地理坐标转换的参数
geotransform = dataset.GetGeoTransform()
if geotransform:
print("Origin = ({}, {})".format(geotransform[0], geotransform[3]))
print("Pixel Size = ({}, {})".format(geotransform[1], geotransform[5]))
# 查看数据类型
band = dataset.GetRasterBand(1)
print("Band Type={}".format(gdal.GetDataTypeName(band.DataType)))
min = band.GetMinimum()
max = band.GetMaximum()
if not min or not max:
(min,max) = band.ComputeRasterMinMax(True)
print("Min={:.3f}, Max={:.3f}".format(min,max))
if band.GetOverviewCount() > 0:
print("Band has {} overviews".format(band.GetOverviewCount()))
if band.GetRasterColorTable():
print("Band has a color table with {} entries".format(band.GetRasterColorTable().GetCount()))
# 将图像读到数组
lc = dataset.ReadAsArray()
print '图像大小（x,y）：'
print lc.shape
# 32767是NaN值，将其替换为0
# 注意这种值如果不替换掉会导致图像显示异常，因为它太大了，可能显示出来全是黑的
lc[lc == 32767] = 0
# 显示 colormap 用 jet
plt.imshow ( lc, interpolation='nearest', vmin=0, cmap=plt.cm.jet)
print '坐标转换的参数:'
print dataset.GetGeoTransform()
print '元数据：'
print dataset.GetMetadata()
print 'space: (Mb)'
print (lc.nbytes/(8*1024*1024.))
print '数组的统计数据：'
print (lc.min(), lc.max(), lc.mean(), lc.std())参考资料：
stackexchange: Reading National Elevation Dataset (ArcGrid/GridFloat/IMG) with Python only tools?
Starting to use Python to work with geospatial data
Choosing Colormaps， 这里可以选colormaps
Welcome to the Python GDAL/OGR Cookbook!
5 Free Global DEM Data Sources – Digital Elevation Models
转载自：https://blog.csdn.net/shanchuan2012/article/details/79039547
