Python + GDAL 栅格数据读取实战教程

在地理信息系统（GIS）领域，栅格数据是最常见的数据类型之一。无论是卫星影像、数字高程模型（DEM）还是气象数据，都以栅格格式存储。今天我们就用 Python + GDAL 来揭开栅格数据的神秘面纱。

一、什么是栅格数据

栅格数据是由规则网格（像素）组成的矩阵，每个像素都包含一个值，代表该位置的地理属性。

栅格 vs 矢量

常见栅格格式

• TIFF/GeoTIFF
  : 最常用的地理栅格格式，支持投影信息
• JPEG/JPEG2000
  : 压缩影像格式
• PNG
  : 无损压缩，支持透明通道
• NetCDF
  : 多维科学数据格式

二、GDAL 环境搭建

2.1 安装方法

Windows 系统：

# 方法1：使用预编译包（推荐）pip install GDAL-3.4.1-cp39-none-win_amd64.whl# 方法2：conda安装conda install -c conda-forge gdal

2.2 导入库

import osgeofrom osgeo import gdal, ogr, osrimport numpy as np

三、核心代码解析

3.1 读取栅格文件

defreadTif():# 文件路径    path = r"G:\data\q.tif"# 启用异常处理    gdal.UseExceptions()# 以只读模式打开文件    data = gdal.Open(path, gdal.GA_ReadOnly)# 检查文件是否成功打开if data isNone:print("无法打开文件！")return

3.2 获取基本信息

# 获取栅格尺寸    x_size = data.RasterXSize  # 列数    y_size = data.RasterYSize  # 行数print(f"栅格尺寸：{x_size} 列 × {y_size} 行")# 获取驱动信息    driver = data.GetDriver()print(f"文件格式：{driver.LongName}")# 获取描述信息    description = data.GetDescription()print(f"文件描述：{description}")

3.3 波段操作

# 获取波段（从1开始计数）    band = data.GetRasterBand(1)# 读取指定区域数据# ReadAsArray(起始列, 起始行, 宽度, 高度)    array = band.ReadAsArray(48, 25, 10, 10).astype(np.int32)print("10×10像素数据：")print(array)

3.4 地理变换参数详解

# 获取地理变换参数    transform = data.GetGeoTransform()# 解析六参数    origin_x = transform[0]  # 左上角X坐标    pixel_width = transform[1]  # 水平分辨率    rotation_x = transform[2]  # 行旋转（通常为0）    origin_y = transform[3]  # 左上角Y坐标    rotation_y = transform[4]  # 列旋转（通常为0）    pixel_height = transform[5]  # 垂直分辨率（通常为负）print(f"\n地理变换参数：")print(f"左上角经度：{origin_x}")print(f"左上角纬度：{origin_y}")print(f"水平分辨率：{pixel_width} 米")print(f"垂直分辨率：{abs(pixel_height)} 米")

坐标转换公式：

地理X坐标 = origin_x + 列号 × pixel_width + 行号 × rotation_x地理Y坐标 = origin_y + 列号 × rotation_y + 行号 × pixel_height

实际应用：

# 计算指定像素的地理坐标    col = 48174    row = 25437    lon = transform[0] + col * transform[1]    lat = transform[3] + row * transform[5]print(f"\n像素({col},{row})对应的地理坐标：")print(f"经度：{lon}")print(f"纬度：{lat}")

3.5 投影信息解析

# 获取投影信息（WKT格式）    projection = data.GetProjection()print(f"\n投影信息：{projection}")# 创建空间参考对象    reference = osr.SpatialReference()    reference.ImportFromWkt(projection)# 获取基准面信息    authority_name = reference.GetAuthorityName("DATUM")    authority_code = reference.GetAuthorityCode("DATUM")print(f"基准面：{authority_name}{authority_code}")# 获取EPSG编码    epsg_code = reference.GetAttrValue('AUTHORITY', 1)print(f"EPSG编码：{epsg_code}")# 获取单位    unit = reference.GetAttrValue('UNIT', 0)print(f"单位：{unit}")

四、实用技巧

4.1 多波段数据处理

# 获取波段数量band_count = data.RasterCountprint(f"波段数：{band_count}")# 遍历所有波段for i inrange(1, band_count + 1):    band = data.GetRasterBand(i)print(f"波段{i}数据类型：{gdal.GetDataTypeName(band.DataType)}")

4.2 数据类型转换

# 常见数据类型映射data_type_mapping = {    gdal.GDT_Byte: "8位无符号整数",    gdal.GDT_UInt16: "16位无符号整数",    gdal.GDT_Int16: "16位有符号整数",    gdal.GDT_UInt32: "32位无符号整数",    gdal.GDT_Int32: "32位有符号整数",    gdal.GDT_Float32: "32位浮点数",    gdal.GDT_Float64: "64位浮点数"}

4.3 保存处理结果

# 创建输出文件driver = gdal.GetDriverByName("GTiff")output_path = "output.tif"output_ds = driver.Create(    output_path,    x_size, y_size, 1,  # 宽度、高度、波段数    gdal.GDT_Float32  # 数据类型)# 设置地理变换和投影output_ds.SetGeoTransform(transform)output_ds.SetProjection(projection)# 写入数据output_band = output_ds.GetRasterBand(1)output_band.WriteArray(array)# 释放资源output_ds = None

五、完整代码

import osgeofrom osgeo import gdal, ogr, osrimport numpy as npdefreadTif():# 1. 打开文件    path = r"G:\项目_数据\黄泥河\q.tif"    gdal.UseExceptions()    data = gdal.Open(path, gdal.GA_ReadOnly)# 2. 获取基本信息    x_size = data.RasterXSize    y_size = data.RasterYSizeprint(f"栅格尺寸：{x_size} × {y_size}")# 3. 读取波段数据    band = data.GetRasterBand(1)    array = band.ReadAsArray(48, 25, 10, 10).astype(np.int32)print("采样数据：\n", array)# 4. 解析地理变换    transform = data.GetGeoTransform()print(f"左上角坐标：({transform[0]:.6f}, {transform[3]:.6f})")# 5. 解析投影信息    projection = data.GetProjection()    reference = osr.SpatialReference()    reference.ImportFromWkt(projection)    epsg_code = reference.GetAttrValue('AUTHORITY', 1)print(f"EPSG编码：{epsg_code}")if __name__ == '__main__':    readTif()

💡 小提示：在实际项目中，记得及时释放 GDAL Dataset 对象，避免内存泄漏哦！