[Python] How can I run a code where I can plot and save multiple hours of the latest GFS model...

I want to plot precip type from the start of the GFS run (hour 0) through hour 240 at 6-hour intervals. (in this code I only try to go to hour 108) Also, at the end of the code when saving the the plots, how do I save them each as separate .png images?

Here is the code:

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)


I tried using this: dt = start_time + time_delta but this only plots the last timedelta which is hour 108 and not all the other timedelta hours.

Plot: Plot

Continue reading...