Plot wind vectors from the ESRL High Resolution Rapid Refresh (HRRR) met model

In [68]:

    

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


    

In [69]:

    

import xray
import datetime
import pytz


    

In [70]:

    

# lon/lat bounding box [lon_min, lon_max, lat_min, lat_max]
ax = [-72,-69.6,40.8, 43.0]
#ax = [-71.1, -70.2, 42.4, 42.8]

# Enter Eastern Time
eastern = pytz.timezone('US/Eastern')
# now....
loc_dt = eastern.localize(datetime.datetime.now())
# or specified time
#loc_dt = eastern.localize(datetime.datetime(2015,7,25,6,0,0))
fmt = '%Y-%m-%d %H:%M:%S %Z%z'
print(loc_dt.strftime(fmt))


    

    




2015-07-26 12:58:46 EDT-0400

In [71]:

    

def get_latest_url(date):  
    # build URL for latest synoptic analysis time

    # keep moving back 1 hours until a valid URL found
    validURL = False; ncount = 0
    while (not validURL and ncount < 24):
        URL = 'http://thredds-jumbo.unidata.ucar.edu/thredds/dodsC/grib/HRRR/CONUS_3km/surface/HRRR_CONUS_3km_surface_%04i%02i%02i%02i%02i.grib2' %\
(date.year,date.month,date.day,1*(date.hour//1),0)
#        URL='http://thredds-jumbo.unidata.ucar.edu/thredds/dodsC/grib/NCEP/HRRR/CONUS_2p5km/HRRR_CONUS_2p5km_%04i%02i%02i_%02i%02i.grib2' %\
#(date.year,date.month,date.day,1*(date.hour//1),0)
#        print(URL)
        try:
            gfs = xray.open_dataset(URL)
            validURL = True
        except RuntimeError:
            date -= datetime.timedelta(hours=1)
            ncount += 1  
    return URL


    

In [72]:

    

url = get_latest_url(loc_dt)
url


    

    
Out[72]:





'http://thredds-jumbo.unidata.ucar.edu/thredds/dodsC/grib/HRRR/CONUS_3km/surface/HRRR_CONUS_3km_surface_201507261200.grib2'

In [73]:

    

nc = xray.open_dataset(url)


    

In [97]:

    

nc


    

    
Out[97]:





<xray.Dataset>
Dimensions:                                                                           (depth_below_surface: 1, depth_below_surface_layer: 9, depth_below_surface_layer_bounds_1: 2, height_above_ground: 1, height_above_ground1: 1, height_above_ground2: 2, height_above_ground3: 1, height_above_ground4: 2, height_above_ground5: 2, height_above_ground6: 6, height_above_ground_layer: 2, height_above_ground_layer1: 2, height_above_ground_layer1_bounds_1: 2, height_above_ground_layer2: 2, height_above_ground_layer2_bounds_1: 2, height_above_ground_layer3: 1, height_above_ground_layer3_bounds_1: 2, height_above_ground_layer_bounds_1: 2, isobaric_layer: 1, isobaric_layer_bounds_1: 2, pressure_difference_layer: 1, pressure_difference_layer1: 1, pressure_difference_layer1_bounds_1: 2, pressure_difference_layer2: 3, pressure_difference_layer2_bounds_1: 2, pressure_difference_layer3: 1, pressure_difference_layer3_bounds_1: 2, pressure_difference_layer_bounds_1: 2, sigma_layer: 1, sigma_layer_bounds_1: 2, time: 10, time1: 9, time1_bounds_1: 2, time2: 46, time3: 17, time3_bounds_1: 2, x: 1799, y: 1059)
Coordinates:
  * x                                                                                 (x) float32 -2697.57 ...
  * y                                                                                 (y) float32 -1587.31 ...
    reftime                                                                           datetime64[ns] ...
  * time                                                                              (time) datetime64[ns] 2015-07-26T12:00:00 ...
  * time1                                                                             (time1) datetime64[ns] 2015-07-26T12:01:00 ...
  * time2                                                                             (time2) datetime64[ns] 2015-07-26T12:00:00 ...
  * time3                                                                             (time3) datetime64[ns] 2015-07-26T12:01:00 ...
  * height_above_ground                                                               (height_above_ground) float32 1000.0 ...
  * height_above_ground1                                                              (height_above_ground1) float32 2.0 ...
  * depth_below_surface                                                               (depth_below_surface) float32 0.0 ...
  * height_above_ground_layer                                                         (height_above_ground_layer) float32 3500.0 ...
  * isobaric_layer                                                                    (isobaric_layer) float32 75000.0 ...
  * height_above_ground_layer1                                                        (height_above_ground_layer1) float32 500.0 ...
  * pressure_difference_layer                                                         (pressure_difference_layer) float32 9000.0 ...
  * height_above_ground_layer2                                                        (height_above_ground_layer2) float32 500.0 ...
  * height_above_ground2                                                              (height_above_ground2) float32 1000.0 ...
  * pressure_difference_layer1                                                        (pressure_difference_layer1) float32 12750.0 ...
  * sigma_layer                                                                       (sigma_layer) float32 0.65 ...
  * depth_below_surface_layer                                                         (depth_below_surface_layer) float32 0.0 ...
  * height_above_ground3                                                              (height_above_ground3) float32 10.0 ...
  * pressure_difference_layer2                                                        (pressure_difference_layer2) float32 4500.0 ...
  * height_above_ground4                                                              (height_above_ground4) float32 1.0 ...
  * pressure_difference_layer3                                                        (pressure_difference_layer3) float32 70000.0 ...
  * height_above_ground5                                                              (height_above_ground5) float32 10.0 ...
  * height_above_ground6                                                              (height_above_ground6) float32 1.0 ...
  * height_above_ground_layer3                                                        (height_above_ground_layer3) float32 3000.0 ...
  * depth_below_surface_layer_bounds_1                                                (depth_below_surface_layer_bounds_1) int64 0 ...
  * height_above_ground_layer1_bounds_1                                               (height_above_ground_layer1_bounds_1) int64 0 ...
  * height_above_ground_layer2_bounds_1                                               (height_above_ground_layer2_bounds_1) int64 0 ...
  * height_above_ground_layer3_bounds_1                                               (height_above_ground_layer3_bounds_1) int64 0 ...
  * height_above_ground_layer_bounds_1                                                (height_above_ground_layer_bounds_1) int64 0 ...
  * isobaric_layer_bounds_1                                                           (isobaric_layer_bounds_1) int64 0 ...
  * pressure_difference_layer1_bounds_1                                               (pressure_difference_layer1_bounds_1) int64 0 ...
  * pressure_difference_layer2_bounds_1                                               (pressure_difference_layer2_bounds_1) int64 0 ...
  * pressure_difference_layer3_bounds_1                                               (pressure_difference_layer3_bounds_1) int64 0 ...
  * pressure_difference_layer_bounds_1                                                (pressure_difference_layer_bounds_1) int64 0 ...
  * sigma_layer_bounds_1                                                              (sigma_layer_bounds_1) int64 0 ...
  * time1_bounds_1                                                                    (time1_bounds_1) int64 0 ...
  * time3_bounds_1                                                                    (time3_bounds_1) int64 0 ...
Data variables:
    LambertConformal_Projection                                                       int64 0 ...
    time1_bounds                                                                      (time1, time1_bounds_1) datetime64[ns] ...
    time3_bounds                                                                      (time3, time3_bounds_1) datetime64[ns] ...
    height_above_ground_layer_bounds                                                  (height_above_ground_layer, height_above_ground_layer_bounds_1) float32 ...
    isobaric_layer_bounds                                                             (isobaric_layer, isobaric_layer_bounds_1) float32 ...
    height_above_ground_layer1_bounds                                                 (height_above_ground_layer1, height_above_ground_layer1_bounds_1) float32 ...
    pressure_difference_layer_bounds                                                  (pressure_difference_layer, pressure_difference_layer_bounds_1) float32 ...
    height_above_ground_layer2_bounds                                                 (height_above_ground_layer2, height_above_ground_layer2_bounds_1) float32 ...
    pressure_difference_layer1_bounds                                                 (pressure_difference_layer1, pressure_difference_layer1_bounds_1) float32 ...
    sigma_layer_bounds                                                                (sigma_layer, sigma_layer_bounds_1) float32 ...
    depth_below_surface_layer_bounds                                                  (depth_below_surface_layer, depth_below_surface_layer_bounds_1) float32 ...
    pressure_difference_layer2_bounds                                                 (pressure_difference_layer2, pressure_difference_layer2_bounds_1) float32 ...
    pressure_difference_layer3_bounds                                                 (pressure_difference_layer3, pressure_difference_layer3_bounds_1) float32 ...
    height_above_ground_layer3_bounds                                                 (height_above_ground_layer3, height_above_ground_layer3_bounds_1) float32 ...
    Best_4-layer_lifted_index_pressure_difference_layer                               (time, pressure_difference_layer, y, x) float32 ...
    Categorical_freezing_rain_surface                                                 (time2, y, x) float32 ...
    Categorical_ice_pellets_surface                                                   (time2, y, x) float32 ...
    Categorical_rain_surface                                                          (time2, y, x) float32 ...
    Categorical_snow_surface                                                          (time2, y, x) float32 ...
    Convective_available_potential_energy_surface                                     (time, y, x) float32 ...
    Convective_available_potential_energy_pressure_difference_layer                   (time, pressure_difference_layer2, y, x) float32 ...
    Convective_inhibition_surface                                                     (time, y, x) float32 ...
    Convective_inhibition_pressure_difference_layer                                   (time, pressure_difference_layer2, y, x) float32 ...
    Dewpoint_temperature_height_above_ground                                          (time2, height_above_ground1, y, x) float32 ...
    Downward_short-wave_radiation_flux_surface                                        (time2, y, x) float32 ...
    Echo_top_cloud_tops                                                               (time2, y, x) float32 ...
    Geopotential_height_surface                                                       (time, y, x) float32 ...
    Geopotential_height_cloud_tops                                                    (time2, y, x) float32 ...
    Geopotential_height_adiabatic_condensation_lifted                                 (time, y, x) float32 ...
    Geopotential_height_UnknownLevelType-215                                          (time2, y, x) float32 ...
    Geopotential_height_UnknownLevelType-247                                          (time, y, x) float32 ...
    High_cloud_cover_UnknownLevelType-234                                             (time, y, x) float32 ...
    Ice_cover_surface                                                                 (time, y, x) float32 ...
    Land_cover_0__sea_1__land_surface                                                 (time, y, x) float32 ...
    Low_cloud_cover_UnknownLevelType-214                                              (time, y, x) float32 ...
    Medium_cloud_cover_UnknownLevelType-224                                           (time, y, x) float32 ...
    Moisture_availability_depth_below_surface                                         (time, depth_below_surface, y, x) float32 ...
    Per_cent_frozen_precipitation_surface                                             (time2, y, x) float32 ...
    Planetary_boundary_layer_height_surface                                           (time, y, x) float32 ...
    Potential_temperature_height_above_ground                                         (time, height_above_ground1, y, x) float32 ...
    Precipitable_water_UnknownLevelType-200                                           (time, y, x) float32 ...
    Precipitation_rate_surface                                                        (time2, y, x) float32 ...
    Pressure_surface                                                                  (time2, y, x) float32 ...
    Pressure_cloud_base                                                               (time, y, x) float32 ...
    Pressure_cloud_tops                                                               (time, y, x) float32 ...
    Pressure_reduced_to_MSL_msl                                                       (time, y, x) float32 ...
    Relative_humidity_height_above_ground                                             (time, height_above_ground1, y, x) float32 ...
    Seconds_prior_to_initial_reference_time_defined_in_Section_1_height_above_ground  (time2, height_above_ground6, y, x) float32 ...
    Snow_cover_surface                                                                (time, y, x) float32 ...
    Snow_depth_surface                                                                (time, y, x) float32 ...
    Soil_temperature_depth_below_surface_layer                                        (time, depth_below_surface_layer, y, x) float32 ...
    Specific_humidity_height_above_ground                                             (time2, height_above_ground1, y, x) float32 ...
    Storm_relative_helicity_height_above_ground_layer                                 (time, height_above_ground_layer2, y, x) float32 ...
    Surface_lifted_index_isobaric_layer                                               (time, isobaric_layer, y, x) float32 ...
    Temperature_surface                                                               (time, y, x) float32 ...
    Temperature_height_above_ground                                                   (time2, height_above_ground1, y, x) float32 ...
    Total_cloud_cover_entire_atmosphere                                               (time, y, x) float32 ...
    Total_column_integrated_graupel_UnknownLevelType-200_1_Minute_Maximum             (time1, y, x) float32 ...
    Total_precipitation_surface_Mixed_intervals_Accumulation                          (time3, y, x) float32 ...
    VAR0-1-242_FROM_59-0--1_entire_atmosphere                                         (time, y, x) float32 ...
    VAR0-16-195_FROM_59-0--1_height_above_ground                                      (time2, height_above_ground2, y, x) float32 ...
    VAR0-16-196_FROM_59-0--1_entire_atmosphere                                        (time2, y, x) float32 ...
    VAR0-16-198_FROM_59-0--1_height_above_ground_1_Minute_Maximum                     (time1, height_above_ground, y, x) float32 ...
    VAR0-16-201_FROM_59-0--1_entire_atmosphere                                        (time2, y, x) float32 ...
    VAR0-17-1_FROM_59-0--1_height_above_ground                                        (time, height_above_ground4, y, x) float32 ...
    VAR0-17-192_FROM_59-0--1_surface                                                  (time, y, x) float32 ...
    VAR0-2-220_FROM_59-0--1_pressure_difference_layer_1_Minute_Maximum                (time1, pressure_difference_layer3, y, x) float32 ...
    VAR0-2-221_FROM_59-0--1_pressure_difference_layer_1_Minute_Maximum                (time1, pressure_difference_layer3, y, x) float32 ...
    VAR0-3-200_FROM_59-0--1_pressure_difference_layer                                 (time, pressure_difference_layer1, y, x) float32 ...
    VAR0-7-197_FROM_59-0--1_height_above_ground_layer                                 (time2, height_above_ground_layer, y, x) float32 ...
    VAR0-7-199_FROM_59-0--1_height_above_ground_layer_1_Minute_Maximum                (time1, height_above_ground_layer, y, x) float32 ...
    VAR2-0-192_FROM_59-0--1_depth_below_surface_layer                                 (time, depth_below_surface_layer, y, x) float32 ...
    VAR2-0-198_FROM_59-0--1_surface                                                   (time, y, x) float32 ...
    VAR3-192-1_FROM_59-0--1_atmosphere_top                                            (time2, y, x) float32 ...
    VAR3-192-2_FROM_59-0--1_atmosphere_top                                            (time2, y, x) float32 ...
    VAR3-192-7_FROM_59-0--1_atmosphere_top                                            (time2, y, x) float32 ...
    VAR3-192-8_FROM_59-0--1_atmosphere_top                                            (time2, y, x) float32 ...
    Upward_long-wave_radiation_flux_atmosphere_top                                    (time2, y, x) float32 ...
    Vertical_u-component_shear_height_above_ground_layer                              (time, height_above_ground_layer1, y, x) float32 ...
    Vertical_v-component_shear_height_above_ground_layer                              (time, height_above_ground_layer1, y, x) float32 ...
    Vertical_velocity_geometric_sigma_layer_1_Minute_Average                          (time1, sigma_layer, y, x) float32 ...
    Vertically_integrated_liquid_water_VIL_entire_atmosphere                          (time2, y, x) float32 ...
    Visibility_surface                                                                (time2, y, x) float32 ...
    Water_equivalent_of_accumulated_snow_depth_surface                                (time, y, x) float32 ...
    Water_equivalent_of_accumulated_snow_depth_surface_Mixed_intervals_Accumulation   (time3, y, x) float32 ...
    Wind_speed_height_above_ground_1_Minute_Maximum                                   (time1, height_above_ground3, y, x) float32 ...
    Wind_speed_gust_surface                                                           (time2, y, x) float32 ...
    u-component_of_wind_height_above_ground                                           (time2, height_above_ground5, y, x) float32 ...
    u-component_storm_motion_height_above_ground_layer                                (time, height_above_ground_layer3, y, x) float32 ...
    v-component_of_wind_height_above_ground                                           (time2, height_above_ground5, y, x) float32 ...
    v-component_storm_motion_height_above_ground_layer                                (time, height_above_ground_layer3, y, x) float32 ...
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
    file_format: GRIB-2
    Conventions: CF-1.6
    history: Read using CDM IOSP GribCollection v3
    featureType: GRID
    _CoordSysBuilder: ucar.nc2.dataset.conv.CF1Convention

In [75]:

    

uvar_name='u-component_of_wind_height_above_ground'
vvar_name='v-component_of_wind_height_above_ground'
uvar = nc[uvar_name]
vvar = nc[vvar_name]


    

In [76]:

    

uvar


    

    
Out[76]:





<xray.DataArray 'u-component_of_wind_height_above_ground' (time2: 46, height_above_ground5: 2, y: 1059, x: 1799)>
[175272972 values with dtype=float32]
Coordinates:
  * time2                 (time2) datetime64[ns] 2015-07-26T12:00:00 ...
  * height_above_ground5  (height_above_ground5) float32 10.0 80.0
    reftime               datetime64[ns] ...
  * y                     (y) float32 -1587.31 -1584.31 -1581.31 -1578.31 ...
  * x                     (x) float32 -2697.57 -2694.57 -2691.57 -2688.57 ...
Attributes:
    long_name: u-component of wind @ Specified height level above ground
    units: m/s
    description: u-component of wind
    grid_mapping: LambertConformal_Projection
    Grib_Variable_Id: VAR_0-2-2_L103
    Grib2_Parameter: [0 2 2]
    Grib2_Parameter_Discipline: Meteorological products
    Grib2_Parameter_Category: Momentum
    Grib2_Parameter_Name: u-component of wind
    Grib2_Level_Type: Specified height level above ground
    Grib2_Generating_Process_Type: Forecast

In [77]:

    

grid = nc[uvar.grid_mapping]
grid


    

    
Out[77]:





<xray.DataArray u'LambertConformal_Projection' ()>
array(0)
Coordinates:
    reftime  datetime64[ns] ...
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 [78]:

    

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


    

In [79]:

    

import cartopy
import cartopy.crs as ccrs


    

In [80]:

    

#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 [81]:

    

#Find indices to read based on lon/lat bounding box
#newcs = ccrs.PlateCarree.transform_points(crs,xx,yy)
dest = ccrs.PlateCarree()
#cartopy wants meters, not km
x2d, y2d = np.meshgrid(uvar.x.data*1000.0, uvar.y.data*1000.0)
lonlat = dest.transform_points(crs, x2d, y2d)


    

In [82]:

    

print(uvar.coords.keys())
timecoord = None
for coord in uvar.coords.keys():
    try:
        if uvar.coords[coord].standard_name == 'time':
            timecoord = coord
    except:
        pass

print(timecoord)


    

    




[u'time2', u'height_above_ground5', u'reftime', u'y', u'x']
time2

In [83]:

    

zmet=10. # want 10 m height
# use xray to select date and height with real numbers, not indices!
kwargs = {timecoord: loc_dt}
uv = uvar.sel(height_above_ground5=zmet ,method='nearest', **kwargs)
vv = vvar.sel(height_above_ground5=zmet ,method='nearest', **kwargs)
tim = uv[timecoord].data


    

In [84]:

    

# determine index range to select data that spans our lon/lat bounding box 
i = np.argwhere((lonlat[:,:,0] >= ax[0]) & 
               (lonlat[:,:,0] <= ax[1]) & 
               (lonlat[:,:,1] >= ax[2]) & 
               (lonlat[:,:,1] <= ax[3]))

j0 = i[:,0].min()
j1 = i[:,0].max()
i0 = i[:,1].min()
i1 = i[:,1].max()


    

In [85]:

    

print(uv.shape)
print(i0,i1,j0,j1)


    

    




(1059, 1799)
(1584, 1670, 714, 810)

In [86]:

    

uv=uv[j0:j1,i0:i1]
vv=vv[j0:j1,i0:i1]


    

In [87]:

    

x = uv.x.data*1000.
y = uv.y.data*1000.
u = uv.data
v = vv.data
spd = np.sqrt(u*u+v*v)


    

In [88]:

    

# create xc,yc for corners of cells

dx = x[1] - x[0]
# subtract 1/2 grid cell from center to get 1st edge
xc = x - dx/2.
# add an extra element to get the last edge
xc = np.append(xc, xc[-1]+dx)

dy = y[1] - y[0]
yc = y + dy/2.
yc = np.append(yc, yc[-1]+dy)


    

In [89]:

    

print(spd.shape)
print(x.shape, y.shape)
print(xc.shape, yc.shape)


    

    




(96, 86)
((86,), (96,))
((87,), (97,))

In [90]:

    

fig = plt.figure(figsize=(14,14))
ax1 = plt.axes(projection=ccrs.PlateCarree())
c = ax1.pcolormesh(xc,yc,spd, transform=crs,zorder=0,vmin=0,vmax=10)
cb = fig.colorbar(c,orientation='vertical',shrink=0.5)
cb.set_label('m/s')
ax1.coastlines(resolution='10m',color='black',zorder=1)
ax1.quiver(x,y,u,v,transform=crs,zorder=2,scale=200)
gl = ax1.gridlines(draw_labels=True)
gl.xlabels_top = False
gl.ylabels_right = False
plt.title(tim);
plt.axis(ax);


    

In [98]:

    

tvar = nc['Temperature_height_above_ground']


    

In [99]:

    

print(tvar.coords.keys())
timecoord = None
for coord in tvar.coords.keys():
    try:
        if tvar.coords[coord].standard_name == 'time':
            timecoord = coord
    except:
        pass

print(timecoord)


    

    




[u'time2', u'height_above_ground1', u'reftime', u'y', u'x']
time2

In [100]:

    

kwargs = {timecoord: loc_dt}
tv = tvar.sel(height_above_ground1=2.,method='nearest',**kwargs)
tv = tv[j0:j1,i0:i1]
tim2 = tv[timecoord].data


    

In [101]:

    

fig = plt.figure(figsize=(14,14))
ax1 = plt.axes(projection=ccrs.PlateCarree())
c = ax1.pcolormesh(xc,yc,((tv.data-273.15)*9./5.+32.), transform=crs,zorder=0)
cb = fig.colorbar(c,orientation='vertical',shrink=0.5)
cb.set_label('deg_F')
ax1.coastlines(resolution='10m',color='black',zorder=1)
gl = ax1.gridlines(draw_labels=True)
gl.xlabels_top = False
gl.ylabels_right = False
plt.title(tim2);
plt.axis(ax);


    

In [102]:

    

uvar.time1.data


    

    




---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-102-6c0b150964cf> in <module>()
----> 1 uvar.time1.data

/home/usgs/miniconda/envs/ioos/lib/python2.7/site-packages/xray/core/common.pyc in __getattr__(self, name)
    134                 pass
    135         raise AttributeError("%r object has no attribute %r" %
--> 136                              (type(self).__name__, name))
    137 
    138     def __dir__(self):

AttributeError: 'DataArray' object has no attribute 'time1'

In [ ]: