计算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 必须使用的 dims 和 coords，它不会为你向量化函数。

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 的更详细说明，此可能会有所帮助。

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

Calculating percentile for each gridpoint in xarray

numpy

probability

multidimensional-array

percentile

python-xarray

1。创建一些数据

2。计算分位数

3。计算异常值共现

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