第04章 - 栅格数据处理基础

第04章：栅格数据处理基础

4.1 栅格数据概述

4.1.1 什么是栅格数据

栅格数据是由规则排列的像元（像素）组成的数据结构，每个像元存储一个或多个数值。在地理空间领域，栅格数据广泛用于表示连续的地理现象。

┌─────────────────────────────────────────────────────────────┐
│                     栅格数据结构                              │
├─────────────────────────────────────────────────────────────┤
│                                                              │
│    ┌───┬───┬───┬───┬───┬───┬───┬───┐                       │
│    │ 45│ 48│ 52│ 55│ 58│ 60│ 62│ 65│  ← 行 0               │
│    ├───┼───┼───┼───┼───┼───┼───┼───┤                       │
│    │ 42│ 46│ 50│ 53│ 56│ 59│ 61│ 63│  ← 行 1               │
│    ├───┼───┼───┼───┼───┼───┼───┼───┤                       │
│    │ 40│ 44│ 48│ 51│ 54│ 57│ 60│ 62│  ← 行 2               │
│    ├───┼───┼───┼───┼───┼───┼───┼───┤                       │
│    │ 38│ 42│ 46│ 49│ 52│ 55│ 58│ 60│  ← 行 3               │
│    └───┴───┴───┴───┴───┴───┴───┴───┘                       │
│      ↑                                                       │
│     列0  列1  列2  列3  列4  列5  列6  列7                   │
│                                                              │
│    每个格子 = 一个像元（Pixel/Cell）                         │
│    格子中的数值 = 像元值（可以是高程、温度、反射率等）         │
│                                                              │
└─────────────────────────────────────────────────────────────┘

4.1.2 常见栅格数据类型

数据类型 描述 示例
卫星影像 遥感传感器获取的图像 Landsat、Sentinel、MODIS
数字高程模型 地表高程数据 SRTM、ASTER GDEM、ALOS
正射影像 经过正射校正的航空/卫星影像 无人机影像、航拍图
土地覆盖 分类后的土地利用数据 GlobeLand30、FROM-GLC
气象数据 气温、降水、风速等 ERA5、GPM
专题图 各种专题分析结果 NDVI、坡度、坡向

4.1.3 常见栅格格式

格式 扩展名 特点
GeoTIFF .tif, .tiff 行业标准，广泛支持
Cloud Optimized GeoTIFF .tif 云优化，支持范围读取
JPEG 2000 .jp2 高压缩率，有损/无损
HDF5 .h5, .hdf5 科学数据格式，多维数组
NetCDF .nc 气象海洋数据标准格式
IMG .img ERDAS IMAGINE 格式
ECW .ecw 高压缩率商业格式

4.2 读取栅格数据

4.2.1 打开数据集

from osgeo import gdal
import numpy as np

# 启用异常处理
gdal.UseExceptions()

# 打开数据集（只读）
ds = gdal.Open('example.tif')

# 打开数据集（读写）
ds = gdal.Open('example.tif', gdal.GA_Update)

# 检查是否成功打开
if ds is None:
    print("无法打开文件")
else:
    print(f"成功打开: {ds.GetDescription()}")

# 获取基本信息
print(f"宽度: {ds.RasterXSize} 像素")
print(f"高度: {ds.RasterYSize} 像素")
print(f"波段数: {ds.RasterCount}")
print(f"驱动: {ds.GetDriver().ShortName}")

# 关闭数据集（重要！）
ds = None

4.2.2 使用上下文管理器

推荐使用上下文管理器来自动管理资源：

from osgeo import gdal
from contextlib import contextmanager

@contextmanager
def open_raster(filepath, mode=gdal.GA_ReadOnly):
    """栅格数据集上下文管理器"""
    ds = gdal.Open(filepath, mode)
    if ds is None:
        raise IOError(f"无法打开文件: {filepath}")
    try:
        yield ds
    finally:
        ds = None

# 使用示例
with open_raster('example.tif') as ds:
    print(f"宽度: {ds.RasterXSize}")
    print(f"高度: {ds.RasterYSize}")
# 自动关闭

4.2.3 获取数据集信息

from osgeo import gdal, osr

def get_raster_info(filepath):
    """获取栅格数据详细信息"""
    ds = gdal.Open(filepath)
    
    info = {
        'filepath': filepath,
        'driver': ds.GetDriver().ShortName,
        'width': ds.RasterXSize,
        'height': ds.RasterYSize,
        'bands': ds.RasterCount,
    }
    
    # 地理变换
    gt = ds.GetGeoTransform()
    if gt:
        info['geotransform'] = {
            'origin_x': gt[0],
            'origin_y': gt[3],
            'pixel_width': gt[1],
            'pixel_height': gt[5],
            'rotation_x': gt[2],
            'rotation_y': gt[4],
        }
        
        # 计算边界
        info['bounds'] = {
            'min_x': gt[0],
            'max_x': gt[0] + ds.RasterXSize * gt[1],
            'min_y': gt[3] + ds.RasterYSize * gt[5],
            'max_y': gt[3],
        }
    
    # 投影信息
    proj = ds.GetProjection()
    if proj:
        srs = osr.SpatialReference()
        srs.ImportFromWkt(proj)
        info['projection'] = {
            'name': srs.GetName(),
            'epsg': srs.GetAuthorityCode(None),
            'is_geographic': srs.IsGeographic(),
            'is_projected': srs.IsProjected(),
        }
    
    # 波段信息
    info['band_info'] = []
    for i in range(1, ds.RasterCount + 1):
        band = ds.GetRasterBand(i)
        band_info = {
            'band_number': i,
            'data_type': gdal.GetDataTypeName(band.DataType),
            'nodata': band.GetNoDataValue(),
            'block_size': band.GetBlockSize(),
            'color_interp': gdal.GetColorInterpretationName(
                band.GetColorInterpretation()
            ),
        }
        
        # 统计信息（如果可用）
        stats = band.GetStatistics(False, False)
        if stats[0] != stats[1]:  # 检查是否有效
            band_info['statistics'] = {
                'min': stats[0],
                'max': stats[1],
                'mean': stats[2],
                'std': stats[3],
            }
        
        info['band_info'].append(band_info)
    
    ds = None
    return info

# 使用示例
info = get_raster_info('example.tif')
import json
print(json.dumps(info, indent=2, ensure_ascii=False))

4.2.4 读取像素数据

from osgeo import gdal
import numpy as np

ds = gdal.Open('example.tif')

# 方法1：读取整个波段为数组
band = ds.GetRasterBand(1)
data = band.ReadAsArray()
print(f"数据形状: {data.shape}")
print(f"数据类型: {data.dtype}")

# 方法2：读取所有波段
all_bands = ds.ReadAsArray()
print(f"所有波段形状: {all_bands.shape}")  # (bands, height, width)

# 方法3：读取指定区域
x_off, y_off = 100, 200  # 起始位置
x_size, y_size = 50, 50  # 读取大小
subset = band.ReadAsArray(x_off, y_off, x_size, y_size)
print(f"子集形状: {subset.shape}")

# 方法4：读取并重采样
target_width = 500
target_height = 500
resampled = band.ReadAsArray(
    buf_xsize=target_width,
    buf_ysize=target_height
)
print(f"重采样后形状: {resampled.shape}")

# 方法5：使用 ReadRaster 读取原始字节
raw_data = band.ReadRaster(0, 0, 100, 100)
arr = np.frombuffer(raw_data, dtype=np.uint8).reshape(100, 100)

ds = None

4.2.5 处理无效值

from osgeo import gdal
import numpy as np

ds = gdal.Open('example.tif')
band = ds.GetRasterBand(1)

# 获取无效值
nodata = band.GetNoDataValue()
print(f"无效值: {nodata}")

# 读取数据
data = band.ReadAsArray()

# 方法1：使用布尔掩码
if nodata is not None:
    valid_mask = data != nodata
    valid_data = data[valid_mask]
    print(f"有效像元数: {valid_data.size}")
    print(f"有效值范围: {valid_data.min()} - {valid_data.max()}")

# 方法2：使用 NumPy 掩码数组
masked_data = np.ma.masked_equal(data, nodata)
print(f"有效像元平均值: {masked_data.mean()}")

# 方法3：替换无效值
data_filled = np.where(data == nodata, np.nan, data)
print(f"均值（忽略NaN）: {np.nanmean(data_filled)}")

ds = None

4.3 创建栅格数据

4.3.1 使用 Create 方法

from osgeo import gdal, osr
import numpy as np

def create_raster(filepath, width, height, bands=1, 
                  dtype=gdal.GDT_Float32, driver='GTiff'):
    """创建新的栅格数据集"""
    
    # 获取驱动
    driver_obj = gdal.GetDriverByName(driver)
    
    # 创建选项
    options = ['COMPRESS=LZW', 'TILED=YES']
    
    # 创建数据集
    ds = driver_obj.Create(filepath, width, height, bands, dtype, options)
    
    return ds

# 创建简单栅格
ds = create_raster('output.tif', 1000, 1000)

# 设置地理变换
# 假设覆盖北京区域
origin_x = 116.0
origin_y = 40.5
pixel_size = 0.001  # 约111米

geotransform = (origin_x, pixel_size, 0, origin_y, 0, -pixel_size)
ds.SetGeoTransform(geotransform)

# 设置投影（WGS84）
srs = osr.SpatialReference()
srs.ImportFromEPSG(4326)
ds.SetProjection(srs.ExportToWkt())

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

# 设置无效值
band.SetNoDataValue(-9999)

# 计算统计信息
band.ComputeStatistics(False)

# 刷新并关闭
ds.FlushCache()
ds = None

4.3.2 使用 CreateCopy 方法

from osgeo import gdal

def copy_raster(src_path, dst_path, driver='GTiff', options=None):
    """复制栅格数据集"""
    
    src_ds = gdal.Open(src_path)
    
    driver_obj = gdal.GetDriverByName(driver)
    
    if options is None:
        options = ['COMPRESS=LZW', 'TILED=YES']
    
    dst_ds = driver_obj.CreateCopy(dst_path, src_ds, 0, options)
    
    src_ds = None
    dst_ds = None

# 简单复制
copy_raster('input.tif', 'output.tif')

# 转换格式（GeoTIFF 转 JPEG）
copy_raster('input.tif', 'output.jpg', driver='JPEG', options=['QUALITY=90'])

# 复制并压缩
copy_raster('input.tif', 'compressed.tif', 
            options=['COMPRESS=DEFLATE', 'PREDICTOR=2', 'ZLEVEL=9'])

4.3.3 创建多波段彩色影像

from osgeo import gdal, osr
import numpy as np

def create_rgb_raster(filepath, width, height, 
                      red_data, green_data, blue_data):
    """创建RGB彩色栅格"""
    
    driver = gdal.GetDriverByName('GTiff')
    
    ds = driver.Create(
        filepath, width, height, 3, gdal.GDT_Byte,
        ['COMPRESS=LZW', 'PHOTOMETRIC=RGB']
    )
    
    # 写入各波段
    for i, (data, interp) in enumerate([
        (red_data, gdal.GCI_RedBand),
        (green_data, gdal.GCI_GreenBand),
        (blue_data, gdal.GCI_BlueBand)
    ], 1):
        band = ds.GetRasterBand(i)
        band.WriteArray(data)
        band.SetColorInterpretation(interp)
    
    return ds

# 创建测试数据
width, height = 500, 500
red = np.random.randint(0, 256, (height, width), dtype=np.uint8)
green = np.random.randint(0, 256, (height, width), dtype=np.uint8)
blue = np.random.randint(0, 256, (height, width), dtype=np.uint8)

ds = create_rgb_raster('rgb_image.tif', width, height, red, green, blue)

# 设置地理信息
ds.SetGeoTransform((116.0, 0.001, 0, 40.0, 0, -0.001))
srs = osr.SpatialReference()
srs.ImportFromEPSG(4326)
ds.SetProjection(srs.ExportToWkt())

ds.FlushCache()
ds = None

4.3.4 创建金字塔

from osgeo import gdal

def build_overviews(filepath, levels=None, resampling='AVERAGE'):
    """创建影像金字塔"""
    
    ds = gdal.Open(filepath, gdal.GA_Update)
    
    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}")
    
    # 创建金字塔
    ds.BuildOverviews(resampling, levels)
    
    ds = None

# 使用示例
build_overviews('large_image.tif', levels=[2, 4, 8, 16, 32])

# 不同的重采样方法
# NEAREST: 最近邻（分类数据）
# AVERAGE: 平均值（连续数据）
# BILINEAR: 双线性
# CUBIC: 三次卷积
# CUBICSPLINE: 三次样条
# LANCZOS: Lanczos
# MODE: 众数（分类数据）

4.4 栅格数据处理

4.4.1 裁剪栅格

from osgeo import gdal

def clip_raster_by_extent(src_path, dst_path, 
                          min_x, min_y, max_x, max_y):
    """按范围裁剪栅格"""
    
    ds = gdal.Warp(
        dst_path,
        src_path,
        outputBounds=(min_x, min_y, max_x, max_y),
        format='GTiff',
        creationOptions=['COMPRESS=LZW']
    )
    
    ds = None

def clip_raster_by_shapefile(src_path, dst_path, shapefile_path):
    """按矢量边界裁剪栅格"""
    
    ds = gdal.Warp(
        dst_path,
        src_path,
        cutlineDSName=shapefile_path,
        cropToCutline=True,
        dstNodata=-9999,
        format='GTiff',
        creationOptions=['COMPRESS=LZW']
    )
    
    ds = None

# 按范围裁剪（北京中心区域）
clip_raster_by_extent(
    'input.tif', 'clipped.tif',
    116.3, 39.8, 116.5, 40.0
)

# 按边界裁剪
clip_raster_by_shapefile(
    'input.tif', 'clipped_by_boundary.tif',
    'boundary.shp'
)

4.4.2 重投影

from osgeo import gdal

def reproject_raster(src_path, dst_path, 
                     dst_crs='EPSG:3857',
                     resampling='bilinear'):
    """重投影栅格"""
    
    # 重采样方法映射
    resample_methods = {
        'nearest': gdal.GRA_NearestNeighbour,
        'bilinear': gdal.GRA_Bilinear,
        'cubic': gdal.GRA_Cubic,
        'cubicspline': gdal.GRA_CubicSpline,
        'lanczos': gdal.GRA_Lanczos,
        'average': gdal.GRA_Average,
        'mode': gdal.GRA_Mode,
    }
    
    ds = gdal.Warp(
        dst_path,
        src_path,
        dstSRS=dst_crs,
        resampleAlg=resample_methods.get(resampling, gdal.GRA_Bilinear),
        format='GTiff',
        creationOptions=['COMPRESS=LZW', 'TILED=YES']
    )
    
    ds = None

# WGS84 转 Web Mercator
reproject_raster('wgs84.tif', 'webmercator.tif', 'EPSG:3857')

# WGS84 转 CGCS2000
reproject_raster('wgs84.tif', 'cgcs2000.tif', 'EPSG:4490')

# 转换到 UTM 投影
reproject_raster('wgs84.tif', 'utm50n.tif', 'EPSG:32650')

4.4.3 重采样

from osgeo import gdal

def resample_raster(src_path, dst_path, 
                    target_width=None, target_height=None,
                    target_resolution=None, resampling='bilinear'):
    """重采样栅格"""
    
    src_ds = gdal.Open(src_path)
    
    # 计算目标尺寸
    if target_resolution:
        gt = src_ds.GetGeoTransform()
        extent_x = src_ds.RasterXSize * abs(gt[1])
        extent_y = src_ds.RasterYSize * abs(gt[5])
        target_width = int(extent_x / target_resolution)
        target_height = int(extent_y / target_resolution)
    
    resample_methods = {
        'nearest': gdal.GRA_NearestNeighbour,
        'bilinear': gdal.GRA_Bilinear,
        'cubic': gdal.GRA_Cubic,
        'average': gdal.GRA_Average,
    }
    
    ds = gdal.Warp(
        dst_path,
        src_ds,
        width=target_width,
        height=target_height,
        resampleAlg=resample_methods.get(resampling, gdal.GRA_Bilinear),
        format='GTiff',
        creationOptions=['COMPRESS=LZW']
    )
    
    src_ds = None
    ds = None

# 按尺寸重采样
resample_raster('input.tif', 'resampled.tif', 
                target_width=500, target_height=500)

# 按分辨率重采样（0.01度）
resample_raster('input.tif', 'resampled_res.tif', 
                target_resolution=0.01)

4.4.4 镶嵌

from osgeo import gdal

def mosaic_rasters(src_files, dst_path, nodata=None):
    """镶嵌多个栅格"""
    
    vrt_options = gdal.BuildVRTOptions(
        resampleAlg='nearest',
        srcNodata=nodata,
        VRTNodata=nodata
    )
    
    # 创建 VRT
    vrt = gdal.BuildVRT('', src_files, options=vrt_options)
    
    # 转换为 GeoTIFF
    ds = gdal.Translate(
        dst_path, vrt,
        format='GTiff',
        creationOptions=['COMPRESS=LZW', 'TILED=YES', 'BIGTIFF=IF_SAFER']
    )
    
    vrt = None
    ds = None

# 镶嵌多个 DEM 文件
dem_files = ['dem_n35e116.tif', 'dem_n35e117.tif', 
             'dem_n36e116.tif', 'dem_n36e117.tif']
mosaic_rasters(dem_files, 'mosaic.tif', nodata=-9999)

4.4.5 波段计算

from osgeo import gdal
import numpy as np

def calculate_ndvi(red_path, nir_path, output_path):
    """计算 NDVI"""
    
    red_ds = gdal.Open(red_path)
    nir_ds = gdal.Open(nir_path)
    
    red = red_ds.GetRasterBand(1).ReadAsArray().astype(np.float32)
    nir = nir_ds.GetRasterBand(1).ReadAsArray().astype(np.float32)
    
    # 计算 NDVI: (NIR - Red) / (NIR + Red)
    # 避免除零
    denominator = nir + red
    ndvi = np.where(denominator > 0, (nir - red) / denominator, 0)
    
    # 创建输出
    driver = gdal.GetDriverByName('GTiff')
    out_ds = driver.Create(
        output_path,
        red_ds.RasterXSize,
        red_ds.RasterYSize,
        1, gdal.GDT_Float32,
        ['COMPRESS=LZW']
    )
    
    out_ds.SetGeoTransform(red_ds.GetGeoTransform())
    out_ds.SetProjection(red_ds.GetProjection())
    
    out_band = out_ds.GetRasterBand(1)
    out_band.WriteArray(ndvi)
    out_band.SetNoDataValue(-9999)
    out_band.ComputeStatistics(False)
    
    red_ds = None
    nir_ds = None
    out_ds = None

# 使用 gdal_calc.py 进行波段计算
def band_calc_with_gdal(expression, output_path, **input_files):
    """使用 gdal_calc 进行波段计算"""
    import subprocess
    
    cmd = ['gdal_calc.py', '-A', input_files.get('A', ''),
           '--outfile', output_path, '--calc', expression,
           '--NoDataValue=-9999']
    
    for key, value in input_files.items():
        if key != 'A':
            cmd.extend([f'-{key}', value])
    
    subprocess.run(cmd)

# 示例：计算 NDVI
# band_calc_with_gdal(
#     '(B-A)/(B+A)', 'ndvi.tif',
#     A='red_band.tif', B='nir_band.tif'
# )

4.5 栅格分析

4.5.1 统计分析

from osgeo import gdal
import numpy as np

def raster_statistics(filepath):
    """计算栅格统计信息"""
    
    ds = gdal.Open(filepath)
    
    results = []
    
    for i in range(1, ds.RasterCount + 1):
        band = ds.GetRasterBand(i)
        
        # 强制计算统计
        stats = band.ComputeStatistics(False)
        
        # 获取直方图
        hist = band.GetHistogram()
        
        band_stats = {
            'band': i,
            'min': stats[0],
            'max': stats[1],
            'mean': stats[2],
            'std': stats[3],
            'nodata': band.GetNoDataValue(),
        }
        
        # 更详细的统计（使用 NumPy）
        data = band.ReadAsArray()
        nodata = band.GetNoDataValue()
        
        if nodata is not None:
            valid_data = data[data != nodata]
        else:
            valid_data = data.flatten()
        
        if len(valid_data) > 0:
            band_stats.update({
                'median': float(np.median(valid_data)),
                'sum': float(np.sum(valid_data)),
                'count': len(valid_data),
                'percentile_25': float(np.percentile(valid_data, 25)),
                'percentile_75': float(np.percentile(valid_data, 75)),
            })
        
        results.append(band_stats)
    
    ds = None
    return results

# 使用示例
stats = raster_statistics('dem.tif')
for s in stats:
    print(f"波段 {s['band']}:")
    print(f"  最小值: {s['min']:.2f}")
    print(f"  最大值: {s['max']:.2f}")
    print(f"  平均值: {s['mean']:.2f}")
    print(f"  标准差: {s['std']:.2f}")

4.5.2 分区统计

from osgeo import gdal, ogr
import numpy as np

def zonal_statistics(raster_path, vector_path, stats=['mean', 'sum', 'count']):
    """计算分区统计"""
    
    raster_ds = gdal.Open(raster_path)
    band = raster_ds.GetRasterBand(1)
    gt = raster_ds.GetGeoTransform()
    nodata = band.GetNoDataValue()
    
    vector_ds = ogr.Open(vector_path)
    layer = vector_ds.GetLayer()
    
    results = []
    
    for feature in layer:
        fid = feature.GetFID()
        geom = feature.GetGeometryRef()
        
        # 获取几何范围
        env = geom.GetEnvelope()  # (minX, maxX, minY, maxY)
        
        # 转换为像素坐标
        x_off = int((env[0] - gt[0]) / gt[1])
        y_off = int((env[3] - gt[3]) / gt[5])
        x_size = int((env[1] - env[0]) / gt[1]) + 1
        y_size = int((env[2] - env[3]) / gt[5]) + 1
        
        # 确保在影像范围内
        x_off = max(0, x_off)
        y_off = max(0, y_off)
        x_size = min(x_size, raster_ds.RasterXSize - x_off)
        y_size = min(y_size, raster_ds.RasterYSize - y_off)
        
        if x_size <= 0 or y_size <= 0:
            continue
        
        # 读取数据
        data = band.ReadAsArray(x_off, y_off, x_size, y_size)
        
        # 掩码处理
        if nodata is not None:
            valid_data = data[data != nodata]
        else:
            valid_data = data.flatten()
        
        # 计算统计
        feature_stats = {'fid': fid}
        
        if len(valid_data) > 0:
            if 'mean' in stats:
                feature_stats['mean'] = float(np.mean(valid_data))
            if 'sum' in stats:
                feature_stats['sum'] = float(np.sum(valid_data))
            if 'count' in stats:
                feature_stats['count'] = len(valid_data)
            if 'min' in stats:
                feature_stats['min'] = float(np.min(valid_data))
            if 'max' in stats:
                feature_stats['max'] = float(np.max(valid_data))
            if 'std' in stats:
                feature_stats['std'] = float(np.std(valid_data))
        
        results.append(feature_stats)
    
    layer.ResetReading()
    raster_ds = None
    vector_ds = None
    
    return results

# 使用示例
# stats = zonal_statistics('dem.tif', 'zones.shp')

4.5.3 地形分析

from osgeo import gdal
import numpy as np

def calculate_slope(dem_path, output_path):
    """计算坡度"""
    
    ds = gdal.DEMProcessing(
        output_path,
        dem_path,
        'slope',
        format='GTiff',
        creationOptions=['COMPRESS=LZW']
    )
    ds = None

def calculate_aspect(dem_path, output_path):
    """计算坡向"""
    
    ds = gdal.DEMProcessing(
        output_path,
        dem_path,
        'aspect',
        format='GTiff',
        creationOptions=['COMPRESS=LZW']
    )
    ds = None

def calculate_hillshade(dem_path, output_path, 
                        azimuth=315, altitude=45):
    """计算山体阴影"""
    
    ds = gdal.DEMProcessing(
        output_path,
        dem_path,
        'hillshade',
        azimuth=azimuth,
        altitude=altitude,
        format='GTiff',
        creationOptions=['COMPRESS=LZW']
    )
    ds = None

def calculate_tpi(dem_path, output_path):
    """计算地形位置指数"""
    
    ds = gdal.DEMProcessing(
        output_path,
        dem_path,
        'TPI',
        format='GTiff',
        creationOptions=['COMPRESS=LZW']
    )
    ds = None

def calculate_tri(dem_path, output_path):
    """计算地形粗糙度指数"""
    
    ds = gdal.DEMProcessing(
        output_path,
        dem_path,
        'TRI',
        format='GTiff',
        creationOptions=['COMPRESS=LZW']
    )
    ds = None

def calculate_roughness(dem_path, output_path):
    """计算地形粗糙度"""
    
    ds = gdal.DEMProcessing(
        output_path,
        dem_path,
        'roughness',
        format='GTiff',
        creationOptions=['COMPRESS=LZW']
    )
    ds = None

# 使用示例
calculate_slope('dem.tif', 'slope.tif')
calculate_aspect('dem.tif', 'aspect.tif')
calculate_hillshade('dem.tif', 'hillshade.tif')

4.5.4 等值线生成

from osgeo import gdal, ogr

def generate_contours(dem_path, output_path, interval=100, base=0):
    """生成等值线"""
    
    dem_ds = gdal.Open(dem_path)
    band = dem_ds.GetRasterBand(1)
    
    # 创建输出
    driver = ogr.GetDriverByName('ESRI Shapefile')
    dst_ds = driver.CreateDataSource(output_path)
    dst_layer = dst_ds.CreateLayer('contour', geom_type=ogr.wkbLineString)
    
    # 添加属性字段
    field_defn = ogr.FieldDefn('elevation', ogr.OFTReal)
    dst_layer.CreateField(field_defn)
    
    # 生成等值线
    gdal.ContourGenerate(
        band,           # 输入波段
        interval,       # 等值线间隔
        base,           # 基础高程
        [],             # 固定高程列表（空）
        0,              # 使用无效值
        band.GetNoDataValue() or 0,  # 无效值
        dst_layer,      # 输出图层
        -1,             # ID 字段索引（-1 不创建）
        0               # 高程字段索引
    )
    
    dem_ds = None
    dst_ds = None

# 使用示例
generate_contours('dem.tif', 'contours.shp', interval=50)

4.6 栅格格式转换

4.6.1 使用 Translate

from osgeo import gdal

def convert_raster(src_path, dst_path, format='GTiff', options=None):
    """转换栅格格式"""
    
    if options is None:
        options = []
    
    translate_options = gdal.TranslateOptions(
        format=format,
        creationOptions=options
    )
    
    ds = gdal.Translate(dst_path, src_path, options=translate_options)
    ds = None

# GeoTIFF 转 JPEG
convert_raster('input.tif', 'output.jpg', 'JPEG', ['QUALITY=90'])

# GeoTIFF 转 PNG
convert_raster('input.tif', 'output.png', 'PNG')

# GeoTIFF 转 Cloud Optimized GeoTIFF
convert_raster('input.tif', 'output_cog.tif', 'COG', 
               ['COMPRESS=LZW', 'OVERVIEW_RESAMPLING=AVERAGE'])

# 转换并更改数据类型
def convert_with_type_change(src_path, dst_path, output_type=gdal.GDT_Byte):
    """转换并更改数据类型"""
    
    options = gdal.TranslateOptions(
        format='GTiff',
        outputType=output_type,
        scaleParams=[[]],  # 自动缩放
        creationOptions=['COMPRESS=LZW']
    )
    
    ds = gdal.Translate(dst_path, src_path, options=options)
    ds = None

# Float32 转 Byte（自动缩放到0-255）
convert_with_type_change('float_dem.tif', 'byte_dem.tif', gdal.GDT_Byte)

4.6.2 批量转换

from osgeo import gdal
import os
from pathlib import Path

def batch_convert(src_dir, dst_dir, src_format='*.tif', 
                  dst_format='GTiff', dst_ext='.tif', options=None):
    """批量转换栅格格式"""
    
    Path(dst_dir).mkdir(parents=True, exist_ok=True)
    
    if options is None:
        options = ['COMPRESS=LZW']
    
    translate_options = gdal.TranslateOptions(
        format=dst_format,
        creationOptions=options
    )
    
    src_path = Path(src_dir)
    converted = 0
    
    for src_file in src_path.glob(src_format):
        dst_file = Path(dst_dir) / (src_file.stem + dst_ext)
        
        print(f"转换: {src_file} -> {dst_file}")
        
        ds = gdal.Translate(str(dst_file), str(src_file), 
                           options=translate_options)
        ds = None
        converted += 1
    
    print(f"共转换 {converted} 个文件")

# 批量转换示例
batch_convert('input_images', 'output_images', 
              src_format='*.img', dst_format='GTiff')

4.7 高级创建选项

4.7.1 GeoTIFF 创建选项

from osgeo import gdal

# 完整的 GeoTIFF 创建选项
geotiff_options = [
    # 压缩选项
    'COMPRESS=LZW',      # 压缩方法: NONE, LZW, JPEG, DEFLATE, PACKBITS, ZSTD, WEBP
    'PREDICTOR=2',       # 预测器: 1(无), 2(水平差分), 3(浮点)
    'ZLEVEL=6',          # DEFLATE/ZSTD 压缩级别: 1-9
    
    # 分块选项
    'TILED=YES',         # 使用分块存储
    'BLOCKXSIZE=256',    # 块宽度
    'BLOCKYSIZE=256',    # 块高度
    
    # 大文件选项
    'BIGTIFF=IF_SAFER',  # BigTIFF: YES, NO, IF_NEEDED, IF_SAFER
    
    # 颜色选项
    'PHOTOMETRIC=RGB',   # 色彩解释: MINISBLACK, RGB, YCBCR, PALETTE
    
    # 其他选项
    'PROFILE=GeoTIFF',   # 配置文件: GeoTIFF, BASELINE, GDALGeoTIFF
    'INTERLEAVE=PIXEL',  # 交织方式: BAND, PIXEL
    'ALPHA=YES',         # 添加 Alpha 通道
    'NBITS=8',           # 位深度（某些数据类型）
]

# Cloud Optimized GeoTIFF 选项
cog_options = [
    'COMPRESS=LZW',
    'TILED=YES',
    'BLOCKSIZE=512',
    'COPY_SRC_OVERVIEWS=YES',
    'OVERVIEWS=AUTO',
]

# 创建高压缩率 GeoTIFF
def create_compressed_tiff(src_path, dst_path):
    """创建高压缩率 GeoTIFF"""
    
    options = gdal.TranslateOptions(
        format='GTiff',
        creationOptions=[
            'COMPRESS=DEFLATE',
            'PREDICTOR=2',
            'ZLEVEL=9',
            'TILED=YES',
            'BLOCKXSIZE=512',
            'BLOCKYSIZE=512',
        ]
    )
    
    ds = gdal.Translate(dst_path, src_path, options=options)
    ds = None

4.7.2 其他格式选项

from osgeo import gdal

# JPEG 选项
jpeg_options = [
    'QUALITY=85',        # 质量: 0-100
    'WORLDFILE=YES',     # 创建 world 文件
]

# PNG 选项
png_options = [
    'ZLEVEL=9',          # 压缩级别: 1-9
    'WORLDFILE=YES',
]

# JPEG 2000 选项
jp2_options = [
    'QUALITY=25',        # 质量因子
    'REVERSIBLE=NO',     # 是否无损
    'NBITS=8',           # 位深度
]

# GeoPackage 选项
gpkg_options = [
    'TILE_FORMAT=PNG',   # 瓦片格式: PNG, PNG8, JPEG, WEBP
    'QUALITY=75',        # JPEG 质量
    'BLOCKSIZE=256',     # 瓦片大小
]

4.8 本章小结

本章详细介绍了 GDAL 栅格数据处理的基础知识：

  1. 栅格数据概述：数据结构、常见类型和格式
  2. 读取栅格数据：
    • 打开数据集和获取信息
    • 读取像素数据（整体/分块/子集）
    • 处理无效值
  3. 创建栅格数据：
    • Create 和 CreateCopy 方法
    • 多波段影像创建
    • 金字塔创建
  4. 栅格处理：
    • 裁剪、重投影、重采样
    • 镶嵌、波段计算
  5. 栅格分析：
    • 统计分析、分区统计
    • 地形分析（坡度、坡向、山体阴影）
    • 等值线生成
  6. 格式转换：
    • Translate 方法
    • 批量转换
  7. 高级选项：
    • GeoTIFF 创建选项
    • 其他格式选项

掌握这些基础操作是进行更复杂栅格数据处理的前提。

4.9 思考与练习

  1. 解释栅格数据的 GeoTransform 参数的作用。
  2. 如何高效读取一个 10GB 的 GeoTIFF 文件？
  3. 比较 Create 和 CreateCopy 两种创建数据集方法的区别。
  4. 编写代码计算两个 DEM 文件的高程差。
  5. 如何将一个 Float32 类型的 DEM 转换为 8 位彩色图？
  6. 实现一个函数，将多景卫星影像镶嵌为一幅完整影像。
  7. 使用 GDAL 计算指定区域的 NDVI 并统计其分布。
posted @ 2025-12-29 11:40  我才是银古  阅读(13)  评论(0)    收藏  举报