科研实战 | 基于CMIP6温度空间趋势图绘制并叠加显著性检验之方法一

1、前言

今天介绍的是基于CMIP6数据，绘制温度空间趋势图，并叠加显著性检验。文末附有源代码和nc文件下载路径，感兴趣的筒子们可以试试。

空间趋势图是研究气象气候领域最常用的研究方法之一，可以从空间上看出不同区域的空间分布差异，而且通过显著性检验可以看出不同区域，那些区域变化比较显著。

2、数据准备

本次测试用的数据是CMIP6格式下的全球温度数据，文件名称为：remapbiled_tas_yearly_EC-Earth3_ssp126_r1i1p1f1_gr_201501-210012.nc。首先解释一下这么文件名的含义（为什么要解释文件名的含义呢？因为从各个不同机构下载的数据，从文件名中，一般我们就能看出其nc文件内容，这是为了良好的可读性，在科研中，当文件非常多的时候，你会发现这么命名文件好处是相当明显的）：remapbiled表示该文件经过了cdo的remapbil命令，比如插值或者掩码处理；tas表示其nc文件中保存气象要素的变量名称，比如当前文件中存储温度数据的变量即为tas；yearly表示该tas变量为年平均数据，当然，你也可以命名为yearmean；EC表示是nc文件的提供机构，这里是nc，还有诸如ERA5或者JRA55等等；EC-Earth3表示是CMIP6中的EC-Earth3这个模式；ssp126表示共享路径下126路径；r1i1p1f1表示模式设置的一些参数；201501-210012表示是时间范围，为2015年1月到2100年12月；nc表示其文件格式，还有其他类型的数据格式，比如grid等等。

如下就是本次测试用数据的内容了。

3、数据提取

path1 = r"remapbiled_tas_yearly_EC-Earth3_ssp126_r1i1p1f1_gr_201501-210012.nc"
cmip6 = xarray.open_dataset(path1).tas
print(cmip6)

读取完了，记得打印出来看一下，多print是个非常好的习惯，这是非常重要的调试手段，数据长啥样，跟预期的是否一致，打印出来一看便知。

这里稍微解释一下time、lon、lat三个参数（为什么要解释呢？因为本人踩过坑，这三个参数的含义不用解释了，我想大家应该都能明白）：首先是time，这个参数呢，不同的文件，其格式可能不一样，但是有一点一定要搞清楚，这是字符串格式，还是整数格式的时间；lon和lat的位置也一定要弄清楚，经过cdo处理以后，这里打印出来的顺序跟nc文件中存储的顺序可能不一致，实际的nc可能是lat和lon的顺序，另外还要注意lat和lon的排序，有些可能是从小到大，有些则可能从大到小排列。以上注意点在loc定位选取数据的时候十分重要，一定不能搞错，因为loc按照字符串去数值和整数取值是不一样的，字符串是需要加引号，而整数不需要。

数据已经取完了，接下来就可以计算趋势了。

4、趋势和斜率计算

trend = np.zeros((cmip6.lat.shape[0],cmip6.lon.shape[0]))
p_value = np.zeros((cmip6.lat.shape[0],cmip6.lon.shape[0]))

for i in range (0,cmip6.lat.shape[0]):
    for j in range (0,cmip6.lon.shape[0]):
        trend[i,j], intercept, r_value, p_value[i,j], std_err=stats.linregress(np.arange(2015,2101),cmip6[:,i,j])

我们画趋势图，主要就是trend（趋势）和p_value（斜率）值，由于这两个值是二位数组，故需要提前定义，而定义的二位数组长度就是cmip6数据的lat格点数乘以lon格点数，cmip6.lat.shape[0]其实就是lat的格点数180，cmip6.lon.shape[0]其实就是lon的格点数360，因此trend和p_value的数组长度是180x360。shape[0]的作用就是取一维数组lat和lon的数组长度。

接下来就是趋势计算，其实就是stats.linregress这个函数，第一个参数是时间，第二个参数是cmip6的数据，其含义就是针对地球上某个特定的点，对其从2015年到2100年计算其变化趋势，外面的两层循环就是遍历全球所有的点。计算完成以后，我们再把trend和p_value打印出来看看。

print(trend)
print(p_value)

上面是trend，下面是p_value，中间有三个省略号，是python工具作了省略，请不要在意（上下两个值区分的应该很明显吧，请仔细看看）。

5、绘制趋势图

绘制趋势图，其实就是绘制trend值对应的图，先上代码，然后再解释

fig = plt.figure(figsize=(6,6),dpi=60)#画布
ax = fig.add_axes([2,2,2,2],projection = ccrs.PlateCarree())#画层
plot = ax.contourf(cmip6.lon, cmip6.lat , trend, zorder=0,  extend = 'both', cmap=cmaps.BlRe, transform=ccrs.PlateCarree())

前两行是画trend数据的准备工作（就比如要画画一样，画板和白纸，你总要准备好吧），准备好了以后，直接就画画了，其实很简单，就是第三行代码，关键的函数就是contourf（这个函数作用就是画画，其参数，大家可以看官网的介绍），这里稍微解释一下个参数的意义：cmip6.lon和cmip6.lat是两个一维数组，分别是经纬度数据，trend就是上面计算出来的趋势数据了，由于我们在上面初始化trend的时候，它只是一个二维数组，不带经纬度数据，因此这三个参数加在一起就表示全球任意一个点上趋势值；zorder可以理解为你在画板上画的第一张图（由于计算机里0才是第一个数字，所以这里用0表示）；cmap表示图片颜色；transform表示地图投影方式，此处为标准平面投影。画出来的效果如下：

先不要管中间的白线，这个问题，后面提供解决办法。这就是trend图，不同的颜色就表示此处温度趋势值不一样。

注解：这里用到了cmaps库，需要通过conda进行安装，该库的作用就是丰富了可用的颜色系。这个库非常强大，后续的文章再详细讲。

6、绘制显著性检验图

画完了趋势图，我们再在此基础上标注显著性检验的结果，同样，我们先上代码：

c1b = ax.contourf(cmip6.lon, cmip6.lat,  p_value,[np.min(p_value),0.05,np.max(p_value)], hatches=['.', None],zorder=1,colors="none", transform=ccrs.PlateCarree())

其实也很简单，还是一行代码解决问题，下面稍微解释一下（前三个参数不解释了，往上看）：[np.min(p_value),0.05,np.max(p_value)] 表示95%的检验，它把p_value的值划分为两个区间，一部分是p_value中的最小值到0.05（即5%），另外一部分是0.05到p_value的最大值，配合hatches参数，就表示，p_value的最小值到5%的检验区间，在底图中用'.'标注出来（也可以用'*'或者其他符号），而另一区间则什么也不画；zorder可以理解为是在上一张画纸已经画完的情况下，我又拿了一张画纸画画，不能覆盖前面的zorder；colors表示当前画纸的底图颜色，如果你把这个参数改成'blue'，你就能知道其代表是那个颜色了（blue的颜色会把zorder=0的画纸全部覆盖了，你就看不到前面contourf画的画了），效果如下：

这个时候，发现趋势图的颜色太深了，打点效果不是很好，然后我们尝试把趋势图的颜色调淡一点，但是为了调整趋势图的显示值的范围，需要先把colorbar调出来

bar=fig.colorbar(plot,ax=ax,orientation="horizontal",shrink=0.8,aspect=30,pad=0.1)  #aspect代表长宽比

由于趋势值主要集中在中间部分，且暖色系居多，所以调整参数范围：

level = np.arange(-0.035,0.025,0.001)
midnorm = mcolors.TwoSlopeNorm(vmin=-0.035, vcenter=0, vmax=0.025)

然后再把这两个参数添加到画趋势图的函数中去，然后之前画趋势图的函数就变成下面这个样子

plot = ax.contourf(cmip6.lon, cmip6.lat , trend, zorder=0,  extend = 'both', 
                     cmap=cmaps.GMT_polar,
                     norm = midnorm,
                     levels = level,
                     transform=ccrs.PlateCarree()

最后，我们再重新打点，这样看上去，貌似好很多。当然你也可以接着调整到自己满意为止。

7、如何去除中间的白线

文章最后，介绍一下如何去除中间的白线。第一行和第二行是解决趋势图上的白条，第三行和第四行是去掉打点图中的白条。

cycle_slope, cycle_lon = add_cyclic_point(trend, coord=lon)
cycle_LON,   cycle_LAT = np.meshgrid(cycle_lon, lat)#解决中间白条
cycle_p, cycle_lon = add_cyclic_point(p_value, coord=lon)
cycle_LON,   cycle_LAT = np.meshgrid(cycle_lon, lat)#解决中间白条

然后趋势图和打点图都需要调整一下：

plot= ax.contourf(cycle_LON,   cycle_LAT , cycle_slope, zorder=0,  extend = 'both', 
                    cmap=cmaps.GMT_polar,
                    norm = midnorm,
                    levels = level, 
                    transform=ccrs.PlateCarree())

上面是趋势图的调整，而下面是打点图的调整

c1b = ax.contourf(cycle_LON,cycle_LAT, p,[np.min(p),0.05,np.max(p)], zorder=1,hatches=['.', None],colors="none", transform=ccrs.PlateCarree())

这是最终效果