我的以下代码实现了绘制一个点位的功能,现在需要将ds文件中,包含的所有点位,都读取出来,一次性的绘制在这一个图层上,请修改代码
def plot_concentration_field_all_time02(ds, pointspec_idx=0, height_idx=0):
"""
绘制同一释放点、同一高度层所有时间步的浓度总和(单图层绘制,颜色由总浓度决定)
:param ds: xarray Dataset(FLEXPART输出)
:param pointspec_idx: 释放点索引(0-13)
:param height_idx: 高度层索引(0-3)
"""
# --------------------------
# 1. 提取固定释放点和高度层的所有时间步数据
# --------------------------
conc_data = ds['spec001_mr'].isel(
nageclass=0,
pointspec=pointspec_idx,
height=height_idx
) # 形状:(time, lat, lon)
# 打印调试信息(关键!)
# print("conc_data 原始形状:", conc_data.shape) # 应输出 (216, 200, 200)
# --------------------------
# 2. 时间维度求和:得到二维总浓度场
# --------------------------
sum_conc = conc_data.sum(dim='time') # 每个(lat, lon)的所有时间步浓度累加
sum_conc_values = sum_conc.values # 提取numpy数组
# print("sum_conc_values 形状:", sum_conc_values.shape) # 应输出 (200, 200)
# 定义阈值(低于100的区域不显示)
threshold = 10
# 生成掩码:小于100的区域标记为True(需要屏蔽)
mask = sum_conc_values < threshold
# 生成masked数组(屏蔽小于100的区域)
masked_sum_conc = np.ma.masked_where(mask, sum_conc_values)
# --------------------------
# 3. 定义总浓度的颜色映射与分档
# --------------------------
cmap = plt.get_cmap('Blues') # 可替换为'viridis'/'hot'等色卡
# 分档边界:覆盖总浓度的最小-最大值
levels = np.linspace(sum_conc.min().values, sum_conc.max().values, 10)
# --------------------------
# 4. 初始化画布和地图投影
# --------------------------
fig = plt.figure(figsize=(16, 12))
ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())
set_dynamic_extent(ax, region_code='370800', buffer=0.05) # 自定义区域范围函数
# --------------------------
# 6. 绘制地理要素(zorder低于浓度场,确保在底层)
# --------------------------
zorderidx = 1
logger.info(f'正在绘制地图上乡镇边界')
# 绘制乡镇街道边界
plot_street_shpfile('', ax, '', '', fig, zorderidx)
zorderidx += 1
logger.info(f'正在绘制地图上区县轮廓')
# 绘制全国县城
plot_china_countys('', ax, '', '', fig, zorderidx)
zorderidx += 1
logger.info(f'正在绘制地图上城市边界')
# 绘制城市边界
plot_city_boundaries('', ax, zorderidx, region_code='370800')
zorderidx += 1
logger.info(f'正在绘制地图上网格')
# 画经纬度网格线
plot_gridline('', ax, 'lon', 'lat', fig, zorderidx)
zorderidx += 1
# --------------------------
# 5. 单次绘制总浓度场(替代循环叠加)
# --------------------------
# 生成经纬度网格(二维)
lon = sum_conc.longitude.values
lat = sum_conc.latitude.values
# lon_2d, lat_2d = np.meshgrid(lon, lat, indexing='ij')
ax.contourf(
lon, lat, masked_sum_conc,
levels=levels, # 总浓度的分档边界
cmap=cmap, # 颜色映射
transform=ccrs.PlateCarree(), # 数据投影
alpha=1, # 固定透明度(避免叠加浑浊)
zorder=0 # 确保浓度场在地理要素上方
)
# --------------------------
# 7. 颜色条(仅显示总浓度范围)
# --------------------------
cbar = plt.colorbar(
plt.cm.ScalarMappable(cmap=cmap, norm=plt.Normalize(vmin=levels[0], vmax=levels[-1])),
ax=ax,
shrink=0.8,
label='总浓度'
)
# --------------------------
# 8. 标题与标注
# --------------------------
relcom = bytes(ds.RELCOM[pointspec_idx].values).decode().strip().replace('city', '')
start_time = datetime.utcfromtimestamp(
ds.time.isel(time=0).values.astype('datetime64[s]').astype(int)
)
end_time = datetime.utcfromtimestamp(
ds.time.isel(time=-1).values.astype('datetime64[s]').astype(int)
)
plt.suptitle('FLEXPART 单释放点污染物扩散总浓度图', fontsize=18, fontweight='bold', y=0.96)
plt.title(f'释放点: {relcom.upper()}, 高度: {ds.height.isel(height=height_idx).values:.0f}m\n'
f'时间范围: {start_time.strftime("%Y-%m-%d %H:%M")} 至 {end_time.strftime("%Y-%m-%d %H:%M")}',
fontsize=14, pad=15)
# --------------------------
# 9. 保存图像
# --------------------------
plt.savefig('single_release_point_total_conc.png', dpi=300, bbox_inches='tight')
plt.close(fig)