1. PythonGIS
  2. GDAL
  3. 读取遥感影像的信息

读取遥感影像的信息

打开已有的GeoTIF文件

下面我们试着读取一个GeoTiff文件的信息。第一步就是打开一个数据集。

In [1]:
from osgeo import gdal
dataset = gdal.Open("/gdata/K52E015007.tif")

既然已经将一个GeoTIFF文件打开为一个GDAL可操作的对象， 下面来看一下都能对其进行怎样的操作。

Python提供了 dir() 内省函数， 可以快速查看一下当前对象可用的操作：

dir() 函数可能是 Python 自省机制中最著名的部分了。它可以返回传递给它的任何对象的属性名称经过排序的列表。 如果不指定对象，则 dir() 返回当前作用域中的名称。

In [2]:
dir(dataset)[:3] + ['... ...'] + dir(dataset)[-3:]
Out[2]:
['__bool__',
 '__class__',
 '__delattr__',
 '... ...',
 '__sizeof__',
 '__str__',
 '__subclasshook__']

下面看一下如何获取文件的一些基本信息，需要用到下面的一些函数与属性。

  • dataset.GetDescription() 获得栅格的描述信息
  • dataset.RasterCount 获得栅格数据集的波段数
  • dataset.RasterXSize 栅格数据的宽度(X方向上的像素个数)
  • dataset.RasterYSize 栅格数据的高度(Y方向上的像素个数)
  • dataset.GetGeoTransform() 栅格数据的六参数。
  • GetProjection() 栅格数据的投影

通过下面小节我们逐一解释一下。

读取影像的元数据

从元数据的作用来看，它更多地是为工程服务的。 客观地说，GDAL对元数据的支持并不好， 它并没有直接的元数据处理接口， 更没有实现元数据的相关标准。 但是，GDAL已经提供了足够方便的函数，可以读取影像的一些信息， 从而方便对数据进行处理。 GDAL一般是以字典的形式对元数据进行组织的， 但是对于不同的栅格数据类型，元数据的类型与键值可能都不一样。

目前， 国际上对空间元数据标准内容进行研究的组织主要有三个，分别是欧洲标准化委员会(CEN/TC 287)、 美国联邦地理数据委员会(FGDC)和国际标准化组织地理信息/地球信息技术委员会(ISO/TC 211)。

如果要进行元数据处理，可以考虑将元数据信息写入到XML文件中。 这个问题扩展开来就不是本书关心的内容的，在此就不再多说了。

我们先来看一下最常用的GeoTIFF文件的元数据信息。 GDAL可以作为数据集级别的元数据来处理下面的基本的TIFF标志。

  • TIFFTAG_DOCUMENTNAME
  • TIFFTAG_IMAGEDESCRIPTION
  • TIFFTAG_SOFTWARE
  • TIFFTAG_DATETIME
  • TIFFTAG_ARTIST
  • TIFFTAG_HOSTCOMPUTER
  • TIFFTAG_COPYRIGHT
  • TIFFTAG_XRESOLUTION
  • TIFFTAG_YRESOLUTION
  • TIFFTAG_RESOLUTIONUNIT
  • TIFFTAG_MINSAMPLEVALUE (read only)
  • TIFFTAG_MAXSAMPLEVALUE (read only)

使用Python来访问一下：

In [3]:
from osgeo import gdal
dataset = gdal.Open("/gdata/geotiff_file.tif")
dataset.GetMetadata()
Out[3]:
{'AREA_OR_POINT': 'Area', 'PyramidResamplingType': 'NEAREST'}

上面的元数据信息，对于每个数据都是不一样的。

比如再打开另外一个文件：

In [4]:
ds = gdal.Open('/gdata/lu75c.tif')
ds.GetMetadata()
Out[4]:
{'AREA_OR_POINT': 'Area',
 'TIFFTAG_XRESOLUTION': '1',
 'TIFFTAG_YRESOLUTION': '1'}

这个文件只对两个TIFF标志进行了定义，还有一个并不是TIFF标志定义的。

GetDescription()获得栅格的描述信息

In [5]:
dataset.GetDescription()
Out[5]:
'/gdata/geotiff_file.tif'

看来这里的图像描述是图像的路径名， 但是这是和各种不同数据集相关的， 不同数据集有不同的描述。

获取栅格数目

栅格数据集是由多个数据构成的，在GDAL中，每一个波段，都是一个数据集； 不仅如此，栅格数据集还可能包含有子数据集，每子数据集又可能包含有波段。

先来看一下刚才打开的数据的RasterCount：

In [6]:
dataset.RasterCount
Out[6]:
3

这是一个由7个波段构成的Landsat遥感影像。

再看一个MODIS L1B数据： mds = gdal.Open(’MOD09A1.A2009193.h28v06.005.2009203125525.hdf’) mds.RasterCount

In [7]:
mds = gdal.Open("/gdata/MOD09A1.A2009193.h28v06.005.2009203125525.hdf")
mds.RasterCount
Out[7]:
0

运行结果居然是0。这意味着当前的数据集 dataset 中的栅格数目是0。 实际上，MODISL1B的数据格式是HDF格式的， 它的数据是以子数据集组织的， 要获取其相关的数据的信息， 需要继续访问其子数据集。 这一部分会在第[subsec:gdalhdf]节再进行介绍。

影像大小

栅格数据的大小指出了影像以像元为单位的宽度与高度。

In [8]:
img_width,img_height = dataset.RasterXSize,dataset.RasterYSize
img_width,img_height
Out[8]:
(1500, 900)

可以看出我们的影像大小是 10572 × 9422 。

获得空间参考

下面看一下如果从栅格数据集中获取其投影与空间参考信息。 更多的关于投影与空间参考的讨论，会在后面章节介绍。

GetGeoTransform() 地理仿射变换参数。

对于遥感影像来说，它需要在地理空间中进行定位。 在GDAL中，这有两种方式，其中一种是使用六个参数坐标转换模型。 这个模型的具体实现在不同的软件中是不一样的。 在GDAL中，这六个参数包括左上角坐标，像元X、Y方向大小，旋转等信息。 要注意，Y方向的像元大小为负值。

In [9]:
ds.GetGeoTransform()
Out[9]:
(1852951.7603168152, 30.0, 0.0, 5309350.360150607, 0.0, -30.0)

获得投影信息

使用 GetProjection() 函数，可以比较容易地获取数据集的投影信息， 但是对于什么是地图投影，以及如何在GDAL中实现，就不是这么容易了。 此处我们只是简单地运行一下，更详细的解释，会在后面章节中展开。

In [10]:
ds.GetProjection()
Out[10]:
'PROJCS["unnamed",GEOGCS["unknown",DATUM["unknown",SPHEROID["unretrievable - using WGS84",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]],PROJECTION["Albers_Conic_Equal_Area"],PARAMETER["standard_parallel_1",25],PARAMETER["standard_parallel_2",47],PARAMETER["latitude_of_center",0],PARAMETER["longitude_of_center",105],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]]]'

使用GDAL获取栅格数据波段信息

上面我们介绍了针对数据集操作的主要函数。 但是如果需要了解栅格数据的更多信息，我们就需要看一下遥感图像处理中更常用到的波段操作的函数。

获取数据集的波段

GetRasterBand() 函数，可以获得栅格数据集的波段。这是函数的参数使用波段的索引值。

In [11]:
from osgeo import gdal
dataset = gdal.Open('/gdata/lu75c.tif')
dataset.RasterCount
Out[11]:
1
In [12]:
band = dataset.GetRasterBand(1)

这里我们获取了第一个波段(红色值组成的表)。 注意！这里的波段获取和通常的C数组获取不一样，开始是1不是0。

查看波段的基本信息

下面我们现在来看看刚才读取出来的 band 有些什么东西可以供我们操作。

与操作数据集一样，GDAL同样提供了查看波段基本信息的函数。

In [13]:
band = dataset.GetRasterBand(1)
dir(band)[:6] + ['... ...'] + dir(band)[-3:]
Out[13]:
['Checksum',
 'ComputeBandStats',
 'ComputeRasterMinMax',
 'ComputeStatistics',
 'CreateMaskBand',
 'DataType',
 '... ...',
 '__weakref__',
 '_s',
 'this']

看一下常用的操作。这些也是用来获取波段的属性信息。

获取波段大小

In [14]:
band.XSize
Out[14]:
6122
In [15]:
band.YSize
Out[15]:
4669
In [16]:
band.DataType
Out[16]:
3

执行以上代码得到了波段图像的宽和高（像元为单位）。 对于我们所使用的影像， 这个与 dataset 中使用 RasterXSize() 与 RasterYSize() 获取的值一致。DataType是图像中实际数值的数据类型，表示8位无符整型。 在[subsec:banddatatype]节有更进一步的解释。

获取波段数据的属性

In [17]:
band.GetNoDataValue()
band.GetMaximum()
band.GetMinimum()
band.ComputeRasterMinMax()
Out[17]:
(-1.0, 66.0)

Maximum 是表示在本波段数值中最大的值，当然 Minimum 就是表示本波段中最小的值。 通过运行结果，我们可以看到在一开始RasterXSize()和RasterYSize()都没有值。 因为对于文件格式不会有固有的最大最小值。 所以我们可以通过函数ComputeRasterMinMax() 计算得到。 注意！这里的最大最小值不包括“无意义值”！ 也就是上面显示的 NoDataValue 。

其他格式

使用GDAL读取ENVI数据格式

ENVI栅格文件格式

ENVI使用的是通用栅格数据格式，包含一个简单的二进制文件（a simple flat binary）和一个相关的ASCII（文本）的头文件。这也保证了单个ENVI栅格文件没有大小上限。

（1）头文件

ENVI头文件包含用于读取图像数据文件的信息，它通常创建于一个数据文件第一次被ENVI读取时。单独的ENVI头文本文件提供关于图像尺寸、 嵌入的头文件（若存在）、数据格式及其它相关信息。所需信息通过交互式输入，或自动地用“文件吸取”创建，并且以后可以编辑修改。 使用者可以在ENVI之外使用一个文本编辑器生成一个ENVI头文件。

（2）数据文件

通用栅格数据都会存储为二进制的字节流，通常它将以BSQ（band sequential，按波段顺序储存）、BIP（band interleaved by pixel，按波段像元交叉储存）或者BIL（band interleaved by line，按波段行交叉储存）的方式进行存储。

ENVI栅格文件必须包含着两个文件，其中头文件的后缀名为：.hdr，数据文件的后缀随意，甚至可以不带后缀名。 这两个文件是通过文件名来关联，即数据文件和头文件名称一致。

GDAL读取HDF数据格式

由于modis卫星数据跟我们经常遇到的geotif数据组织方式不一样，读取的时候一定要特别注意。geotif数据，一般是一个文件，包含了多个波段的数据； 而modis呢，一个文件包含了多各SUBDATASETSGDAL， 每个SUBDATASETS又包含多个波段数据。 另外默认编译的GDAL并不包含对MODIS数据支持， 需要单独下载针对HDF4，HDF5的源码，再修改下make.opt文件， 这时再编译GDAL，就支持modis数据的读写了。

开始使用GDAL

从GDAL提供的实用程序来看，很多程序的后缀都是 .py ，这充分地说明了Python语言在GDAL的开发中得到了广泛的应用。

导入GDAL

在Python中使用GDAL，只需要导入 gdal 模块。 在早期的版本（1.5以前）中，GDAL是使用下面的语句导入的：

In [18]:
import gdal

但是后来GDAL成为OSGEO的子项目后，对代码进行了重新组织。 在 GDAL RFC 17号文件 中， 实现了Python的新的名称空间osgeo， 并将gdal与ogr都包含在这个名称空间之下。

RFC(Request For Comments)，意即“请求评议”，一系列以编号排定的文件。 当某家机构或团体开发出了一套标准或提出对某种标准的设想，想要征询外界的意见时， 就会在Internet上发放一份RFC，对这一问题感兴趣的人可以阅读该RFC并提出自己的意见。

1.6以后，推荐使用下面的语句导入：

In [19]:
from osgeo import gdal

当然早期版本也是支持的，但是在有些版本使用时候会产生一个弃用警告：

In [20]:
import gdal

usr(/lib/python2.6/dist-packages/osgeo/gdal.py:99: DeprecationWarning: gdal.py was placed in a namespace, it is now available DeprecationWarning

```

而在最新的版本中，这个弃用警告又消失了。

为了保持兼容性，可以使用下面的语句来导入：

In [21]:
try:
    import gdal
except:
    from osgeo import gdal

除了gdal包，还有一个gdalconst包最好也导入。 gdalconst也是osgeo的一个包， 它只是在代码中对GDAL中用到的一些常量进行了绑定。 其中gdalconst中的常量都加了前缀，力图与其他模块冲突最小。 所以对gdalconst你可以直接这样导入：

In [22]:
from osgeo.gdalconst import *

驱动

要读取某种类型的数据，必须要先载入数据驱动（Driver），也就是初始化一个对象，让它“知道”某种数据结构。 可以使用下面语句一次性注册所有的数据驱动，但是只能读不能写：

In [23]:
gdal.AllRegister()

单独注册某一类型的数据驱动，可以读也可以写，可以新建数据集 （这最终还要取决于GDAL是否已经进行了实现）。 一些不同类型的数据驱动的编码，在表[tab:gdaL~f~ormat]中已经介绍过了。

下面的语句注册了Erdas的栅格数据类型。

In [24]:
driver = gdal.GetDriverByName('HFA')
driver.Register()
Out[24]:
5

可以使用下面的语句判断driver是否注册成功。

In [25]:
driver = gdal.GetDriverByName('GTiff')
driver == None
Out[25]:
False

上面的注册就失败了，因为不存在名称为“GTiff”的数据格式（正确的格式为“GeoTiff”）。

查看系统支持的数据格式

GDAL除了使用 GetDriverByName()获得驱动，还可以使用 GetDriver() 。下面的代码获取了系统支持的所有的驱动的名称。

In [26]:
drv_count = gdal.GetDriverCount()
drv_count
Out[26]:
202

对于不同的Linux发行版，以及安装的GDAL的版本与编译选项的不同， 上面程序的结果是不一样的。 所以一般情况下要避免使用 gdal.GetDriver()， 而要使用gdal.GetDriverByName() 这个函数来获取驱动。

我使用的系统是Debian Squeeze，返回的驱动的个数是88； 在Debian Wheezy中，返回的个数是114；在Debian Jessie中，则增加到120个。这个一方面说明GDAL的发展还是很快的， 另一方面说明GDAL目前已经比较成熟了。

In [27]:
for idx in range(10):
    driver = gdal.GetDriver(idx)
    print( "%10s: %s" % (driver.ShortName, driver.LongName))

       VRT: Virtual Raster
     GTiff: GeoTIFF
      NITF: National Imagery Transmission Format
    RPFTOC: Raster Product Format TOC format
   ECRGTOC: ECRG TOC format
       HFA: Erdas Imagine Images (.img)
  SAR_CEOS: CEOS SAR Image
      CEOS: CEOS Image
JAXAPALSAR: JAXA PALSAR Product Reader (Level 1.1/1.5)
       GFF: Ground-based SAR Applications Testbed File Format (.gff)

上面第4行，直接使用了索引值来获得驱动，而在第5行则打印了驱动的名称。注意到驱动有 ShortName 与 LongName。ShortName与栅格数据格式在GDAL中定义的编码是一致的， 而LongName则可以看成是描述性的文字。

版权归 OSGeo中国中心（http://www.osgeo.cn/ ） 所有
邮箱： bukun@iga.ac.cn