从GDAL提供的实用程序来看,很多程序的后缀都是 .py
,这充分地说明了Python语言在GDAL的开发中得到了广泛的应用。
from osgeo import gdal
dataset = gdal.Open("/data/gdata/geotiff_file.tif")
既然已经将一个GeoTIFF文件打开为一个GDAL可操作的对象, 下面来看一下都能对其进行怎样的操作。
Python提供了 dir()
内省函数, 可以快速查看一下当前对象可用的操作:
dir()
函数可能是 Python 自省机制中最著名的部分了。它可以返回传递给它的任何对象的属性名称经过排序的列表。 如果不指定对象,则dir()
返回当前作用域中的名称。
dir(dataset)[:10]
['AbortSQL', 'AddBand', 'AddFieldDomain', 'AddRelationship', 'AdviseRead', 'BeginAsyncReader', 'BuildOverviews', 'ClearStatistics', 'Close', 'CommitTransaction']
下面看一下如何获取文件的一些基本信息,需要用到下面的一些函数与属性。
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 打开数据来读取一下数据的信息:
dataset = gdal.Open("/data/gdata/geotiff_file.tif")
dataset.GetMetadata()
{'AREA_OR_POINT': 'Area', 'PyramidResamplingType': 'NEAREST'}
GetMetadata()
方法可以访问数据的元数据信息,元数据信息对于每个数据都是不一样的。
比如再打开另外一个文件:
ds = gdal.Open('/data/gdata/lu75c.tif')
ds.GetMetadata()
{'AREA_OR_POINT': 'Area', 'TIFFTAG_XRESOLUTION': '1', 'TIFFTAG_YRESOLUTION': '1'}
dataset.GetDescription()
'/data/gdata/geotiff_file.tif'
dataset.RasterCount
3
这是一个由7个波段构成的Landsat遥感影像。
再看一个MODIS L1B数据:
mds = gdal.Open("/data/gdata/MOD09A1.A2009193.h28v06.005.2009203125525.hdf")
mds.RasterCount
0
运行结果居然是 0
。这意味着当前的数据集 dataset
中的栅格数目是 0
。
实际上,MODISL1B的数据格式是HDF格式的,它的数据是以子数据集组织的,要获取其相关的数据的信息,需要继续访问其子数据集。
dataset.RasterXSize,dataset.RasterYSize
(1500, 900)
ds.GetGeoTransform()
(1852951.7603168152, 30.0, 0.0, 5309350.360150607, 0.0, -30.0)
ds.GetProjection()
'PROJCS["unnamed",GEOGCS["unknown",DATUM["unnamed",SPHEROID["unretrievable - using WGS84",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]]],PROJECTION["Albers_Conic_Equal_Area"],PARAMETER["latitude_of_center",0],PARAMETER["longitude_of_center",105],PARAMETER["standard_parallel_1",25],PARAMETER["standard_parallel_2",47],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH]]'
from osgeo import gdal
dataset = gdal.Open('/data/gdata/lu75c.tif')
dataset.RasterCount
1
band = dataset.GetRasterBand(1)
band = dataset.GetRasterBand(1)
dir(band)[:10]
['AdviseRead', 'AsMDArray', 'Checksum', 'ComputeBandStats', 'ComputeRasterMinMax', 'ComputeStatistics', 'CreateMaskBand', 'DataType', 'DeleteNoDataValue', 'Fill']
band.XSize, band.YSize, band.DataType
(6122, 4669, 3)
band.GetNoDataValue()
band.GetMaximum()
band.GetMinimum()
band.ComputeRasterMinMax()
(-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数据的读写了。