numpy - 计算 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.ological_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]),
)

输出是一个 bool 数据数组:

<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