在多年的数据中获取 SON、DJF、MAM 变量的 95 个百分点

Get 95 percentile of the variables for SON, DJF, MAM over multiple years' data

我有一个名为 ds 的 45 年数据,格式为 netCDF(.nc)。它包含三个坐标:timelatitudelongitude.

print(ds)

<xarray.Dataset>
Dimensions:    (latitude: 106, longitude: 193, time: 403248)
Coordinates:
  * latitude   (latitude) float32 -39.2 -39.149525 ... -33.950478 -33.9
  * longitude  (longitude) float32 140.8 140.84792 140.89584 ... 149.95209 150.0
  * time       (time) datetime64[ns] 1972-01-01 ... 2017-12-31T23:00:00
Data variables:
    FFDI       (time, latitude, longitude) float32 dask.array<shape=(403248, 106, 193), chunksize=(744, 106, 193)>
Attributes:
    creationTime:        1525925611
    creationTimeString:  Wed May  9 21:13:31 PDT 2018
    Conventions:         COARDS

我需要按季节计算 FFDI 的 95 个百分位数,即 SON(9 月、10 月、11 月)、DJF(12 月、1 月、2 月)、MAM(3 月、4 月、5 月)、JJA(6 月、7 月、八月)。

da_ffdi_95th = ds['FFDI'].reduce(np.percentile, dim='time', q=95)

这创建了一个带有百分位数变量的新 DataArray 对象,但删除了时间维度。

groupby 如何与 np.percentile 函数一起使用?

信不信由你,我认为您已经完成了大部分工作!有关详细信息,请参阅 DataArrayGroupBy.reduce

da_ffdi_95th = ds['FFDI'].groupby('time.season').reduce(
    np.percentile, dim='time', q=95)

然而,由于我们使用的是 NumPy 函数,数据将被急切加载。为了使这个 dask 兼容,我们传递给 reduce 的函数必须能够在 NumPy 或 dask 数组上运行。虽然 dask 实现了一个函数,dask.array.percentile, it only operates on 1D arrays, and is not a perfect match to the NumPy function

幸运的是,有了 dask.array.map_blocks,编写我们自己的就足够容易了。这使用 percentile 的 NumPy 实现并将其应用于 dask 数组的每个块;我们唯一需要注意的是确保我们应用它的数组没有沿着我们想要计算百分位数的维度分块。

import dask.array as dask_array

def dask_percentile(arr, axis=0, q=95):
    if len(arr.chunks[axis]) > 1:
        msg = ('Input array cannot be chunked along the percentile '
               'dimension.')
        raise ValueError(msg)
    return dask_array.map_blocks(np.percentile, arr, axis=axis, q=q,
                                 drop_axis=axis)

然后我们可以编写一个包装函数,根据输入数组的类型(NumPy 或 dask)调用适当的 percentile 实现:

def percentile(arr, axis=0, q=95):
    if isinstance(arr, dask_array.Array):
        return dask_percentile(arr, axis=axis, q=q)
    else:
        return np.percentile(arr, axis=axis, q=q)

现在如果我们调用 reduce,确保添加 allow_lazy=True 参数,这个操作 returns 一个 dask 数组(如果底层数据存储在一个 dask 数组中并且是适当分块):

da_ffdi_95th = ds['FFDI'].groupby('time.season').reduce(
    percentile, dim='time', q=95, allow_lazy=True)