超大规模逐日NETCDF数据逐值百分位计算方案咨询

问题背景 #

我有存储逐日实测数据的NetCDF文件，数据为三维格式（time, lat, lon），规模极大：涵盖全球150年、0.1°空间分辨率，包含10个变量×25个集合×54750天×3600个经度格点×1500个纬度格点。计划按变量和集合逐一处理，计算每个数据值对应的百分位并替换原值后写入NetCDF，但尝试四种方法均失败：

尝试1：xarray直接计算分位数导致维度异常 #

xarrayelems = []
startdatestr = '1953-01-03'
enddatestr = '2014-12-31'
startdate = datetime.datetime.strptime(startdatestr, '%Y-%m-%d')
enddate = datetime.datetime.strptime(enddatestr, '%Y-%m-%d')
x =0
while startdate<=enddate:
    try:
        vv = xr.open_mfdataset(inpath)
        xx = vv.isel(ensemble=1).isel(time=0)['Evap_tavg']
        xarrayelems.append(xx)
        startdate = startdate + datetime.timedelta(days=1)
    except:
        startdate = startdate + datetime.timedelta(days=1)
        continue
    
combined = xr.concat(xarrayelems, dim='time')

dates = pd.date_range(start='1/3/1953', end='12/31/2014')
print(dates)
combined['time'] = dates
combined['east_west'] = lon[:360]
combined['north_south'] = lat[:300]

quantiles = np.arange(0.01,1,0.01)
print(quantiles)
percentiles = combined.chunk(dict(time=len(combined.time))).quantile(q=quantiles, dim='time')

newds = xr.where(combined.notnull(), percentiles, combined)

测试小数据集时，结果维度错误，新增99个百分位维度，不符合预期。

尝试2：xarray.apply_ufunc调用自定义函数异常 #

def replace_with_percentile(value, percentiles):
    return np.interp(value, percentiles, np.arange(0, 101))

# Apply the function to each element of the original dataset
ds_percentiles = xr.apply_ufunc(
    replace_with_percentile, combined, percentiles,
    input_core_dims=[['time'], []],
    dask='allowed',
    output_dtypes=[float]
)

出现无法理解的异常。

尝试3：直接调用scipy.stats.percentileofscore结果全NaN #

from scipy import stats
percentiles =  stats.percentileofscore(combined, combined, kind='rank')

运行耗时极长，输出全为NaN（包含无NaN的格点）。

尝试4：嵌套循环效率极低且维度不匹配 #

percentiles = np.empty_like(combined)
print("Starting . . .")
for i in range(combined.shape[1]):
    for j in range(combined.shape[2]):
        print(j)
        percentiles[:, i, j] = stats.percentileofscore(combined[:, i, j], combined[:, i, j], kind='rank')

# Create a new DataArray with percentiles
percentile_data = xr.DataArray(
    percentiles,
    dims=combined.dims,
    coords=combined.coords,
    name="percentile_temperature"
)

# Create a new dataset with only the percentile data
ds_percentile = xr.Dataset({"percentile_temperature": percentile_data})

嵌套循环效率极低，无法完成计算；后续xr.apply_ufunc调用出现维度不匹配问题。

咨询问题 #

1. 处理此类超大规模数据的最高效方案（可使用pandas、xarray等主流库）；
2. 若使用stats.percentileofscore，如何忽略NaN值？能否通过xarray.apply_ufunc实现？

解决方案 #

1. 超大规模数据的最高效处理方案 #

针对这种**(time, lat, lon)的三维格点数据，核心思路是利用Dask+xarray的分块并行计算**，避免一次性加载全量数据到内存，同时针对每个空间格点（lat, lon）计算时间序列的百分位排名，保持原数据维度不变。

核心步骤：

• 分块加载数据：用xr.open_mfdataset时指定chunks参数，将数据按time、lat、lon合理分块（比如time分块为365天，lat/lon分块为100×100），让Dask自动处理并行和内存管理。
• 自定义百分位计算函数：针对每个格点的时间序列，计算每个值在序列中的百分位排名，忽略NaN。
• 用xr.apply_ufunc批量并行计算：指定input_core_dims和vectorize=True，让函数自动遍历每个空间格点的时间序列。

示例代码：

import xarray as xr
import numpy as np
from scipy import stats

# 1. 分块加载单个变量+集合的数据
# 假设inpath是对应单个变量、单个集合的文件路径（按变量+集合逐一处理）
ds = xr.open_mfdataset(inpath, chunks={'time': 365, 'lat': 100, 'lon': 100})
var_data = ds['Evap_tavg']  # 替换为目标变量名

# 2. 定义忽略NaN的百分位计算函数
def calc_percentile_rank(ts):
    # ts是单个格点的时间序列（一维数组）
    mask = ~np.isnan(ts)
    if not np.any(mask):
        return np.full_like(ts, np.nan)
    # 计算每个有效值在非NaN序列中的百分位排名
    ranks = stats.percentileofscore(ts[mask], ts[mask], kind='rank')
    # 还原到原序列长度，NaN位置保持NaN
    result = np.full_like(ts, np.nan)
    result[mask] = ranks
    return result

# 3. 用apply_ufunc并行计算
percentile_da = xr.apply_ufunc(
    calc_percentile_rank,
    var_data,
    input_core_dims=[['time']],  # 按time维度逐格点处理
    output_core_dims=[['time']],  # 输出保持time维度
    vectorize=True,  # 自动遍历lat和lon维度
    dask='parallelized',  # 启用Dask并行
    output_dtypes=[float]
)

# 4. 保存结果到NetCDF
# 保留原数据集的坐标和属性
output_ds = ds.copy()
output_ds['Evap_tavg_percentile'] = percentile_da
output_ds.to_netcdf('output_percentile.nc', compute=True)

优化点：

• 分块策略：根据你的内存大小调整chunks，确保每个分块能被内存容纳；time分块不宜过大，lat/lon分块尽量均衡。
• 逐一处理变量和集合：按变量+集合循环处理，避免同时加载过多数据。
• 利用集群资源：如果有计算集群，可以配置Dask分布式集群，进一步提升计算速度。

2. stats.percentileofscore忽略NaN及apply_ufunc实现 #

• 忽略NaN的方法：在自定义函数中先过滤NaN值，只对非NaN序列计算百分位排名，再将结果映射回原序列的非NaN位置。
• apply_ufunc实现：如上述示例代码所示，通过input_core_dims=[['time']]指定按时间维度处理每个格点，vectorize=True自动遍历空间维度，同时结合Dask实现并行计算。

关键说明：

• stats.percentileofscore的第一个参数需要是非NaN的一维数组，第二个参数是要计算排名的数值（这里是同一数组的每个元素）。
• 函数中先判断序列是否全NaN，避免报错；然后用mask提取有效值计算排名，再将结果填充回原序列的对应位置。

内容的提问来源于stack exchange，提问作者nishanuw