Python调用GDAL方法用法介绍|极客教程

Python调用GDAL方法用法介绍

1. 简介

GDAL (Geospatial Data Abstraction Library) 是一个用于处理地理信息数据的开源库，它包含了大量的功能，可以读取和写入各种栅格和矢量数据格式，进行数据变换和投影转换，以及执行空间分析等操作。GDAL 是许多 GIS 软件的核心组件，如 QGIS、ArcGIS 等。

Python 是一种简单易用且功能强大的编程语言，其拥有丰富的第三方库支持，便于进行地理空间数据处理和分析。通过调用 GDAL 库的 Python 绑定，我们可以使用 Python 进行高效的地理空间数据处理。

本文将详细介绍如何使用 Python 调用 GDAL 方法，包括数据读取、数据写入、数据转换处理等常用功能的示例与说明。

2. 安装

要在 Python 中使用 GDAL，首先需要安装 GDAL 库。GDAL 有两个版本，分别是 GDAL 1.x 和 GDAL 2.x。如果你使用的是 Python 2.x 版本，可以选择安装 GDAL 1.x，而如果你使用的是 Python 3.x 版本，建议安装 GDAL 2.x。

可以通过 pip 命令来快速安装 GDAL：

pip install GDAL

在安装过程中可能会因为缺少依赖库而导致失败，可以根据提示手动安装依赖库，然后重新运行安装命令。

3. 数据读取

使用 GDAL 可以读取各种栅格和矢量数据格式，如 TIFF、JPEG、Shapefile 等。下面是一个读取 TIFF 栅格数据的示例：

from osgeo import gdal

# 打开栅格数据集
dataset = gdal.Open('input.tif')

# 获取栅格数据的宽度和高度
width = dataset.RasterXSize
height = dataset.RasterYSize

# 获取栅格数据的通道数
bands = dataset.RasterCount

# 获取栅格数据的投影信息
projection = dataset.GetProjection()

# 获取栅格数据的地理变换参数
geotransform = dataset.GetGeoTransform()

# 读取栅格数据的像素值
band = dataset.GetRasterBand(1)
data = band.ReadAsArray(0, 0, width, height)

# 关闭栅格数据集
dataset = None

# 打印读取的数据信息
print('Width:', width)
print('Height:', height)
print('Bands:', bands)
print('Projection:', projection)
print('GeoTransform:', geotransform)
print('Data:', data)

上述代码首先使用 gdal.Open() 方法打开 TIFF 栅格数据集，然后通过 RasterXSize 和 RasterYSize 可以获取栅格数据的宽度和高度。使用 RasterCount 可以获取栅格数据的通道数。

使用 GetProjection() 和 GetGeoTransform() 可以获取栅格数据的投影信息和地理变换参数。

使用 GetRasterBand() 方法可以获取栅格数据的通道，然后使用 ReadAsArray() 方法读取栅格数据的像素值。

最后通过 print() 方法打印读取的数据信息。

4. 数据写入

GDAL 不仅可以读取地理空间数据，还可以将处理后的数据写入到文件中。下面是一个将 NumPy 数组保存为 GeoTIFF 文件的示例：

from osgeo import gdal

# 创建栅格数据
data = numpy.zeros((height, width), dtype=numpy.uint8)
data[100:200, 200:300] = 255

# 创建栅格驱动对象
driver = gdal.GetDriverByName('GTiff')

# 创建栅格数据集
dataset = driver.Create('output.tif', width, height, 1, gdal.GDT_Byte)

# 设置栅格数据的地理变换参数
dataset.SetGeoTransform(geotransform)

# 设置栅格数据的投影信息
dataset.SetProjection(projection)

# 写入栅格数据的像素值
band = dataset.GetRasterBand(1)
band.WriteArray(data, 0, 0)

# 关闭栅格数据集
dataset = None

上述代码首先创建了一个全零的 height x width 的 NumPy 数组，然后将其中的一部分像素值设置为 255。

使用 gdal.GetDriverByName() 方法获取栅格驱动对象，这里选择了 GeoTIFF 格式。

使用 driver.Create() 方法创建栅格数据集，指定了输出文件名、宽度、高度、通道数、像素数据类型等参数。

使用 SetGeoTransform() 和 SetProjection() 方法设置栅格数据的地理变换参数和投影信息。

使用 GetRasterBand() 方法获取栅格数据的通道，然后使用 WriteArray() 方法将像素值写入栅格数据。

最后通过将数据集对象设为 None 关闭栅格数据集。

5. 数据转换处理

GDAL 提供了丰富的数据转换和处理方法，可以进行栅格数据的投影转换、重采样、剪裁、融合、光栅代数运算等操作。下面是一个实现栅格数据剪裁和融合的示例：

from osgeo import gdal, osr

# 打开待剪裁的栅格数据集
input_dataset = gdal.Open('input.tif')

# 创建剪裁后栅格数据集
output_dataset = gdal.Warp('output.tif', input_dataset, outputBounds=(100, 200, 300, 400))

# 创建融合后栅格数据集
merge_dataset = gdal.Warp('merge.tif', [input_dataset, output_dataset], options=gdal.WarpOptions(format='GTiff', outputBounds=(0, 0, 500, 500), outputResolution=(10, 10)))

# 关闭栅格数据集
input_dataset = None
output_dataset = None
merge_dataset = None

上述代码首先使用 gdal.Open() 方法打开待剪裁的栅格数据集，然后使用 gdal.Warp() 方法创建剪裁后的栅格数据集，通过设定 outputBounds 参数来指定剪裁的范围。

使用 gdal.Warp() 方法还可以实现多个栅格数据集的融合，通过将待融合的数据集以列表形式传入函数中，并设定 outputResolution 参数、outputBounds 参数等来指定融合后数据的分辨率和范围。

最后通过将数据集对象设为 None 关闭栅格数据集。