Python气象数据处理与绘图：相关性分析之散点图
1、前言
应粉丝要求更新一篇散点图相关分析的文章，这个图是否叫这个名字我也不太确定，考虑到这种图的画法大部分是使用散点的形式进行展示，那本文叫相关分析图吧。
很多文章中会提到2个变量之间有什么潜在联系的时候经常用的一张图就是散点图，比如说某地降水和温度之间是否存在联系，此时建立笛卡尔坐标系，横轴纵轴分别为降水和温度在同一位置随时间变化的量，这个时候在图中如果能看出来两者存在一些线性关系或者非线性的话，这个图就有价值了。说了这么多，我找几张图给大家看一下。
--降水和温度与动力之间的关系--
该图来自论文：Chen Z, Zhou T, Zhang L, et al. Global land monsoon precipitation changes in CMIP6 projections[J]. Geophysical Research Letters, 2020, 47(14): e2019GL086902.
图中还可以叠加相关系数的显著性检验，如果有多条模式或者不同时段进行对比的话可以在图中设置不同颜色的字体，比如上图使用的是不同时段的SSPs路径下CMIP6 MME结果来进行对比，目前学术上认为2021-2040为近期，2041-2060为中期，2080-2100为长期。这样的目的是为了看在未来不同时段下，CMIP6模式的结果是否呈现不同的特征，可以有一个横向的对比。
废话不多说了，开始吧~
2、读取数据
第一步，读入数据，代码如下：
data=xarray.open_dataset('/mnt/d/science/large_ensemble/picture/Fig15.tas_gardient_sfc/scatter_tas_sfcwind_add_members.nc')
tas_gradient_NH=data.tas_gradient_nh.values
tas_gradient_SH=data.tas_gradient_sh.values
sfc_NH= (data.sfc_nh.values)
sfc_SH= (data.sfc_sh.values) #NH 的梯度需要反转
sfc_NH_full= (data.sfc_nh_full.values)
sfc_SH_full= (data.sfc_sh_full.values) #NH 的梯度需要反转
tas_gradient_NH_full= (data.tas_gradient_nh_full[:,1995-1950:].values)
tas_gradient_SH_full= (data.tas_gradient_sh_full[:,1995-1950:].values) 这一步的话，不用太深究每个数据是什么意思。我在本篇文章中使用的数据是南北半球某区域的变量A和温度梯度在1950年到1995年的时间序列。也就是说A和温度对应的长度都是46，必须保证长度一致才能画图。这里可以这么理解: 平时时间序列图大家都画过了，x是Time，y是画的变量。这时候x和y的长度一致才能画图，散点图的话就是把x改为另一个长度相同的变量而已。
3、封装画图函数
封装函数是程序走向简洁性、易用性的重大一步，同时也是批量出图的一个关键点。后期会专门开一篇讲一下封装函数的一些注意点。
首先我们定义一个函数用来绘制图中的拟合一元线性方程：
def f(x,intercept,trend):# 计算出拟合方程
y=trend*x+intercept
return y然后我们定义散点图的函数：
函数的步骤分别是计算散点图的线性趋势，通过linregress函数返回的数据计算拟合方程。然后通过personr函数计算出相关系数和进行显著性检验。最后把显著性检验的结果放到legend上面进行展示。
def draw_scatter(ax,x,y,label,color,linecolor,loc,bbox_to_anchor,):
trend=np.zeros(2)
intercept=np.zeros(2)
p_value=np.zeros(2)
trend[0], intercept[0], rsq1, p_value[0], std_err=stats.linregress(x,y)
print("each 0.5 ℃ change of LTG over could cause a XX±XX m s-1=",0.5*trend[0],'sterr=',std_err)
# 绘出图像
CC=pearsonr(x, y)
p=[]
#if CC[1]<0.001:
# p="P<0.001"
if CC[1]<0.1:
p=r"$\it{P<0.01}$"
if 0.1<=CC[1]<0.5:
p=r"$\it{P<0.5}$"
if CC[1]>=0.5:
p=" "
print("p=",p)
ax.scatter(x,y,c = color,marker = '.', s=100)
ax.plot(x, f(x,intercept[0],trend[0]), color=linecolor, linewidth = 3,
label=("$R$= %2.2f" %(CC[0])+", " +p )
)
ax.grid()
leg=ax.legend(loc=loc,prop={'size':20},frameon=False,
handlelength=0,
bbox_to_anchor=bbox_to_anchor
)
ax.tick_params(labelsize=17)
ax.minorticks_on()
ax.tick_params(which='major', width=1.5)
ax.tick_params(which='major', width=1.0)
ax.tick_params(which='major', length=9)
ax.tick_params(which='minor', length=6)
ax.set_ylabel('Precipitation anomaly ', loc='center' , fontsize =20,labelpad=12)
ax.yaxis.set_major_formatter(FormatStrFormatter('%.2f'))# 保留2位小数点
return ax
这个图传递的参数比较多，ax是创建的图层，x和y是散点图的x和y轴对应的数据，后面2个color是散点和线的颜色，loc是legend的位置，bbox是自定义位置相对于loc而已。有同学会说怎么这么复杂呢。其实函数复杂一点的话对很多图自定义来说反而简单了，这样可以增加后期修改图时的麻烦，也可以高度自定义化，形成自己的画图风格。我建议大家对函数要真正理解清楚每一句话的意义，这样才能达到随心所欲的境界。
4、画图
画图的代码就是添加fig和ax，不使用subplot，这一步可以参考：
好了，代码奉上：
fig_31 =plt.figure(figsize=(6, 5),dpi=60)
ax0 =fig_31.add_axes([1.4, 6.3, 0.9, 0.9])
ax1 =fig_31.add_axes([2.6, 6.3, 0.9, 0.9])
bbox_to_anchor1=(0.8,0.95)
bbox_to_anchor2=(0.2,0.95)
draw_scatter(ax0,tas_gradient_NH,sfc_NH," ",'orange',"black",
'upper right',(1.01,1))
draw_scatter(ax1,tas_gradient_SH,sfc_SH," ", 'lime',"black",
'upper left',(0.00,1))
ax0.invert_xaxis()
ax0.set_title('Northern Hemisphere ',loc='right',fontsize =22)
ax1.set_title('Southern Hemisphere ',loc='right',fontsize =22)
ax0.set_title('(c)',loc='left',fontsize =24)
ax1.set_title('(d)',loc='left',fontsize =24)
ax1.set_xlabel('Air temperature gradient increase :$^\circ$C ', loc='center' , fontsize =20,labelpad=12)
ax0.set_xlabel('Air temperature gradient decrease:$^\circ$C ', loc='center' , fontsize =20,labelpad=12)
出图效果如下：
图片画出来后就是看图说话，比如说a图中，我们可以看到北半球的变量A和温度梯度之间呈现明显的正相关，由于温度梯度在南北半球的方向相反，所以我把a图的温度梯度进行了反转，这样更加有利于理解。其他的同学要具体问题具体分析，这里要学会修改代码。a图中相关系数为0.88，P<0.01，这个结果显示出很强的相关性，并且通过线性回归方程的表现来看，基本呈现线性相关。从b图中也可以看到差不多的结果，但是b图中的相关系数更高，所以从图的整体可以得到结论：A变量与温度梯度的相关性很高，其中南半球A与温度梯度变化更为密切相关。
腾讯云开发者
