Как запустить код, позволяющий построить график и сохранить несколько часов работы последней модели GFS? ⇐ Python

[Anonymous]

Я хочу построить график типа осадков от начала прогона GFS (0-й час) до 240-го часа с 6-часовыми интервалами. (в этом коде я пытаюсь перейти только к часу 108). Кроме того, в конце кода при сохранении графиков, как мне сохранить их каждый как отдельные изображения PNG?
Здесь это код:

Код: Выделить всё

from datetime import datetime
from datetime import timedelta
import pandas as pd
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt
import matplotlib.colors as mcolors
import xarray as xr
import numpy as np
import metpy.calc as mpcalc
from metpy.plots import USCOUNTIES
import netCDF4
from netCDF4 import Dataset
from netCDF4 import num2date
from metpy.units import units
from scipy.ndimage import gaussian_filter
import scipy.ndimage as ndimage
from siphon.catalog import TDSCatalog

start_time = datetime(2025, 1, 7, 12, 0, 0)
time_deltas = [timedelta(hours=6), timedelta(hours=12), timedelta(hours=18), timedelta(hours=24), timedelta(hours=30), timedelta(hours=36),
timedelta(hours=42), timedelta(hours=48), timedelta(hours=54), timedelta(hours=60), timedelta(hours=66), timedelta(hours=72),
timedelta(hours=78), timedelta(hours=84), timedelta(hours=90), timedelta(hours=96), timedelta(hours=102), timedelta(hours=108)]
for time_delta in time_deltas:
dt = start_time + time_delta
#dt = datetime(2025,1,4,12)
best_gfs = TDSCatalog('https://thredds.ucar.edu/thredds/catalog/grib/NCEP/GFS/Global_0p25deg/catalog.xml?dataset=grib/NCEP/GFS/Global_0p25deg/Best')
best_ds = best_gfs.datasets[0]
ncss = best_ds.subset()
query = ncss.query()
query.accept('netcdf')
query.lonlat_box(north=75, south=15, east=320, west=185)
query.time(dt)
query.variables('Geopotential_height_isobaric', 'Pressure_reduced_to_MSL_msl', 'Precipitation_rate_surface', 'Snow_depth_surface', 'Categorical_Snow_surface','Categorical_Freezing_Rain_surface', 'Categorical_Ice_Pellets_surface')

data = ncss.get_data(query)
print(list(data.variables))

plev = list(data.variables['isobaric'][:])

lat = data.variables['latitude'][:].squeeze()
lon = data.variables['longitude'][:].squeeze()
time1 = data['time']
vtime = num2date(time1[:].squeeze(), units=time1.units)
emsl_var = data.variables['Pressure_reduced_to_MSL_msl']
preciprate = data.variables['Precipitation_rate_surface'][:].squeeze()
snowdepth = data.variables['Snow_depth_surface'][:].squeeze()
catsnow = data.variables['Categorical_Snow_surface'][:].squeeze()
catice = data.variables['Categorical_Freezing_Rain_surface'][:].squeeze()
catsleet = data.variables['Categorical_Ice_Pellets_surface'][:].squeeze()
EMSL = units.Quantity(emsl_var[:], emsl_var.units).to('hPa')
mslp = gaussian_filter(EMSL[0], sigma=3.0)
hght_1000 = data.variables['Geopotential_height_isobaric'][0, plev.index(100000)]
hght_500 = data.variables['Geopotential_height_isobaric'][0, plev.index(50000)]
thickness_1000_500 = gaussian_filter((hght_500 - hght_1000)/10, sigma=3.0)
lon_2d, lat_2d = np.meshgrid(lon, lat)

precip_inch_hour = preciprate * 141.73228346457
precip2 = mpcalc.smooth_n_point(precip_inch_hour, 5, 1)

precip_colors = [
"#bde9bf",  # 0.01 - 0.02 inches 1
"#adddb0",  # 0.02 - 0.03 inches 2
"#9ed0a0",  # 0.03 - 0.04 inches 3
"#8ec491",  # 0.04 - 0.05 inches 4
"#7fb882",  # 0.05 - 0.06 inches 5
"#70ac74",  # 0.06 - 0.07 inches 6
"#60a065",  # 0.07 - 0.08 inches 7
"#519457",  # 0.08 - 0.09 inches 8
"#418849",  # 0.09 - 0.10 inches 9
"#307c3c",  # 0.10 - 0.12 inches 10
"#1c712e",  # 0.12 - 0.14 inches 11
"#f7f370",  # 0.14 - 0.16 inches 12
"#fbdf65",  # 0.16 - 0.18 inches 13
"#fecb5a",  # 0.18 - 0.2 inches 14
"#ffb650",  # 0.2 - 0.3 inches 15
"#ffa146",  # 0.3 - 0.4 inches 16
"#ff8b3c",   # 0.4 - 0.5 inches 17
"#f94609",   # 0.5 - 0.6 inches 18
]

precip_colormap = mcolors.ListedColormap(precip_colors)

clev_precip =  np.concatenate((np.arange(0.01, 0.1, .01), np.arange(.1, .2, .02), np.arange(.2, .61, .1)))
norm = mcolors.BoundaryNorm(clev_precip, 18)

datacrs = ccrs.PlateCarree()
plotcrs = ccrs.LambertConformal(central_latitude=35, central_longitude=-100,standard_parallels=(30, 60))
bounds = ([-105, -90, 30, 40])
fig = plt.figure(figsize=(14,12))
ax = fig.add_subplot(1,1,1, projection=plotcrs)
ax.set_extent(bounds, crs=ccrs.PlateCarree())
ax.add_feature(cfeature.COASTLINE.with_scale('50m'), linewidth = 0.75)
ax.add_feature(cfeature.STATES, linewidth = 1)
ax.add_feature(USCOUNTIES, edgecolor='grey',  linewidth = .5)
clevs = (np.arange(0, 540, 6),
np.array([540]),
np.arange(546, 700, 6))
colors = ('tab:blue', 'b', 'tab:red')
kw_clabels = {'fontsize': 11, 'inline': True, 'inline_spacing': 5, 'fmt': '%i',
'rightside_up': True, 'use_clabeltext': True}

# Plot MSLP
clevmslp = np.arange(800., 1120., 2)
cs2 = ax.contour(lon_2d, lat_2d, mslp, clevmslp, colors='k', linewidths=1.25,
linestyles='solid', transform=ccrs.PlateCarree())

cf = ax.contourf(lon_2d, lat_2d, precip2, clev_precip, cmap=precip_colormap, norm=norm, extend='max', transform=ccrs.PlateCarree())

ax.set_title('GFS Precip Type, Rate(in/hr), MSLP (hPa), & 1000-500mb Thickness (dam)', loc='left', fontsize=10, weight = 'bold')
ax.set_title('Valid Time: {}z'.format(vtime), loc = 'right', fontsize=8)

Я пробовал использовать это: dt = start_time + time_delta, но это отображает только последнюю дельту времени, которая равна 108 часам, а не все остальные часы дельты времени.
Сюжет: Сюжет


Подробнее здесь: https://stackoverflow.com/questions/793 ... latest-gfs