Siphon XARRAY Cartopy HRRR
Extract HRRR data using Unidata's Siphon package and Xarray
Unidata Python Workshop
In [1]:
import matplotlib.pyplot as plt
import numpy as np
%matplotlib inline
In [2]:
# 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 [3]:
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 [4]:
import xarray as xr
In [5]:
nc = xr.open_dataset('test.nc')
In [6]:
nc
Out[6]:
In [7]:
var='Temperature_height_above_ground'
ncvar = nc[var]
ncvar
Out[7]:
In [8]:
grid = nc[ncvar.grid_mapping]
grid
Out[8]:
In [9]:
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 [10]:
import cartopy
import cartopy.crs as ccrs
In [11]:
#cartopy wants meters, not km
x = ncvar.x.data*1000.
y = ncvar.y.data*1000.
In [12]:
#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 [13]:
print(ncvar.x.data.shape)
print(ncvar.y.data.shape)
print(ncvar.data.shape)
In [14]:
# find the correct time dimension name
for d in ncvar.dims:
if "time" in d:
timevar = d
nc[timevar].data[6]
Out[14]:
In [15]:
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]);
