基于python的nc文件偏相关度量

前言

偏相关是一种统计度量，用于衡量两个变量之间的相关性，同时控制一个或多个其他变量的影响。在多元变量分析中，偏相关可以帮助我们了解变量间的直接关系，而不受其他变量的干扰。

一、Pingouin.partial_corr 参数说明

data: Pandas DataFrame 数据集。
x, y: 字符串，指定需要计算偏相关的两个变量的列名。
covar: 字符串或列表，指定一个或多个协变量的列名，这些协变量在计算偏相关时会被控制。
x_covar: 字符串或列表，指定仅从 x 变量中去除影响的协变量。使用此参数时，即计算半偏相关。
y_covar: 字符串或列表，指定仅从 y 变量中去除影响的协变量。使用此参数时，即计算半偏相关。
alternative: 字符串，定义备择假设的尾部，可以是 “two-sided”（默认）、“greater” 或 “less”。
method: 字符串，指定相关性类型，可以是 “pearson”（皮尔逊相关系数）或 “spearman”（斯皮尔曼相关系数）。

返回一个包含以下内容的 Pandas DataFrame：

‘n’: 样本大小(去除缺测值后)。
‘r’: 偏相关系数。
‘CI95’: 95% 置信区间。
‘p-val’: p 值。
注意事项
函数会自动去除含有缺失值的行。
偏相关系数的取值范围是 -1 到 1，表示完美的负相关到完美的正相关。
半偏相关与偏相关类似，但控制变量集仅从 x 或 y 中去除，而不是两者都去除

二、代码

# 导入必要的库
import xarray as xr  # 用于处理多维数组的数据结构
import numpy as np  # 提供支持大规模矩阵运算和数学函数
from matplotlib import pyplot as plt  # 用于绘图
import cartopy.crs as ccrs  # 用于地图投影和地理数据可视化
from tqdm import tqdm  # 用于显示循环进度条
import pandas as pd  # 用于数据处理和分析
import pingouin as pg  # 用于统计分析
import warnings  # 用于控制警告信息

# 忽略 RuntimeWarning 警告
warnings.filterwarnings('ignore', category=RuntimeWarning)

# 定义一个创建图形的函数
def createFigure(figsize=(12, 8), dpi=300, subplotAdj=None, **kwargs):
    """
    创建一个 Matplotlib 图形对象，并可选择性地调整子图参数。
    
    参数:
    figsize (tuple): 图形的宽和高
    dpi (int): 每英寸的点数
    subplotAdj (dict): 子图参数调整
    **kwargs: 其他 Matplotlib 图形参数
    
    返回:
    figure: 创建的图形对象
    """
    # 设置图形的尺寸和分辨率
    figure = plt.figure(figsize=figsize, dpi=dpi, **kwargs)
    # 如果提供了子图调整参数，则应用这些调整
    if subplotAdj is not None:
        plt.subplots_adjust(**subplotAdj)
    return figure

# 打开并读取 "apo.nc" 文件中的数据，并提取 'apo0' 变量的值
apo = xr.open_dataset("apo.nc")['apo0'].values

# 定义要分析的月份
month = [6, 7, 8]

# 读取 "NINO3.4_index.nc" 文件中的数据，并选择指定月份的时间范围
nino = xr.open_dataset("NINO3.4_index.nc").where(
    lambda m: m.time.dt.month.isin(month), drop=True).groupby(
        "time.year").mean(
            "time")['NINO3.4'].sel(year=slice(1980, 2014)).values
            
# 读取 "sst.mnmean.nc" 文件中的数据，并选择指定月份的时间范围
sst = xr.open_dataset("sst.mnmean.nc").where(
    lambda m: m.time.dt.month.isin(month)).groupby(
        "time.year").mean(
        "time")["sst"].sel(year=slice(1980, 2014))

# 提取经度和纬度的值
lat, lon = sst.lat.values, sst.lon.values
# 提取 'sst' 数据的值
sst = sst.values
# 初始化相关系数矩阵和 p 值矩阵
r_matrix, p_matrix = [np.zeros_like(sst[0]) for _ in range(2)]

# 对每个纬度和经度的网格点计算部分相关系数
for i in tqdm(range(len(lat)), desc='Calc apo and sst partial corr...'):
    for j in range(len(lon)):
        # 检查是否存在缺失值
        if np.isnan(sst[:, i, j]).any():
            continue
        else:
            # 创建包含 'apo'、'Nino3_4' 和 'sst' 的 DataFrame
            df = pd.DataFrame({
                'apo': apo,
                'Nino3_4': nino,
                'sst': sst[:, i, j]
            })
            # 计算 'sst' 和 'apo' 的部分相关系数，控制变量为 'Nino3_4'
            pc = pg.partial_corr(data=df, x='sst', y='apo', covar='Nino3_4', method='pearson')
            # 将部分相关系数和 p 值存储到矩阵中
            r_matrix[i, j] = pc['r'][0]  # 部分相关系数
            p_matrix[i, j] = pc['p-val'][0]  # p 值

# 将 p 值为 0 的元素替换为 NaN
p_matrix[p_matrix == 0] = np.nan

# 创建图形对象，设置图形大小、分辨率以及子图调整参数
fig = createFigure(figsize=(12, 8), dpi=300, 
             subplotAdj=dict(left=0.04, right=0.98,
                             top=0.9, bottom=0.05, 
                             wspace=0.05, hspace=0.1))

# 创建一个地图投影对象，使用 Robinson 投影
ax = plt.subplot(111, projection=ccrs.Robinson(central_longitude=180))
# 设置全局显示
ax.set_global()
# 绘制海岸线
ax.coastlines()

# 绘制部分相关系数的填充等高线图
CF = plt.contourf(lon, lat, 
                  r_matrix, 
                  cmap="RdBu_r",  # 颜色映射
                  levels=np.linspace(-0.5, 0.5, 21),  # 等高线级别
                  extend='both',  # 颜色条扩展
                  transform=ccrs.PlateCarree())  # 坐标转换

# 绘制 p 值的填充等高线图，表示显著性水平
c1b = ax.contourf(lon, lat, 
                  p_matrix, 
                  levels=[0,0.05,0.5], 
                  zorder=1,  # 绘制顺序
                  hatches=['..', None],  # 图案填充
                  colors="none", transform=ccrs.PlateCarree())

# 设置图标题
plt.title('Partial Correlation between APO and SST', 
           fontweight='bold', 
           fontsize=20)

# 创建颜色条
position = fig.add_axes([0.31, 0.05,  0.4, 0.025])
fig.colorbar(CF, cax=position, orientation='horizontal', format='%.1f',)

# 显示图形
plt.show()

出图效果

参考文献

赵平,代玮 & 肖子牛.(2011).亚洲—太平洋涛动研究进展.气象科技进展(02),6-10.