Basemap在3D轴加载地形图层时报无法设置aspect错误求助

报错原因 #

Basemap的bluemarble()、shadedrelief()、etopo()这类栅格底图绘制方法，内部会强制调用坐标轴的aspect比例设置逻辑，而matplotlib 2.2.3版本的3D坐标轴不支持手动调整aspect，即使你传入了fix_aspect=False也无法绕过方法内部的aspect设置调用，因此会触发报错。

可行解决方案 #

绕开Basemap自带的栅格底图绘制接口，手动将底图导出为纹理图片，再贴到3D坐标系的z=0地表平面上，修改后的代码如下：

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from mpl_toolkits.basemap import Basemap
import numpy as np
from matplotlib.image import AxesImage

# 读取射线数据部分保留原有逻辑
f = open("best-working4D.ray_paths", 'r') 
data = f.read()
lines = data.split('\n        ')

fig = plt.figure()
ax = plt.axes(projection='3d')

extent = [300, 360, 50, 75]
bm = Basemap(llcrnrlon=extent[0], llcrnrlat=extent[2],
             urcrnrlon=extent[1], urcrnrlat=extent[3], resolution='l', fix_aspect=False)

# ----------------------修改部分：手动获取底图纹理----------------------
# 新建临时2D画布绘制底图，导出为RGBA数组
fig_tmp, ax_tmp = plt.subplots()
bm.bluemarble(ax=ax_tmp)
# 从临时画布提取底图数组
for im in ax_tmp.get_children():
    if isinstance(im, AxesImage):
        bg_img = im.get_array()
        break
plt.close(fig_tmp) # 关闭临时画布

# 生成地表平面网格，z设置为0（可根据自身坐标体系调整地表高度）
xx, yy = np.meshgrid(np.linspace(extent[0], extent[1], bg_img.shape[1]),
                     np.linspace(extent[2], extent[3], bg_img.shape[0]))
zz = np.zeros_like(xx)

# 将底图纹理贴到3D地表平面
ax.plot_surface(xx, yy, zz, rstride=1, cstride=1, facecolors=bg_img/255, shade=False)
# --------------------------------------------------------------------

# 射线绘制部分保留原有逻辑
for i in range(1, 119):
    wave = lines[i]
    j = wave.split('\n')
    k = []
    for item in j:
        k.append(item.split())
    x=[]
    y=[]
    z=[]
    for m in k[1:]:
        x.append(m[0])
        y.append(m[1])
        z.append(m[2])
    x= np.array(x).astype('float32')
    y= np.array(y).astype('float32')
    z= np.array(z).astype('float32')
    ax.plot3D(x,y,z, color='red')

# 对流层顶绘制保留
T_hi = 20
xx_trop, yy_trop = np.meshgrid(range(300,360), range(50,75))
zz_trop = yy_trop*0 + T_hi
ax.plot_surface(xx_trop, yy_trop, zz_trop, alpha=0.15)

ax.set_xlabel("Latitude [deg]")
ax.set_ylabel("Longitude [deg]")
ax.set_zlabel("Altitude [km]")
ax.add_collection3d(bm.drawcoastlines(linewidth=0.25))

plt.show()

补充说明 #

• 如果你需要加载shaded relief或者topo图层，只需要把临时画布中的bm.bluemarble(ax=ax_tmp)替换为bm.shadedrelief(ax=ax_tmp)或者bm.etopo(ax=ax_tmp)即可
• 该方法完全绕开了Basemap底图绘制对3D轴aspect的设置逻辑，不会触发对应报错，同时能保留底图的显示效果

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