第03章 - GDAL核心架构与数据模型
第03章：GDAL核心架构与数据模型
3.1 GDAL 整体架构
3.1.1 架构概览
GDAL 采用分层架构设计，将复杂的地理空间数据处理分解为多个清晰的层次：
┌─────────────────────────────────────────────────────────────────────────┐
│ 应用程序接口层 │
│ C/C++ API │ Python API │ Java API │ C# API │ 命令行工具 │
└─────────────────────────────────────────────────────────────────────────┘
│
▼
┌─────────────────────────────────────────────────────────────────────────┐
│ 核心抽象层 │
│ ┌───────────────┐ ┌───────────────┐ ┌───────────────┐ │
│ │ GDALDataset │ │ OGRDataSource│ │OGRSpatialRef │ │
│ │ (数据集) │ │ (数据源) │ │(空间参考) │ │
│ └───────┬───────┘ └───────┬───────┘ └───────────────┘ │
│ │ │ │
│ ┌───────▼───────┐ ┌───────▼───────┐ │
│ │GDALRasterBand │ │ OGRLayer │ │
│ │ (栅格波段) │ │ (图层) │ │
│ └───────────────┘ └───────┬───────┘ │
│ │ │
│ ┌───────▼───────┐ │
│ │ OGRFeature │ │
│ │ (要素) │ │
│ └───────────────┘ │
└─────────────────────────────────────────────────────────────────────────┘
│
▼
┌─────────────────────────────────────────────────────────────────────────┐
│ 驱动层 │
│ ┌─────────────┐ ┌─────────────┐ ┌─────────────┐ ┌─────────────┐ │
│ │ GeoTIFF │ │ Shapefile │ │ PostGIS │ │ NetCDF │ │
│ │ Driver │ │ Driver │ │ Driver │ │ Driver │ │
│ └─────────────┘ └─────────────┘ └─────────────┘ └─────────────┘ │
│ ...更多驱动... │
└─────────────────────────────────────────────────────────────────────────┘
│
▼
┌─────────────────────────────────────────────────────────────────────────┐
│ 数据源层 │
│ 文件系统 │ 数据库 │ 网络服务 │ 内存 │ 云存储 │
└─────────────────────────────────────────────────────────────────────────┘
3.1.2 核心设计原则
GDAL 的架构设计遵循以下原则：
| 原则 | 描述 | 实现方式 |
|---|---|---|
| 抽象化 | 隐藏底层格式差异 | 统一的 Dataset/Layer 接口 |
| 可扩展 | 支持添加新格式 | 驱动插件机制 |
| 高性能 | 处理大数据集 | 分块读取、缓存、延迟加载 |
| 跨平台 | 多操作系统支持 | 标准 C/C++ 实现 |
| 向后兼容 | API 稳定性 | 版本控制、弃用策略 |
3.1.3 模块组成
GDAL 库
├── gcore/ # 栅格核心模块
│ ├── gdal.h # 主头文件
│ ├── gdal_priv.h # 私有头文件
│ ├── gdaldataset.cpp # 数据集实现
│ └── gdalrasterband.cpp # 波段实现
│
├── ogr/ # 矢量核心模块
│ ├── ogr_feature.h # 要素定义
│ ├── ogr_geometry.h # 几何定义
│ ├── ogrsf_frmts/ # 矢量格式驱动
│ └── ogrutils.cpp # 工具函数
│
├── frmts/ # 栅格格式驱动
│ ├── gtiff/ # GeoTIFF 驱动
│ ├── hdf5/ # HDF5 驱动
│ ├── netcdf/ # NetCDF 驱动
│ └── ...
│
├── alg/ # 算法模块
│ ├── gdalwarpkernel.cpp # Warp 核心
│ ├── gdal_rpc.cpp # RPC 校正
│ └── contour.cpp # 等值线生成
│
└── apps/ # 命令行应用
├── gdalinfo.cpp
├── gdal_translate.cpp
└── ...
3.2 驱动模型
3.2.1 驱动概念
GDAL 驱动是连接抽象层和具体数据格式的桥梁。每个驱动负责读写特定的数据格式。
┌─────────────────────────────────────────────────────────────┐
│ GDALDriverManager │
│ (驱动管理器) │
├─────────────────────────────────────────────────────────────┤
│ │
│ ┌──────────┐ ┌──────────┐ ┌──────────┐ ┌──────────┐ │
│ │ GeoTIFF │ │ PNG │ │ HDF5 │ │ Shapefile│ │
│ │ Driver │ │ Driver │ │ Driver │ │ Driver │ │
│ └────┬─────┘ └────┬─────┘ └────┬─────┘ └────┬─────┘ │
│ │ │ │ │ │
│ ▼ ▼ ▼ ▼ │
│ GDALDriver GDALDriver GDALDriver OGRSFDriver │
│ │
└─────────────────────────────────────────────────────────────┘
3.2.2 驱动类型
| 类型 | 描述 | 示例 |
|---|---|---|
| 栅格驱动 | 处理栅格/影像数据 | GeoTIFF, PNG, HDF5 |
| 矢量驱动 | 处理矢量数据 | Shapefile, GeoJSON, PostgreSQL |
| 混合驱动 | 同时支持栅格和矢量 | GeoPackage, PDF |
| 虚拟驱动 | 不对应实际格式 | VRT, MEM |
| 网络驱动 | 访问网络资源 | WMS, WCS, WFS |
3.2.3 驱动能力标识
每个驱动声明其支持的能力：
from osgeo import gdal
# 获取驱动
driver = gdal.GetDriverByName('GTiff')
# 检查驱动能力
metadata = driver.GetMetadata()
# 常见能力标识
capabilities = {
'DCAP_CREATE': '支持创建新数据集',
'DCAP_CREATECOPY': '支持复制创建',
'DCAP_VIRTUALIO': '支持虚拟文件系统',
'DCAP_RASTER': '支持栅格数据',
'DCAP_VECTOR': '支持矢量数据',
'DMD_CREATIONDATATYPES': '支持的数据类型',
'DMD_CREATIONOPTIONLIST': '创建选项列表',
}
for cap, desc in capabilities.items():
value = metadata.get(cap, 'N/A')
print(f"{desc}: {value}")
3.2.4 驱动注册与管理
from osgeo import gdal, ogr
# 注册所有驱动（通常自动完成）
gdal.AllRegister()
# 获取驱动数量
raster_count = gdal.GetDriverCount()
vector_count = ogr.GetDriverCount()
print(f"栅格驱动数: {raster_count}")
print(f"矢量驱动数: {vector_count}")
# 按名称获取驱动
gtiff_driver = gdal.GetDriverByName('GTiff')
shp_driver = ogr.GetDriverByName('ESRI Shapefile')
# 按索引获取驱动
for i in range(gdal.GetDriverCount()):
driver = gdal.GetDriver(i)
print(f"{driver.ShortName}: {driver.LongName}")
# 注销驱动（很少使用）
gdal.GetDriverByName('GTiff').Deregister()
# 重新注册
gdal.GetDriverByName('GTiff').Register()
3.2.5 常用驱动详解
GeoTIFF 驱动
from osgeo import gdal
driver = gdal.GetDriverByName('GTiff')
# 创建选项
creation_options = [
'COMPRESS=LZW', # 压缩方式
'TILED=YES', # 分块存储
'BLOCKXSIZE=256', # 块宽度
'BLOCKYSIZE=256', # 块高度
'BIGTIFF=IF_SAFER', # 大文件支持
'PREDICTOR=2', # 预测器（整数）
'PHOTOMETRIC=RGB', # 色彩解释
]
# 创建数据集
ds = driver.Create(
'output.tif',
width=1000,
height=1000,
bands=3,
eType=gdal.GDT_Byte,
options=creation_options
)
Shapefile 驱动
from osgeo import ogr
driver = ogr.GetDriverByName('ESRI Shapefile')
# 创建数据源
ds = driver.CreateDataSource('output.shp')
# Shapefile 限制
limitations = {
'字段名长度': '最多10个字符',
'字段类型': '有限的类型支持',
'文件大小': '最大2GB',
'几何类型': '每个文件只能有一种',
}
3.3 栅格数据模型
3.3.1 数据集（GDALDataset）
GDALDataset 是栅格数据的顶层容器：
┌─────────────────────────────────────────────────────────────┐
│ GDALDataset │
├─────────────────────────────────────────────────────────────┤
│ 属性： │
│ ├── RasterXSize : 宽度（像素） │
│ ├── RasterYSize : 高度（像素） │
│ ├── RasterCount : 波段数 │
│ ├── GeoTransform : 地理变换参数 │
│ ├── Projection : 投影信息 │
│ └── Metadata : 元数据 │
├─────────────────────────────────────────────────────────────┤
│ 方法： │
│ ├── GetRasterBand(n) : 获取第n个波段 │
│ ├── ReadAsArray() : 读取所有波段为数组 │
│ ├── GetGeoTransform() : 获取地理变换 │
│ ├── SetGeoTransform() : 设置地理变换 │
│ ├── GetProjection() : 获取投影 │
│ ├── SetProjection() : 设置投影 │
│ └── FlushCache() : 刷新缓存到磁盘 │
└─────────────────────────────────────────────────────────────┘
│
│ 包含多个
▼
┌─────────────────────────────────────────────────────────────┐
│ GDALRasterBand │
│ (波段1) │
└─────────────────────────────────────────────────────────────┘
3.3.2 地理变换（GeoTransform）
地理变换是一个包含6个参数的数组，定义了像素坐标到地理坐标的转换关系：
GeoTransform = [originX, pixelWidth, rotation1, originY, rotation2, pixelHeight]
其中：
- originX : 左上角X坐标（地理坐标）
- pixelWidth : 像素宽度（通常为正值）
- rotation1 : X方向旋转（通常为0）
- originY : 左上角Y坐标（地理坐标）
- rotation2 : Y方向旋转（通常为0）
- pixelHeight : 像素高度（通常为负值）
坐标转换公式：
地理X = originX + 像素列 * pixelWidth + 像素行 * rotation1
地理Y = originY + 像素列 * rotation2 + 像素行 * pixelHeight
Python 示例：
from osgeo import gdal
ds = gdal.Open('example.tif')
gt = ds.GetGeoTransform()
print(f"左上角X: {gt[0]}")
print(f"像素宽度: {gt[1]}")
print(f"X旋转: {gt[2]}")
print(f"左上角Y: {gt[3]}")
print(f"Y旋转: {gt[4]}")
print(f"像素高度: {gt[5]}")
# 像素坐标转地理坐标
def pixel_to_geo(gt, col, row):
x = gt[0] + col * gt[1] + row * gt[2]
y = gt[3] + col * gt[4] + row * gt[5]
return x, y
# 地理坐标转像素坐标
def geo_to_pixel(gt, x, y):
# 简化版（无旋转时）
col = int((x - gt[0]) / gt[1])
row = int((y - gt[3]) / gt[5])
return col, row
# 测试
geo_x, geo_y = pixel_to_geo(gt, 100, 200)
print(f"像素(100, 200) -> 地理({geo_x}, {geo_y})")
3.3.3 波段（GDALRasterBand）
波段是栅格数据的基本组成单元：
┌─────────────────────────────────────────────────────────────┐
│ GDALRasterBand │
├─────────────────────────────────────────────────────────────┤
│ 属性： │
│ ├── XSize : 波段宽度 │
│ ├── YSize : 波段高度 │
│ ├── DataType : 数据类型 │
│ ├── BlockXSize : 块宽度 │
│ ├── BlockYSize : 块高度 │
│ ├── NoDataValue : 无效值 │
│ ├── Minimum : 最小值 │
│ ├── Maximum : 最大值 │
│ └── ColorInterp : 颜色解释 │
├─────────────────────────────────────────────────────────────┤
│ 方法： │
│ ├── ReadAsArray() : 读取为NumPy数组 │
│ ├── WriteArray() : 写入NumPy数组 │
│ ├── ReadRaster() : 读取原始字节 │
│ ├── WriteRaster() : 写入原始字节 │
│ ├── GetStatistics() : 获取统计信息 │
│ ├── ComputeStatistics() : 计算统计信息 │
│ ├── GetHistogram() : 获取直方图 │
│ └── GetOverview() : 获取金字塔层 │
└─────────────────────────────────────────────────────────────┘
3.3.4 数据类型
GDAL 支持多种数据类型：
| 常量 | 描述 | 范围 | 字节数 |
|---|---|---|---|
| GDT_Byte | 无符号8位整数 | 0-255 | 1 |
| GDT_UInt16 | 无符号16位整数 | 0-65535 | 2 |
| GDT_Int16 | 有符号16位整数 | -32768 - 32767 | 2 |
| GDT_UInt32 | 无符号32位整数 | 0 - 4294967295 | 4 |
| GDT_Int32 | 有符号32位整数 | -2147483648 - 2147483647 | 4 |
| GDT_Float32 | 32位浮点数 | ±3.4e38 | 4 |
| GDT_Float64 | 64位浮点数 | ±1.7e308 | 8 |
| GDT_CInt16 | 复数（16位整数） | - | 4 |
| GDT_CInt32 | 复数（32位整数） | - | 8 |
| GDT_CFloat32 | 复数（32位浮点） | - | 8 |
| GDT_CFloat64 | 复数（64位浮点） | - | 16 |
from osgeo import gdal
import numpy as np
# 数据类型与 NumPy 类型对应
gdal_to_numpy = {
gdal.GDT_Byte: np.uint8,
gdal.GDT_UInt16: np.uint16,
gdal.GDT_Int16: np.int16,
gdal.GDT_UInt32: np.uint32,
gdal.GDT_Int32: np.int32,
gdal.GDT_Float32: np.float32,
gdal.GDT_Float64: np.float64,
}
# 获取数据类型名称
ds = gdal.Open('example.tif')
band = ds.GetRasterBand(1)
dtype = band.DataType
dtype_name = gdal.GetDataTypeName(dtype)
print(f"数据类型: {dtype_name}")
3.3.5 分块读取
对于大文件，应使用分块读取以节省内存：
from osgeo import gdal
import numpy as np
def read_raster_in_blocks(filepath):
"""分块读取栅格数据"""
ds = gdal.Open(filepath)
band = ds.GetRasterBand(1)
# 获取块大小
block_x, block_y = band.GetBlockSize()
x_size = band.XSize
y_size = band.YSize
print(f"影像大小: {x_size} x {y_size}")
print(f"块大小: {block_x} x {block_y}")
# 分块读取
for y in range(0, y_size, block_y):
for x in range(0, x_size, block_x):
# 计算实际读取大小
cols = min(block_x, x_size - x)
rows = min(block_y, y_size - y)
# 读取块数据
data = band.ReadAsArray(x, y, cols, rows)
# 处理数据...
print(f"读取块: ({x}, {y}), 大小: ({cols}, {rows})")
ds = None
# 自定义块大小读取
def read_with_custom_block(filepath, block_size=512):
"""使用自定义块大小读取"""
ds = gdal.Open(filepath)
band = ds.GetRasterBand(1)
x_size = band.XSize
y_size = band.YSize
for y in range(0, y_size, block_size):
for x in range(0, x_size, block_size):
cols = min(block_size, x_size - x)
rows = min(block_size, y_size - y)
data = band.ReadAsArray(x, y, cols, rows)
yield x, y, data
ds = None
3.4 矢量数据模型
3.4.1 数据源（DataSource）
在 GDAL 3.x 中，矢量数据源也使用 GDALDataset 类：
┌─────────────────────────────────────────────────────────────┐
│ GDALDataset │
│ (矢量数据源) │
├─────────────────────────────────────────────────────────────┤
│ 属性： │
│ ├── GetLayerCount() : 图层数量 │
│ ├── GetDescription() : 数据源描述 │
│ └── GetDriver() : 获取驱动 │
├─────────────────────────────────────────────────────────────┤
│ 方法： │
│ ├── GetLayer(n) : 按索引获取图层 │
│ ├── GetLayerByName() : 按名称获取图层 │
│ ├── CreateLayer() : 创建新图层 │
│ ├── DeleteLayer() : 删除图层 │
│ ├── ExecuteSQL() : 执行SQL查询 │
│ └── StartTransaction() : 开始事务 │
└─────────────────────────────────────────────────────────────┘
│
│ 包含多个
▼
┌─────────────────────────────────────────────────────────────┐
│ OGRLayer │
│ (图层) │
└─────────────────────────────────────────────────────────────┘
3.4.2 图层（OGRLayer）
图层是矢量要素的容器：
┌─────────────────────────────────────────────────────────────┐
│ OGRLayer │
├─────────────────────────────────────────────────────────────┤
│ 属性： │
│ ├── GetName() : 图层名称 │
│ ├── GetGeomType() : 几何类型 │
│ ├── GetFeatureCount() : 要素数量 │
│ ├── GetExtent() : 范围 │
│ ├── GetSpatialRef() : 空间参考 │
│ └── GetLayerDefn() : 图层定义 │
├─────────────────────────────────────────────────────────────┤
│ 遍历方法： │
│ ├── GetNextFeature() : 获取下一个要素 │
│ ├── ResetReading() : 重置读取位置 │
│ ├── GetFeature(fid) : 按FID获取要素 │
│ └── SetNextByIndex(n) : 设置读取位置 │
├─────────────────────────────────────────────────────────────┤
│ 筛选方法： │
│ ├── SetSpatialFilter() : 设置空间过滤器 │
│ ├── SetAttributeFilter(): 设置属性过滤器 │
│ └── SetSpatialFilterRect(): 设置矩形过滤器 │
├─────────────────────────────────────────────────────────────┤
│ 编辑方法： │
│ ├── CreateFeature() : 创建要素 │
│ ├── SetFeature() : 更新要素 │
│ └── DeleteFeature() : 删除要素 │
└─────────────────────────────────────────────────────────────┘
3.4.3 要素（OGRFeature）
要素是矢量数据的基本单元，包含几何和属性：
┌─────────────────────────────────────────────────────────────┐
│ OGRFeature │
├─────────────────────────────────────────────────────────────┤
│ │
│ ┌─────────────────────────────────────────────────────┐ │
│ │ OGRGeometry │ │
│ │ (几何对象) │ │
│ │ │ │
│ │ Point / LineString / Polygon / MultiPoint / ... │ │
│ └─────────────────────────────────────────────────────┘ │
│ │
│ ┌─────────────────────────────────────────────────────┐ │
│ │ 属性字段 │ │
│ │ │ │
│ │ 字段1: 值1 │ │
│ │ 字段2: 值2 │ │
│ │ 字段3: 值3 │ │
│ │ ... │ │
│ └─────────────────────────────────────────────────────┘ │
│ │
│ 其他属性： │
│ ├── FID : 要素ID │
│ └── StyleString : 样式字符串 │
│ │
└─────────────────────────────────────────────────────────────┘
3.4.4 几何类型
OGR 支持 OGC Simple Features 规范定义的所有几何类型：
| 类型常量 | 描述 | WKT 示例 |
|---|---|---|
| wkbPoint | 点 | POINT (30 10) |
| wkbLineString | 线 | LINESTRING (30 10, 10 30, 40 40) |
| wkbPolygon | 多边形 | POLYGON ((30 10, 40 40, 20 40, 10 20, 30 10)) |
| wkbMultiPoint | 多点 | MULTIPOINT ((10 40), (40 30)) |
| wkbMultiLineString | 多线 | MULTILINESTRING ((10 10, 20 20)) |
| wkbMultiPolygon | 多多边形 | MULTIPOLYGON (((30 20, 45 40, 10 40, 30 20))) |
| wkbGeometryCollection | 几何集合 | GEOMETRYCOLLECTION (POINT(4 6), LINESTRING(4 6, 7 10)) |
带 Z/M 值的几何类型：
| 类型 | 描述 |
|---|---|
| wkbPoint25D | 带高程的点 |
| wkbPointM | 带测量值的点 |
| wkbPointZM | 带高程和测量值的点 |
from osgeo import ogr
# 创建点几何
point = ogr.Geometry(ogr.wkbPoint)
point.AddPoint(116.4, 39.9)
# 创建线几何
line = ogr.Geometry(ogr.wkbLineString)
line.AddPoint(116.4, 39.9)
line.AddPoint(121.5, 31.2)
line.AddPoint(113.3, 23.1)
# 创建多边形
ring = ogr.Geometry(ogr.wkbLinearRing)
ring.AddPoint(116.0, 39.0)
ring.AddPoint(117.0, 39.0)
ring.AddPoint(117.0, 40.0)
ring.AddPoint(116.0, 40.0)
ring.AddPoint(116.0, 39.0) # 闭合
polygon = ogr.Geometry(ogr.wkbPolygon)
polygon.AddGeometry(ring)
# 从 WKT 创建
wkt = "POINT (116.4 39.9)"
point_from_wkt = ogr.CreateGeometryFromWkt(wkt)
# 从 GeoJSON 创建
geojson = '{"type":"Point","coordinates":[116.4,39.9]}'
point_from_json = ogr.CreateGeometryFromJson(geojson)
3.4.5 字段定义
属性字段由 OGRFieldDefn 定义：
| 字段类型 | 描述 | Python 对应 |
|---|---|---|
| OFTInteger | 32位整数 | int |
| OFTInteger64 | 64位整数 | int |
| OFTReal | 浮点数 | float |
| OFTString | 字符串 | str |
| OFTDate | 日期 | datetime.date |
| OFTTime | 时间 | datetime.time |
| OFTDateTime | 日期时间 | datetime.datetime |
| OFTBinary | 二进制 | bytes |
| OFTIntegerList | 整数数组 | list |
| OFTRealList | 浮点数组 | list |
| OFTStringList | 字符串数组 | list |
from osgeo import ogr
# 创建字段定义
field_defn = ogr.FieldDefn('name', ogr.OFTString)
field_defn.SetWidth(100) # 字符串宽度
# 整数字段
int_field = ogr.FieldDefn('population', ogr.OFTInteger64)
# 浮点字段
float_field = ogr.FieldDefn('area', ogr.OFTReal)
float_field.SetWidth(12)
float_field.SetPrecision(2) # 小数位数
# 日期字段
date_field = ogr.FieldDefn('create_date', ogr.OFTDate)
3.5 空间参考系统
3.5.1 OGRSpatialReference
空间参考系统（Spatial Reference System, SRS）定义了坐标的含义：
┌─────────────────────────────────────────────────────────────┐
│ OGRSpatialReference │
├─────────────────────────────────────────────────────────────┤
│ 定义来源： │
│ ├── ImportFromEPSG() : 从EPSG代码导入 │
│ ├── ImportFromWkt() : 从WKT导入 │
│ ├── ImportFromProj4() : 从PROJ字符串导入 │
│ └── ImportFromESRI() : 从ESRI WKT导入 │
├─────────────────────────────────────────────────────────────┤
│ 导出格式： │
│ ├── ExportToWkt() : 导出为WKT │
│ ├── ExportToPrettyWkt(): 导出为格式化WKT │
│ ├── ExportToProj4() : 导出为PROJ字符串 │
│ └── ExportToXML() : 导出为XML │
├─────────────────────────────────────────────────────────────┤
│ 属性查询： │
│ ├── GetName() : 获取名称 │
│ ├── GetAuthorityCode() : 获取权威代码(如EPSG) │
│ ├── IsGeographic() : 是否为地理坐标系 │
│ ├── IsProjected() : 是否为投影坐标系 │
│ └── GetLinearUnits() : 获取线性单位 │
└─────────────────────────────────────────────────────────────┘
3.5.2 常用坐标系
from osgeo import osr
# WGS 84 地理坐标系
wgs84 = osr.SpatialReference()
wgs84.ImportFromEPSG(4326)
print(f"WGS84: {wgs84.GetName()}")
# Web Mercator 投影
web_mercator = osr.SpatialReference()
web_mercator.ImportFromEPSG(3857)
print(f"Web Mercator: {web_mercator.GetName()}")
# 中国 CGCS2000
cgcs2000 = osr.SpatialReference()
cgcs2000.ImportFromEPSG(4490)
print(f"CGCS2000: {cgcs2000.GetName()}")
# UTM 投影（北半球50带）
utm50n = osr.SpatialReference()
utm50n.ImportFromEPSG(32650)
print(f"UTM 50N: {utm50n.GetName()}")
# 检查坐标系类型
print(f"\nWGS84 是地理坐标系: {wgs84.IsGeographic()}")
print(f"Web Mercator 是投影坐标系: {web_mercator.IsProjected()}")
3.5.3 坐标转换
from osgeo import osr, ogr
def transform_point(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[0], result[1]
def transform_geometry(geom, 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)
# 克隆几何对象并转换
new_geom = geom.Clone()
new_geom.Transform(transform)
return new_geom
# 示例：WGS84 转 Web Mercator
x, y = transform_point(116.4, 39.9, 4326, 3857)
print(f"WGS84 (116.4, 39.9) -> Web Mercator ({x:.2f}, {y:.2f})")
# 几何对象转换
point = ogr.Geometry(ogr.wkbPoint)
point.AddPoint(116.4, 39.9)
transformed_point = transform_geometry(point, 4326, 3857)
print(f"转换后: {transformed_point.ExportToWkt()}")
3.6 虚拟文件系统
3.6.1 VFS 概念
GDAL 虚拟文件系统（Virtual File System）允许访问非本地文件：
┌─────────────────────────────────────────────────────────────┐
│ GDAL 虚拟文件系统 │
├─────────────────────────────────────────────────────────────┤
│ │
│ /vsimem/ 内存文件系统 │
│ /vsizip/ ZIP 压缩包 │
│ /vsigzip/ GZIP 压缩文件 │
│ /vsitar/ TAR 归档文件 │
│ /vsicurl/ HTTP/HTTPS URL │
│ /vsis3/ Amazon S3 │
│ /vsigs/ Google Cloud Storage │
│ /vsiaz/ Azure Blob Storage │
│ /vsiadls/ Azure Data Lake Storage │
│ /vsioss/ 阿里云 OSS │
│ /vsiswift/ OpenStack Swift │
│ /vsihdfs/ Hadoop HDFS │
│ │
└─────────────────────────────────────────────────────────────┘
3.6.2 内存文件系统
from osgeo import gdal
import numpy as np
# 创建内存文件
mem_path = '/vsimem/temp_file.tif'
# 创建驱动
driver = gdal.GetDriverByName('GTiff')
# 创建数据集
ds = driver.Create(mem_path, 100, 100, 1, gdal.GDT_Byte)
band = ds.GetRasterBand(1)
# 写入数据
data = np.random.randint(0, 255, (100, 100), dtype=np.uint8)
band.WriteArray(data)
ds.FlushCache()
ds = None
# 读取内存文件
ds = gdal.Open(mem_path)
read_data = ds.GetRasterBand(1).ReadAsArray()
print(f"读取数据形状: {read_data.shape}")
ds = None
# 删除内存文件
gdal.Unlink(mem_path)
3.6.3 访问压缩文件
from osgeo import gdal, ogr
# 读取 ZIP 文件中的 Shapefile
shp_in_zip = '/vsizip/data.zip/layer.shp'
ds = ogr.Open(shp_in_zip)
if ds:
layer = ds.GetLayer()
print(f"要素数: {layer.GetFeatureCount()}")
ds = None
# 读取 GZIP 压缩的 GeoTIFF
gzip_tif = '/vsigzip/image.tif.gz'
ds = gdal.Open(gzip_tif)
if ds:
print(f"影像大小: {ds.RasterXSize} x {ds.RasterYSize}")
ds = None
3.6.4 访问网络资源
from osgeo import gdal
# 配置网络访问
gdal.SetConfigOption('GDAL_HTTP_TIMEOUT', '30')
gdal.SetConfigOption('GDAL_HTTP_MAX_RETRY', '3')
gdal.SetConfigOption('CPL_VSIL_CURL_ALLOWED_EXTENSIONS', '.tif,.tiff,.vrt')
# 访问 HTTP URL
url = 'https://example.com/data.tif'
vsicurl_path = f'/vsicurl/{url}'
ds = gdal.Open(vsicurl_path)
if ds:
print(f"远程影像大小: {ds.RasterXSize} x {ds.RasterYSize}")
ds = None
# 访问 S3 存储
gdal.SetConfigOption('AWS_ACCESS_KEY_ID', 'your_key')
gdal.SetConfigOption('AWS_SECRET_ACCESS_KEY', 'your_secret')
gdal.SetConfigOption('AWS_REGION', 'us-east-1')
s3_path = '/vsis3/bucket-name/path/to/file.tif'
ds = gdal.Open(s3_path)
3.7 元数据系统
3.7.1 元数据域
GDAL 支持多个元数据域：
| 域名 | 描述 |
|---|---|
| "" (默认) | 默认元数据域 |
| IMAGE_STRUCTURE | 影像结构信息 |
| SUBDATASETS | 子数据集列表 |
| GEOLOCATION | 地理定位信息 |
| RPC | 有理多项式系数 |
| DERIVED_SUBDATASETS | 派生子数据集 |
| xml:* | XML 格式元数据 |
3.7.2 读写元数据
from osgeo import gdal
ds = gdal.Open('example.tif')
# 获取默认域元数据
metadata = ds.GetMetadata()
for key, value in metadata.items():
print(f"{key}: {value}")
# 获取特定域元数据
image_structure = ds.GetMetadata('IMAGE_STRUCTURE')
print(f"\n影像结构:")
for key, value in image_structure.items():
print(f" {key}: {value}")
# 获取单个元数据项
value = ds.GetMetadataItem('AREA_OR_POINT')
print(f"\nAREA_OR_POINT: {value}")
# 波段元数据
band = ds.GetRasterBand(1)
band_meta = band.GetMetadata()
print(f"\n波段元数据: {band_meta}")
ds = None
# 写入元数据
ds = gdal.Open('example.tif', gdal.GA_Update)
ds.SetMetadataItem('AUTHOR', 'Your Name')
ds.SetMetadataItem('DATE', '2024-01-01')
ds.SetMetadata({'KEY1': 'Value1', 'KEY2': 'Value2'})
ds.FlushCache()
ds = None
3.8 本章小结
本章深入介绍了 GDAL 的核心架构和数据模型：
- 整体架构：分层设计（应用层、抽象层、驱动层、数据源层）
- 驱动模型：驱动注册、能力标识、常用驱动配置
- 栅格数据模型：
- GDALDataset：数据集容器
- GDALRasterBand：波段操作
- GeoTransform：地理变换
- 数据类型和分块读取
- 矢量数据模型：
- DataSource/Layer/Feature 层次结构
- 几何类型（点、线、面等）
- 字段定义
- 空间参考系统：坐标系定义和转换
- 虚拟文件系统：内存文件、压缩文件、网络资源访问
- 元数据系统：元数据域和读写操作
理解这些核心概念是有效使用 GDAL 的基础。在后续章节中，我们将详细学习如何运用这些知识处理实际的地理空间数据。
3.9 思考与练习
- 解释 GDAL 的分层架构设计有什么优势？
- GeoTransform 的6个参数分别代表什么含义？
- 如何判断一个 GDAL 驱动是否支持创建新数据集？
- OGR 的 Layer/Feature/Geometry 三者之间是什么关系？
- 虚拟文件系统 /vsicurl/ 和 /vsis3/ 分别用于什么场景？
- 编写代码创建一个包含3个波段的内存栅格数据集。
- 编写代码将 WGS84 坐标转换为 Web Mercator 坐标。
浙公网安备 33010602011771号