使用magpylib绘制全球磁场等值线图遇输入错误的解决方法

magpylib的getB()方法要求输入的观测点坐标数组必须满足最后一维长度为3（对应x/y/z三个坐标分量），也就是数组形状格式为(n1, n2, ..., 3)。
你当前代码里直接把三个meshgrid生成的x、y、z坐标数组以元组形式传入，numpy会自动将其转换为形状(3, 50, 50)的数组——把长度为3的坐标维度放在了最前面，完全不符合接口的维度顺序要求，因此触发输入校验报错。

只需要调整坐标数组的维度拼接方式，同步匹配磁场返回值的计算维度即可，不需要改动其他绘图逻辑：

• 用np.stack()沿最后一个维度拼接x、y、z三个网格数组，生成形状为(50, 50, 3)的合规观测点数组，再传入getB()
• 调整磁场模长的计算逻辑：接口返回的磁场数组形状和输入观测点形状完全一致，也就是(50, 50, 3)，最后一维对应Bx/By/Bz三个分量，直接沿最后一维计算模长即可，不需要手动写循环遍历分量。

import magpylib as magpy
import numpy as np
from matplotlib import pyplot as plt

# 常量与偶极子定义
R = 6371000000  # 地球半径，单位mm
mx = -0.40e27  # 磁偶极矩，单位mT·mm3 
my =  -2.39e27
mz = -7.04e27
rx = -0.01e9 # 偶极子位置，单位mm
ry = -0.28e9
rz = 0.40e9
x1 = magpy.misc.Dipole((mx, my, mz), (rx, ry, rz))

# 经纬度网格生成
xlist = np.linspace(-180, 180) # 经度
ylist = np.linspace(-90, 90) # 纬度
X, Y = np.meshgrid(xlist, ylist)

# 球面经纬度坐标转三维直角坐标
x = R * np.cos(np.deg2rad(Y)) * np.cos(np.deg2rad(X))
y = R * np.cos(np.deg2rad(Y)) * np.sin(np.deg2rad(X))
z = R * np.sin(np.deg2rad(Y))

# 维度适配修复
observers = np.stack([x, y, z], axis=-1)
B = x1.getB(observers)
b_fin = np.linalg.norm(B, axis=-1)

b_fin *= 1000
Z = b_fin

# 等值线绘制
fig = plt.figure(figsize=(6,5))
left, bottom, width, height = 0.1, 0.1, 0.8, 0.8
ax = fig.add_axes([left, bottom, width, height]) 
cp = ax.contour(X, Y, Z)
ax.clabel(cp, inline=True, fontsize=10)
ax.set_title('Geomagnetic Field Intensity Contour')
ax.set_xlabel('Longitude (°)')
ax.set_ylabel('Latitude (°)')
plt.show()

注：代码中移除了原文件里未被使用的scipy旋转模块、time、Collection导入，减少不必要的依赖加载。

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