用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的指导

1
2
3
4
5
6
# 导入包
from osgeo import gdal
import matplotlib.pyplot as plt
# 使用gdal打开img文件
dataset = gdal.Open('ASTGTM_S34E141F.img')

查看driver,可以理解为看它的格式:

1
2
print("Driver: {}/{}".format(dataset.GetDriver().ShortName,
dataset.GetDriver().LongName))

输出:

1
Driver: HFA/Erdas Imagine Images (.img)

这里附带提一下,我们知道ArcGrid格式是ESRI的二进制格式,GridFloat是USGS的DEM数据格式,IMG是Erdas的。

1
2
3
4
# 查看大小、波段
print("Size is {} x {} x {}".format(dataset.RasterXSize,
dataset.RasterYSize,
dataset.RasterCount))

输出:

1
2
Size is 3115 x 3712 x 1
# 行数 列数 波段数

查看投影方式:

1
print("Projection is {}".format(dataset.GetProjection()))

地理坐标转换的参数:

1
2
3
4
geotransform = dataset.GetGeoTransform()
if geotransform:
print("Origin = ({}, {})".format(geotransform[0], geotransform[3]))
print("Pixel Size = ({}, {})".format(geotransform[1], geotransform[5]))

输出:

1
2
Origin = (499987.025558, 6348728.45327)
Pixel Size = (30.0, -30.0)

查看数据类型:

1
2
band = dataset.GetRasterBand(1)
print("Band Type={}".format(gdal.GetDataTypeName(band.DataType)))

输出:

1
Band Type=Int16

还有一些:

1
2
3
4
5
6
7
8
print '元数据:'
print dataset.GetMetadata()
print 'space: (Mb)'
print (lc.nbytes/(8*1024*1024.))
print '数组的统计数据:'
print (lc.min(), lc.max(), lc.mean(), lc.std())

完整代码:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
# -*- 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())

参考资料:

GDAL API Tutorial

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