#超前滞后相关
import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from scipy.stats import pearsonr, t
# 设置中文字体和负号显示
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False
# ================================
# 1. 数据加载与预处理
# ================================
start_date = '1981-01-01'
end_date = '2022-12-31'
pre = pre_index
sst = sst_summer
# 检查并确保时间对齐
if not np.array_equal(pre.time.values, sst.time.values):
print("时间坐标未对齐,执行时间匹配...")
common_time = np.intersect1d(pre.time.values, sst.time.values)
pre = pre.sel(time=common_time)
sst = sst.sel(time=common_time)
# 数据质量检查
print(f"降水数据时间范围: {pre.time.min().values} 到 {pre.time.max().values}")
print(f"海温数据时间范围: {sst.time.min().values} 到 {sst.time.max().values}")
print(f"共同时间点数量: {len(pre.time)}")
# ================================
# 2. 超前滞后相关分析
# ================================
# 设置滞后参数
lags = range(-60, 1, 5) # 从-60天到0天,步长为5天
n_lags = len(lags)
# 初始化结果数组 - 基于海温空间结构
corr_map = xr.DataArray(
np.full((len(lags), len(sst.lat), len(sst.lon)), np.nan),
dims=['lag', 'lat', 'lon'],
coords={'lag': lags, 'lat': sst.lat, 'lon': sst.lon}
)
p_value_map = corr_map.copy() # 显著性检验结果
# 循环计算每个滞后值的相关系数
for i, lag in enumerate(lags):
print(f"处理滞后 {lag} 天...")
# 应用滞后处理
if lag < 0:
# 海温提前 |lag| 天
sst_lagged = sst.shift(time=lag)
# 截断对齐时间序列
aligned_pre = pre.isel(time=slice(-lag, None))
aligned_sst = sst_lagged.isel(time=slice(-lag, None))
else:
# 正滞后或零滞后
aligned_pre = pre
aligned_sst = sst
# 确保降水数据为一维
aligned_pre = aligned_pre.values.ravel() # 转换为扁平数组
# 计算每个网格点的相关系数
for lat_idx in range(len(sst.lat)):
for lon_idx in range(len(sst.lon)):
# 提取海温网格点时间序列
sst_ts = aligned_sst.isel(lat=lat_idx, lon=lon_idx).values
# 移除可能的缺失值(虽然降水无缺失,但海温可能有)
valid_mask = ~np.isnan(sst_ts) & ~np.isnan(aligned_pre)
sst_valid = sst_ts[valid_mask]
pre_valid = aligned_pre[valid_mask]
# 确保有足够的数据点
if len(pre_valid) > 10:
# 计算相关系数和p值
corr, p_value = pearsonr(pre_valid, sst_valid)
corr_map[i, lat_idx, lon_idx] = corr
p_value_map[i, lat_idx, lon_idx] = p_value
# ================================
# 3. 显著性检验修正
# ================================
# 计算自由度 (n-2)
dof = len(aligned_pre) - 2
# 计算t统计量的临界值
t_critical = t.ppf(1 - 0.025, dof)
# 计算相关系数的临界值
critical_r = t_critical / np.sqrt(t_critical**2 + dof)
# 创建显著性掩码 (p<0.05)
significant = p_value_map < 0.05
# ================================
# 4. 可视化结果(修正后)
# ================================
# 设置投影
proj = ccrs.PlateCarree()
# 创建图形
fig = plt.figure(figsize=(16, 10))
for i, lag in enumerate(lags):
ax = plt.subplot(3, 5, i+1, projection=proj)
# 添加地理特征
ax.add_feature(cfeature.COASTLINE, linewidth=0.5)
ax.add_feature(cfeature.BORDERS, linestyle='-', linewidth=0.5)
# 绘制相关系数
corr_data = corr_map.sel(lag=lag)
im = ax.pcolormesh(
corr_data.lon, corr_data.lat, corr_data,
transform=proj,
cmap='coolwarm',
vmin=-1,
vmax=1
)
# 叠加显著性点
sig_data = significant.sel(lag=lag)
ax.scatter(
corr_data.lon.values[sig_data.values],
corr_data.lat.values[sig_data.values],
s=1,
color='black',
alpha=0.5,
transform=proj
)
# 设置标题
lead_days = -lag
ax.set_title(f"超前 {lead_days} 天", fontsize=10)
ax.set_global()
# 添加整体标题和颜色条
plt.suptitle('海温超前降水相关性分析 (1981-2022)', fontsize=16, y=0.98)
cbar_ax = fig.add_axes([0.15, 0.05, 0.7, 0.02])
cbar = fig.colorbar(im, cax=cbar_ax, orientation='horizontal')
cbar.set_label('相关系数', fontsize=12)
#plt.savefig('corrected_correlation_map.png', dpi=300, bbox_inches='tight')
plt.show()
处理完滞后,在画图的时候报错:
---------------------------------------------------------------------------
IndexError Traceback (most recent call last)
Cell In[5], line 137
134 # 叠加显著性点
135 sig_data = significant.sel(lag=lag)
136 ax.scatter(
--> 137 corr_data.lon.values[sig_data.values],
138 corr_data.lat.values[sig_data.values],
139 s=1,
140 color='black',
141 alpha=0.5,
142 transform=proj
143 )
145 # 设置标题
146 lead_days = -lag
IndexError: too many indices for array: array is 1-dimensional, but 2 were indexed