Greetings! Apologies for the delayed response. In terms of plotting, you'll probably be better off using a combination of xarray and MetPy, which can help you out providing some standardized access and handling parsing projection metadata. Here's how I've redone the script you provided: from datetime import datetime, timedelta import cartopy.crs as ccrs import cartopy.feature as cfeature import matplotlib.pyplot as plt from metpy.units import units from metpy.plots import USCOUNTIES, ctables import numpy as np from siphon.catalog import TDSCatalog from siphon.ncss import NCSS import xarray as xr from xarray.backends import NetCDF4DataStore catalog_url = 'https://thredds.ucar.edu/thredds/catalog/grib/NCEP/NBM/CONUS/catalog.xml' cat = TDSCatalog(catalog_url) dataset = cat.datasets['Best National Blend of Models CONUS Grids Time Series'] ncss = dataset.subset() if 'Total_snowfall_surface_1_Hour_Accumulation' in ncss.metadata.variables: field_name = 'Total_snowfall_surface_1_Hour_Accumulation' else: field_name = 'Total_snowfall_surface_Mixed_intervals_Accumulation' query = ncss.query() now = datetime.utcnow() start_time = now end_time = now + timedelta(hours=264) north = 46 south = 42 east = -71 west = -77 query.lonlat_box(north=north, south=south, east=east, west=west) query.accept('netcdf') query.strides(1) query.add_lonlat() query.time_range(start_time, end_time) query.variables(field_name) nbm_ds = xr.open_dataset(NetCDF4DataStore(ncss.get_data(query))) # Parse the CF metadata for the projection information field = nbm_ds.metpy.parse_cf(field_name) snow_cumulative = field.cumsum(dim=field.metpy.time.name) proj = field.metpy.cartopy_crs fig = plt.figure(figsize=(10,10), dpi=100) ax = fig.add_subplot(1, 1, 1, projection=proj) max_color_value = snow_cumulative.max() norm, cmap = ctables.registry.get_with_range('rainbow_r', 0.0, max_color_value) # This color table was chosen mostly to take white out of the color table ax.set_extent((west - 2, east + 2, south - 2, north + 2)) # Bigger extent so I could see put the actual extent into larger geographic context ax.add_feature(cfeature.STATES.with_scale('50m'), zorder=1, edgecolor='black', linewidth=2) artists =  # This naturally iterates over each of the time slices in the accumulation array for time_data in snow_cumulative: display_time = time_data.metpy.time.dt.strftime("%m/%d/%Y at %I:%M:%S %p").item() text = ax.text(0.5, 1.02, "Forecast datetime: " + display_time, ha='center', fontsize='x-large', transform=ax.transAxes) mesh = ax.pcolormesh(field.metpy.x, field.metpy.y, snow_cumulative, transform=proj, zorder=0) artists.append([text, mesh]) #the following 5 lines are just to give me a pause at the end of the animation # because the anim.save parameter to do that does not work artists.append([text, mesh]) artists.append([text, mesh]) artists.append([text, mesh]) artists.append([text, mesh]) artists.append([text, mesh]) from matplotlib.animation import ArtistAnimation anim = ArtistAnimation(fig, artists, interval=200) anim.save('nbmsnowloop.gif', writer='pillow') Note how above I used ncss.metadata.variables to fine the set of variables that are available in the dataset to frame the query. Regarding the times: yes, some of that is complicated. This response might help give you some insight into how to handle the GRIB mixed time intervals: https://www.unidata.ucar.edu/support/help/MailArchives/thredds/msg02568.html Cheers! Ryan > Ryan has been very kind already to help me understand the nbm snow forecast > datasets available in thredds. I still have a couple questions. The dataset > has been a tough nut for me to crack, but I'm almost there... > > I have attached my code, the specific dataset I was using, and the > resulting gif animation I was using to understand the data in the dataset. > > Two questions: > 1. What am I doing wrong in simply trying to see the plot from the x,y > coordinates in the dataset? I thred grabbing the lats/longs also, but it > gives me 2-d arrays for each and pcolormesh cannot handle that. The simple > line I thought SHOULD HAVE worked is commented out in the code. The rough > approximation in the current code was just so I could see something. > 2. I'd like to use this code in an automated process. Short of exception > handling in the code, how can the code know exactly which snow variable > will be in the dataset. That question involves both the ncep schedule and > when it's processed by unidata. I'm not even 100% sure which ncep nbm snow > variable is represented by which thredds dataset variable. The ncep > standard has Snow01 and Snow06 variables for different hourly runs. I'm > guessing, Snow01 becomes 'Total_snowfall_surface_1_Hour_Accumulation' and > Snow06 becomes 'Total_snowfall_surface_Mixed_intervals_Accumulation' in the > Unidata dataset. Ticket Details =================== Ticket ID: UTR-100670 Department: Support THREDDS Priority: Critical Status: Closed =================== NOTE: All email exchanges with Unidata User Support are recorded in the Unidata inquiry tracking system and then made publicly available through the web. If you do not want to have your interactions made available in this way, you must let us know in each email you send to us.
NOTE: All email exchanges with Unidata User Support are recorded in the Unidata inquiry tracking system and then made publicly available through the web. If you do not want to have your interactions made available in this way, you must let us know in each email you send to us.