[THREDDS #UTR-100670]: nbm thredds data

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[0], 
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.