基于Python对三维数据做空间相关分析

前言

读取两个 .nc 格式的气候数据文件 (SPEI.nc 和 SRI.nc)，并进行相关性分析。具体来说，代码从数据文件中提取特定变量，计算这些变量与另一个变量之间的相关系数，并进行显著性检验。最终，相关系数和显著性结果以全球地图的形式可视化展示。

代码

import numpy as np
import xarray as xr
from sklearn.feature_selection import f_regression
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
import cmaps

# 使用上下文管理器读取 .nc 数据，并提取数据中的变量
# 可以提前用 NASA 的 Panoply 软件查看 .nc 文件信息
with xr.open_dataset(r'G:\2024codes\Codes\Copulas\Copulacodestest\SPEI.nc') as f1:
    pre = f1['spei_gamma_01'][:-1, :, :]  # 提取三维数据 (时间，纬度，经度)
    lat, lon = f1['lat'], f1['lon']  # 提取经纬度

# 将三维数据转换为二维数据，便于计算相关系数
pre2d = pre.values.reshape(pre.shape[0], pre.shape[1] * pre.shape[2])

# 读取另一个 .nc 文件的数据
with xr.open_dataset(r'G:\2024codes\Codes\Copulas\Copulacodestest\SRI.nc') as f2:
    pc = f2['spi_gamma_01'][0, :].values.flatten()

# 计算相关系数矩阵
pre_cor = np.corrcoef(pre2d.T, pc)[:-1, -1].reshape(len(lat), len(lon))

# 做显著性检验
pre_cor_sig = f_regression(np.nan_to_num(pre2d), pc)[1].reshape(len(lat), len(lon))  # 用 0 代替 NaN
significant_area = np.where(pre_cor_sig < 0.05)  # 提取显著区域

# 创建格网化的经纬度网格
nx, ny = np.meshgrid(lon, lat)

# 创建绘图画布
plt.figure(figsize=(16, 8))
plt.rcParams['font.family'] = 'Times New Roman'
plt.rcParams['font.size'] = 16

# 设置投影为 PlateCarree 并将中心经度设置为 180 度
ax2 = plt.subplot(projection=ccrs.PlateCarree(central_longitude=180))
ax2.coastlines(lw=0.4)  # 添加海岸线
ax2.set_global()  # 设置显示全球范围

# 绘制相关系数的等值线图
c2 = ax2.contourf(nx, ny, pre_cor, extend='both', cmap=cmaps.nrl_sirkes, transform=ccrs.PlateCarree())
plt.colorbar(c2, fraction=0.05, orientation='horizontal', shrink=0.4, pad=0.06)  # 添加颜色条

# 在显著性区域打点
ax2.scatter(nx[significant_area], ny[significant_area], marker='+', s=1, c='k', alpha=0.6, transform=ccrs.PlateCarree())

# 设置标题和字体
plt.title('Correlation Analysis', fontdict={'family': 'Times New Roman', 'size': 16})

# 设置经度和纬度的刻度
ax2.set_xticks(np.arange(0, 361, 30), crs=ccrs.PlateCarree())
ax2.set_yticks(np.arange(-90, 90, 15), crs=ccrs.PlateCarree())
plt.xticks(fontproperties='Times New Roman', size=16)
plt.yticks(fontproperties='Times New Roman', size=16)

# 设置经纬度格式
ax2.xaxis.set_major_formatter(LongitudeFormatter(zero_direction_label=False))
ax2.yaxis.set_major_formatter(LatitudeFormatter())

# 设置地图显示的空间范围
ax2.set_extent([-180, 180, -90, 90], crs=ccrs.PlateCarree())

# 添加网格线
ax2.gridlines(draw_labels=True, color='k', linestyle='--')

# 显示图像
plt.show()