目录
一、栅格数据驱动
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 同步更新