Как я могу запустить код, где я могу построить и сэкономить несколько часов последнего запуска модели GFS? ⇐ Python

[Anonymous]

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

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

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 Но это только запланирует последнюю Timedelta , который составляет час 108, а не все другие Timeedelta Часы>

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