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

Posted

技术标签:

【中文标题】使用 xarray 重新采样非标准 CFTimeIndex 日历(360 天,无闰年)以供 pandas 使用的方法【英文标题】:Ways to resample non-standard CFTimeIndex calendars (360-day, no-leap-year) with xarray for pandas usage 【发布时间】:2021-05-17 05:03:38 【问题描述】:

#60198708 带我打开这个问题,因为我还没有找到好的解决方案。

问题

我已经从 EURO-CORDEX 集合下载了几个用于每日降水通量的气候模型。虽然有些型号使用标准日历,兼容 Pandas datetime,但其他型号,特别是 MOHC HadGem2 ES,使用 360 天 CFTimeIndex

主要问题是,如何使用这些日历有效地对月度数据进行重新采样,以便能够对其进行协调并在以后生成整体统计数据。

降水通量数据(2011-2015 节选)可能如下所示 可以here下载。

<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: &lt;MonthEnd: n=1&gt; is a non-fixed frequency 上失败,我没有找到纠正月末的方法 xarray.groupby(time='time.month') 可以处理非标准日历,但是,它的用例是沿不同的轴分组,这是不希望的

我肯定错过了什么,因为这是一个复杂的问题。任何帮助表示赞赏。

【问题讨论】:

如果 1) 您提供示例数据并且 2) 提供预期的输出,您更有可能获得帮助。插值是可以实现的,我们可以提供帮助,但通常不想经历欺骗数据的麻烦,所以请提供给我们:) @anon01 我不得不稍微重写一下我的问题,因为它困扰了我好几天,我想更清楚地说明这一点。所以请再读一遍,主要是到目前为止的解决方案。我会在它们加载后立即附加数据。 @anon01 链接到问题中附加的数据。 【参考方案1】:

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

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

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

【讨论】:

月开始的东西真的很方便。是的,这是我部分希望的。如果我需要每日粒度,我仍然可以坚持我的解决方案:)

以上是关于使用 xarray 重新采样非标准 CFTimeIndex 日历(360 天,无闰年)以供 pandas 使用的方法的主要内容,如果未能解决你的问题,请参考以下文章

为啥 xarray 重新采样均值计算会产生额外的时间?

使用 groupby 重新采样非定时 df

重新采样一个 numpy 数组

1.7 非平衡数据的处理方法大全

基于列标准的熊猫数据框重采样

如何将 pandas Dataframe 时间序列数据从 8hz 重新采样到 16hz?