做相关系数的显著性检验

Bone

计算1979-2019年的EU遥相关指数与同期环流场（500hPa高度场或海平面气压场）的相关系数，并做相关系数的显著性检验

Bone

ds=xr.open_dataset('D:\dqqh\dq2022\dq2022\shixi1\hgt.mon.mean.nc')
hgt=ds.hgt
print(hgt.shape)
hgt_Jan_year=hgt.groupby('time.month')[1]
print(hgt_Jan_year.shape)
hgt_Jan=hgt_Jan_year[31:72,5,:,:]
print(hgt_Jan.shape)

lat=ds.variables['lat'][:]
lon=ds.variables['lon'][:]

##################################EU指数#####################################
hgt_1=hgt_Jan.loc[:,55,20]
hgt_2=hgt_Jan.loc[:,55,75]
hgt_3=hgt_Jan.loc[:,40,145]

I=hgt_2*(0.5)-hgt_1*(0.25)-hgt_3*(0.25)
print(I)

I_mean=np.mean(I) #标准化处理
I_sd=math.sqrt(sum((I-I_mean)**2)/41.0)
I_std=(I-I_mean)/I_sd
print(I_std.shape)

fig=plt.figure(figsize=(10,8))
ax1=fig.add_subplot()

ax1.plot(I_std,marker='o',linewidth=1.5,ms=5)
ax1.set_title('Std EU Index')
ax1.set_xlabel('Years')
plt.xticks(np.arange(0, 41, 5), np.arange(1979,2020,5))  #或者ax.set_xticks(np.arange(0, 41,5))//ax.set_xticklabels(np.arange(1979,2020,5))
plt.savefig(r'D:\dqqh\dq2022\dq2022\shixi3-4\EU-Index.jpg',dpi=300)
plt.show()
plt.close()

##################################与500hPa位势高度的相关系数############################
r,p = np.zeros((73,144)),np.zeros((73,144))
for i in range(73):
    for j in range(144):
        r[i,j],p[i,j]=pearsonr(hgt_Jan[:,i,j],I_std)
print(r.shape)

fig=plt.figure(figsize=(10,10))
ax1  = fig.add_subplot(111)
cmap = plt.get_cmap('RdBu_r')  #选择coloarmap

ax1 = plt.subplot(111, projection=ccrs.PlateCarree(central_longitude=180))
ax1.set_global() #使得轴域（Axes即两条坐标轴围城的区域）适应地图的大小
ax1.coastlines() #画出海岸线

ax1.set_xticks([0, 60, 120, 180, 240, 300, 360], crs=ccrs.PlateCarree())
ax1.set_yticks([-90, -60, -30, 0, 30, 60, 90], crs=ccrs.PlateCarree())
#zero_direction_label用来设置经度的0度加不加E和W
lon_formatter = LongitudeFormatter(zero_direction_label=False)
lat_formatter = LatitudeFormatter()
ax1.xaxis.set_major_formatter(lon_formatter)
ax1.yaxis.set_major_formatter(lat_formatter)

plt.contourf(lon,lat,r,cmap='RdBu_r',levels=np.arange(-1,1.1,0.1),transform=ccrs.PlateCarree(),zorder=0)
plt.colorbar(shrink=0.87,aspect=20)
#ax1.set_title(f'2019-500hPa Jan_anomaly gpm')
plt.savefig(r'D:\dqqh\dq2022\dq2022\shixi3-4\EU_hgt_r_plate.jpg',dpi=300)
plt.show()

#相关系数的t检验

代码如上

GXYYYYY

看样子同是信带大三学子

Bone

GXYYYYY 发表于 2022-5-21 09:01
看样子同是信带大三学子

哈哈，是的！

希诺漾

请问楼主如果缺测数量偏多的话该怎么处理呢

wei来守虎者

GXYYYYY 发表于 2022-5-21 09:01
看样子同是信带大三学子

我们老师是zy你们呢。。不会还是一个班的吧。。

GXYYYYY

wei来守虎者 发表于 2022-6-13 16:14
我们老师是zy你们呢。。不会还是一个班的吧。。

看来不是哈哈哈哈哈 我们短期是华老师