Python库：gdal

---

GDAL（Geospatial Data Abstraction Library）是处理地理空间数据的开源库，广泛用于地理信息系统（GIS）和遥感领域。GDAL 提供了读取、写入、转换和处理栅格数据（如卫星图像、DEM等）和矢量数据（如Shapefile）的功能。

Python 通过 gdal 库提供了对 GDAL 的绑定，使得开发者可以在 Python 环境中使用 GDAL 的功能。

1. 安装 gdal 库

在安装 gdal 库之前，需要确保已经安装了 GDAL 库本身。然后可以通过 pip 安装 Python 的 gdal 绑定：

pip install gdal

2. 基本概念

在使用 gdal 库之前，了解一些基本概念是很有帮助的：

• 栅格数据：栅格数据是由像素组成的图像，每个像素都有一个值，通常表示某种测量值（如高度、温度、反射率等）。
• 矢量数据：矢量数据由点、线、面等几何对象组成，通常用于表示地理特征（如道路、建筑物、河流等）。
• 投影：地理数据通常需要投影到二维平面上，常见的投影方式包括 UTM、WGS84 等。
• 坐标系：地理数据通常使用特定的坐标系来表示位置，如 WGS84（经纬度）、UTM 等。

3. 读取栅格数据

gdal 库可以读取多种栅格数据格式，如 GeoTIFF、JPEG、PNG 等。以下是一个简单的示例，展示如何读取一个 GeoTIFF 文件并获取其基本信息。

from osgeo import gdal

# 打开栅格数据文件
dataset = gdal.Open('example.tif')

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

# 获取栅格数据的波段数量
bands = dataset.RasterCount

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

# 获取栅格数据的地理变换信息（仿射变换参数）
geotransform = dataset.GetGeoTransform()

# 获取第一个波段的数据
band = dataset.GetRasterBand(1)
data = band.ReadAsArray()

# 打印基本信息
print(f"Width: {width}, Height: {height}, Bands: {bands}")
print(f"Projection: {projection}")
print(f"GeoTransform: {geotransform}")
print(f"First band data shape: {data.shape}")

4. 读取矢量数据

gdal 库也可以读取矢量数据，如 Shapefile。以下是一个简单的示例，展示如何读取一个 Shapefile 并获取其几何对象和属性。

from osgeo import ogr

# 打开矢量数据文件
driver = ogr.GetDriverByName('ESRI Shapefile')
dataSource = driver.Open('example.shp', 0)

# 获取图层
layer = dataSource.GetLayer()

# 遍历图层中的每个要素
for feature in layer:
    # 获取几何对象
    geometry = feature.GetGeometryRef()
    print(f"Geometry: {geometry.GetGeometryName()}")

    # 获取属性
    attributes = feature.items()
    print(f"Attributes: {attributes}")

5. 写入栅格数据

gdal 库不仅可以读取数据，还可以创建和写入栅格数据。以下是一个简单的示例，展示如何创建一个新的 GeoTIFF 文件并写入数据。

from osgeo import gdal, osr
import numpy as np

# 创建一个新的栅格数据集
driver = gdal.GetDriverByName('GTiff')
dataset = driver.Create('output.tif', 512, 512, 1, gdal.GDT_Float32)

# 设置地理变换参数
geotransform = (0, 1, 0, 0, 0, -1)
dataset.SetGeoTransform(geotransform)

# 设置投影信息
srs = osr.SpatialReference()
srs.ImportFromEPSG(4326)  # WGS84
dataset.SetProjection(srs.ExportToWkt())

# 创建一个波段并写入数据
band = dataset.GetRasterBand(1)
data = np.random.rand(512, 512).astype(np.float32)
band.WriteArray(data)

# 关闭数据集
dataset = None

6. 写入矢量数据

同样，gdal 库也可以创建和写入矢量数据。以下是一个简单的示例，展示如何创建一个新的 Shapefile 并写入几何对象和属性。

from osgeo import ogr, osr

# 创建一个新的矢量数据集
driver = ogr.GetDriverByName('ESRI Shapefile')
dataSource = driver.CreateDataSource('output.shp')

# 创建一个图层
srs = osr.SpatialReference()
srs.ImportFromEPSG(4326)  # WGS84
layer = dataSource.CreateLayer('output', srs, ogr.wkbPoint)

# 创建一个字段
field_defn = ogr.FieldDefn('Name', ogr.OFTString)
field_defn.SetWidth(24)
layer.CreateField(field_defn)

# 创建一个要素并写入几何对象和属性
feature_defn = layer.GetLayerDefn()
feature = ogr.Feature(feature_defn)
feature.SetField('Name', 'Point 1')

# 创建几何对象
point = ogr.Geometry(ogr.wkbPoint)
point.AddPoint(10, 20)
feature.SetGeometry(point)

# 将要素写入图层
layer.CreateFeature(feature)

# 关闭数据源
dataSource = None

7. 常见操作

7.1 重采样（Resampling）

重采样是指改变栅格数据的分辨率，即改变每个像素的大小。常见的重采样方法包括最近邻插值、双线性插值、三次卷积插值等。

from osgeo import gdal, gdalconst

# 打开原始数据集
src_ds = gdal.Open('input.tif')

# 设置目标分辨率
x_res = src_ds.RasterXSize * 2  # 例如，将宽度增加一倍
y_res = src_ds.RasterYSize * 2  # 例如，将高度增加一倍

# 创建目标数据集
dst_ds = gdal.GetDriverByName('GTiff').Create('output_resampled.tif', x_res, y_res, src_ds.RasterCount, gdal.GDT_Float32)

# 设置地理变换参数
geotransform = list(src_ds.GetGeoTransform())
geotransform[1] /= 2  # 调整分辨率
geotransform[5] /= 2  # 调整分辨率
dst_ds.SetGeoTransform(geotransform)

# 设置投影信息
dst_ds.SetProjection(src_ds.GetProjection())

# 执行重采样
gdal.ReprojectImage(src_ds, dst_ds, src_ds.GetProjection(), src_ds.GetProjection(), gdalconst.GRA_Bilinear)

# 关闭数据集
src_ds = None
dst_ds = None

7.2 裁剪（Clipping）

裁剪是指根据几何对象（如多边形）裁剪栅格数据，只保留几何对象内的部分。

from osgeo import gdal, ogr, osr

# 打开原始数据集
src_ds = gdal.Open('input.tif')

# 打开几何对象数据集（例如，一个Shapefile）
shp_ds = ogr.Open('clip_polygon.shp')
layer = shp_ds.GetLayer()

# 创建目标数据集
dst_ds = gdal.GetDriverByName('GTiff').Create('output_clipped.tif', src_ds.RasterXSize, src_ds.RasterYSize, src_ds.RasterCount, gdal.GDT_Float32)

# 设置地理变换参数和投影信息
dst_ds.SetGeoTransform(src_ds.GetGeoTransform())
dst_ds.SetProjection(src_ds.GetProjection())

# 执行裁剪
gdal.Warp(dst_ds, src_ds, cutlineDSName='clip_polygon.shp', cropToCutline=True)

# 关闭数据集
src_ds = None
dst_ds = None
shp_ds = None

7.3 投影转换（Reprojection）

投影转换是指将数据从一个坐标系转换到另一个坐标系。常见的投影包括 WGS84、UTM 等。

from osgeo import gdal, osr

# 打开原始数据集
src_ds = gdal.Open('input.tif')

# 设置目标投影（例如，UTM Zone 18N）
dst_srs = osr.SpatialReference()
dst_srs.ImportFromEPSG(32618)  # UTM Zone 18N

# 创建目标数据集
dst_ds = gdal.GetDriverByName('GTiff').Create('output_reprojected.tif', src_ds.RasterXSize, src_ds.RasterYSize, src_ds.RasterCount, gdal.GDT_Float32)

# 执行投影转换
gdal.Warp(dst_ds, src_ds, dstSRS=dst_srs.ExportToWkt())

# 关闭数据集
src_ds = None
dst_ds = None

7.4 数据格式转换（Format Conversion）

数据格式转换是指将数据从一种格式转换为另一种格式。常见的栅格数据格式包括 GeoTIFF、JPEG、PNG 等。

from osgeo import gdal

# 打开原始数据集
src_ds = gdal.Open('input.tif')

# 创建目标数据集（例如，转换为JPEG格式）
dst_ds = gdal.GetDriverByName('JPEG').CreateCopy('output_converted.jpg', src_ds)

# 关闭数据集
src_ds = None
dst_ds = None

8. 总结

gdal 库是一个功能强大的工具，适用于处理地理空间数据。通过 Python 的 gdal 绑定，开发者可以在 Python 环境中轻松地读取、写入、转换和处理栅格和矢量数据。无论是进行简单的数据读取，还是复杂的图像处理和分析，gdal 库都能提供强大的支持。