基于.tiff文件绘制3D表面图为何显示异常？

我尝试用matplotlib将两个.tiff格式的DEM文件（Nisqually_1987.tif和Nisqually_2017.tif）可视化为3D表面图，但生成的图像显示异常。实现代码如下：

# import relevant moduls
from osgeo import gdal
import numpy as np
import matplotlib.pyplot as plt

# Specify input files
dem_old = '../data/Nisqually_1987.tif'
dem_new = '../data/Nisqually_2017.tif'

# Creating first 3D figure
# set up a figure twice as wide as it is tall
fig = plt.figure(figsize=plt.figaspect(0.5))

# Old DEM plot
ax = fig.add_subplot(1, 2, 1, projection='3d')
dem = gdal.Open(dem_old)
dem_array = dem.ReadAsArray()

lin_x = np.linspace(0,1,dem_array.shape[0],endpoint=False)
lin_y = np.linspace(0,1,dem_array.shape[1],endpoint=False)
y,x = np.meshgrid(lin_y,lin_x)
z = dem_array
surf = ax.plot_surface(x,y,z,cmap='terrain', edgecolor='none')
fig.colorbar(surf, shrink=0.5, aspect=15)
ax.set_title('3D Nisqually-Gletscher Oberfläche (1987)')
plt.xticks([])
plt.yticks([])

# New DEM 3D plot
ax = fig.add_subplot(1, 2, 2, projection='3d')
dem_2 = gdal.Open(dem_new)
dem_array_2 = dem_2.ReadAsArray()

lin_x_2 = np.linspace(0,1,dem_array_2.shape[0],endpoint=False)
lin_y_2 = np.linspace(0,1,dem_array_2.shape[1],endpoint=False)
y,x = np.meshgrid(lin_y_2,lin_x_2)
z = dem_array_2
surf = ax.plot_surface(x,y,z,cmap='terrain', edgecolor='none')
fig.colorbar(surf, shrink=0.5, aspect=10)
ax.set_title('3D Nisqually-Gletscher Oberfläche (2017)')
plt.xticks([])
plt.yticks([])

# show plot
plt.savefig ('glaciers_3D.png')
plt.show()

异常效果参考：

数据来源：ScienceBase平台的Nisqually冰川DEM数据集

异常原因分析 #

从代码和异常图来看，核心问题有三个：

• 网格维度不匹配：plot_surface要求x、y、z的维度严格对应，但你把DEM的行数（y轴方向）对应给了x轴的线性空间，列数（x轴方向）对应给了y轴，再加上meshgrid的参数顺序错误，直接导致网格和高程数据错位，地形扭曲。
• 未处理NoData值：DEM数据里通常存在无效的NoData标记（比如-9999），这些值会被当成正常高程绘制，产生突兀的尖刺或凹陷，破坏地形形态。
• 坐标与视角问题：用0-1的归一化坐标没有体现DEM的真实地理范围，再加上默认3D视角的比例失衡，让地形起伏看起来变形严重。

修正后的代码 #

from osgeo import gdal
import numpy as np
import matplotlib.pyplot as plt

# 读取DEM并处理NoData、生成地理坐标网格
def read_dem(file_path):
    dem = gdal.Open(file_path)
    band = dem.GetRasterBand(1)
    dem_array = band.ReadAsArray()
    # 替换NoData为NaN，matplotlib会自动忽略这些值
    nodata = band.GetNoDataValue()
    if nodata is not None:
        dem_array[dem_array == nodata] = np.nan
    # 从地理变换参数获取真实坐标范围
    geotrans = dem.GetGeoTransform()
    x_min, x_res, _, y_max, _, y_res = geotrans
    x_size, y_size = dem.RasterXSize, dem.RasterYSize
    # 创建匹配DEM维度的地理坐标网格
    x = np.linspace(x_min, x_min + x_res * x_size, x_size)
    y = np.linspace(y_max, y_max + y_res * y_size, y_size)
    y, x = np.meshgrid(y, x)
    return x, y, dem_array

# 输入文件路径
dem_old = '../data/Nisqually_1987.tif'
dem_new = '../data/Nisqually_2017.tif'

# 创建3D画布
fig = plt.figure(figsize=plt.figaspect(0.5))

# 1987年DEM绘图
ax = fig.add_subplot(1, 2, 1, projection='3d')
x_old, y_old, z_old = read_dem(dem_old)
# rstride和cstride控制网格密度，减少绘制压力同时让地形更平滑
surf = ax.plot_surface(x_old, y_old, z_old, cmap='terrain', edgecolor='none', rstride=10, cstride=10)
fig.colorbar(surf, shrink=0.5, aspect=15)
ax.set_title('1987年尼斯奎利冰川表面3D图')
ax.set_xlabel('X坐标')
ax.set_ylabel('Y坐标')
ax.set_zlabel('高程')
# 设置统一视角，方便对比两年冰川形态
ax.view_init(elev=30, azim=45)

# 2017年DEM绘图
ax = fig.add_subplot(1, 2, 2, projection='3d')
x_new, y_new, z_new = read_dem(dem_new)
surf = ax.plot_surface(x_new, y_new, z_new, cmap='terrain', edgecolor='none', rstride=10, cstride=10)
fig.colorbar(surf, shrink=0.5, aspect=10)
ax.set_title('2017年尼斯奎利冰川表面3D图')
ax.set_xlabel('X坐标')
ax.set_ylabel('Y坐标')
ax.set_zlabel('高程')
ax.view_init(elev=30, azim=45)

plt.tight_layout()
plt.savefig('glaciers_3D.png')
plt.show()

关键修正说明 #

• 匹配网格维度：基于DEM的地理变换参数生成真实坐标网格，确保x、y、z的维度完全对应，解决地形扭曲问题。
• 清理无效值：把NoData替换为np.nan，避免异常凸起破坏地形形态。
• 优化绘图参数：添加网格步长参数提升绘制效率和平滑度，设置统一视角方便对比，补充坐标轴标签让可视化更直观。

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