Python绘制南极NOAA-NSIDC V4海冰浓度数据集的错误解决及投影适配方法

我来帮你解决这个问题！你遇到的TypeError: 'GeometryCollection' object is not subscriptable本质是投影不匹配导致的，咱们一步步拆解问题并修复：

错误原因分析 #

你用ccrs.SouthPolarStereo()创建了南极极射投影的坐标轴，但在pcolormesh里指定了transform=ccrs.PlateCarree()——可你的xgrid和ygrid根本不是经纬度坐标（从属性能看到是米单位的投影坐标），这就导致Cartopy在尝试把错误的坐标转换到目标投影时，生成了无法处理的几何对象，最终抛出报错。

另外后续的NameError都是因为im对象创建失败（前面的pcolormesh报错中断了），所以只要解决投影问题，这些错误自然就消失了。

解决方案步骤 #

1. 匹配数据集的原生投影 #

从你提供的数据集属性来看，它用的是NSIDC标准的南极极射投影，参数和Cartopy的SouthPolarStereo默认参数略有差异，我们需要手动指定椭球参数来对齐：

# 创建和数据集完全匹配的南极极射投影
data_proj = ccrs.SouthPolarStereo(
    true_scale_latitude=-70,
    globe=ccrs.Globe(
        semimajor_axis=6378273.0,
        semiminor_axis=6356889.449
    )
)

2. 处理无效数据值 #

数据集里的flag_values（251-255）代表无效值（极点空洞、陆地、缺失数据等），需要先把这些值替换为NaN，避免绘图异常：

sea_ice_extent = ds.cdr_seaice_conc_monthly.sel(tdim=ds.tdim[0]).values
# 替换无效值为NaN
invalid_flags = [251,252,253,254,255]
sea_ice_extent[np.isin(sea_ice_extent, invalid_flags)] = np.nan

3. 修正绘图的transform参数 #

因为xgrid和ygrid已经是data_proj投影下的坐标，所以pcolormesh不需要指定transform=ccrs.PlateCarree()，直接用坐标轴的投影即可（或者显式指定transform=data_proj）。

完整修正代码 #

import xarray as xr
import numpy as np
import cmocean
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import cartopy.feature as cfeature
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER

# 1. 读取数据集
ds = xr.open_dataset('seaice_conc_monthly_sh_202212_f17_v04r00.nc')

# 2. 创建匹配数据集的南极极射投影
data_proj = ccrs.SouthPolarStereo(
    true_scale_latitude=-70,
    globe=ccrs.Globe(
        semimajor_axis=6378273.0,
        semiminor_axis=6356889.449
    )
)

# 3. 提取并预处理海冰浓度数据
sea_ice_data = ds.cdr_seaice_conc_monthly.sel(tdim=ds.tdim[0]).values
# 替换无效标记值为NaN
invalid_flags = [251,252,253,254,255]
sea_ice_data[np.isin(sea_ice_data, invalid_flags)] = np.nan

# 提取投影坐标
xgrid = ds.xgrid.values
ygrid = ds.ygrid.values

# 4. 创建绘图
fig = plt.figure(figsize=(10, 10))
ax = fig.add_subplot(1, 1, 1, projection=data_proj)

# 设置显示范围（用投影坐标或者经纬度范围都可以）
ax.set_extent([-3950000, 3950000, -3950000, 4350000], crs=data_proj)
# 或者用经纬度范围：ax.set_extent([-180, 180, -90, -60], ccrs.PlateCarree())

# 绘制海冰浓度，transform指定为数据的投影（和ax一致，也可以省略）
im = ax.pcolormesh(xgrid, ygrid, sea_ice_data, transform=data_proj, cmap='cmo.ice', vmin=0, vmax=100)

# 添加地图元素
ax.add_feature(cfeature.LAND, color='lightgray')
ax.coastlines(linewidth=0.8)

# 添加网格线和标签
gl = ax.gridlines(draw_labels=True, linewidth=2, color='gray', alpha=0.5, linestyle='--')
gl.top_labels = False
gl.right_labels = False
gl.xformatter = LONGITUDE_FORMATTER
gl.yformatter = LATITUDE_FORMATTER

# 添加色标
cbar = fig.colorbar(im, ax=ax, shrink=0.4)
cbar.ax.set_ylabel('Sea Ice Concentration (%)', fontsize=16)

plt.show()

额外说明 #

• 如果你想用默认的ccrs.SouthPolarStereo()，也可以尝试，但手动匹配数据集的椭球参数能避免坐标偏移问题；
• 预处理无效值很重要，不然极点空洞、陆地这些区域会干扰绘图效果；
• 添加vmin=0, vmax=100可以让色标更贴合海冰浓度的实际范围（0-100%）。

备注：内容来源于stack exchange，提问作者asheef_ik