如何将STILT模拟的49层footprints时间维度累加生成单张2D地图？

我有一个包含lon、lat、time和foot（代表footprints）的nc文件，foot变量维度为(lon, lat, time)（注：从代码索引逻辑看实际可能为(time, lon, lat)）。该文件来自48小时前向轨迹STILT模拟，time维度包含49个切片。我希望对时间维度进行累加以生成一张2D地图，但当前代码会生成49张地图，不知如何解决，以下是我的代码：

footprints = "myfile.nc"
print(footprints)
 
with Dataset(footprints) as root:
    time = root.variables['time'][:]
    dates = num2date(time, root.variables['time'].units)
    print(dates[0].strftime('%Y-%m-%d-%H-%M-%S'))
    lats= root.variables['lat'][:]
    lons= root.variables['lon'][:]
    foot= root.variables['foot'][:, :, :]
 
Foot = np.arange(0,49)
for i in Foot:
   ax = plt.axes(projection=ccrs.PlateCarree())
   plt.contourf(lons[:], lats[:], np.squeeze(foot[i,:,:]), cmap = 'YlOrBr', transform=ccrs.PlateCarree())
   ax.set_extent([-178.546, -130.946, 50, 75])   # regional map, x=longitude, y=latitude (xmin, xmax, ymin, ymax)  
   ax.coastlines()
   ax.gridlines(draw_labels=True)
   plt.show()

我不想仅绘制单个切片，而是将所有切片合并为一张地图，现在却得到了49张地图，请帮忙解决。

解决方案 #

核心思路是对foot变量的时间维度进行求和累加，而非循环每个时间步单独绘图。以下提供两种实现方式：

方式1：基于现有netCDF4代码修改 #

从代码中foot[i,:,:]的索引逻辑可以看出，time是foot的第0个维度（数据形状为(49, lon, lat)），直接对第0轴求和即可得到累加后的2D足迹数据：

import numpy as np
import matplotlib.pyplot as plt
from netCDF4 import Dataset, num2date
import cartopy.crs as ccrs

footprints = "myfile.nc"
print(footprints)
 
with Dataset(footprints) as root:
    time = root.variables['time'][:]
    dates = num2date(time, root.variables['time'].units)
    print(dates[0].strftime('%Y-%m-%d-%H-%M-%S'))
    lats = root.variables['lat'][:]
    lons = root.variables['lon'][:]
    foot = root.variables['foot'][:, :, :]

# 对时间维度（第0轴）求和，得到累加后的2D数据
foot_total = np.sum(foot, axis=0)

# 仅绘制一次累加结果
ax = plt.axes(projection=ccrs.PlateCarree())
plt.contourf(lons, lats, foot_total, cmap='YlOrBr', transform=ccrs.PlateCarree())
ax.set_extent([-178.546, -130.946, 50, 75])
ax.coastlines()
ax.gridlines(draw_labels=True)
plt.colorbar(label='Total Footprints')  # 可选：添加颜色条说明
plt.show()

如果实际foot的维度是(lon, lat, time)，则将求和轴改为axis=2即可：

foot_total = np.sum(foot, axis=2)

方式2：使用xarray（更简洁高效） #

xarray能自动识别维度名称，无需手动判断轴索引，代码可读性更强：

import xarray as xr
import matplotlib.pyplot as plt
import cartopy.crs as ccrs

ds = xr.open_dataset("myfile.nc")
# 对time维度求和
foot_total = ds['foot'].sum(dim='time')

# 绘图
ax = plt.axes(projection=ccrs.PlateCarree())
foot_total.plot.contourf(ax=ax, cmap='YlOrBr', transform=ccrs.PlateCarree())
ax.set_extent([-178.546, -130.946, 50, 75])
ax.coastlines()
ax.gridlines(draw_labels=True)
plt.show()

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