使用xarray重新采样非标准CFTimeIndex日历(360天,无闰年)以供pandas使用的方法



#60198708让我开始思考这个问题,因为我还没有找到好的解决方案。

问题我从EURO-CORDEX集合中下载了几个daily的气候模型。前次通量。虽然一些模型使用标准日历,兼容Pandasdatetime,但其他模型,特别是MOHC HadGem2 ES,使用360天CFTimeIndex

主要的问题是,如何有效地用这些日历重新采样每月的数据,以便能够协调它,并在以后产生综合统计。

降水通量数据(2011-2015年节选)如下图你可以在这里下载。

<xarray.Dataset>
Dimensions:       (bnds: 2, rlat: 412, rlon: 424, time: 1800)
Coordinates:
lat           (rlat, rlon) float64 ...
lon           (rlat, rlon) float64 ...
* rlat          (rlat) float64 -23.38 -23.26 -23.16 ... 21.61 21.73 21.83
* rlon          (rlon) float64 -28.38 -28.26 -28.16 ... 17.93 18.05 18.16
* time          (time) object 2011-01-01 12:00:00 ... 2015-12-30 12:00:00
Dimensions without coordinates: bnds
Data variables:
pr            (time, rlat, rlon) float32 ...
rotated_pole  |S1 ...
time_bnds     (time, bnds) object ...
Attributes:
CDI:                            Climate Data Interface version 1.3.2
Conventions:                    CF-1.6
NCO:                            4.4.2
CDO:                            Climate Data Operators version 1.3.2 (htt...
contact:                        Fredrik Boberg, Danish Meteorological Ins...
creation_date:                  2019-11-16 14:39:25
experiment:                     Scenario experiment using HadGEM as drivi...
experiment_id:                  rcp45
driving_experiment:             MOHC-HadGEM2-ES,rcp45,r1i1p1
driving_model_id:               MOHC-HadGEM2-ES
driving_model_ensemble_member:  r1i1p1
driving_experiment_name:        rcp45
frequency:                      day
institution:                    Danish Meteorological Institute
institute_id:                   DMI
model_id:                       DMI-HIRHAM5
rcm_version_id:                 v2
project_id:                     CORDEX
CORDEX_domain:                  EUR-11
product:                        output
tracking_id:                    hdl:21.14103/158e462e-499c-4d6e-8462-ac3e...
c3s_disclaimer:                 This data has been produced in the contex...

可以看到,数据集的时间维为cftime.Datetime360Day。所有月份都是30天,这有时对气候预测很好,但对pandas不是。

<xarray.DataArray 'time' (time: 1800)>
array([cftime.Datetime360Day(2011-01-01 12:00:00),
cftime.Datetime360Day(2011-01-02 12:00:00),
cftime.Datetime360Day(2011-01-03 12:00:00), ...,
cftime.Datetime360Day(2015-12-28 12:00:00),
cftime.Datetime360Day(2015-12-29 12:00:00),
cftime.Datetime360Day(2015-12-30 12:00:00)], dtype=object)
Coordinates:
* time     (time) object 2011-01-01 12:00:00 ... 2015-12-30 12:00:00
Attributes:
standard_name:  time
long_name:      time
bounds:         time_bnds

我所做的一切

我走了肮脏的方式,将CFTimeIndex转换为字符串,放入pandas.DataFrame并将时间转换为pd.to_datetimeerrors=coerce

ds = xarray.open_dataset('data/mohc_hadgem2_es.nc')
def cft_to_string(cfttime_obj):
month = str(cfttime_obj.month)
day = str(cfttime_obj.day)
# This is awful but there were no two-digit months/days by default
month = '0'+month if len(month)==1 else month
day = '0'+day if len(day)==1 else day
return f'{cfttime_obj.year}-{month}-{day}'
# Apply above function
ds_time_strings = list(map(cft_to_string, ds['time']))
# Get precipitation values only (to use in pandas dataframe)
# Suppose the data are from multiple pixels (for whole of Europe)
# - that's why the mean(axis=(1,2))
precipitation = ds['pr'].values.mean(axis=(1,2))
# To dataframe
df = pd.DataFrame(index=ds_time_strings, data={'precipitation': precipitation})
# Coerce erroneous dates
df.index = pd.to_datetime(df.index, errors='coerce') # Now, dates such as 2011-02-30 are omitted

这给出了一个非标准日期作为NaT的数据帧,并且缺少一些日期(31天)。我不介意,因为我做了90年的预测。

precipitation
2011-01-01  0.000049
2011-01-02  0.000042
2011-01-03  0.000031
2011-01-04  0.000030
2011-01-05  0.000038
... ...
2011-02-28  0.000041
NaT         0.000055
NaT         0.000046
2011-03-01  0.000031
... ...
2015-12-26  0.000028
2015-12-27  0.000034
2015-12-28  0.000028
2015-12-29  0.000025
2015-12-30  0.000024
1800 rows × 1 columns

现在我可以使用pandas轻松地重新采样到每月数据。

虽然这似乎工作,是否有一个更干净的方式与xarray/pandas ?可能不是基于字符串的?

  • ds.indexes['time'].to_datetimeindex()在非标准日历上失败
  • ds.resample(time='M')将进行重新采样,但是,它产生非标准月末。我没有找到正确月末的方法,因为ds['time'].dt.floor('M')ValueError: <MonthEnd: n=1> is a non-fixed frequency
  • 上失败了。
  • xarray.groupby(time='time.month')可以处理非标准日历,然而,它的用例是沿着不同的轴分组,这是不希望的

我肯定错过了什么,因为这是一个复杂的问题。谢谢你的帮助。

感谢详细的示例!如果你的分析可以接受月均值的时间序列,我认为最干净的方法是重新抽样到"月开始"。频率,然后协调日期类型,例如,对于CFTimeIndex索引的数据集,类似于:

resampled = ds.resample(time="MS").mean()
resampled["time"] = resampled.indexes["time"].to_datetimeindex()

这基本上是你的第二个要点,但有一个小的变化。重新采样到月开始频率可以解决360天日历中包含标准日历中不存在的月结束的问题,例如2月30日。

相关内容

最新更新