带缺值的相关系数及其检验

leoaries

数据：一维时间序列（含缺值），空间三维场（无缺值）

1：Scipy.stats.pearsonr(x,y) 官方说明没有解释输入的数据维数要求，实测应该是只限于两条一维序列

这样计算相关场必须写循环遍历三维场中的空间二维

优势是自动跳过两个序列的缺值，显著性p值直接输出，缺点是对三维以上而言，代码要循环

但如果一维序列和三维场都有缺值，那就只有遍历所有格点这个笨办法了。

2. Xarray.corr 默认计算pearson相关， 无法直接处理NaN值， 需要先去除Nan值（在时间序列和空间场同时去除）

# 去除温度序列中的NaN值，并保持与高度场的时间维度一致

temp_aug_da_valid = temp_aug_da.where(~np.isnan(temp_aug_da), drop=True)

hgt_500_aug_valid = hgt_500_aug.sel(time=temp_aug_da_valid['time'])

# 计算相关系数

#ccr = xr.corr(hgt_500_aug_valid, temp_aug_da_valid, dim='time')

3. Geocat.comp.pearson_r 方案 优势是自动适配时间维，适用多维数据场，自动略过缺值

from geocat.comp import pearson_r

ccr = pearson_r(temp_aug_da, hgt_500_aug, skipna=True) # 优势是自动略过缺值



显著性计算：

上述3方案都是pearson相关系数，对应的显著性计算如下：

from scipy.stats import t # 应用在后面的 't.sf'

#2.显著性计算，对1.1和1.2都适用，因为都是pearson相关

# 计算有效样本数量，即去除NaN后的样本数

valid_count = (~np.isnan(temp_aug_da)).sum().item()

print(valid_count)

# 计算t值

t_values = ccr * np.sqrt((valid_count - 2) / (1 - ccr**2))

# 计算p值 (双尾检验)

p_values = 2 * t.sf(np.abs(t_values), df=valid_count - 2)



完整程序如下：

import xarray as xr

import pandas as pd

import numpy as np

import matplotlib.pyplot as plt

import cartopy.crs as ccrs

import cartopy.feature as cfeature

from scipy.stats import t

#from geocat.comp import pearson_r

#rom scipy import stats



# 打开NetCDF数据集

hgt_ds = xr.open_dataset('/mnt/d/Data/ERA5/pressure/t&hgt_100&500&1000hpa_1979-2023.nc')

hgt_ds = hgt_ds.sel(time=slice('1980-01', '2023-12'))



# 提取500hPa高度场的8月数据

#hgt_500_aug = hgt_ds['z'].sel(level=500).sel(time=hgt_ds['time.month'] == 8)

hgt_500_aug = hgt_ds['z'].sel(level=500).sel(time=hgt_ds['time.month'] == 8)



# 定义南极投影

antarctic_proj = ccrs.SouthPolarStereo()



# 读取站点数据并计算相关性

file_path = '/mnt/d/Data/READER/temp/'

sta_names = ['Byrd', 'Ferrell', 'Mcmurdo', 'Marble_point']

# 站点的经纬度位置，确保顺序为（经度，纬度）

sta_coords = [(-120, -80), (170.8, -77.9), (166.7, -77.9), (163.7, -77.4)]



# 创建一个图形和子图（2x2 排列）

fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(15, 10), subplot_kw={'projection': antarctic_proj})

axes = axes.flatten()



# 定义一个函数来绘制每个站点的相关性图

def plot_correlation(ax, sta_name, hgt_500_aug, file_path, sta_coord):

df = pd.read_csv(f'{file_path}{sta_name}.All.temperature.txt', sep=r'\s+', skiprows=1, header=None, na_values='-')

df.columns = ['Year'] + [f'Month_{j+1}' for j in range(12)]

df = df[df['Year'] != 2023]

df_aug = df[['Month_8']]



# 将温度数据转为 xarray.DataArray，并确保时间对齐

temp_aug_da = xr.DataArray(df_aug['Month_8'].values, dims='time', coords={'time': hgt_500_aug['time']})

#1.2 Geocat.comp 方案 pearson_r

# 计算相关系数

ccr = pearson_r(temp_aug_da, hgt_500_aug, skipna=True) # 优势是自动略过缺值

#1.1 Xarray.corr 方案

# 去除温度序列中的NaN值，并保持与高度场的时间维度一致

#temp_aug_da_valid = temp_aug_da.where(~np.isnan(temp_aug_da), drop=True)

#hgt_500_aug_valid = hgt_500_aug.sel(time=temp_aug_da_valid['time'])

# 计算相关系数

#ccr = xr.corr(hgt_500_aug_valid, temp_aug_da_valid, dim='time')

#2.显著性计算，对1.1和1.2都适用，因为都是pearson相关

# 计算有效样本数量，即去除NaN后的样本数

valid_count = (~np.isnan(temp_aug_da)).sum().item()

print(valid_count)

# 计算t值

t_values = ccr * np.sqrt((valid_count - 2) / (1 - ccr**2))

# 计算p值 (双尾检验)

p_values = 2 * t.sf(np.abs(t_values), df=valid_count - 2)

#3.绘图

# 定义等值线和填色图的级别 (-0.6到0.6，以0.1为间隔)

levels = np.arange(-0.6, 0.7, 0.1)

# 绘制填色图 (correlation)

#ccr.plot.contourf(ax=ax, transform=ccrs.PlateCarree(), cmap='coolwarm', levels=levels, add_colorbar=True)

ccr.plot.contourf(ax=ax, transform=ccrs.PlateCarree(), cmap='coolwarm', levels=levels, extend='both', add_colorbar=True)



# 绘制等值线 (contour lines)

#ccr.plot.contour(ax=ax, transform=ccrs.PlateCarree(), colors='black', levels=levels)



# 设置地图特征

ax.add_feature(cfeature.LAND)

ax.add_feature(cfeature.OCEAN)

ax.add_feature(cfeature.COASTLINE)

ax.add_feature(cfeature.BORDERS, linestyle=':')

ax.set_extent([-180, 180, -90, -60], crs=ccrs.PlateCarree())

# 标记站点位置

ax.plot(sta_coord[0], sta_coord[1], 'k*', transform=ccrs.PlateCarree(), markersize=8,color='yellow')



# 设置标题

#ax.set_title(f'Corr:{sta_name} SAT(Aug) vs Hgt500(Aug)')

ax.set_title(f'Corr:{sta_name} SAT(Aug) vs Hgt500(Jul)')

# 绘制显著性等值线

contour = ax.contour(hgt_500_aug.longitude, hgt_500_aug.latitude, p_values, levels=[0.01], colors='red', transform=ccrs.PlateCarree())

ax.clabel(contour, inline=True, fontsize=8, fmt='%.2f')



# 循环绘制每个站点的相关性图

for ax, sta_name, sta_coord in zip(axes, sta_names, sta_coords):

plot_correlation(ax, sta_name, hgt_500_aug, file_path, sta_coord)



# 调整子图间距

plt.tight_layout()

plt.savefig('ccr_obs.sat(aug).vs.hgt500(jul).pdf',dpi=300)

plt.show()