第11章 - 栅格数据高级处理
第11章：栅格数据高级处理
11.1 影像金字塔与概览
11.1.1 金字塔原理
影像金字塔是一种多分辨率影像表示方法，通过预先计算不同分辨率的影像层，加速显示和访问。
原始影像 (Level 0): 4000 x 4000
│
▼
Level 1: 2000 x 2000 (1/2)
│
▼
Level 2: 1000 x 1000 (1/4)
│
▼
Level 3: 500 x 500 (1/8)
│
▼
Level 4: 250 x 250 (1/16)
11.1.2 创建金字塔
from osgeo import gdal
def build_overviews(filepath, levels=None, resampling='AVERAGE'):
"""创建影像金字塔"""
ds = gdal.Open(filepath, gdal.GA_Update)
if ds is None:
raise Exception(f"无法打开: {filepath}")
if levels is None:
# 自动计算级别
max_dim = max(ds.RasterXSize, ds.RasterYSize)
levels = []
level = 2
while max_dim // level >= 256:
levels.append(level)
level *= 2
print(f"创建金字塔级别: {levels}")
# 可选的重采样方法
resample_methods = {
'NEAREST': gdal.GRA_NearestNeighbour,
'AVERAGE': gdal.GRA_Average,
'BILINEAR': gdal.GRA_Bilinear,
'CUBIC': gdal.GRA_Cubic,
'CUBICSPLINE': gdal.GRA_CubicSpline,
'LANCZOS': gdal.GRA_Lanczos,
'MODE': gdal.GRA_Mode,
}
ds.BuildOverviews(
resampling,
levels,
callback=gdal.TermProgress_nocb
)
ds.FlushCache()
ds = None
print("金字塔创建完成")
def build_external_overviews(filepath, levels=None):
"""创建外部金字塔(.ovr文件)"""
gdal.SetConfigOption('COMPRESS_OVERVIEW', 'LZW')
gdal.SetConfigOption('INTERLEAVE_OVERVIEW', 'PIXEL')
ds = gdal.Open(filepath, gdal.GA_ReadOnly)
if levels is None:
levels = [2, 4, 8, 16, 32]
ds.BuildOverviews('AVERAGE', levels)
ds = None
# 使用示例
# build_overviews('large_image.tif')
11.2 VRT 虚拟数据集
11.2.1 VRT 概念
VRT (Virtual Raster) 是一种 XML 格式的虚拟数据集，可以将多个栅格文件组合成一个逻辑数据集，而无需实际复制数据。
11.2.2 创建 VRT
from osgeo import gdal
def create_vrt_mosaic(input_files, output_vrt, nodata=None):
"""创建 VRT 镶嵌"""
vrt_options = gdal.BuildVRTOptions(
resampleAlg='nearest',
srcNodata=nodata,
VRTNodata=nodata,
resolution='average',
)
vrt = gdal.BuildVRT(output_vrt, input_files, options=vrt_options)
vrt = None
print(f"VRT 创建完成: {output_vrt}")
def create_vrt_stack(input_files, output_vrt):
"""创建波段堆叠 VRT"""
vrt_options = gdal.BuildVRTOptions(
separate=True, # 分离波段
resampleAlg='nearest',
)
vrt = gdal.BuildVRT(output_vrt, input_files, options=vrt_options)
vrt = None
def vrt_to_tiff(vrt_path, tiff_path):
"""将 VRT 转换为 GeoTIFF"""
translate_options = gdal.TranslateOptions(
format='GTiff',
creationOptions=[
'COMPRESS=LZW',
'TILED=YES',
'BIGTIFF=IF_SAFER'
]
)
ds = gdal.Translate(tiff_path, vrt_path, options=translate_options)
ds = None
11.2.3 VRT XML 结构
<VRTDataset rasterXSize="2000" rasterYSize="2000">
<SRS>EPSG:4326</SRS>
<GeoTransform>116.0, 0.001, 0, 40.0, 0, -0.001</GeoTransform>
<VRTRasterBand dataType="Float32" band="1">
<NoDataValue>-9999</NoDataValue>
<ComplexSource>
<SourceFilename relativeToVRT="1">tile1.tif</SourceFilename>
<SourceBand>1</SourceBand>
<SrcRect xOff="0" yOff="0" xSize="1000" ySize="1000"/>
<DstRect xOff="0" yOff="0" xSize="1000" ySize="1000"/>
</ComplexSource>
<ComplexSource>
<SourceFilename relativeToVRT="1">tile2.tif</SourceFilename>
<SourceBand>1</SourceBand>
<SrcRect xOff="0" yOff="0" xSize="1000" ySize="1000"/>
<DstRect xOff="1000" yOff="0" xSize="1000" ySize="1000"/>
</ComplexSource>
</VRTRasterBand>
</VRTDataset>
11.3 Cloud Optimized GeoTIFF (COG)
11.3.1 COG 概念
Cloud Optimized GeoTIFF 是一种特殊的 GeoTIFF 格式，经过优化以便在云存储中高效访问。
11.3.2 创建 COG
from osgeo import gdal
def create_cog(input_path, output_path, compression='LZW'):
"""创建 Cloud Optimized GeoTIFF"""
translate_options = gdal.TranslateOptions(
format='COG',
creationOptions=[
f'COMPRESS={compression}',
'OVERVIEWS=AUTO',
'BLOCKSIZE=512',
]
)
ds = gdal.Translate(output_path, input_path, options=translate_options)
ds = None
print(f"COG 创建完成: {output_path}")
def validate_cog(filepath):
"""验证是否为有效的 COG"""
ds = gdal.Open(filepath)
if ds is None:
return False, "无法打开文件"
errors = []
# 检查驱动
driver = ds.GetDriver().ShortName
if driver != 'GTiff':
errors.append(f"不是 GeoTIFF 格式: {driver}")
# 检查分块
band = ds.GetRasterBand(1)
block_size = band.GetBlockSize()
if block_size[0] == ds.RasterXSize:
errors.append("未使用分块存储")
# 检查金字塔
if band.GetOverviewCount() == 0:
errors.append("没有金字塔")
ds = None
if errors:
return False, errors
return True, "有效的 COG"
# 使用示例
# create_cog('input.tif', 'output_cog.tif')
# valid, msg = validate_cog('output_cog.tif')
11.4 高级裁剪与镶嵌
11.4.1 复杂裁剪
from osgeo import gdal, ogr
def clip_raster_by_geometry(raster_path, geometry, output_path, nodata=-9999):
"""使用几何对象裁剪栅格"""
# 创建临时矢量数据源
mem_driver = ogr.GetDriverByName('Memory')
mem_ds = mem_driver.CreateDataSource('temp')
mem_layer = mem_ds.CreateLayer('clip', geom_type=ogr.wkbPolygon)
feature = ogr.Feature(mem_layer.GetLayerDefn())
feature.SetGeometry(geometry)
mem_layer.CreateFeature(feature)
# 裁剪
warp_options = gdal.WarpOptions(
cutlineDSName=mem_ds,
cutlineLayer='clip',
cropToCutline=True,
dstNodata=nodata,
format='GTiff',
creationOptions=['COMPRESS=LZW']
)
ds = gdal.Warp(output_path, raster_path, options=warp_options)
ds = None
mem_ds = None
def clip_raster_by_extent(raster_path, extent, output_path):
"""按范围裁剪"""
min_x, min_y, max_x, max_y = extent
warp_options = gdal.WarpOptions(
outputBounds=(min_x, min_y, max_x, max_y),
format='GTiff',
creationOptions=['COMPRESS=LZW']
)
ds = gdal.Warp(output_path, raster_path, options=warp_options)
ds = None
11.4.2 高级镶嵌
from osgeo import gdal
import numpy as np
def mosaic_with_feathering(input_files, output_path, overlap=10):
"""带羽化的镶嵌"""
# 首先创建基础镶嵌
vrt = gdal.BuildVRT('', input_files)
# 读取数据
data = vrt.ReadAsArray()
# 这里可以添加羽化处理逻辑
# ...
# 保存结果
driver = gdal.GetDriverByName('GTiff')
out_ds = driver.Create(
output_path,
vrt.RasterXSize,
vrt.RasterYSize,
1,
gdal.GDT_Float32
)
out_ds.SetGeoTransform(vrt.GetGeoTransform())
out_ds.SetProjection(vrt.GetProjection())
out_ds.GetRasterBand(1).WriteArray(data)
out_ds = None
vrt = None
def weighted_mosaic(input_files, output_path, weights=None):
"""加权镶嵌"""
if weights is None:
weights = [1.0] * len(input_files)
datasets = [gdal.Open(f) for f in input_files]
# 获取参考信息
ref_ds = datasets[0]
width = ref_ds.RasterXSize
height = ref_ds.RasterYSize
# 加权平均
result = np.zeros((height, width), dtype=np.float32)
weight_sum = np.zeros((height, width), dtype=np.float32)
for ds, w in zip(datasets, weights):
data = ds.GetRasterBand(1).ReadAsArray()
nodata = ds.GetRasterBand(1).GetNoDataValue()
if nodata is not None:
mask = data != nodata
result = np.where(mask, result + data * w, result)
weight_sum = np.where(mask, weight_sum + w, weight_sum)
# 归一化
result = np.where(weight_sum > 0, result / weight_sum, -9999)
# 保存
driver = gdal.GetDriverByName('GTiff')
out_ds = driver.Create(output_path, width, height, 1, gdal.GDT_Float32)
out_ds.SetGeoTransform(ref_ds.GetGeoTransform())
out_ds.SetProjection(ref_ds.GetProjection())
out_ds.GetRasterBand(1).WriteArray(result)
out_ds.GetRasterBand(1).SetNoDataValue(-9999)
for ds in datasets:
ds = None
out_ds = None
11.5 影像配准与校正
11.5.1 基于控制点的配准
from osgeo import gdal
def georeference_with_gcps(input_path, output_path, gcps, srs='EPSG:4326'):
"""使用控制点进行地理配准"""
ds = gdal.Open(input_path)
# 创建 GCP 对象列表
gcp_list = []
for i, (pixel_x, pixel_y, geo_x, geo_y) in enumerate(gcps):
gcp = gdal.GCP(geo_x, geo_y, 0, pixel_x, pixel_y)
gcp_list.append(gcp)
# 设置 GCP
ds.SetGCPs(gcp_list, srs)
# 使用 Warp 进行校正
warp_options = gdal.WarpOptions(
format='GTiff',
dstSRS=srs,
tps=True, # 使用薄板样条
creationOptions=['COMPRESS=LZW']
)
out_ds = gdal.Warp(output_path, ds, options=warp_options)
ds = None
out_ds = None
# 使用示例
# gcps = [
# (0, 0, 116.0, 40.0), # 左上角
# (1000, 0, 117.0, 40.0), # 右上角
# (0, 1000, 116.0, 39.0), # 左下角
# (1000, 1000, 117.0, 39.0) # 右下角
# ]
# georeference_with_gcps('input.tif', 'output.tif', gcps)
11.5.2 RPC 校正
from osgeo import gdal
def rpc_orthorectify(input_path, output_path, dem_path=None):
"""RPC 正射校正"""
ds = gdal.Open(input_path)
# 检查 RPC 元数据
rpc = ds.GetMetadata('RPC')
if not rpc:
raise Exception("没有 RPC 元数据")
# Warp 选项
options = ['-t_srs', 'EPSG:4326']
if dem_path:
options.extend(['-rpc', '-to', f'RPC_DEM={dem_path}'])
else:
options.extend(['-rpc', '-to', 'RPC_HEIGHT=0'])
warp_options = gdal.WarpOptions(
format='GTiff',
creationOptions=['COMPRESS=LZW', 'TILED=YES'],
options=options
)
out_ds = gdal.Warp(output_path, ds, options=warp_options)
ds = None
out_ds = None
11.6 影像增强与处理
11.6.1 直方图均衡化
from osgeo import gdal
import numpy as np
def histogram_equalization(input_path, output_path):
"""直方图均衡化"""
ds = gdal.Open(input_path)
bands = ds.RasterCount
width = ds.RasterXSize
height = ds.RasterYSize
driver = gdal.GetDriverByName('GTiff')
out_ds = driver.Create(output_path, width, height, bands, gdal.GDT_Byte)
out_ds.SetGeoTransform(ds.GetGeoTransform())
out_ds.SetProjection(ds.GetProjection())
for i in range(1, bands + 1):
band = ds.GetRasterBand(i)
data = band.ReadAsArray()
# 计算直方图
hist, bins = np.histogram(data.flatten(), 256, [0, 256])
cdf = hist.cumsum()
# 归一化
cdf_normalized = cdf * 255 / cdf[-1]
# 应用
equalized = np.interp(data.flatten(), bins[:-1], cdf_normalized)
equalized = equalized.reshape(height, width).astype(np.uint8)
out_ds.GetRasterBand(i).WriteArray(equalized)
ds = None
out_ds = None
def linear_stretch(input_path, output_path, percent=2):
"""线性拉伸"""
ds = gdal.Open(input_path)
bands = ds.RasterCount
width = ds.RasterXSize
height = ds.RasterYSize
driver = gdal.GetDriverByName('GTiff')
out_ds = driver.Create(output_path, width, height, bands, gdal.GDT_Byte)
out_ds.SetGeoTransform(ds.GetGeoTransform())
out_ds.SetProjection(ds.GetProjection())
for i in range(1, bands + 1):
band = ds.GetRasterBand(i)
data = band.ReadAsArray().astype(np.float32)
# 计算百分比截断
min_val = np.percentile(data, percent)
max_val = np.percentile(data, 100 - percent)
# 拉伸
stretched = (data - min_val) / (max_val - min_val) * 255
stretched = np.clip(stretched, 0, 255).astype(np.uint8)
out_ds.GetRasterBand(i).WriteArray(stretched)
ds = None
out_ds = None
11.6.2 滤波处理
from osgeo import gdal
import numpy as np
from scipy import ndimage
def apply_filter(input_path, output_path, filter_type='mean', size=3):
"""应用滤波"""
ds = gdal.Open(input_path)
band = ds.GetRasterBand(1)
data = band.ReadAsArray().astype(np.float32)
if filter_type == 'mean':
filtered = ndimage.uniform_filter(data, size=size)
elif filter_type == 'median':
filtered = ndimage.median_filter(data, size=size)
elif filter_type == 'gaussian':
filtered = ndimage.gaussian_filter(data, sigma=size/3)
elif filter_type == 'sobel':
sx = ndimage.sobel(data, axis=0)
sy = ndimage.sobel(data, axis=1)
filtered = np.hypot(sx, sy)
else:
raise ValueError(f"未知的滤波类型: {filter_type}")
# 保存
driver = gdal.GetDriverByName('GTiff')
out_ds = driver.Create(
output_path,
ds.RasterXSize,
ds.RasterYSize,
1,
gdal.GDT_Float32
)
out_ds.SetGeoTransform(ds.GetGeoTransform())
out_ds.SetProjection(ds.GetProjection())
out_ds.GetRasterBand(1).WriteArray(filtered)
ds = None
out_ds = None
11.7 多维数组数据处理
11.7.1 NetCDF 数据
from osgeo import gdal
def read_netcdf(filepath, variable=None):
"""读取 NetCDF 数据"""
ds = gdal.Open(filepath)
if ds is None:
raise Exception(f"无法打开: {filepath}")
# 获取子数据集
subdatasets = ds.GetSubDatasets()
print(f"子数据集数量: {len(subdatasets)}")
for i, (path, desc) in enumerate(subdatasets):
print(f" {i}: {desc}")
if variable:
# 打开指定变量
for path, desc in subdatasets:
if variable in desc:
sub_ds = gdal.Open(path)
return sub_ds
ds = None
return None
def netcdf_to_geotiff(nc_path, variable, output_path, time_index=0):
"""将 NetCDF 变量导出为 GeoTIFF"""
ds = gdal.Open(nc_path)
subdatasets = ds.GetSubDatasets()
# 查找变量
var_path = None
for path, desc in subdatasets:
if variable in path:
var_path = path
break
if var_path is None:
raise Exception(f"未找到变量: {variable}")
# 打开子数据集
sub_ds = gdal.Open(var_path)
# 如果是多时相数据，选择特定时间
band = sub_ds.GetRasterBand(time_index + 1)
data = band.ReadAsArray()
# 创建输出
driver = gdal.GetDriverByName('GTiff')
out_ds = driver.Create(
output_path,
sub_ds.RasterXSize,
sub_ds.RasterYSize,
1,
band.DataType
)
out_ds.SetGeoTransform(sub_ds.GetGeoTransform())
out_ds.SetProjection(sub_ds.GetProjection())
out_ds.GetRasterBand(1).WriteArray(data)
sub_ds = None
out_ds = None
ds = None
11.8 本章小结
本章介绍了栅格数据高级处理技术：
- 影像金字塔：创建和管理多分辨率概览
- VRT 虚拟数据集：灵活组合多个数据集
- COG 格式：云优化 GeoTIFF
- 高级裁剪镶嵌：复杂的几何裁剪和加权镶嵌
- 影像配准：GCP 和 RPC 校正
- 影像增强：直方图均衡化、滤波
- 多维数据：NetCDF 等格式处理
11.9 思考与练习
- 为什么使用 COG 格式可以提高云端访问效率？
- VRT 和实际镶嵌有什么区别，各适用于什么场景？
- 实现一个支持多种重采样方法的影像缩放函数。
- 编写代码处理多时相的 NetCDF 气象数据。
- 如何优化大规模影像镶嵌的内存使用？
浙公网安备 33010602011771号