Python绘制两TIFF高程差报错：TypeError: 图像数据形状无效

问题原因 #

• 你创建的dh_list是一个列表，里面包含了一个2D高程差数组。当把这个列表传给ax.imshow()时，它会被识别为3D数组（shape为(1, 2355, 2676)），而imshow()对于单通道图像仅接受2D数组（(高度, 宽度)），因此触发形状不匹配错误。

修复方法 #

直接使用计算得到的高程差数组，无需将其包装在列表中，修改两处代码即可：

修改后的核心代码片段 #

# 计算高程差，去掉列表包装
dh = dem_new - dem_old

# 绘图时传入dh这个2D数组
im = ax.imshow(dh, origin='upper', interpolation='None', cmap='seismic_r', vmin = -40, vmax = 40)

完整修复代码 #

# import relevant moduls
from osgeo import gdal
import numpy as np
import matplotlib.pyplot as plt
from pygeotools.lib import iolib, warplib, geolib, timelib, malib

# Specify input files
dem_old = '../data/2007.tif'
dem_new = '../data/2018.tif'
# generate a list of the files
dem_fn_list = [dem_old, dem_new]

# combine dems with warplib and set to extent of intersection
ds_list = warplib.memwarp_multi_fn(dem_fn_list, extent='intersection', res='min', t_srs=dem_new)

# Load datasets to NumPy masked arrays
dem_old, dem_new = [iolib.ds_getma(i) for i in ds_list]

# Calculate elevation difference for each time period
dh = dem_new - dem_old

# plot elevation difference
print('Plotting DEM difference...')
plt.rcParams.update({'font.size': 15})
fig = plt.figure(figsize=(12, 12))
ax = fig.add_subplot(1,1,1)
im = ax.imshow(dh, origin='upper', interpolation='None', cmap='seismic_r', vmin = -40, vmax = 40)
# add colorbar
fig.colorbar(im)
# add title
plt.title("DSM elevation differences")
# disable xticks
plt.xticks([]) 
plt.yticks([])
# save figure
plt.savefig('glacier_change.png')
# show figure
plt.show()

额外说明 #

如果后续需要处理多组高程差（比如多个时间对），可以保留列表，但绘图时需取出列表中的单个2D数组，示例：

dh_list = [dem_new - dem_old, dem_2020 - dem_2018]
# 绘制第一组差值
ax.imshow(dh_list[0], ...)

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