notebook.community

Edit and run



In [97]:

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



In [98]:

    
import xray
import datetime



In [99]:

    
# lon/lat bounding box [lon_min, lon_max, lat_min, lat_max]
ax = [-72,-69.6,40.6, 43.5]

# Enter time in EDT
date = datetime.datetime.now()
date = datetime.datetime(2015,7,25,1,0,0)
# convert from EDT to UTC
date += datetime.timedelta(hours=+5)
date









    Out[99]:





datetime.datetime(2015, 7, 25, 6, 0)



In [100]:

    
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 [101]:

    
url = get_latest_url(date)
url









    Out[101]:





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

Try reading extracted data with Xray



In [102]:

    
url









    Out[102]:





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



In [103]:

    
nc = xray.open_dataset(url)



In [104]:

    
#nc



In [105]:

    
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 [106]:

    
uvar









    Out[106]:





<xray.DataArray 'u-component_of_wind_height_above_ground' (time3: 80, height_above_ground5: 2, y: 1059, x: 1799)>
[304822560 values with dtype=float32]
Coordinates:
  * time3                 (time3) datetime64[ns] 2015-07-24T17: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 [107]:

    
grid = nc[uvar.grid_mapping]
grid









    Out[107]:





<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 [108]:

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

Find points of LambertConformal grid withing lon/lat bounding box using Cartopy



In [109]:

    
import cartopy
import cartopy.crs as ccrs



In [110]:

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

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

    
uvar.coords









    Out[112]:





Coordinates:
  * time3                 (time3) datetime64[ns] 2015-07-24T17: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 ...



In [114]:

    
nc.time3









    Out[114]:





<xray.DataArray 'time3' (time3: 80)>
array(['2015-07-24T13:00:00.000000000-0400',
       '2015-07-24T13:01:00.000000000-0400',
       '2015-07-24T13:02:00.000000000-0400',
       '2015-07-24T13:03:00.000000000-0400',
       '2015-07-24T13:04:00.000000000-0400',
       '2015-07-24T13:05:00.000000000-0400',
       '2015-07-24T13:06:00.000000000-0400',
       '2015-07-24T13:07:00.000000000-0400',
       '2015-07-24T13:08:00.000000000-0400',
       '2015-07-24T13:09:00.000000000-0400',
       '2015-07-24T13:10:00.000000000-0400',
       '2015-07-24T13:11:00.000000000-0400',
       '2015-07-24T13:12:00.000000000-0400',
       '2015-07-24T13:13:00.000000000-0400',
       '2015-07-24T13:14:00.000000000-0400',
       '2015-07-24T13:15:00.000000000-0400',
       '2015-07-24T13:16:00.000000000-0400',
       '2015-07-24T13:17:00.000000000-0400',
       '2015-07-24T13:18:00.000000000-0400',
       '2015-07-24T13:19:00.000000000-0400',
       '2015-07-24T13:20:00.000000000-0400',
       '2015-07-24T13:21:00.000000000-0400',
       '2015-07-24T13:22:00.000000000-0400',
       '2015-07-24T13:23:00.000000000-0400',
       '2015-07-24T13:24:00.000000000-0400',
       '2015-07-24T13:30:00.000000000-0400',
       '2015-07-24T13:45:00.000000000-0400',
       '2015-07-24T14:00:00.000000000-0400',
       '2015-07-24T14:15:00.000000000-0400',
       '2015-07-24T14:30:00.000000000-0400',
       '2015-07-24T14:45:00.000000000-0400',
       '2015-07-24T15:00:00.000000000-0400',
       '2015-07-24T15:15:00.000000000-0400',
       '2015-07-24T15:30:00.000000000-0400',
       '2015-07-24T15:45:00.000000000-0400',
       '2015-07-24T16:00:00.000000000-0400',
       '2015-07-24T16:15:00.000000000-0400',
       '2015-07-24T16:30:00.000000000-0400',
       '2015-07-24T16:45:00.000000000-0400',
       '2015-07-24T17:00:00.000000000-0400',
       '2015-07-24T17:15:00.000000000-0400',
       '2015-07-24T17:30:00.000000000-0400',
       '2015-07-24T17:45:00.000000000-0400',
       '2015-07-24T18:00:00.000000000-0400',
       '2015-07-24T18:15:00.000000000-0400',
       '2015-07-24T18:30:00.000000000-0400',
       '2015-07-24T18:45:00.000000000-0400',
       '2015-07-24T19:00:00.000000000-0400',
       '2015-07-24T19:15:00.000000000-0400',
       '2015-07-24T19:30:00.000000000-0400',
       '2015-07-24T19:45:00.000000000-0400',
       '2015-07-24T20:00:00.000000000-0400',
       '2015-07-24T20:15:00.000000000-0400',
       '2015-07-24T20:30:00.000000000-0400',
       '2015-07-24T20:45:00.000000000-0400',
       '2015-07-24T21:00:00.000000000-0400',
       '2015-07-24T21:15:00.000000000-0400',
       '2015-07-24T21:30:00.000000000-0400',
       '2015-07-24T21:45:00.000000000-0400',
       '2015-07-24T22:00:00.000000000-0400',
       '2015-07-24T22:15:00.000000000-0400',
       '2015-07-24T22:30:00.000000000-0400',
       '2015-07-24T22:45:00.000000000-0400',
       '2015-07-24T23:00:00.000000000-0400',
       '2015-07-24T23:15:00.000000000-0400',
       '2015-07-24T23:30:00.000000000-0400',
       '2015-07-24T23:45:00.000000000-0400',
       '2015-07-25T00:00:00.000000000-0400',
       '2015-07-25T01:15:00.000000000-0400',
       '2015-07-25T01:30:00.000000000-0400',
       '2015-07-25T01:45:00.000000000-0400',
       '2015-07-25T02:00:00.000000000-0400',
       '2015-07-25T02:15:00.000000000-0400',
       '2015-07-25T02:30:00.000000000-0400',
       '2015-07-25T02:45:00.000000000-0400',
       '2015-07-25T03:00:00.000000000-0400',
       '2015-07-25T03:15:00.000000000-0400',
       '2015-07-25T03:30:00.000000000-0400',
       '2015-07-25T03:45:00.000000000-0400',
       '2015-07-25T04:00:00.000000000-0400'], dtype='datetime64[ns]')
Coordinates:
  * time3    (time3) datetime64[ns] 2015-07-24T17:00:00 2015-07-24T17:01:00 ...
    reftime  datetime64[ns] ...
Attributes:
    standard_name: time
    long_name: GRIB forecast or observation time
    _CoordinateAxisType: Time



In [113]:

    
zmet=10. # want 10 m height
# use xray to select date and height with real numbers, not indices!
uv = uvar.sel(time3=date, height_above_ground5=zmet ,method='nearest')
vv = vvar.sel(time3=date, height_above_ground5=zmet ,method='nearest')
tim = uv.time3.data









    



---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
<ipython-input-113-30596b5eaa9d> in <module>()
      1 zmet=10. # want 10 m height
      2 # use xray to select date and height with real numbers, not indices!
----> 3 uv = uvar.sel(time=date, height_above_ground5=zmet ,method='nearest')
      4 vv = vvar.sel(time=date, height_above_ground5=zmet ,method='nearest')
      5 tim = uv.time.data

/home/usgs/miniconda/envs/ioos/lib/python2.7/site-packages/xray/core/dataarray.pyc in sel(self, method, **indexers)
    548         """
    549         return self.isel(**indexing.remap_label_indexers(self, indexers,
--> 550                                                          method=method))
    551 
    552     def reindex_like(self, other, method=None, copy=True):

/home/usgs/miniconda/envs/ioos/lib/python2.7/site-packages/xray/core/indexing.pyc in remap_label_indexers(data_obj, indexers, method)
    161     return dict((dim, convert_label_indexer(data_obj[dim].to_index(), label,
    162                                             dim, method))
--> 163                 for dim, label in iteritems(indexers))
    164 
    165 

/home/usgs/miniconda/envs/ioos/lib/python2.7/site-packages/xray/core/indexing.pyc in <genexpr>((dim, label))
    161     return dict((dim, convert_label_indexer(data_obj[dim].to_index(), label,
    162                                             dim, method))
--> 163                 for dim, label in iteritems(indexers))
    164 
    165 

/home/usgs/miniconda/envs/ioos/lib/python2.7/site-packages/xray/core/dataarray.pyc in __getitem__(self, key)
    365     def __getitem__(self, key):
    366         if isinstance(key, basestring):
--> 367             return self.coords[key]
    368         else:
    369             # orthogonal array indexing

/home/usgs/miniconda/envs/ioos/lib/python2.7/site-packages/xray/core/coordinates.pyc in __getitem__(self, key)
     48             return self._dataset[key]
     49         else:
---> 50             raise KeyError(key)
     51 
     52     def __setitem__(self, key, value):

KeyError: 'time'



In [ ]:

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

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



In [ ]:

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



In [ ]:

    
uv

Plot Lambert Conformal Conic data in PlateCarree coordinates using Cartopy



In [ ]:

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



In [ ]:

    
print(spd.shape)
print(x.shape)
print(y.shape)



In [ ]:

    
spd2 = spd[1:,1:]
print(spd2.shape)



In [ ]:

    
fig = plt.figure(figsize=(14,14))
ax1 = plt.axes(projection=ccrs.PlateCarree())
c = ax1.pcolormesh(x,y,spd2, 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 [ ]:



In [ ]:



In [ ]:



In [ ]: