基于Python和Matplotlib优化/平滑DEM地形3D可视化效果

问题背景 #

• 使用Matplotlib绘制DEM的3D地形时，存在内存占用过高、地形表面过于尖锐的问题
• 已完成基础工作：
  • 成功读取DEM GeoTIFF数据，2D可视化效果正常
  • 参考在线代码实现3D绘制，但生成的3D表面呈方格状，地形尖锐
  • 尝试调整Z轴比例未解决问题，可接受牺牲少量精度来平滑地形

优化方案及实现 #

通过数据平滑处理、优化绘图参数等方式，实现了更贴近真实地形的3D可视化，同时控制内存占用。以下是优化后的代码及辅助类：

优化后的核心代码 #

# 导入必要库
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
from osgeo import gdal
from scipy.ndimage import gaussian_filter

# 读取DEM数据
def read_dem(file_path):
    ds = gdal.Open(file_path)
    band = ds.GetRasterBand(1)
    dem_data = band.ReadAsArray()
    # 获取地理范围
    transform = ds.GetGeoTransform()
    x = np.arange(transform[0], transform[0] + ds.RasterXSize * transform[1], transform[1])
    y = np.arange(transform[3], transform[3] + ds.RasterYSize * transform[5], transform[5])
    x, y = np.meshgrid(x, y)
    return x, y, dem_data

# 平滑DEM数据
def smooth_dem(dem_data, sigma=2):
    # 使用高斯滤波平滑地形，sigma控制平滑程度
    return gaussian_filter(dem_data, sigma=sigma)

# 3D可视化函数
def plot_3d_dem(x, y, dem_data, z_scale=100):
    fig = plt.figure(figsize=(12, 8))
    ax = fig.add_subplot(111, projection='3d')
    
    # 降低采样密度减少内存占用（可选，根据数据大小调整）
    step = 2
    x_sampled = x[::step, ::step]
    y_sampled = y[::step, ::step]
    dem_sampled = dem_data[::step, ::step]
    
    # 保持地理投影比例
    ax.set_box_aspect([np.ptp(x), np.ptp(y), np.ptp(dem_data)/z_scale])
    
    # 绘制3D表面
    surf = ax.plot_surface(x_sampled, y_sampled, dem_sampled, cmap=cm.terrain,
                           linewidth=0, antialiased=True)
    
    # 添加颜色条
    fig.colorbar(surf, shrink=0.5, aspect=5)
    ax.set_xlabel('Longitude')
    ax.set_ylabel('Latitude')
    ax.set_zlabel('Elevation (m)')
    plt.title('Smoothed 3D DEM Visualization')
    plt.show()

# 主程序
if __name__ == "__main__":
    dem_file = "your_dem_file.tif"
    x, y, dem_data = read_dem(dem_file)
    smoothed_dem = smooth_dem(dem_data)
    plot_3d_dem(x, y, smoothed_dem)

3D箭头绘制辅助类 #

class Arrow3D:
    def __init__(self, xs, ys, zs, *args, **kwargs):
        self.xs = xs
        self.ys = ys
        self.zs = zs
        self.args = args
        self.kwargs = kwargs

    def draw(self, renderer):
        from matplotlib.patches import FancyArrowPatch
        arrow = FancyArrowPatch((self.xs[0], self.ys[0]), (self.xs[1], self.ys[1]),
                               *self.args, **self.kwargs)
        arrow.set_zorder(10)
        arrow.draw(renderer)

# 使用示例
# arrow = Arrow3D([x0, x1], [y0, y1], [z0, z1], mutation_scale=20, lw=1, arrowstyle="-|>", color="red")
# ax.add_artist(arrow)

待改进方向建议 #

• 内存优化：可尝试使用dask.array处理超大规模DEM数据，实现分块加载与计算，避免一次性读入内存
• 平滑效果精细化：可结合地形特征（如山脊、山谷）进行自适应平滑，避免过度平坦化关键地形
• 渲染效率提升：对于超精细地形，可考虑切换到plotly或pyvista等更适合3D地理数据的可视化库，提升交互性与渲染速度
• 投影一致性：若需严格保持地理投影精度，可结合cartopy库实现带投影的3D可视化，确保地形的空间位置准确

内容的提问来源于stack exchange，提问作者lilaaffe