基于Python统计NetCDF降雨超阈值网格数并排序日期

原代码问题修正 #

原代码中标记极端事件的部分存在错误：所有extreme_Xd变量均错误使用了accum_2d_90p作为阈值，应替换为对应时长的90分位数阈值。同时需移除drop=True，否则会直接删除不满足条件的网格维度，无法保留时间序列完整结构用于后续统计。修正后的代码段：

#=== 修正后的极端事件标记 ===
data['extreme_2d'] = data['precip'].where(data['precip'] > data['accum_2d_90p'])
data['extreme_3d'] = data['precip'].where(data['precip'] > data['accum_3d_90p'])
data['extreme_5d'] = data['precip'].where(data['precip'] > data['accum_5d_90p'])
data['extreme_7d'] = data['precip'].where(data['precip'] > data['accum_7d_90p'])

统计每日极端网格数并导出 #

完整实现代码 #

import numpy as np
import pandas as pd
import xarray as xr

#=== 读取数据 ===
data_path = '/home/wilson/Documents/PH_D/GPCC/GPCC/GPCC_daily_1982-2020.nc'
data = xr.open_dataset(data_path)

#=== 计算多日累积降雨量 ===
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)

#=== 计算各网格点的90分位数 ===
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)

#=== 标记极端事件 ===
data['extreme_2d'] = data['precip'].where(data['precip'] > data['accum_2d_90p'])
data['extreme_3d'] = data['precip'].where(data['precip'] > data['accum_3d_90p'])
data['extreme_5d'] = data['precip'].where(data['precip'] > data['accum_5d_90p'])
data['extreme_7d'] = data['precip'].where(data['precip'] > data['accum_7d_90p'])

#=== 统计每日极端网格数（以2日累积阈值为例，可替换为其他类型） ===
# 统计每个时间点的非NaN网格数
extreme_counts = data['extreme_2d'].count(dim=['lon', 'lat']).to_dataframe(name='Number of grid cells/point')
# 重置索引并格式化日期
extreme_counts.reset_index(inplace=True)
extreme_counts['time'] = pd.to_datetime(extreme_counts['time']).dt.strftime('%Y-%m-%d')
# 按网格数降序排序
extreme_counts_sorted = extreme_counts.sort_values(by='Number of grid cells/point', ascending=False)
# 重命名日期列
extreme_counts_sorted.rename(columns={'time': 'Date'}, inplace=True)

#=== 导出为TXT文件 ===
output_path = '/home/wilson/Documents/PH_D/GPCC/extreme_event_counts.txt'
with open(output_path, 'w') as f:
    # 写入对齐格式的表头
    f.write(f"{'Date':^20} {'Number of grid cells/point':^30}\n")
    f.write('-' * 50 + '\n')
    # 逐行写入数据
    for _, row in extreme_counts_sorted.iterrows():
        f.write(f"{row['Date']:^20} {row['Number of grid cells/point']:^30}\n")

print(f"结果已保存至：{output_path}")

导出结果示例 #

生成的TXT文件内容格式如下：

Date                 Number of grid cells/point        
--------------------------------------------------
      1992-07-01                      432                      
      1983-09-23                      407                      
      2009-08-12                      388                      

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