Python进行基本地形分析

Python 进行地形分析的基本步骤

地形分析是通过数字高程模型（DEM）数据计算地形特征（如坡度、坡向、高程等）的过程。Python 中常用的地形分析工具包括 rasterio（用于读取 DEM 数据）和 numpy（用于数组操作）。以下是地形分析的基本步骤：

地形分析的基本步骤

读取 DEM 数据：使用 rasterio 读取 DEM 文件。
计算坡度：通过 DEM 数据的梯度计算每个像元的坡度（单位：度）。
计算坡向：通过 DEM 数据的梯度计算每个像元的坡向（单位：度）。
可视化结果：使用 matplotlib 绘制高程图、坡度图和坡向图。
保存结果：将可视化结果保存为 PNG 文件。
统计信息：计算并打印 DEM 和坡度的统计信息（如最小值、最大值、平均值、标准差）。

示例代码注释说明

import rasterio  # 导入 rasterio 库，用于读取 DEM 数据
import numpy as np  # 导入 numpy 库，用于数组操作
import matplotlib.pyplot as plt  # 导入 matplotlib 库，用于绘图
from matplotlib.colors import LightSource  # 导入 LightSource 类，用于光照效果
import os  # 导入 os 库，用于文件路径操作

def calculate_slope_aspect(dem, cell_size):
    """
    计算坡度和坡向
    :param dem: 数字高程模型（DEM）数据
    :param cell_size: 像元大小（单位：米）
    :return: 坡度（slope）和坡向（aspect）
    """
    # 计算 DEM 的梯度
    dy, dx = np.gradient(dem, cell_size)
    
    # 计算坡度（单位：度）
    slope = np.degrees(np.arctan(np.sqrt(dx**2 + dy**2)))
    
    # 计算坡向（单位：度）
    aspect = np.degrees(np.arctan2(dy, dx))
    aspect = np.where(aspect < 0, aspect + 360, aspect)  # 将负值转换为 0-360 度
    
    return slope, aspect

def analyze_dem(dem_file):
    """
    分析 DEM 数据
    :param dem_file: DEM 文件路径
    """
    print("开始DEM数据分析...")
    
    # 打开 DEM 文件
    with rasterio.open(dem_file) as src:
        dem = src.read(1)  # 读取第一个波段的数据
        transform = src.transform  # 获取变换矩阵
        cell_size = abs(transform[0])  # 计算像元大小
    
    # 计算坡度和坡向
    slope, aspect = calculate_slope_aspect(dem, cell_size)
    
    # 创建绘图窗口
    fig = plt.figure(figsize=(15, 5))
    
    # 绘制高程图
    ax1 = fig.add_subplot(131)
    ls = LightSource(azdeg=315, altdeg=45)  # 设置光源方向
    dem_plot = ax1.imshow(dem, cmap=plt.cm.terrain)  # 绘制 DEM
    rgb = ls.shade(dem, cmap=plt.cm.terrain, blend_mode='soft')  # 添加光照效果
    ax1.imshow(rgb)
    ax1.set_title('高程图')
    plt.colorbar(dem_plot, ax=ax1, label='高程 (m)')
    ax1.axis('off')
    
    # 绘制坡度图
    ax2 = fig.add_subplot(132)
    im2 = ax2.imshow(slope, cmap='YlOrRd')  # 使用黄-橙-红色彩映射
    ax2.set_title('坡度图')
    plt.colorbar(im2, ax=ax2, label='坡度 (°)')
    ax2.axis('off')
    
    # 绘制坡向图
    ax3 = fig.add_subplot(133)
    im3 = ax3.imshow(aspect, cmap='hsv')  # 使用 HSV 色彩映射
    ax3.set_title('坡向图')
    plt.colorbar(im3, ax=ax3, label='坡向 (°)')
    ax3.axis('off')
    
    # 调整布局
    plt.tight_layout()
    
    # 保存图表为 PNG 文件
    output_dir = r'第六章\3\结果'
    os.makedirs(output_dir, exist_ok=True)  # 创建输出目录（如果不存在）
    output_file = os.path.join(output_dir, 'DEM分析结果.png')
    plt.savefig(output_file, dpi=300, bbox_inches='tight')
    print(f"分析结果已保存至：{output_file}")
    
    # 显示图表
    plt.show()
    
    # 计算 DEM 和坡度的统计信息
    stats = {
        '高程': {
            '最小值': np.min(dem),
            '最大值': np.max(dem),
            '平均值': np.mean(dem),
            '标准差': np.std(dem)
        },
        '坡度': {
            '最小值': np.min(slope),
            '最大值': np.max(slope),
            '平均值': np.mean(slope),
            '标准差': np.std(slope)
        }
    }
    
    # 打印统计信息
    print("\nDEM统计信息：")
    for key, value in stats.items():
        print(f"\n{key}统计：")
        for stat_name, stat_value in value.items():
            print(f"{stat_name}: {stat_value:.2f}")

if __name__ == "__main__":
    # 设置 matplotlib 的中文字体和符号显示
    plt.rcParams['font.sans-serif'] = ['SimHei']
    plt.rcParams['axes.unicode_minus'] = False
    
    # 定义 DEM 文件路径
    dem_file = r'<DEM的文件路径>'
    
    try:
        # 调用函数分析 DEM 数据
        analyze_dem(dem_file)
    except Exception as e:
        # 捕获并打印错误信息
        print(f"处理过程中出现错误: {str(e)}")