在使用Python的netcdf文件中，当日降雨量大于阈值时，根据网格单元的数量对日期进行计数和排序

Question

我有维度的每日网格降雨量数据(时间：14245，经度：40，纬度：20(。我计算了数据域中每个网格点的2、3、5和7天累积降雨量及其各自的第90个百分位数。我已经使用DataArray设置了我的条件。其中(condition，drop=True(以了解日降雨量何时超过阈值，如下代码所示。我当前的工作代码在这里：

import numpy as np
import pandas as pd
import xarray as xr
#=== reading in the data ===
data_path = '/home/wilson/Documents/PH_D/GPCC/GPCC/GPCC_daily_1982-2020.nc'
data = xr.open_dataset(data_path)
#=== computing 2, 3, 5 and 7-days acummulated rainfall amount ===
data[['precip_2d']] = np.around(data.precip.rolling(time=2).sum(),decimals=2)
data[['precip_3d']] = np.around(data.precip.rolling(time=3).sum(),decimals=2)
data[['precip_5d']] = np.around(data.precip.rolling(time=5).sum(),decimals=2)
data[['precip_7d']] = np.around(data.precip.rolling(time=7).sum(),decimals=2)
#=== Computing 10% largest at each grid point (per grid cel) this is 90th percentile ===
data[['accum_2d_90p']] = np.around(data.precip_2d.quantile(0.9, dim='time'), decimals=2)
data[['accum_3d_90p']] = np.around(data.precip_3d.quantile(0.9, dim='time'), decimals=2)
data[['accum_5d_90p']] = np.around(data.precip_5d.quantile(0.9, dim='time'), decimals=2)
data[['accum_7d_90p']] = np.around(data.precip_7d.quantile(0.9, dim='time'), decimals=2)
#=== locating extreme events, i.e., when daily precip greater than 90th percentile of each of the accumulated rainfall amount ===
data[['extreme_2d']] = data['precip'].where(data['precip'] > data['accum_2d_90p'], drop=True)
data[['extreme_3d']] = data['precip'].where(data['precip'] > data['accum_2d_90p'], drop=True)
data[['extreme_5d']] = data['precip'].where(data['precip'] > data['accum_2d_90p'], drop=True)
data[['extreme_7d']] = data['precip'].where(data['precip'] > data['accum_2d_90p'], drop=True)

我现在的问题是如何计算域中特定日期条件为true的网格单元/点的数量，并使用计数结果按降序排列日期。预期的结果应该看起来像一个可以保存为txt文件的表。例如：cells_count是一个包含所需结果的变量，当print(cells_count(给出时

日期1992-07-014321983-09-234072009-08-12388

Answer 1

好吧，根据您的评论，听起来您想要的是基于3D(时间、纬度、经度(数据的全局汇总统计，在time坐标中获得日期列表。我将以条件(data['precip'] > data['accum_2d_90p'])为例。

我绝对建议在处理日期之前，先降低病情的维度，因为处理粗糙的3D日期时间数组真的很痛苦。因此，既然你提到想要满足标准的像素数，你可以简单地这样做：

global_pixel_count = (
(data['precip'] > data['accum_2d_90p']).sum(dim=("lat", "lon"))
)

现在，您有一个1-D数组global_pixel_count，它仅由time索引。您可以使用xr.DataArray.argsort:从中获得日期的排序顺序

sorted_date_positions = global_pixel_count.argsort()
sorted_dates = global_pixel_count.time.isel(
time=sorted_date_positions.values
)

这将返回一个日期DataArray，按升序对数组进行排序。您可以使用sorted_date_positions.values[::-1]将其反转以按降序获取日期，也可以使用按降序选择所有像素计数

global_pixel_count.isel(time=sorted_date_positions.values[::-1])

您还可以使用以下索引器对整个数组进行索引：

data.isel(time=sorted_date_positions.values[::-1])