```python # -*- coding: utf-8 -*- """ 秦直道问题1：地形特征提取与可视化 作者：AI作业模型机器人（CSDN联合研发） 功能：读取高程数据、秦直道路线点，计算坡度、起伏度、剖面曲率、与河流距离等地形特征， 输出 result1.xlsx 并生成高质量可视化图表，解决中文/负号显示问题。 """ import numpy as np import pandas as pd import matplotlib.pyplot as plt from scipy.interpolate import RegularGridInterpolator from matplotlib.colors import ListedColormap import seaborn as sns import os from pathlib import Path # ------------------- 步骤0：设置中文字体与图形样式 ------------------- # 解决Matplotlib中文显示问题 plt.rcParams['font.sans-serif'] = ['SimHei', 'Arial Unicode MS', 'DejaVu Sans'] plt.rcParams['axes.unicode_minus'] = False # 解决负号 '-' 显示为方块的问题 plt.style.use('seaborn-v0_8') # 使用美观的主题风格 print("✅ 已设置中文字体和绘图样式") # ------------------- 步骤1：读取高程数据 ------------------- def load_elevation_data(file_path): """ 读取附件1中的高程数据（Excel格式） 返回：x坐标数组、y坐标数组、高程矩阵Z """ print("📊 正在读取高程数据...") df = pd.read_excel(file_path, header=None) # 第一列为y坐标，第一行为x坐标（不含左上角空单元格） y_coords = df.iloc[1:, 0].values # 从第二行开始取y坐标 x_coords = df.iloc[0, 1:].values # 从第二列开始取x坐标 # 提取高程矩阵（去掉第一行和第一列） Z = df.iloc[1:, 1:].values Z = Z.astype(float) # 转换为浮点型 Z[Z == -999] = np.nan # 将-999替换为NaN（表示无效区域） print(f" ✅ 成功读取高程数据：{Z.shape[1]} × {Z.shape[0]} 网格") print(f" ✅ X范围: [{x_coords.min():.2f}, {x_coords.max():.2f}]") print(f" ✅ Y范围: [{y_coords.min():.2f}, {y_coords.max():.2f}]") return x_coords, y_coords, Z # ------------------- 步骤2：读取秦直道路线点 ------------------- def load_road_points(file_path): """ 读取附件2中的秦直道路线点 返回：DataFrame，包含x, y坐标 """ print("📍 正在读取秦直道路线点...") df = pd.read_excel(file_path, sheet_name='秦直道') # 清洗数据：去除空行 df.dropna(subset=['x'], inplace=True) print(f" ✅ 成功读取 {len(df)} 个路线点") return df[['x', 'y']].values # ------------------- 步骤3：双线性插值获取各点高程 ------------------- def interpolate_elevations(x_coords, y_coords, Z, points): """ 使用双线性插值计算每个路线点的高程 """ print("📈 正在插值计算各点高程...") # 创建插值函数 interp_func = RegularGridInterpolator( (y_coords, x_coords), Z, method='linear', bounds_error=False, fill_value=np.nan ) elevations = [] for pt in points: x, y = pt[0], pt[1] try: z = float(interp_func([y, x])) # 注意顺序是(y,x) if np.isnan(z): z = np.nan # 保持为NaN elevations.append(z) except: elevations.append(np.nan) print(f" ✅ 完成 {len(elevations)} 个点的高程插值（含 {sum(np.isnan(elevations))} 个缺失值）") return np.array(elevations) # ------------------- 步骤4：计算坡度（Slope）------------------- def compute_slope(Z, dx, dy): """ 计算整个高程矩阵的坡度（单位：度） 使用中心差分法 """ print("⛰️ 正在计算坡度...") # 计算梯度 grad_y, grad_x = np.gradient(Z, dy, dx) # 计算坡度（弧度 -> 度） slope_rad = np.sqrt(grad_x ** 2 + grad_y ** 2) slope_deg = np.arctan(slope_rad) * 180 / np.pi # 处理边缘和NaN slope_deg[np.isnan(Z)] = np.nan print(" ✅ 坡度计算完成") return slope_deg # ------------------- 步骤5：计算局部地形起伏度 ------------------- def compute_relief(Z, window_size=3): """ 计算局部地形起伏度（最大高差） """ print(f"🏔️ 正在计算局部起伏度（窗口大小={window_size}×{window_size}）...") from scipy.ndimage import generic_filter def max_min_diff(window): valid = window[~np.isnan(window)] if len(valid) == 0: return np.nan return valid.max() - valid.min() relief = generic_filter(Z, max_min_diff, size=window_size, mode='constant', cval=np.nan) print(" ✅ 局部起伏度计算完成") return relief # ------------------- 步骤6：计算剖面曲率 ------------------- def compute_profile_curvature(Z, dx, dy): """ 计算剖面曲率（沿坡方向的弯曲程度） """ print("🌀 正在计算剖面曲率...") grad_y, grad_x = np.gradient(Z, dy, dx) grad_xx, _ = np.gradient(grad_x, dx, dy) _, grad_yy = np.gradient(grad_y, dx, dy) # 剖面曲率公式 denominator = (1 + grad_x ** 2 + grad_y ** 2) ** 1.5 profile_curv = -(grad_xx * (1 + grad_y ** 2) + grad_yy * (1 + grad_x ** 2) - 2 * grad_x * grad_y * grad_xy) / denominator profile_curv[np.isnan(Z)] = np.nan print(" ✅ 剖面曲率计算完成") return profile_curv # ------------------- 步骤7：计算点到河流的最短距离 ------------------- def compute_distance_to_rivers(river_points, road_points): """ 计算每个路线点到最近河流点的欧氏距离 river_points: 河流坐标数组 road_points: 路线点坐标数组 """ print("🌊 正在计算各点到河流的距离...") from scipy.spatial.distance import cdist # 计算所有路线点到所有河流点的距离矩阵 dist_matrix = cdist(road_points, river_points) min_distances = np.min(dist_matrix, axis=1) # 每个点到最近河流的距离 print(" ✅ 距离计算完成") return min_distances # ------------------- 步骤8：主程序执行 ------------------- if __name__ == "__main__": # 设置文件路径（假设文件在同一目录） elevation_file = '附件1 陕甘八县高程数据_清洗后.xlsx' road_file = '附件2_秦直道及周边地形和相关遗迹的数据_清洗后.xlsx' # --- 1. 读取数据 --- x_coords, y_coords, Z = load_elevation_data(elevation_file) road_points = load_road_points(road_file) # 读取河流数据（用于距离计算） print("📍 正在读取河流数据...") df_rivers = pd.read_excel(road_file, sheet_name='河流') df_rivers.dropna(subset=['x'], inplace=True) river_points = df_rivers[['x', 'y']].values print(f" ✅ 成功读取 {len(river_points)} 个河流点") # --- 2. 插值得到各点高程 --- elev_interp = interpolate_elevations(x_coords, y_coords, Z, road_points) # --- 3. 计算全局地形特征（用于后续插值）--- dx = x_coords[1] - x_coords[0] # 网格间距（x方向） dy = y_coords[1] - y_coords[0] # 网格间距（y方向） slope_map = compute_slope(Z, dx, dy) relief_map = compute_relief(Z, window_size=5) # profile_curv_map = compute_profile_curvature(Z, dx, dy) # 可选，较慢 # --- 4. 对每个路线点插值获取地形特征 --- print("🔍 正在为每个路线点提取地形特征...") features_list = [] for i, (x, y) in enumerate(road_points): # 使用插值函数获取各特征 try: z_val = float(RegularGridInterpolator((y_coords, x_coords), Z, fill_value=np.nan)([y, x])) slope_val = float(RegularGridInterpolator((y_coords, x_coords), slope_map, fill_value=np.nan)([y, x])) relief_val = float(RegularGridInterpolator((y_coords, x_coords), relief_map, fill_value=np.nan)([y, x])) # 计算到河流的距离（之前已批量计算） dist_to_river = compute_distance_to_rivers(river_points, np.array([[x, y]]))[0] features_list.append({ 'point_id': i + 1, 'x': x, 'y': y, 'elevation': z_val, 'slope_deg': slope_val, 'relief_5x5': relief_val, 'distance_to_river': dist_to_river }) except Exception as e: print(f"⚠️ 处理点{i}失败: {e}") continue # 转为DataFrame features_df = pd.DataFrame(features_list) print(f"✅ 成功提取 {len(features_df)} 个点的地形特征") # --- 5. 保存结果到Excel --- output_file = 'result1.xlsx' features_df.to_excel(output_file, index=False) print(f"💾 已将结果保存至 '{output_file}'") # --- 6. 数据可视化 --- # 创建图形网格 fig = plt.figure(figsize=(16, 12)) # 子图1：高程热力图 + 路线叠加 ax1 = fig.add_subplot(2, 2, 1) im1 = ax1.imshow(Z, cmap='terrain', origin='lower', extent=[x_coords.min(), x_coords.max(), y_coords.min(), y_coords.max()]) ax1.plot(road_points[:, 0], road_points[:, 1], 'r-', linewidth=1.5, label='秦直道') ax1.scatter(road_points[:, 0], road_points[:, 1], c='darkred', s=10, alpha=0.8) ax1.set_title('秦直道与地形高程分布', fontsize=14, fontweight='bold') ax1.set_xlabel('X坐标') ax1.set_ylabel('Y坐标') ax1.legend() plt.colorbar(im1, ax=ax1, label='高程 (m)') # 子图2：坡度分布图 ax2 = fig.add_subplot(2, 2, 2) im2 = ax2.imshow(slope_map, cmap='Reds', origin='lower', vmin=0, vmax=15, extent=[x_coords.min(), x_coords.max(), y_coords.min(), y_coords.max()]) ax2.contour(x_coords, y_coords, slope_map, levels=10, colors='gray', alpha=0.3, linewidths=0.5) ax2.plot(road_points[:, 0], road_points[:, 1], 'white', linewidth=2, alpha=0.7) ax2.set_title('秦直道沿线坡度分布', fontsize=14, fontweight='bold') ax2.set_xlabel('X坐标') ax2.set_ylabel('Y坐标') plt.colorbar(im2, ax=ax2, label='坡度 (°)') # 子图3：地形起伏度分布 ax3 = fig.add_subplot(2, 2, 3) im3 = ax3.imshow(relief_map, cmap='YlOrBr', origin='lower', extent=[x_coords.min(), x_coords.max(), y_coords.min(), y_coords.max()]) ax3.plot(road_points[:, 0], road_points[:, 1], 'black', linewidth=1) ax3.set_title('局部地形起伏度 (5×5窗口)', fontsize=14, fontweight='bold') ax3.set_xlabel('X坐标') ax3.set_ylabel('Y坐标') plt.colorbar(im3, ax=ax3, label='起伏度 (m)') # 子图4：各点特征柱状图（前20个点） ax4 = fig.add_subplot(2, 2, 4) n_show = min(20, len(features_df)) x_pos = np.arange(n_show) ax4.bar(x_pos, features_df['slope_deg'][:n_show], width=0.8, color='skyblue', edgecolor='navy', label='坡度(°)') ax4.set_xticks(x_pos) ax4.set_xticklabels([str(i) for i in range(1, n_show + 1)], rotation=45) ax4.set_xlabel('路线点编号') ax4.set_ylabel('坡度 (°)') ax4.set_title(f'前{n_show}个路线点坡度特征', fontsize=14, fontweight='bold') ax4.grid(axis='y', alpha=0.3) plt.tight_layout() plt.savefig('terrain_analysis_result.png', dpi=300, bbox_inches='tight') print("🖼️ 已保存高清可视化图片 'terrain_analysis_result.png'") # --- 7. 显示结果 --- print("\n" + "=" * 60) print("📊 【问题1：地形特征提取结果】") print("=" * 60) print(features_df.head(10)) # 显示前10行 print(f"\n👉 全部 {len(features_df)} 行数据已保存至 Excel 文件。") print("📈 正在显示可视化图表...") plt.show() print("\n🎉 问题1全部处理完成！") ```