Python滞后相关系数(Lagged correlation)代码分享,气象相关

本文分享了在气象研究中使用Python计算滞后相关系数的方法,适用于南海海温等相关数据分析。涉及的库包括catorpy、numpy和scipy。

最近研究南海海温时用到了这个,向大家进行一个代码分享

需要:catorpy,numpy,scipy

import matplotlib.pyplot as plt
import numpy as np
from scipy.stats import pearsonr

def LCC(data1, data2, n, path):
#   作滞后相关,n>0时,data2先于data1发生
#   如data2为海温,data1为降水,n为1——>LCC为海温关于当年降水和来年降水的相关
    a = n * 2 + 1
    b = len(data1)
    x = np.arange(-n, n + 1, 1)
    r = np.zeros((a))
    p = np.zeros((a))

    for i in range(a):
        if i < (n):
            r[n - i - 1], p[n - i - 1] = pearsonr(data1[:(b - i - 1)], data2[i + 1:])
        else :
            r[i], p[i] = pearsonr(data1[x[i]:], data2[:b - x[i]])
#   画图
    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.axhline(0,linewidth=3, color='k')   #零刻度线
    cs = ax.bar(x[n:], r[n:], align='center',edgecolor='k',linewidth=3,color = 'white')
#            从无滞后相关开始绘图,对齐方式为中心对齐
    count = n
    for bar, height in zip(cs, r):         #开始显著性检验
        if p[count] < 0.01:
            bar.set(hatch='///')         
#绘制超前滞后相关 import numpy as np import matplotlib.pyplot as plt import cartopy.crs as ccrs import cartopy.feature as cfeature import xarray as xr # 坐标转换(保持不变) results = results.assign_coords(lon=(((results.lon + 180) % 360) - 180)) results = results.sortby('lon') # 创建图形(增加尺寸和调整布局) fig = plt.figure(figsize=(18, 15)) nrows, ncols = 4, 4 cbar_ax = fig.add_axes([0.15, 0.07, 0.7, 0.02]) vmin, vmax = -0.4, 0.4 sig_level = 0.05 # 显著性水平 dot_size = 1.5 sparse_stride = 3 # 稀疏采样步长,值越大点越稀疏[^1] # 创建投影 proj = ccrs.PlateCarree(central_longitude=180) # 循环绘制子图 for i, lag in enumerate(results.lag.values): ax = fig.add_subplot(nrows, ncols, i+1, projection=proj) # 设置地图范围 ax.set_global() # 绘制相关系数场 corr_map = results.correlation.sel(lag=lag) im = corr_map.plot.pcolormesh( ax=ax, cmap="RdBu_r", vmin=vmin, vmax=vmax, add_colorbar=False, transform=ccrs.PlateCarree() ) # 创建显著性掩膜 sig_mask = results.p_value.sel(lag=lag) < sig_level # 稀疏采样显著点 sparse_mask = np.zeros_like(sig_mask, dtype=bool) sparse_mask[::sparse_stride, ::sparse_stride] = True # 组合条件:既显著又在采样点上 combined_mask = np.logical_and(sig_mask, sparse_mask) # 获取稀疏点的经纬度 sig_points = np.where(combined_mask) sig_lat = corr_map.lat.values[sig_points[0]] sig_lon = corr_map.lon.values[sig_points[1]] # 添加显著性点(统一使用星形标记) ax.scatter( sig_lon, sig_lat, s=dot_size , # 稍大一些的点 edgecolors='none', # 移除边缘线条 marker='o', color='k', transform=ccrs.PlateCarree() ) # 添加地图特征 ax.add_feature(cfeature.COASTLINE, linewidth=0.5) # 添加网格线和标签 gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True, linewidth=0.5, color='gray', alpha=0.5) gl.top_labels = False gl.right_labels = False # 设置标题 ax.set_title(f"Lag = {lag}d", fontsize=12) # 添加共享色标 cbar = fig.colorbar(im, cax=cbar_ax, orientation='horizontal', extend='both') cbar.set_label('Correlation Coefficient', fontsize=12) # 调整布局 plt.tight_layout(rect=[0, 0.1, 1, 0.98]) plt.subplots_adjust(bottom=0.15, hspace=0.2, wspace=0.2) plt.suptitle("Snow_Cover-Precipitation Lagged Correlation Analysis", fontsize=16, y=0.99) # 保存和显示 #plt.savefig("/online1/hjwang_group/zhoufang/Students/ningziyi/olr-pre_lag=75png", dpi=800, bbox_inches='tight') plt.show() 我的积雪数据纬度为20-70N,调整上述代码
最新发布
11-28
#超前滞后相关 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
11-25
评论 3
成就一亿技术人!
拼手气红包6.0元
还能输入1000个字符
 
红包 添加红包
表情包 插入表情
 条评论被折叠 查看
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值