GDAL Python 绑定使用指南

简介

GDAL (Geospatial Data Abstraction Library) 是一个开源的地理空间数据转换库。Python 绑定提供了简洁易用的 API，是地理空间数据处理的事实标准。

安装步骤

1. 使用 conda 安装（推荐，最省事）

# conda-forge 同时提供 libgdal 原生库和 Python 绑定，一步到位
conda install -c conda-forge gdal

说明: pip install gdal 只安装 Python 绑定层，需要系统里已存在匹配版本的 GDAL 原生库（libgdal）。在 Windows/macOS 上直接 pip install gdal 通常无法一步完成。下面按平台给出可行方案。

2. Windows 安装

方法一：使用 conda（最简便，推荐）

conda install -c conda-forge gdal

方法二：先装 GDAL 运行时，再装 Python 绑定

# 1. 从 https://trac.osgeo.org/gdal/reports/query?status=closed&component=Release+Notes
#    或 GitHub OSGeo/gdal releases 下载 GDAL 的 standalone MSI/ZIP（含 DLL）
# 2. 将 GDAL bin 目录加入 PATH
# 3. 安装与原生库版本一致的 Python 绑定
pip install gdal==<与系统 GDAL 相同的版本号>

注意：早期文档常指向 Christoph Gohlke 的预编译 wheel（gohlke pythonlibs），该站点已停止维护并下线，不要再使用。

方法三：使用 OSGeo4W

# 下载并运行 OSGeo4W 安装程序
# 选择安装 GDAL 和 python3-gdal

3. Linux (Ubuntu/Debian)

# 安装系统 GDAL 及开发头文件
sudo apt-get update
sudo apt-get install gdal-bin libgdal-dev python3-dev

# 安装与系统版本一致的 Python 绑定
pip install gdal==$(gdal-config --version)

4. macOS

# 使用 Homebrew
brew install gdal
pip install gdal

# 或使用 conda
conda install -c conda-forge gdal

5. 验证安装

from osgeo import gdal
print(gdal.__version__)
print(gdal.VersionInfo())

基本用法

读取栅格数据

from osgeo import gdal

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

if dataset is None:
    raise Exception("无法打开文件")

# 获取基本信息
print(f"尺寸: {dataset.RasterXSize} x {dataset.RasterYSize}")
print(f"波段数: {dataset.RasterCount}")

# 获取地理变换参数
geotransform = dataset.GetGeoTransform()
print(f"左上角 X: {geotransform[0]}")
print(f"左上角 Y: {geotransform[3]}")
print(f"像素宽度: {geotransform[1]}")
print(f"像素高度: {geotransform[5]}")

# 读取波段数据
band = dataset.GetRasterBand(1)
data = band.ReadAsArray()
print(f"数据类型: {data.dtype}")
print(f"最小值: {data.min()}, 最大值: {data.max()}")

# 关闭数据集
dataset = None

写入栅格数据

from osgeo import gdal
import numpy as np

# 创建新数据集
driver = gdal.GetDriverByName('GTiff')
width, height = 1000, 1000
dataset = driver.Create('output.tif', width, height, 1, gdal.GDT_Float32)

# 设置地理变换：minX, 像素宽, 0, maxY, 0, -像素高（最后一项必须为负，否则图像上下颠倒）
geotransform = (0, 1, 0, height, 0, -1)
dataset.SetGeoTransform(geotransform)

# 写入数据
band = dataset.GetRasterBand(1)
data = np.random.rand(height, width).astype(np.float32)
band.WriteArray(data)

# 设置 NoData 值
band.SetNoDataValue(-9999)

# 刷新缓存并关闭
dataset.FlushCache()
dataset = None

读取矢量数据

from osgeo import ogr

# 打开数据源
datasource = ogr.Open('input.shp')

if datasource is None:
    raise Exception("无法打开文件")

# 获取图层
layer = datasource.GetLayer(0)

# 遍历要素
for feature in layer:
    # 获取几何对象
    geometry = feature.GetGeometryRef()
    
    # 获取属性
    name = feature.GetField('name')
    
    # 处理几何
    if geometry.GetGeometryType() == ogr.wkbPoint:
        x, y, z = geometry.GetPoint()
        print(f"点坐标: ({x}, {y}, {z})")
    
    # 重要：释放要素
    feature = None

datasource = None

创建矢量数据

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('points', 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', 'Test Point')

# 创建几何
point = ogr.Geometry(ogr.wkbPoint)
point.AddPoint(116.4, 39.9)
feature.SetGeometry(point)

# 添加要素到图层
layer.CreateFeature(feature)

# 清理
feature = None
datasource = None

坐标转换

from osgeo import osr

def transform_coordinates(x, y, src_epsg, dst_epsg):
    """坐标转换"""
    # 定义源坐标系
    src_srs = osr.SpatialReference()
    src_srs.ImportFromEPSG(src_epsg)
    
    # 定义目标坐标系
    dst_srs = osr.SpatialReference()
    dst_srs.ImportFromEPSG(dst_epsg)
    
    # 创建坐标转换器
    transform = osr.CoordinateTransformation(src_srs, dst_srs)
    
    # 执行转换
    result = transform.TransformPoint(x, y)
    return result[:2]  # 返回 (x, y)

# 示例：WGS84 转 Web Mercator
lon, lat = 116.4, 39.9
x, y = transform_coordinates(lon, lat, 4326, 3857)
print(f"Web Mercator 坐标：({x}, {y})")

重投影栅格数据

from osgeo import gdal

def reproject_raster(input_file, output_file, dst_epsg=3857):
    """重投影栅格数据"""
    input_ds = gdal.Open(input_file)
    
    # 计算目标分辨率和范围
    warp_options = gdal.WarpOptions(
        format='GTiff',
        dstSRS=f'EPSG:{dst_epsg}',
        resampleAlg=gdal.GRA_Bilinear,
        creationOptions=['COMPRESS=LZW']
    )
    
    output_ds = gdal.Warp(output_file, input_ds, options=warp_options)
    
    output_ds = None
    input_ds = None

# 使用示例
reproject_raster('input.tif', 'output_reprojected.tif')

裁剪栅格数据

from osgeo import gdal

def clip_raster(input_file, output_file, clip_shp):
    """使用矢量边界裁剪栅格"""
    input_ds = gdal.Open(input_file)
    
    warp_options = gdal.WarpOptions(
        format='GTiff',
        cutlineDSName=clip_shp,
        cropToCutline=True,
        creationOptions=['COMPRESS=LZW']
    )
    
    output_ds = gdal.Warp(output_file, input_ds, options=warp_options)
    
    output_ds = None
    input_ds = None

# 使用示例
clip_raster('input.tif', 'clipped.tif', 'boundary.shp')

合并栅格数据

from osgeo import gdal

def merge_rasters(input_files, output_file):
    """合并多个栅格文件"""
    merge_options = gdal.WarpOptions(
        format='GTiff',
        resampleAlg=gdal.GRA_Bilinear,
        creationOptions=['COMPRESS=LZW', 'TILED=YES']
    )
    
    output_ds = gdal.Warp(output_file, input_files, options=merge_options)
    output_ds = None

# 使用示例
input_files = ['tile1.tif', 'tile2.tif', 'tile3.tif']
merge_rasters(input_files, 'mosaic.tif')

生成金字塔

from osgeo import gdal

def build_pyramids(raster_file):
    """为栅格文件生成金字塔（概览图 overviews）"""
    # 需要以可写方式打开
    dataset = gdal.Open(raster_file, gdal.GA_Update)

    # resampling 算法可选 NEAREST / AVERAGE / CUBIC 等；decimation 因子按需求调整
    dataset.BuildOverviews('NEAREST', [2, 4, 8, 16, 32])

    dataset.FlushCache()
    dataset = None

# 使用示例
build_pyramids('large_image.tif')

常用工具函数

获取 GDAL 支持的文件格式

from osgeo import gdal

def list_supported_formats():
    """列出所有支持的栅格格式"""
    for i in range(gdal.GetDriverCount()):
        driver = gdal.GetDriver(i)
        print(f"{driver.ShortName}: {driver.LongName}")

list_supported_formats()

读取栅格统计信息

from osgeo import gdal

def get_raster_stats(raster_file):
    """获取栅格统计信息"""
    dataset = gdal.Open(raster_file)
    band = dataset.GetRasterBand(1)
    
    stats = band.GetStatistics(True, True)
    print(f"最小值: {stats[0]}")
    print(f"最大值: {stats[1]}")
    print(f"平均值: {stats[2]}")
    print(f"标准差: {stats[3]}")
    
    dataset = None

get_raster_stats('input.tif')

批量格式转换

from osgeo import gdal
import glob

def batch_convert(input_pattern, output_format='GTiff'):
    """批量转换栅格格式"""
    files = glob.glob(input_pattern)
    
    for input_file in files:
        output_file = input_file.rsplit('.', 1)[0] + '.tif'
        
        input_ds = gdal.Open(input_file)
        translate_options = gdal.TranslateOptions(format=output_format)
        
        output_ds = gdal.Translate(output_file, input_ds, options=translate_options)
        
        output_ds = None
        input_ds = None
        
        print(f"已转换：{input_file} -> {output_file}")

# 使用示例
batch_convert('*.jp2')

性能优化建议

1. 分块读取大文件

def read_in_chunks(raster_file, chunk_size=1024):
    """分块读取大栅格文件"""
    dataset = gdal.Open(raster_file)
    band = dataset.GetRasterBand(1)
    
    width = dataset.RasterXSize
    height = dataset.RasterYSize
    
    for y in range(0, height, chunk_size):
        for x in range(0, width, chunk_size):
            xsize = min(chunk_size, width - x)
            ysize = min(chunk_size, height - y)
            
            chunk = band.ReadAsArray(x, y, xsize, ysize)
            # 处理数据块
            
    dataset = None

2. 使用压缩选项

# 创建压缩的 GeoTIFF
creation_options = [
    'COMPRESS=LZW',      # LZW 压缩
    'TILED=YES',         # 启用分块
    'BIGTIFF=YES',       # 支持大文件
    'NUM_THREADS=ALL_CPUS'  # 多线程
]

3. 使用内存优化模式

# 使用虚拟数据集减少内存占用
vrt = gdal.BuildVRT('combined.vrt', ['file1.tif', 'file2.tif'])
data = vrt.ReadAsArray()
vrt = None

常见问题解决

1. 导入错误

错误: ModuleNotFoundError: No module named 'osgeo'

解决方案:

# 确保正确安装
pip install gdal

# 或者使用 conda
conda install -c conda-forge gdal

2. 版本不匹配 / DLL 加载失败

错误: ImportError: DLL load failed 或 Cannot link to GDAL library

原因: Python 绑定版本与系统原生 GDAL（libgdal）版本不一致，或原生库不在搜索路径上。

解决方案:

# Linux/macOS：查看系统 GDAL 版本
gdal-config --version
# 然后安装相同版本的绑定
pip install gdal==<上面的版本号>

# Windows：确认 GDAL 的 bin 目录（含 gdal*.dll）已加入 PATH，
# 或在代码开头显式添加后再导入
import os
os.add_dll_directory(r"C:\Program Files\GDAL\bin")
from osgeo import gdal

3. 中文路径问题

Windows 下部分旧版 GDAL 对非 ASCII 路径支持不佳。可行做法：

import os
# 1. 尽量使用绝对路径，避免相对路径中的中文
path = os.path.abspath("中文路径/文件.tif")

# 2. 若仍报错，可将文件临时复制到纯英文路径再打开
import shutil, tempfile
tmp = os.path.join(tempfile.gettempdir(), "data.tif")
shutil.copy(path, tmp)
dataset = gdal.Open(tmp)

说明：str.encode('utf-8').decode('utf-8') 是恒等变换，对解决中文路径无效，请勿使用。

4. 内存不足

# 使用 VRT 虚拟数据集
vrt_dataset = gdal.BuildVRT('', source_datasets)

# 或使用低内存读取模式
band.ReadAsArray(xoff, yoff, xsize, ysize, buf_xsize, buf_ysize)

参考资源

• GDAL 官方文档
• GDAL Python API
• GDAL GitHub
• Conda-forge GDAL