计算 xarray 中每个网格点的百分位数

Question

我目前正在使用 xarray 来制作概率图。我想使用像“计数”练习这样的统计评估。这意味着，对于 NEU 中的所有数据点，计算两个变量共同超过其阈值的次数。这意味着降水数据的第 1 个百分位和温度数据的第 99 个百分位。那么连接发生的概率 (P) 就是联合超出数除以数据集中的数据点数。

<xarray.Dataset>
Dimensions:    (latitude: 88, longitude: 200, time: 6348)
Coordinates:
  * latitude   (latitude) float64 49.62 49.88 50.12 50.38 ... 70.88 71.12 71.38
  * longitude  (longitude) float64 -9.875 -9.625 -9.375 ... 39.38 39.62 39.88
  * time       (time) datetime64[ns] 1950-06-01 1950-06-02 ... 2018-08-31
Data variables:
    rr         (time, latitude, longitude) float32 dask.array<chunksize=(6348, 88, 200), meta=np.ndarray>
    tx         (time, latitude, longitude) float32 dask.array<chunksize=(6348, 88, 200), meta=np.ndarray>
    Ellipsis   float64 0.0

我想计算每个网格点的降水和温度的百分位数，这基本上意味着我想为每个网格点重复下面的函数。

Neu_Precentile=np.nanpercentile(NEU.rr[:,0,0],1)

谁能帮我解决这个问题。我也尝试使用 xr.apply_ufunc 但不幸的是效果不佳。

Answer 1

我不确定你想如何处理分位数，但这里有一个你可以适应的版本。

此外，我选择在计算分位数时保留数据集结构，因为它显示了如何检索异常值的值（如果这是相关的）（它距离检索有效数据点的值仅一步之遥，这可能相关）。

1。创建一些数据

coords = ("time", "latitude", "longitude")
sizes = (500, 80, 120)

ds = xr.Dataset(
    coords={c: np.arange(s) for c, s in zip(coords, sizes)},
    data_vars=dict(
        precipitation=(coords, np.random.randn(*sizes)),
        temperature=(coords, np.random.randn(*sizes)),
    ),
)

查看数据：

<xarray.Dataset>
Dimensions:        (latitude: 80, longitude: 120, time: 500)
Coordinates:
  * time           (time) int64 0 1 2 3 ... 496 497 498 499
  * latitude       (latitude) int64 0 1 2 3 ... 76 77 78 79
  * longitude      (longitude) int64 0 1 2 3 ... 117 118 119
Data variables:
    precipitation  (time, latitude, longitude) float64 -1.673 ... -0.3323
    temperature    (time, latitude, longitude) float64 -0.331 ... -0.03728

2。计算分位数

qt_dims = ("latitude", "longitude")
qt_values = (0.1, 0.9)

ds_qt = ds.quantile(qt_values, dim=qt_dims)

它是一个数据集，分析维度（“纬度”、“经度”）丢失，并具有新的“分位数”维度：

<xarray.Dataset>
Dimensions:        (quantile: 2, time: 500)
Coordinates:
  * time           (time) int64 0 1 2 3 ... 496 497 498 499
  * quantile       (quantile) float64 0.1 0.9
Data variables:
    precipitation  (quantile, time) float64 -1.305 ... 1.264
    temperature    (quantile, time) float64 -1.267 ... 1.254

3。计算异常值共现

对于异常值的位置： （编辑：使用np.logical_and，比& 操作符更易读）

da_outliers_loc = np.logical_and(
    ds.precipitation > ds_qt.precipitation.sel(quantile=qt_values[0]),
    ds.temperature > ds_qt.temperature.sel(quantile=qt_values[1]),
)

输出是一个布尔数据数组：

<xarray.DataArray (time: 500, latitude: 80, longitude: 120)>
array([[[False, ...]]])
Coordinates:
  * time       (time) int64 0 1 2 3 4 ... 496 497 498 499
  * latitude   (latitude) int64 0 1 2 3 4 ... 75 76 77 78 79
  * longitude  (longitude) int64 0 1 2 3 ... 116 117 118 119

如果这些值是相关的：

ds_outliers = ds.where(
    (ds.precipitation > ds_qt.precipitation.sel(quantile=qt_values[0]))
    & (ds.temperature > ds_qt.temperature.sel(quantile=qt_values[1]))
)

4。计算每个时间步的异常值

outliers_count = da_outliers_loc.sum(dim=qt_dims)

最后，这是一个只有时间维度的 DataArray，其值为每个时间戳的异常值数量。

<xarray.DataArray (time: 500)>
array([857, ...])
Coordinates:
  * time     (time) int64 0 1 2 3 4 ... 495 496 497 498 499

Answer 2

np.nanpercentile 默认情况下适用于扁平数组，但是，在这种情况下，目标是仅减少第一个维度，生成包含每个网格点结果的二维数组。为此，可以使用nanpercentile 的axis 参数：

np.nanpercentile(NEU.rr, 1, axis=0)

但是，这将删除标记的尺寸和坐标。这是为了保留 apply_ufunc 必须使用的暗淡和坐标，它不会为您向量化函数。

xr.apply_ufunc(
    lambda x: np.nanpercentile(x, 1, axis=-1), NEU.rr, input_core_dims=[["time"]]
)

注意现在轴是-1，我们使用input_core_dims，它告诉apply_ufunc，这个维度将被缩小，并将其移动到最后一个位置（因此-1）。有关apply_ufunc 的更详细说明，此other answer 可能会有所帮助。