import matplotlib.pyplot as plt
import numpy as np
%matplotlib inline
# 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()
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)
import xarray
nc = xarray.open_dataset('test.nc')
nc
<xarray.Dataset> Dimensions: (height_above_ground1: 1, time2: 16, x: 365, y: 240) Coordinates: * time2 (time2) datetime64[ns] 2016-02-04T15:00:00 ... * height_above_ground1 (height_above_ground1) float32 2.0 * y (y) float32 1907.66 1910.2 1912.73 ... * x (x) float32 1508.56 1511.1 1513.64 ... Data variables: Temperature_height_above_ground (time2, height_above_ground1, y, x) float32 277.922 ... LambertConformal_Projection int32 0 Attributes: Originating_or_generating_Center: US National Weather Service, National Centres for Environmental Prediction (NCEP) Originating_or_generating_Subcenter: 0 GRIB_table_version: 2,1 Type_of_generating_process: Forecast Conventions: CF-1.6 history: Read using CDM IOSP GribCollection v3 featureType: GRID History: Translated to CF-1.0 Conventions by Netcdf-Java CDM (CFGridWriter2) Original Dataset = /data/ldm/pub/native/grid/NCEP/HRRR/CONUS_2p5km/HRRR_CONUS_2p5km_20160204_1500.grib2.ncx3#LambertConformal_1377X2145-38p22N-95p43W; Translation Date = 2016-02-04T16:46:28.014Z geospatial_lat_min: 39.7308379134 geospatial_lat_max: 46.1872175578 geospatial_lon_min: -77.7129121541 geospatial_lon_max: -65.8516561614
var='Temperature_height_above_ground'
ncvar = nc[var]
ncvar
<xarray.DataArray 'Temperature_height_above_ground' (time2: 16, height_above_ground1: 1, y: 240, x: 365)> [1401600 values with dtype=float32] Coordinates: * time2 (time2) datetime64[ns] 2016-02-04T15:00:00 ... * height_above_ground1 (height_above_ground1) float32 2.0 * y (y) float32 1907.66 1910.2 1912.73 1915.27 1917.81 ... * x (x) float32 1508.56 1511.1 1513.64 1516.18 1518.72 ... Attributes: long_name: Temperature @ Specified height level above ground units: K abbreviation: TMP grid_mapping: LambertConformal_Projection Grib_Variable_Id: VAR_0-0-0_L103 Grib2_Parameter: [0 0 0] Grib2_Parameter_Discipline: Meteorological products Grib2_Parameter_Category: Temperature Grib2_Parameter_Name: Temperature Grib2_Level_Type: Specified height level above ground Grib2_Generating_Process_Type: Forecast
grid = nc[ncvar.grid_mapping]
grid
<xarray.DataArray 'LambertConformal_Projection' ()> array(0) Attributes: grid_mapping_name: lambert_conformal_conic latitude_of_projection_origin: 25.0 longitude_of_central_meridian: 265.0 standard_parallel: 25.0 earth_radius: 6371229.0 _CoordinateTransformType: Projection _CoordinateAxisTypes: GeoX GeoY
lon0 = grid.longitude_of_central_meridian
lat0 = grid.latitude_of_projection_origin
lat1 = grid.standard_parallel
earth_radius = grid.earth_radius
import cartopy
import cartopy.crs as ccrs
#cartopy wants meters, not km
x = ncvar.x.data*1000.
y = ncvar.y.data*1000.
#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)
print(ncvar.x.data.shape)
print(ncvar.y.data.shape)
print(ncvar.data.shape)
(365,) (240,) (16, 1, 240, 365)
# find the correct time dimension name
for d in ncvar.dims:
if "time" in d:
timevar = d
nc[timevar].data[6]
numpy.datetime64('2016-02-04T14:00:00.000000000-0700')
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)
/Users/lesserwhirls/miniconda3/envs/devel/lib/python3.4/site-packages/matplotlib/artist.py:221: MatplotlibDeprecationWarning: This has been deprecated in mpl 1.5, please use the axes property. A removal date has not been set. warnings.warn(_get_axes_msg, mplDeprecation, stacklevel=1)