为什么 xarray reftime 键的末尾突然有一个 1?

Why does the xarray reftime key suddenly have a 1 at the end?

我正在使用 siphon 从 Unidata Thredds 服务器下载 GFS 数据,因此我可以使用 MetPy 绘制它。我写了一个脚本来执行此操作,并且昨天运行良好:

#Get data using siphon
best_gfs = TDSCatalog('http://thredds.ucar.edu/thredds/catalog/grib/NCEP/GFS/Global_0p25deg/catalog.xml?dataset=grib/NCEP/GFS/Global_0p25deg/Best')
best_ds = best_gfs.datasets[0]
ncss = best_ds.subset()
query = ncss.query()
query.lonlat_box(north=55, south=20, east=-60, west=-120).time(datetime.utcnow())
query.accept('netcdf4')
query.variables('Geopotential_height_isobaric')

data = ncss.get_data(query)

#Parse data using MetPy
ds = xr.open_dataset(NetCDF4DataStore(data))
data = ds.metpy.parse_cf()

time_of_run = data['reftime'][0].dt.strftime('%Y%m%d_%H%MZ').values
print(time_of_run)

当我在美国东部时间下午 2 点左右 运行 它时,此代码输出 2020-03-29 12:00Z 并且一切正常。

今天早上我 运行 它时,出现错误:

Traceback (most recent call last):
  File "C:\Users\jacks\Anaconda3\envs\metpy_test\lib\site-packages\xarray\core\dataset.py", line 1155, in _construct_dataarray
    variable = self._variables[name]
KeyError: 'reftime'

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "h5_rh_wind_gph_temp.py", line 51, in <module>
    time_of_run = data['reftime'][0].dt.strftime('%Y%m%d_%H%MZ').values
  File "C:\Users\jacks\Anaconda3\envs\metpy_test\lib\site-packages\xarray\core\dataset.py", line 1245, in __getitem__
    return self._construct_dataarray(key)
  File "C:\Users\jacks\Anaconda3\envs\metpy_test\lib\site-packages\xarray\core\dataset.py", line 1158, in _construct_dataarray
    self._variables, name, self._level_coords, self.dims
  File "C:\Users\jacks\Anaconda3\envs\metpy_test\lib\site-packages\xarray\core\dataset.py", line 165, in _get_virtual_variable
    ref_var = variables[ref_name]
KeyError: 'reftime'

这表明对 'reftime' 键的引用无效。为了调查,我打印了 'data' xarray:

<xarray.Dataset>
Dimensions:                       (isobaric6: 34, lat: 141, lon: 241, time1: 1)
Coordinates:
    reftime1                      (time1) datetime64[ns] ...
  * time1                         (time1) datetime64[ns] 2020-03-30T12:00:00
  * isobaric6                     (isobaric6) float32 40.0 100.0 ... 100000.0
  * lat                           (lat) float32 55.0 54.75 54.5 ... 20.25 20.0
  * lon                           (lon) float32 240.0 240.25 ... 299.75 300.0
    crs                           object Projection: latitude_longitude
Data variables:
    Geopotential_height_isobaric  (time1, isobaric6, lat, lon) float32 ...
    LatLon_Projection             int32 ...
Attributes:
    Originating_or_generating_Center:                                        ...
    Originating_or_generating_Subcenter:                                     ...
    GRIB_table_version:                                                      ...
    Type_of_generating_process:                                              ...
    Analysis_or_forecast_generating_process_identifier_defined_by_originating...
    Conventions:                                                             ...
    history:                                                                 ...
    featureType:                                                             ...
    History:                                                                 ...
    geospatial_lat_min:                                                      ...
    geospatial_lat_max:                                                      ...
    geospatial_lon_min:                                                      ...
    geospatial_lon_max:                                                      ...

表示我想要的信息(模型的运行时间)现在存储为'reftime1'。为什么这个1突然出现在reftime键的末尾?它的出现运行是否有规律地发生?我希望最终 运行 这个脚本作为 cron 作业来自动生成绘图,所以最好能找到一种方法来预测名称的这种变化或完全绕过键名。

reftimereftime1 的变化来自 THREDDS 和 netCDF-java 如何处理基础 GRIB 数据的 netCDF 表示。 GRIB 从根本上作为单独的 2D 数据切片到达。为了创建时间、reftime 和各种垂直维度,netCDF-java 正在查看可用于给定字段(例如 Geopotential_height_isobaric)的 GRIB 消息集。如果字段具有不同的 times/vertical 维度集,则会创建具有唯一名称的单独维度,例如reftimereftime1reftime2。哪个字段以哪个维度名称结尾取决于 netCDF-java 在集合中遇到特定 GRIB 消息的顺序。

使用此方法的方法是避免依赖于名称,而是使用元数据来确定您要查找的内容。 MetPy 可以通过为各种 dimensions/coordinates:

提供适当的别名来做到这一点
# Will point to appropriate time1, time2, etc.
time = data.Geopotential_height_isobaric.metpy.time

这适用于特定变量的坐标。在 reftime 的情况下,因为它不是变量的坐标,您还可以通过查找气候和预测 (CF) 元数据标准名称 forecast_reference_time:

来查找它
filtered_ds = data_.filter_by_attrs(standard_name='forecast_reference_time')

这仍然会留下一个 xarray Dataset,您需要找到一些方法来提取内部唯一的变量——我不确定从那里开始的最佳方法是什么。