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

Question

我有 45 年的数据，名为 ds，格式为 netCDF(.nc)。它包含三个坐标:时间、纬度和经度。

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 月、8 月) .

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

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

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

Answer 1

信不信由你，我认为您已经完成了大部分工作!参见 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 , 它只在一维数组上运行，并且 is not a perfect match to the NumPy function .

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

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 参数，此操作返回一个 dask 数组(如果底层数据存储在 dask 数组中并适本地分块):

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