In [1]:

    

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


    

In [2]:

    

# Resolve the latest HRRR dataset
from siphon.catalog import TDSCatalog
latest_hrrr = TDSCatalog('http://thredds-jumbo.unidata.ucar.edu/thredds/catalog/grib/HRRR/CONUS_3km/surface/latest.xml')
hrrr_ds = list(latest_hrrr.datasets.values())[0]

# Set up access via NCSS
from siphon.ncss import NCSS
ncss = NCSS(hrrr_ds.access_urls['NetcdfSubset'])

# 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_surface')
query.lonlat_box(45, 41., -63, -71.5)

# 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)


    

In [4]:

    

import xray


    

In [5]:

    

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


    

In [6]:

    

nc


    

    
Out[6]:





<xray.Dataset>
Dimensions:                      (time1: 23, x: 223, y: 216)
Coordinates:
  * time1                        (time1) datetime64[ns] 2015-07-23T16:00:00 ...
  * y                            (y) float32 584.694 587.694 590.694 593.694 ...
  * x                            (x) float32 2030.43 2033.43 2036.43 2039.43 ...
Data variables:
    Temperature_surface          (time1, y, x) float32 306.89 305.953 305.14 ...
    LambertConformal_Projection  int32 0
Attributes:
    Originating_or_generating_Center: The NOAA Forecast Systems Laboratory, Boulder, CO, United States
    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/NOAA_GSD/HRRR/CONUS_3km/surface/HRRR_CONUS_3km_surface_201507231600.grib2.ncx3#LambertConformal_1059X1799-38p51N-97p48W; Translation Date = 2015-07-23T18:18:05.150Z
    geospatial_lat_min: 39.4740842168
    geospatial_lat_max: 46.8690460758
    geospatial_lon_min: -72.9539117154
    geospatial_lon_max: -62.6364751161

In [7]:

    

var='Temperature_surface'
ncvar = nc[var]
ncvar


    

    
Out[7]:





<xray.DataArray 'Temperature_surface' (time1: 23, y: 216, x: 223)>
[1107864 values with dtype=float32]
Coordinates:
  * y        (y) float32 584.694 587.694 590.694 593.694 596.694 599.694 ...
  * x        (x) float32 2030.43 2033.43 2036.43 2039.43 2042.43 2045.43 ...
  * time1    (time1) datetime64[ns] 2015-07-23T16:00:00 2015-07-23T16:01:00 ...
Attributes:
    long_name: Temperature @ Ground or water surface
    units: K
    description: Temperature
    grid_mapping: LambertConformal_Projection
    Grib_Variable_Id: VAR_0-0-0_L1
    Grib2_Parameter: [0 0 0]
    Grib2_Parameter_Discipline: Meteorological products
    Grib2_Parameter_Category: Temperature
    Grib2_Parameter_Name: Temperature
    Grib2_Level_Type: Ground or water surface
    Grib2_Generating_Process_Type: Forecast

In [8]:

    

grid = nc[ncvar.grid_mapping]
grid


    

    
Out[8]:





<xray.DataArray u'LambertConformal_Projection' ()>
array(0)
Attributes:
    grid_mapping_name: lambert_conformal_conic
    latitude_of_projection_origin: 38.5
    longitude_of_central_meridian: 262.5
    standard_parallel: 38.5
    earth_radius: 6371229.0
    _CoordinateTransformType: Projection
    _CoordinateAxisTypes: GeoX GeoY

In [9]:

    

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


    

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)


    

    




(223,)
(216,)
(23, 216, 223)

In [14]:

    

ncvar[6,:,:].time1.data


    

    
Out[14]:





numpy.datetime64('2015-07-23T12:06:00.000000000-0400')

In [15]:

    

istep =6
fig = plt.figure(figsize=(12,8))
ax = plt.axes(projection=ccrs.PlateCarree())
ax.pcolormesh(x,y,ncvar[istep,:,:].data, 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(ncvar[istep].time1.data);


    

In [15]: