GDAL / OGR 学习手册 [02] ：栅格数据读取

目录

一、栅格数据驱动

二、gdal.Open

三、gdal.Dataset

四、获取影像的基本信息

1. 获取影像元数据

2. 获取影像基本信息

---

一、栅格数据驱动

GDAL 通过数据驱动来识别各种类型的栅格数据，目前已经支持GeoTIFF、ERDAS IMAGINE、HDF、netCDF等百余种数据格式。

 gdal 支持的栅格数据格式：https://gdal.org/drivers/raster/index.html#raster-drivers

早期的版本中，在读取栅格数据时，需要先载入数据驱动，通过 GetDriverByName 方法，传入驱动名（ShortName）获取到对应的数据驱动，然后使用 Register 方法完成注册。

driver = gdal.GetDriverByName('HFA')
driver.Register()

除此之外，还可以使用 GetDrive 方法来获取驱动（传入驱动索引值），如下代码可以查询所有驱动的索引值对应关系以及 ShortName 和 LongName：

from osgeo import gdal
drive_count = gdal.GetDriverCount()
for drive_index in range(drive_count):
    driver = gdal.GetDriver(drive_index)
    print( "%10s: %s" % (driver.ShortName, driver.LongName))

在最近的版本中，读写数据已经不需要注册驱动，只有在创建新数据时，需要通过数据驱动来完成数据创建。

二、gdal.Open

我们试着读取一个 GeoTiff 格式的影像文件。

根据 gdal 官方提供的 API 文档，使用 gdal.Open 方法，可将栅格类型的影像文件打开为一个 GDAL 可操作的对象：gdal 数据集（gdal.Dataset）。

 gdal.Open 两个参数：

     ① 文件路径 [ type: str ]

     ② 访问模式 [ type: int 可选 ]       0：<只读>   1：<可读写>

首先按照最简单的方式读取（不传入第二个参数），此时默认以只读模式打开数据集，在只读模式无法对栅格数据进行写入操作。

from osgeo import gdal
​
dataset = gdal.Open('testGeoTiff.tif')
print(dataset)
# <osgeo.gdal.Dataset; proxy of <Swig Object of type 'GDALDatasetShadow *' at 0x000001F2E54F8FC0>

当需要对栅格数据进行写入操作时，可以在打开数据集时指定第二个参数为 1（可读写）。

from osgeo import gdal
​
dataset = gdal.Open('testGeoTiff.tif',1)

对于访问模式的指定，推荐引入 gdalconst 包，它对 gdal 中用到的一些常量进行了绑定，添加了特定的前缀，可以减少与其他模块的冲突。不仅如此，gdalconst 中常量的命名与 API 文档一致，能够帮助调用者快速定位传入的参数。在 gdalconst 包中 GA_Update 被定义为1，GA_ReadOnly 被定义为0，使用如下写法引入 GA_Update ，也能起到同样的效果。

from osgeo import gdal
from osgeo.gdalconst import GA_Update
​
dataset = gdal.Open('testGeoTiff.tif',GA_Update)

三、gdal.Dataset

GDAL 数据集（gdal.Dataset）是基于 OGC 格式定义的数据集。它是相关栅格波段和一些相关信息的集合，不仅包含所有波段的地理参考和坐标系统定义，且本身也携带相关的元数据。

使用 Python 内置的 dir 函数可以查看 gdal.Dataset 携带的属性和方法。

from osgeo import gdal
from osgeo.gdalconst import GA_Update
​
# 打开Dataset
dataset = gdal.Open('testGeoTiff.tif',GA_Update)
# 遍历Dataset的属性和方法：
Mtds_and_Attrs = dir(dataset)
print('gdal.dataset的属性和方法：')
for name in Mtds_and_Attrs:
    print(' => ',name)

 gdal.Dataset 相关API：https://gdal.org/api/python/osgeo.gdal.html#osgeo.gdal.Dataset

四、获取影像的基本信息

1. 获取影像元数据

gdal.Dataset 可以调用继承其基类 gdal.MajorObject 的一系列方法来完成元数据的读写。

例如使用 GetMetadata 方法获取影像文件的元数据，元数据信息会以键值对（dict 格式）返回。

from osgeo import gdal
from osgeo.gdalconst import GA_Update
​
# 打开Dataset
dataset = gdal.Open('testGeoTiff.tif',GA_Update)
# 获取影像文件的元数据
metaData = dataset.GetMetadata()
# 打印获取到的元数据
print(metaData)
​
# {'AREA_OR_POINT': 'Area',
# 'TIFFTAG_XRESOLUTION': '1',
# 'TIFFTAG_YRESOLUTION': '1'}

使用 GetDescription 方法可以获取影像文件的描述信息：

from osgeo import gdal
from osgeo.gdalconst import GA_Update
# 打开Dataset
dataset = gdal.Open('testGeoTiff.tif',GA_Update)
# 获取影像文件的描述信息
dataDescription = dataset.GetDescription()
# 打印获取到的描述信息
print(dataDescription)
​
# testGeoTiff.tif 
# 该影像的描述信息是图像的路径名，这与数据集相关

2. 获取影像基本信息

（1）影像大小

利用 RasterXSize 和 RasterYSize 属性可以获取影像以像元（pixel）为单位的宽和高。

from osgeo import gdal
from osgeo.gdalconst import GA_Update
# 打开Dataset
dataset=gdal.Open('testGeoTiff.tif',GA_Update)
# 获取影像的宽和高
Width=dataset.RasterXSize
Height=dataset.RasterYSize
print("Width:{0},Height:{1}".format(Width,Height))
​
# Width:7721,Height:786

（2）空间参考和投影信息

利用 GetGeoTransform 和 GetProjection 可以分别获取影像的空间参考和投影信息。

from osgeo import gdal
from osgeo.gdalconst import GA_Update
# 打开Dataset
dataset = gdal.Open('testGeoTiff.tif',GA_Update)
# 打印影像的空间参考和投影信息
print(dataset.GetGeoTransform())
print(dataset.GetProjection())
​
# (244485.0, 30.0, 0.0, 3471015.0, 0.0, -30.0)
# PROJCS["WGS 84 / UTM zone 51N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]]......

gdal.Dataset 在表示栅格位置（像素）和地理参考坐标之间的关系时，使用的是仿射变换。GetGeoTransform 方法返回的是仿射变换的六个参数，gdal 通过六个参数构建坐标转换模型，将像素坐标转换到地理参考空间。

GetProjection 方法返回的是影像文件的投影信息，具体地图投影的相关内容，会在后面的章节中详细展开。

（3）波段数和波段

利用 RasterCount 属性，可以获取到影像的波段数。利用 GetRasterBand 方法，传入波段的索引值 [ type: int ]，可以获取到对应的波段。

from osgeo import gdal
from osgeo.gdalconst import GA_Update
​
# 打开Dataset
dataset = gdal.Open('testGeoTiff.tif',GA_Update)
# 获取影像的波段数
bandCount = dataset.RasterCount
print(bandCount)
# 获取影像的第一个波段
red_Band = dataset.GetRasterBand(1)
print(red_Band)
​
# 4
# <osgeo.gdal.Band; proxy of <Swig Object of type 'GDALRasterBandShadow *' at 0x0000016902F49ED0> >

P.S. 这里的波段索引和通常的数组索引不一样，是从1开始而不是0。

通过 GetRasterBand 方法获取到的波段为 gdal.Band 类型，后续波段数据的读写也需要借助这个类来实现。下一节将以 gdal.Band 为核心详细展开波段数据的处理方法。

【往期内容】

GDAL / OGR 学习手册 [01] ：Python环境配置https://blog.csdn.net/eternity1412/article/details/125111165GDAL / OGR 学习手册 | 前言https://blog.csdn.net/eternity1412/article/details/124561608

微信公众号 GISea 同步更新