I need to open a grib2 file with xarray. To do so, I'm using python 2.7 and pynio as engine in xarray:
grbs = xr.open_dataset('hrrr.t06z.wrfsubhf02.grib2'], engine = 'pynio')
Output:
<xarray.Dataset>
Dimensions: (forecast_time0: 4, lv_HTGL0: 2, lv_HTGL1: 2, xgrid_0: 1799, ygrid_0: 1059)
Coordinates:
* forecast_time0 (forecast_time0) timedelta64[ns] 5 days 15:00:00 ...
* lv_HTGL1 (lv_HTGL1) float32 1000.0 4000.0
* lv_HTGL0 (lv_HTGL0) float32 10.0 80.0
gridlat_0 (ygrid_0, xgrid_0) float32 21.1381 21.1451 ...
gridlon_0 (ygrid_0, xgrid_0) float32 -122.72 -122.693 ...
Dimensions without coordinates: xgrid_0, ygrid_0
Data variables:
ULWRF_P0_L1_GLC0 (forecast_time0, ygrid_0, xgrid_0) float64 446.3 ...
WIND_P8_L103_GLC0_avg5min (forecast_time0, ygrid_0, xgrid_0) float64 5.31 ...
SBT124_P0_L8_GLC0 (forecast_time0, ygrid_0, xgrid_0) float64 288.7 ...
VDDSF_P0_L1_GLC0 (forecast_time0, ygrid_0, xgrid_0) float64 0.0 ...
VIS_P0_L1_GLC0 (forecast_time0, ygrid_0, xgrid_0) float64 8.4e+03 ...
DSWRF_P0_L1_GLC0 (forecast_time0, ygrid_0, xgrid_0) float64 0.0 ...
CICEP_P0_L1_GLC0 (forecast_time0, ygrid_0, xgrid_0) float64 0.0 ...
DLWRF_P0_L1_GLC0 (forecast_time0, ygrid_0, xgrid_0) float64 429.7 ...
USWRF_P0_L1_GLC0 (forecast_time0, ygrid_0, xgrid_0) float64 0.0 ...
ULWRF_P0_L8_GLC0 (forecast_time0, ygrid_0, xgrid_0) float64 291.0 ...
HGT_P0_L3_GLC0 (forecast_time0, ygrid_0, xgrid_0) float64 1.4e+03 ...
VGRD_P8_L103_GLC0_avg5min (forecast_time0, ygrid_0, xgrid_0) float64 -4.92 ...
gridrot_0 (ygrid_0, xgrid_0) float32 -0.274008 ...
VIL_P0_L10_GLC0 (forecast_time0, ygrid_0, xgrid_0) float64 0.0048 ...
CSNOW_P0_L1_GLC0 (forecast_time0, ygrid_0, xgrid_0) float64 0.0 ...
SBT123_P0_L8_GLC0 (forecast_time0, ygrid_0, xgrid_0) float64 251.5 ...
GUST_P0_L1_GLC0 (forecast_time0, ygrid_0, xgrid_0) float64 6.469 ...
SBT114_P0_L8_GLC0 (forecast_time0, ygrid_0, xgrid_0) float64 289.4 ...
DPT_P0_L103_GLC0 (forecast_time0, ygrid_0, xgrid_0) float64 295.8 ...
UGRD_P8_L103_GLC0_avg5min (forecast_time0, ygrid_0, xgrid_0) float64 -2.02 ...
RETOP_P0_L3_GLC0 (forecast_time0, ygrid_0, xgrid_0) float64 -999.0 ...
REFD_P0_L103_GLC0 (forecast_time0, lv_HTGL1, ygrid_0, xgrid_0) float64 -10.0 ...
TMP_P0_L103_GLC0 (forecast_time0, ygrid_0, xgrid_0) float64 297.0 ...
UGRD_P0_L103_GLC0 (forecast_time0, lv_HTGL0, ygrid_0, xgrid_0) float64 -1.998 ...
HGT_P0_L215_GLC0 (forecast_time0, ygrid_0, xgrid_0) float64 461.2 ...
UPHL_P0_2L103_GLC0 (forecast_time0, ygrid_0, xgrid_0) float64 0.0 ...
SBT113_P0_L8_GLC0 (forecast_time0, ygrid_0, xgrid_0) float64 253.9 ...
VBDSF_P0_L1_GLC0 (forecast_time0, ygrid_0, xgrid_0) float64 0.0 ...
PRATE_P0_L1_GLC0 (forecast_time0, ygrid_0, xgrid_0) float64 0.0 ...
CFRZR_P0_L1_GLC0 (forecast_time0, ygrid_0, xgrid_0) float64 0.0 ...
CPOFP_P0_L1_GLC0 (forecast_time0, ygrid_0, xgrid_0) float64 -50.0 ...
CRAIN_P0_L1_GLC0 (forecast_time0, ygrid_0, xgrid_0) float64 0.0 ...
REFC_P0_L10_GLC0 (forecast_time0, ygrid_0, xgrid_0) float64 1.0 ...
DSWRF_P8_L1_GLC0_avg15min (forecast_time0, ygrid_0, xgrid_0) float64 0.0 ...
VBDSF_P8_L1_GLC0_avg15min (forecast_time0, ygrid_0, xgrid_0) float64 0.0 ...
PRES_P0_L1_GLC0 (forecast_time0, ygrid_0, xgrid_0) float64 1.014e+05 ...
SPFH_P0_L103_GLC0 (forecast_time0, ygrid_0, xgrid_0) float64 0.01705 ...
VGRD_P0_L103_GLC0 (forecast_time0, lv_HTGL0, ygrid_0, xgrid_0) float64 -4.939 ...
HGT_P0_L2_GLC0 (forecast_time0, ygrid_0, xgrid_0) float64 869.5 ...
I can load the file and see it's contents, I can also read the data of a given variable with slice method:
data = grbs['DSWRF_P8_L1_GLC0_avg15min'].isel(**{'forecast_time0': 0}).sel(**{'ygrid_0': slice(22,24), 'xgrid_0': slice(-115,-110)})
Output:
0}).sel(**{'ygrid_0': slice(22,24), 'xgrid_0': slice(-115,-110)})
<xarray.DataArray 'DSWRF_P8_L1_GLC0_avg15min' (ygrid_0: 2, xgrid_0: 5)>
array([[ 0., 0., 0., 0., 0.],
[ 0., 0., 0., 0., 0.]])
Coordinates:
forecast_time0 timedelta64[ns] 5 days 15:00:00
gridlat_0 (ygrid_0, xgrid_0) float32 22.4459 22.4396 22.4334 ...
gridlon_0 (ygrid_0, xgrid_0) float32 -75.2096 -75.1823 -75.155 ...
Dimensions without coordinates: ygrid_0, xgrid_0
Attributes:
production_status: Operational products
center: US National Weather Servi...
level: [ 0.]
type_of_statistical_processing: Average
long_name: Downward short-wave radia...
parameter_template_discipline_category_number: [8 0 4 7]
initial_time: 09/06/2017 (06:00)
grid_type: Lambert Conformal can be ...
units: W m-2
statistical_process_duration: 15 minutes (ending at for...
level_type: Ground or water surface
parameter_discipline_and_category: Meteorological products, ...
Now I would like to use the interpolation='nearest' method to retrieve the value of a given variable near a given latitude/longitude:
data = grbs['DSWRF_P8_L1_GLC0_avg15min'].sel(gridlon_0=-75.2096, gridlat_0=22.4396, method='nearest')
Output:
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
<ipython-input-37-0c24d2bdb040> in <module>()
7 #data = grbs['DSWRF_P8_L1_GLC0_avg15min'].sel(gridlon_0=-75.2096, gridlat_0=22.4396, method='nearest') # ValueError: dimensions or multi-index levels ['gridlat_0', 'gridlon_0'] do not exist
8 #data = grbs['DSWRF_P8_L1_GLC0_avg15min'].sel(lat=gridlat_0, lon = gridlon_0) # ValueError: dimensions or multi-index levels ['gridlat_0', 'gridlon_0'] do not exist
----> 9 data = grbs['DSWRF_P8_L1_GLC0_avg15min'].sel(gridlon_0=-75.2096, gridlat_0=22.4396, method='nearest')
/Users/maurice/anaconda3/envs/py27/lib/python2.7/site-packages/xarray/core/dataarray.pyc in sel(self, method, tolerance, drop, **indexers)
690 """
691 pos_indexers, new_indexes = indexing.remap_label_indexers(
--> 692 self, indexers, method=method, tolerance=tolerance
693 )
694 result = self.isel(drop=drop, **pos_indexers)
/Users/maurice/anaconda3/envs/py27/lib/python2.7/site-packages/xarray/core/indexing.pyc in remap_label_indexers(data_obj, indexers, method, tolerance)
275 new_indexes = {}
276
--> 277 dim_indexers = get_dim_indexers(data_obj, indexers)
278 for dim, label in iteritems(dim_indexers):
279 try:
/Users/maurice/anaconda3/envs/py27/lib/python2.7/site-packages/xarray/core/indexing.pyc in get_dim_indexers(data_obj, indexers)
243 if invalid:
244 raise ValueError("dimensions or multi-index levels %r do not exist"
--> 245 % invalid)
246
247 level_indexers = defaultdict(dict)
ValueError: dimensions or multi-index levels ['gridlat_0', 'gridlon_0'] do not exist
Any help will be appreciated!
Note: The files can be found at http://nomads.ncep.noaa.gov/pub/data/nccf/com/hrrr/prod/ then the links with format hrrr.YYYMMDD/hrrr.tHHz.wrfsubhf**.grib2)
Your lat/lon coordinate variables are 2-dimensional. A nearest neighbor lookup like this is not currently possible with xarray. The xarray developers are actively discussing how this may be included in the package but nothing concrete has materialized yet. This github issue is quite relevant: https://github.com/pydata/xarray/issues/475.