使用 xarray 读取 grib2 文件
Read grib2 file with xarray
我需要用 xarray 打开一个 grib2 文件。为此，我在 xarray:
中使用 python 2.7 和 pynio 作为引擎
grbs = xr.open_dataset('hrrr.t06z.wrfsubhf02.grib2'], engine = 'pynio')
输出：
<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 ...
我可以加载文件并查看它的内容，我还可以使用 slice 方法读取给定变量的数据：
data = grbs['DSWRF_P8_L1_GLC0_avg15min'].isel(**{'forecast_time0': 0}).sel(**{'ygrid_0': slice(22,24), 'xgrid_0': slice(-115,-110)})
输出：
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, ...
现在我想使用 interpolation='nearest' 方法来检索给定变量附近给定 latitude/longitude:
的值
data = grbs['DSWRF_P8_L1_GLC0_avg15min'].sel(gridlon_0=-75.2096, gridlat_0=22.4396, method='nearest')
输出：
---------------------------------------------------------------------------
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
任何帮助将不胜感激！
注意：文件可以在 http://nomads.ncep.noaa.gov/pub/data/nccf/com/hrrr/prod/ 找到，然后是格式为 hrrr.YYYMMDD/hrrr.tHHz.wrfsubhf**.grib2)
的链接
您的 lat/lon 坐标变量是二维的。 xarray 目前无法实现像这样的最近邻查找。 xarray 开发人员正在积极讨论如何将其包含在包中，但尚未具体化。这个 github 问题非常相关：https://github.com/pydata/xarray/issues/475.
我需要用 xarray 打开一个 grib2 文件。为此，我在 xarray:
中使用 python 2.7 和 pynio 作为引擎grbs = xr.open_dataset('hrrr.t06z.wrfsubhf02.grib2'], engine = 'pynio')
输出：
<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 ...
我可以加载文件并查看它的内容，我还可以使用 slice 方法读取给定变量的数据：
data = grbs['DSWRF_P8_L1_GLC0_avg15min'].isel(**{'forecast_time0': 0}).sel(**{'ygrid_0': slice(22,24), 'xgrid_0': slice(-115,-110)})
输出：
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, ...
现在我想使用 interpolation='nearest' 方法来检索给定变量附近给定 latitude/longitude:
的值 data = grbs['DSWRF_P8_L1_GLC0_avg15min'].sel(gridlon_0=-75.2096, gridlat_0=22.4396, method='nearest')
输出：
---------------------------------------------------------------------------
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
任何帮助将不胜感激！
注意：文件可以在 http://nomads.ncep.noaa.gov/pub/data/nccf/com/hrrr/prod/ 找到，然后是格式为 hrrr.YYYMMDD/hrrr.tHHz.wrfsubhf**.grib2)
的链接您的 lat/lon 坐标变量是二维的。 xarray 目前无法实现像这样的最近邻查找。 xarray 开发人员正在积极讨论如何将其包含在包中，但尚未具体化。这个 github 问题非常相关：https://github.com/pydata/xarray/issues/475.