第15章 - 实战案例与项目应用
第15章：实战案例与项目应用
15.1 概述
本章通过实际项目案例，展示如何综合运用 GDAL 的各项功能解决实际问题。这些案例涵盖了遥感影像处理、地理数据转换、空间分析等多个领域。
15.2 案例一：批量影像预处理流水线
15.2.1 需求分析
对一批卫星影像进行标准化预处理：
- 坐标系统统一到 EPSG:4326
- 裁剪到指定区域
- 重采样到统一分辨率
- 输出为压缩的 GeoTIFF 格式
15.2.2 Python 实现
from osgeo import gdal, osr
import os
from pathlib import Path
class ImageProcessor:
"""影像批量处理器"""
def __init__(self, output_srs='EPSG:4326', resolution=0.0001):
self.output_srs = output_srs
self.resolution = resolution
gdal.SetCacheMax(1024 * 1024 * 1024) # 1GB 缓存
def process_image(self, input_path, output_path, bbox=None):
"""处理单个影像"""
print(f"处理: {input_path}")
src_ds = gdal.Open(input_path)
if src_ds is None:
raise RuntimeError(f"无法打开文件: {input_path}")
try:
warp_options = {
'format': 'GTiff',
'dstSRS': self.output_srs,
'xRes': self.resolution,
'yRes': self.resolution,
'resampleAlg': gdal.GRA_Bilinear,
'creationOptions': ['COMPRESS=LZW', 'TILED=YES']
}
if bbox:
warp_options['outputBounds'] = bbox
ds = gdal.Warp(output_path, src_ds, **warp_options)
if ds:
print(f"✓ 完成: {output_path}")
ds = None
return True
else:
return False
finally:
src_ds = None
def batch_process(self, input_dir, output_dir, bbox=None):
"""批量处理目录中的所有影像"""
input_path = Path(input_dir)
output_path = Path(output_dir)
output_path.mkdir(parents=True, exist_ok=True)
extensions = ['.tif', '.tiff', '.img', '.jp2']
success_count = 0
fail_count = 0
for img_file in input_path.iterdir():
if img_file.suffix.lower() in extensions:
out_file = output_path / f"{img_file.stem}_processed.tif"
try:
if self.process_image(str(img_file), str(out_file), bbox):
success_count += 1
else:
fail_count += 1
except Exception as e:
print(f"✗ 错误: {img_file.name} - {e}")
fail_count += 1
print(f"\n处理完成: 成功 {success_count}, 失败 {fail_count}")
# 使用示例
processor = ImageProcessor(output_srs='EPSG:4326', resolution=0.00025)
beijing_bbox = (115.7, 39.4, 117.4, 41.6)
processor.batch_process('./raw_images', './processed_images', bbox=beijing_bbox)
15.3 案例二：地理数据 ETL 系统
15.3.1 完整实现
from osgeo import ogr, osr
import logging
class GeoDataETL:
"""地理数据 ETL 处理器"""
def __init__(self):
self.logger = logging.getLogger('GeoDataETL')
self.logger.setLevel(logging.INFO)
def extract(self, source_path, layer_name=None):
"""提取数据"""
self.logger.info(f"提取数据: {source_path}")
ds = ogr.Open(source_path, 0)
if ds is None:
raise RuntimeError(f"无法打开数据源: {source_path}")
layer = ds.GetLayerByName(layer_name) if layer_name else ds.GetLayer(0)
if layer is None:
raise RuntimeError(f"无法获取图层: {layer_name}")
return layer, ds
def transform(self, layer, target_srs='EPSG:4326', attribute_filter=None):
"""转换数据"""
if attribute_filter:
layer.SetAttributeFilter(attribute_filter)
source_srs = layer.GetSpatialRef()
target_srs_obj = osr.SpatialReference()
target_srs_obj.SetFromUserInput(target_srs)
transform = None
if not source_srs.IsSame(target_srs_obj):
transform = osr.CoordinateTransformation(source_srs, target_srs_obj)
return transform, target_srs_obj
def load(self, layer, output_path, format_name='ESRI Shapefile',
transform=None, target_srs=None):
"""加载数据到目标"""
driver = ogr.GetDriverByName(format_name)
if os.path.exists(output_path):
driver.DeleteDataSource(output_path)
out_ds = driver.CreateDataSource(output_path)
out_layer = out_ds.CreateLayer(
'output',
srs=target_srs if target_srs else layer.GetSpatialRef(),
geom_type=layer.GetGeomType()
)
# 复制字段定义
layer_defn = layer.GetLayerDefn()
for i in range(layer_defn.GetFieldCount()):
out_layer.CreateField(layer_defn.GetFieldDefn(i))
# 复制要素
count = 0
for feature in layer:
out_feature = ogr.Feature(out_layer.GetLayerDefn())
for i in range(layer_defn.GetFieldCount()):
out_feature.SetField(i, feature.GetField(i))
geom = feature.GetGeometryRef()
if geom:
geom_copy = geom.Clone()
if transform:
geom_copy.Transform(transform)
out_feature.SetGeometry(geom_copy)
out_layer.CreateFeature(out_feature)
count += 1
self.logger.info(f"✓ 成功加载 {count} 个要素")
out_ds = None
# 使用示例
etl = GeoDataETL()
etl.run(source_path='input.shp', output_path='output.geojson',
format='GeoJSON', target_srs='EPSG:4326')
15.4 案例三：DEM 地形分析
from osgeo import gdal
class DEMAnalyzer:
"""DEM 地形分析器"""
def __init__(self, dem_path):
self.ds = gdal.Open(dem_path)
if self.ds is None:
raise RuntimeError(f"无法打开 DEM: {dem_path}")
def calculate_slope(self, output_path):
"""计算坡度"""
gdal.DEMProcessing(output_path, self.ds, 'slope',
slopeFormat='percent', computeEdges=True)
print(f"✓ 坡度已保存到: {output_path}")
def calculate_aspect(self, output_path):
"""计算坡向"""
gdal.DEMProcessing(output_path, self.ds, 'aspect', computeEdges=True)
print(f"✓ 坡向已保存到: {output_path}")
def calculate_hillshade(self, output_path, azimuth=315, altitude=45):
"""计算山体阴影"""
gdal.DEMProcessing(output_path, self.ds, 'hillshade',
azimuth=azimuth, altitude=altitude, computeEdges=True)
print(f"✓ 山体阴影已保存到: {output_path}")
# 使用示例
analyzer = DEMAnalyzer('dem.tif')
analyzer.calculate_slope('slope.tif')
analyzer.calculate_aspect('aspect.tif')
analyzer.calculate_hillshade('hillshade.tif')
15.5 性能监控
import time
import psutil
from contextlib import contextmanager
@contextmanager
def performance_monitor(operation_name):
"""性能监控上下文管理器"""
process = psutil.Process()
start_time = time.time()
start_memory = process.memory_info().rss / 1024 / 1024
print(f"\n{'='*60}")
print(f"开始: {operation_name}")
print(f"{'='*60}")
try:
yield
finally:
end_time = time.time()
end_memory = process.memory_info().rss / 1024 / 1024
print(f"\n{'='*60}")
print(f"完成: {operation_name}")
print(f"耗时: {end_time - start_time:.2f} 秒")
print(f"内存使用: {end_memory - start_memory:.2f} MB")
print(f"{'='*60}\n")
# 使用示例
with performance_monitor("影像重采样"):
gdal.Warp('output.tif', 'input.tif', xRes=30, yRes=30)
15.6 Docker 容器化部署
# Dockerfile
FROM osgeo/gdal:ubuntu-full-latest
WORKDIR /app
COPY requirements.txt .
RUN pip install -r requirements.txt
COPY . .
ENV GDAL_CACHEMAX=1024
ENV GDAL_NUM_THREADS=ALL_CPUS
CMD ["python", "main.py"]
15.7 本章小结
本章通过三个完整的实战案例，展示了 GDAL 在实际项目中的应用：
- 影像预处理流水线：自动化批量处理遥感影像
- 地理数据 ETL 系统：构建完整的数据转换流程
- DEM 地形分析：地形指标计算和分析
关键要点：
- 合理使用缓存和并行处理提升性能
- 完善的错误处理和日志记录
- 生产环境的配置和部署
- 性能监控和优化
相关资源：
浙公网安备 33010602011771号