想用滑动阈值法定义极端降水事件，但是结果都是0，不知道哪里出了问题，求各位...

饼饼

import numpy as np
import netCDF4 as nc
import csv
from datetime import datetime

def process_extreme_precip(file_path):
    dataset = nc.Dataset(file_path)

    precipitation = dataset.variables['tp'][:]
    times = dataset.variables['time'][:]

    extreme_precip_days_sum = np.zeros((4140 // 30))  # 存储每年的极端降水事件总数

    def calculate_90th_percentile(data):
        data_above_01 = data[data > 0.1]
        if len(data_above_01) > 0:
            sorted_data = np.sort(data_above_01)
            percentile_index = int(len(sorted_data) * 0.9)
            return sorted_data[percentile_index]
        else:
            return np.nan

    window_size = 11
    for t in range(window_size//2, precipitation.shape[0] - window_size//2):
        # 计算当前年份
        year = datetime.strptime(str(times[t])[:4], '%Y').year
        year_index = year - 1979  # 计算在数组中的索引
        for i in range(precipitation.shape[1]):
            for j in range(precipitation.shape[2]):
                window_data = precipitation[t - window_size//2:t + window_size//2 + 1, i, j]
                window_data = window_data[window_data > 0.1]
                if len(window_data) > 0:
                    threshold = calculate_90th_percentile(window_data)
                    if precipitation[t, i, j] > threshold:
                        extreme_precip_days_sum[year_index] += 1

    csv_file_path = r'C:\Users\Lenovo\Desktop\extreme_precip_days_sum1.csv'
    with open(csv_file_path, 'w', newline='') as csvfile:
        writer = csv.writer(csvfile)
        writer.writerow(['year', 'extreme_precip_total'])  # 修改标题行

        for year, total in enumerate(extreme_precip_days_sum, start=1979):
            writer.writerow([year, total])

    dataset.close()

if __name__ == "__main__":
    file_path = r'C:\Users\Lenovo\Desktop\extreme preciptation\ERA5_China_1979-2023_summer_dailysum_precipitation.nc'
    process_extreme_precip(file_path)


另外也有一个代码是拉长了时间窗口的但结果一样
import numpy as np
import netCDF4 as nc
import csv
from datetime import datetime

# 读取nc文件
file_path = r'C:\Users\Lenovo\Desktop\extreme preciptation\ERA5_China_1979-2023_summer_dailysum_precipitation.nc'
dataset = nc.Dataset(file_path)

# 获取降水数据，假设变量名为'precipitation'
precipitation = dataset.variables['tp'][:]

# 获取时间维度数据，假设变量名为'time'
times = dataset.variables['time'][:]

# 初始化极端降水天数总和数组
extreme_precip_days_sum = np.zeros((precipitation.shape[1], precipitation.shape[2], 4140 // 30))

# 定义一个函数来计算滑动窗口内大于0.1的降水数据的第90个百分位数
def calculate_90th_percentile(data):
    # 只对大于0.1的数据进行排序
    data_above_01 = data[data > 0.1]
    if len(data_above_01) > 0:
        sorted_data = np.sort(data_above_01)
        percentile_index = int(len(sorted_data) * 0.9)
        return sorted_data[percentile_index]
    else:
        return np.nan

# 计算滑动阈值
window_size = 11
for i in range(precipitation.shape[1]):
    for j in range(precipitation.shape[2]):
        for t in range(window_size//2, precipitation.shape[0] - window_size//2):
            # 计算窗口内的数据
            window_data = precipitation[t-window_size//2:t+window_size//2+1, i, j]
            window_data = window_data[window_data > 0.1]  # 过滤出大于0.1的数据
            if len(window_data) > 0:
                # 计算第90个百分位数
                threshold = calculate_90th_percentile(window_data)
                # 检查是否为极端降水事件
                if precipitation[t, i, j] > threshold:
                    extreme_precip_days_sum[i, j, t//30] += 1

# 输出结果到CSV文件
csv_file_path = r'C:\Users\Lenovo\Desktop\extreme_precip_days_sum.csv'
with open(csv_file_path, 'w', newline='') as csvfile:
    writer = csv.writer(csvfile)
    # 写入标题行
    writer.writerow(['lat', 'lon', 'year', 'extreme_precip_days'])

    # 遍历每个格点
    for i in range(precipitation.shape[1]):
        for j in range(precipitation.shape[2]):
            for t in range(precipitation.shape[0] // 30):
                # 获取当前年份
                year = datetime.strptime(str(times[t*30])[:4], '%Y').year
                # 写入每个格点的数据
                writer.writerow([i, j, year, extreme_precip_days_sum[i, j, t]])

# 关闭文件
dataset.close()

llyxxn.

码一下，看看后续

生气的粉色条纹

llyxxn. 发表于 2024-7-17 15:14
码一下，看看后续

谢谢，已经解决了，主要是浮点数的问题