Extract HRRR data using Unidata's Siphon package and Xarray

Unidata Python Workshop



In [ ]:

    
import matplotlib.pyplot as plt
import numpy as np
%matplotlib inline



In [ ]:

    
# Resolve the latest HRRR dataset
from siphon.catalog import get_latest_access_url

hrrr_catalog = "http://thredds.ucar.edu/thredds/catalog/grib/NCEP/HRRR/CONUS_2p5km/catalog.xml"
latest_hrrr_ncss = get_latest_access_url(hrrr_catalog, "NetcdfSubset")

# Set up access via NCSS
from siphon.ncss import NCSS
ncss = NCSS(latest_hrrr_ncss)

# Create a query to ask for all times in netcdf4 format for
# the Temperature_surface variable, with a bounding box
query = ncss.query()



In [ ]:

    
query.all_times().accept('netcdf4').variables('Temperature_height_above_ground')
query.lonlat_box(north=45, south=41, east=-67, west=-77)

# Get the raw bytes and write to a file.
data = ncss.get_data_raw(query)
with open('test.nc', 'wb') as outf:
    outf.write(data)

Try reading extracted data with Xarray



In [ ]:

    
import xarray



In [ ]:

    
nc = xarray.open_dataset('test.nc')



In [ ]:

    
nc



In [ ]:

    
var='Temperature_height_above_ground'
ncvar = nc[var]
ncvar



In [ ]:

    
grid = nc[ncvar.grid_mapping]
grid



In [ ]:

    
lon0 = grid.longitude_of_central_meridian
lat0 = grid.latitude_of_projection_origin
lat1 = grid.standard_parallel
earth_radius = grid.earth_radius

Try plotting the LambertConformal data with Cartopy



In [ ]:

    
import cartopy
import cartopy.crs as ccrs



In [ ]:

    
#cartopy wants meters, not km
x = ncvar.x.data*1000.
y = ncvar.y.data*1000.



In [ ]:

    
#globe = ccrs.Globe(ellipse='WGS84') #default
globe = ccrs.Globe(ellipse='sphere', semimajor_axis=grid.earth_radius)

crs = ccrs.LambertConformal(central_longitude=lon0, central_latitude=lat0, 
                            standard_parallels=(lat0,lat1), globe=globe)



In [ ]:

    
print(ncvar.x.data.shape)
print(ncvar.y.data.shape)
print(ncvar.data.shape)



In [ ]:

    
# find the correct time dimension name
for d in ncvar.dims:
    if "time" in d: 
        timevar = d
nc[timevar].data[6]



In [ ]:

    
istep = 6
fig = plt.figure(figsize=(12,8))
ax = plt.axes(projection=ccrs.PlateCarree())
mesh = ax.pcolormesh(x,y,ncvar[istep,::].data.squeeze(), transform=crs,zorder=0)
ax.coastlines(resolution='10m',color='black',zorder=1)
gl = ax.gridlines(draw_labels=True)
gl.xlabels_top = False
gl.ylabels_right = False
plt.title(nc[timevar].data[istep]);

# TODO: This is to workaround a bug in cartopy 0.14
from matplotlib.transforms import Bbox
if hasattr(mesh, '_corners'): mesh._corners = Bbox(mesh._corners)