使用Basemap绘制NC文件时遭遇IndexError问题求助

问题背景 #

处理IASI卫星大气N₂O气体辐射数据的NC文件，使用netCDF4、numpy、matplotlib结合Basemap绘图时，调用m.contour触发IndexError，提示数组维度不匹配。当前X（经度）、Y（纬度）是含235586个元素的1D数组，Z（n2o_retrieval）是235586×14的2D数组，调整维度后仍无法解决。

报错信息 #

Traceback (most recent call last):
File "C:\Users\Lucas\Documents\Projet N2O\Plot.py", line 70, in 
cs = m.contour(latitude,longitude,n2O_r,1000,linewidths=1.5,cmap=viridis, 
File "C:\Users\Lucas\AppData\Local\Programs\Python\Python38\lib\site-packages\mpl_toolkits\basemap\__init__.py", line 549, in with_transform
return plotfunc(self,x,y,data,*args,**kwargs)
File "C:\Users\Lucas\AppData\Local\Programs\Python\Python38\lib\site-packages\mpl_toolkits\basemap\__init__.py", line 3570, in contour
xx = x[x.shape[0]//2,:]
IndexError: too many indices for array: array is 1-dimensional, but 2 were indexed

关联代码 #

import netCDF4 as nc
import numpy as np
import matplotlib.colors as mcolors
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap, cm, shiftgrid, addcyclic # mapping

# 打开NC文件
filePath = 'C:\\Users\\Lucas\\Documents\\ProjetN2O\\iasi_ret_L2_n2o_ch4_2019_05_02.nc'
ds = nc.Dataset(filePath)
print(ds)

# 打印数据集的维度和变量
for dimension in ds.dimensions.values():
    print(dimension)
for var in ds.variables.values():
    print(var)

# 提取变量为numpy数组
temps = np.array(ds.variables["time"])
latitude = np.array(ds.variables["lat"])
longitude = np.array(ds.variables["lon"])
n2O_r = np.array(ds.variables["n2o_retrieval"])
n2O_avk = np.array(ds.variables["n2o_AVK"])
n2O_ap = np.array(ds.variables["n2o_apriori"])
SurfP = np.array(ds.variables["Surf_P"])
SurfPid = np.array(ds.variables["Surf_P_id"])

# 选取n2o_retrieval的第一层数据
A = n2O_r[0:235587, 0:1] # 取第一列所有值
B = A.ravel() # 转为1D数组

# 绘图参数设置
lat0 = 0; lon0 = 0 
lon_max = longitude.max(); lon_min = longitude.min()
lat_max =  latitude.max(); lat_min = latitude.min()

# 初始化Basemap
m = Basemap(projection='cyl',lat_0=lat0,lon_0=lon0,resolution='c',
llcrnrlat=lat_min,urcrnrlat=lat_max,
llcrnrlon=lon_min,urcrnrlon=lon_max)
m.drawcoastlines(linewidth=1.2, linestyle='solid', color='k', antialiased=1, zorder=2)
m.drawcountries()
m.drawlsmask(land_color='none', ocean_color='aqua', zorder=1)
viridis =plt.get_cmap('viridis', 12)

# 触发报错的代码行
cs = m.contour(latitude,longitude,B,1000,linewidths=1.5,cmap=viridis, colors='b', alpha=0.3)

问题原因 #

Basemap的contour函数要求X、Y必须是2D网格数组（例如通过np.meshgrid生成的N×M形状数组），且Z需与X/Y的网格维度完全匹配。但IASI的L2数据是离散的观测点，经纬度是1D数组，对应的N₂O数据也是1D的离散值，直接传入contour会触发维度不匹配的错误。

解决方案 #

方案1：使用散点图绘制离散观测点（推荐） #

卫星L2数据本身是离散的观测记录，用散点图能更准确展示数据分布，替换contour为scatter：

# 替换原有报错的contour代码
# 转换经纬度到地图投影坐标
x, y = m(longitude, latitude)
# 绘制散点图，s控制点大小，c指定颜色映射变量
sc = m.scatter(x, y, c=B, cmap=viridis, alpha=0.3, s=1)
# 添加颜色条
plt.colorbar(sc, label='N₂O 反演值')
plt.title('IASI N₂O 反演结果（第一层）')
plt.show()

方案2：插值到规则网格后绘制等值线 #

如果需要生成平滑的等值线图，需先将离散点插值到规则网格上，使用scipy.interpolate.griddata实现：

# 导入插值模块
from scipy.interpolate import griddata

# 生成规则网格的经纬度
lon_grid = np.linspace(lon_min, lon_max, 100)  # 100为网格密度，可调整
lat_grid = np.linspace(lat_min, lat_max, 100)
lon_mesh, lat_mesh = np.meshgrid(lon_grid, lat_grid)

# 将离散数据插值到规则网格（method可选linear/nearest/cubic）
z_grid = griddata((longitude, latitude), B, (lon_mesh, lat_mesh), method='linear')

# 转换网格坐标到地图投影
x_mesh, y_mesh = m(lon_mesh, lat_mesh)

# 绘制等值线
cs = m.contour(x_mesh, y_mesh, z_grid, 10, linewidths=1.5, cmap=viridis, alpha=0.3)
plt.colorbar(cs, label='N₂O 反演值')
plt.title('IASI N₂O 反演等值线（第一层）')
plt.show()

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