|
- import numpy as np
-
- # 加载数据
- # data = np.load('data.npz')['arr_0']
- data = sst
-
- n_days, n_lat, n_lon = data.shape
-
- # 按照MHW定义识别热浪事件并计算矩阵信息
- window_size = 11
- ma_window_size = 5
- threshold = np.zeros((n_days, n_lat, n_lon))
- mask = np.zeros((n_days, n_lat, n_lon), dtype=bool)
- events = []
- label, cnt = np.zeros(n_days), 0
-
- for d in range(n_days):
- # 计算11个时间窗口期内的90th百分位
- start_date = max(0, d - window_size // 2)
- end_date = min(n_days, d + window_size // 2 + 1)
- pct = np.percentile(data[start_date:end_date], 90, axis=0)
- print(pct.shape)
- # 计算当前日期的滑动平均
- ma_start_date = max(0, d - ma_window_size + 1)
- ma = np.mean(data[ma_start_date:d + 1], axis=0)
-
- # 当前日期的阈值为历史上该日期90th百分位的平均值,排序后取平均值
- sorted_pct = np.sort(pct, axis=0)
- threshold[d] = np.mean(sorted_pct[-ma_window_size:], axis=0)
- print(threshold[d])
- # 当前日期是否为热浪事件的判断
- mask[d] = data[d] > threshold[d]
-
- if np.all(mask[d]):
- cnt += 1
- label[d] = cnt
- elif cnt >= 5:
- # 符合MHW定义,记录热浪事件信息
- indices = np.nonzero(label)[0]
- start, end = indices[0], indices[-1]
- duration = end - start + 1
- if duration >= 5 and f"热浪事件({start}, {end})" not in events:
- max_intensity = np.max(data[start:end + 1])
- mean_intensity = np.mean(data[start:end + 1])
- cum_intensity = np.sum(data[start:end + 1] - threshold[start:end + 1])
- growth_rate = max_intensity - data[start]
- area = np.sum(mask[start:end + 1])
-
- # 计算当前日期的滑动平均90th百分位
- ma_pct = np.percentile(ma, 90)
-
- event_info = f"热浪事件({start}, {end})信息:持续时间{duration}天,最大强度{max_intensity},平均强度{mean_intensity},积累强度{cum_intensity},增长率{growth_rate},面积{area},滑动平均90th百分位强度{ma_pct}"
-
- events.append(f"热浪事件({start}, {end})")
- print(event_info)
- label, cnt = np.zeros(n_days), 0
- else:
- label, cnt = np.zeros(n_days), 0
- print(f"总共发生了{len(events)}个热浪事件")
|