【附代码】时间序列、空间序列的相关、显著性检验打点
在气象科研与业务经常使用的相关有：时间序列与时间序列的相关、时间序列与空间场的相关、空间场与空间场的相关。其中最常使用的就是皮尔逊相关系数。
什么是皮尔逊相关系数
该相关系数是由卡尔·皮尔逊在前人的研究基础上所提出的相关统计量，可以用来度量两个变量之间的简单线性关系。它的计算公式如下：
通过该公式计算得到的相关系数r，取值范围为[-1,1]。
• 当r=0时，表明两个变量X和Y之间无线性关系(注意，r=0并不代表X和Y一定相互独立，可能存在非线性等其他关系，具体可以自行带入 进行体会)；
• 当0<r<1时，表明两个变量X和Y之间存在正相关关系，即当X的值增大(减小)时，Y的值也增大(减小)；
• 当-1<r<0时，表明两个变量X和Y之间存在负相关关系，即当X的值增大(减小)时，Y的值减小(增大)。
由于Pearson相关系数的形式简单，因此其应用十分广泛。但其也存在缺点，即：
• 该相关系数只能识别简单的线性相关关系，无法处理非线性相关关系；
• 对异常值(或离群点)和样本容量较为敏感；
• 要求研究的变量是数值变量，且变量符合或较为接近正态分布。
气象实例
时间序列与时间序列的相关系数计算：
#导入库
import xarray as xr #读取、处理nc数据的包
import numpy as np #进行数学处理的包
from scipy.stats import pearsonr #进行Pearson相关系数计算的包
from scipy.stats import normaltest #检验数据是否符合正态分布的包
import cartopy.crs as ccrs #进行画图的cartopy和matplotlib包
import cartopy.feature as cfeature
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
import matplotlib.pyplot as plt
#原始的T2和RAIN气象变量是具有时间一维、空间二维的三维变量，为了将其变为仅有时间维度的一维时间序列，我们分别对这两个变量用 mean() 方法沿着 south_north 和 south_north 两个空间维度求平均，并赋值给新变量 T2_series, RAIN_series
T2_series=data['T2'].mean('south_north').mean('west_east')
RAIN_series=data['RAIN'].mean('south_north').mean('west_east')
#用 normaltest() 方法检验两个气象变量是否符合或接近正态分布(当pvalue>0.05时，即表示数据接近正态分布)，判断这两个变量是否适合使用Pearson相关系数
normaltest(T2_series)
normaltest(RAIN_series)
# 由于Pearson相关系数对于异常值和缺省值比较敏感，所以一般需要用 np.isnan 来检测数据是否存在缺省值(存在为True，不存在为False)，并通过绘制散点图等方式观察是否存在显著的离群点。如果存在需要剔除，比如可以用均值进行替代或直接去掉。
True in np.isnan(T2_series) #即不存在异常值nan值
True in np.isnan(RAIN_series) #即不存在异常值nan值
plt.scatter(T2_series,RAIN_series)#未观察到显著的离群点
#对两个时间序列 T2_series 和 RAIN_series 使用 pearsonr() 方法计算相关系数，返回值为 (r, p)。其中 r 为相关系数，p 为显著性检验结果(小于0.05即为显著)
r,p=pearsonr(T2_series,RAIN_series)
print('r=',np.round(r,3),'p=',np.round(p,3))#np.round(x,3)表示将x保留3位小数
空间场与空间场的相关系数计算：
计算场与场之间相关系数的思路是：将场中的每一个格点都看作为一条时间序列，对两个场的对应格点分别做序列与序列的相关，再将计算结果赋给该格点即可。这样得到的是一个相关场(2维的)。
如果想得到一个相关序列，则可以将时间作为循环，将每一个时刻的两个空间场reshape成一个1维的空间序列，再对这两个序列做相关性计算。
p.s. 可以看到，计算场与场之间相关系数最终还是回到了对序列和序列的关系进行处理。
相关场(空间2D)：
#定义两个空数组 r2 和 p2，并将数组的大小设置为 (south_north, west_east)，r2 和 p2 会用来存放每个格点对应的 r 值(Pearson相关系数) 及 p 值(显著性检验结果)
r2=np.nan*np.zeros((len(data.south_north),len(data.west_east)))
p2=np.nan*np.zeros((len(data.south_north),len(data.west_east)))
#通过 i * j 对二维空间场中的每个格点进行循环，分别计算 T2 和 RAIN 两个场对应格点上时间序列的相关系数，并存储在 r2 和 p2 中
for i in range(len(data.south_north)):
for j in range(len(data.west_east)):
r2,p2=pearsonr(data['T2'][:,i,j],data['RAIN'][:,i,j])
相关序列(时间1D)：
#定义两个空数组 r2 和 p2
r2=np.nan*np.zeros((len(data.Time)))
p2=np.nan*np.zeros((len(data.Time)))
#通过 i * j 对二维空间场中的每个格点进行循环，分别计算 T2 和 RAIN 两个场对应格点上时间序列的相关系数，并存储在 r2 和 p2 中
for i in range(len(data.Time)):
r2<i>,p2<i>=pearsonr(np.reshape(data['T2'],len(data.south_north)*len(data.west_east)),np.reshape(data['RAIN'],len(data.south_north)*len(data.west_east)))</i></i>
显著性检验打点
在完成场与场之间相关的计算以后，我们就可以根据显著性检验的结果了——使用Matplotlib中的 scatter() 方法绘制打点图。
打点图可以呈现出：哪些区域的相关性是通过显著性检验的，而哪些区域是没有通过显著性检验的。
#先根据显著性检验结果 p2，使用Numpy中的 where() 方法选择出 p值小于0.05的格点 所对应的索引。
p2 = np.where(p2<0.05)
#使用 Cartopy 绘制地图底图，并用 Matplotlib 中的 scatter() 方法将显著检验结果p2可视化。
fig = plt.figure(figsize=(4, 4), dpi=200) #设置画布
ax = plt.axes(projection=ccrs.PlateCarree(central_longitude=180)) #设置投影方式
ax.add_feature(cfeature.LAND) #添加陆地
ax.add_feature(cfeature.OCEAN) #添加海洋
lon=np.linspace(-179.5,179.5,360)
lat=np.linspace(-67.5,76.5,145)
mesh_lon,mesh_lat=np.meshgrid(lon,lat) #生成经纬度网格
ax.scatter(mesh_lon,mesh_lat,color='k',s=0.2,linewidths=0.1,transform=ccrs.PlateCarree()) #绘制散点图(即打点)
ax.set_extent(,crs = ccrs.PlateCarree()) #为了方便大家看的更清楚，我们限制显示的区域为70°E-140°E，纬度为0°-55°N
时间序列与空间场的相关系数计算
要想计算计算温度时间序列数据 T2_series 与降水场数据 RAIN 的相关系数，就是将降水场 RAIN 中的每个格点看作为一条时间序列，计算每个格点的降水时间序列与温度时间序列 T2_series 之间的相关系数。
import xarray as xr #读取、处理nc数据的包
import numpy as np #进行数学处理的包
from scipy.stats import pearsonr #进行Pearson相关系数计算的包
from scipy.stats import normaltest #检验数据是否符合正态分布的包
import cartopy.crs as ccrs #进行画图的cartopy和matplotlib包
import cartopy.feature as cfeature
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
import matplotlib.pyplot as plt
#计算相关系数
r3=np.nan*np.zeros((len(data.south_north),len(data.west_east)))
p3=np.nan*np.zeros((len(data.south_north),len(data.west_east)))
for i in range(len(data.south_north)):
for j in range(len(data.west_east)):
r3,p3=pearsonr(T2_series,data['RAIN'][:,i,j])
#画图
fig = plt.figure(figsize=(6, 4), dpi=200)
ax = plt.axes(projection=ccrs.PlateCarree()) #设置投影方式
ax.coastlines(lw=0.3) #设置海岸线
ax.set_xticks([ -120,-60,0,60,120], crs=ccrs.PlateCarree()) #设置x轴坐标刻度
ax.set_yticks([ -60, -30, 0, 30, 60], crs=ccrs.PlateCarree())#设置y轴坐标刻度
lon_formatter = LongitudeFormatter(zero_direction_label=False)#给经度加上E和W，0°不标
lat_formatter = LatitudeFormatter() #给经度加上S和N，0°不标
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)
c1= ax.contourf(data.XLONG,data.XLAT,r3,cmap='rainbow',transform=ccrs.PlateCarree())
plt.colorbar(c1,shrink=0.45)
<div></div>
文章来源于微信公众号：自学气象人
有没有更高效的实现方式？我在寻找优化方法。
页:
[1]