登录后查看更多精彩内容~
您需要 登录 才可以下载或查看,没有帐号?立即注册
x
本帖最后由 不想去气象局 于 2020-7-3 16:29 编辑
写在前面:说实话昨天的文章被推到气象家园公众号了之后属实有些惊喜,今天下午服务器后台在处理CMIP6的数据,就突然闲下来不知道干嘛了,所以打算整理一下我这学期的上机课程的代码。就先从短期气候预测开始吧,下次看看气象统计方法有没有要整理的,如果有的话会和今天留下的一点小尾巴放在一起。 下面开始进入正题: 实习一就是简单的环流平均图,距平图,纬偏图。因为真的很简单,我就不多说了,但是需要提一句的是,Cartopy在画完整360°经度的图时,需要对数据多一步处理,否则会出现180°或者0°经线上的数据是缺省的情况,会出现一条白线。出过这种错误的人肯定知道我在说啥,不知道我在说什么的就留个印象就行,以后180°或者0°经线出问题的时候回来看看就行,关键代码贴在下面,所以记得关注我!代码应该是参考气象学家公众号的,但是哪一篇文章我忘了,反正不好找。所以我的文章主要记录那些零碎的点,大家如果要系统地学习建议还是多看看其他文章。这个函数的作用相当于就是把0和360连起来了吧。 代码放上来出错了,只好改用分割线了。。
from cartopy.util import add_cyclic_point cycle_hgt, cycle_lon =add_cyclic_point(hgt, coord=lon) mapp=f1_ax1.contourf(cycle_lon,lat,cycle_hgt,levels=np.arange(1100,1625,40), zorder=0, extend ='both',transform=ccrs.PlateCarree(), cmap=plt.cm.jet)
丑丑的图我也贴上了,毕竟只是交个作业而且报告还是黑白打印的哈哈哈。 实习二是对亚澳季风区环流的eof分解,先给大家贴两张图,不知道大家能不能发现两者的区别。同样是eof分解,究竟哪个是正确的,哪个是错误的呢? 先不管空间两张图空间分布的特征,我们可以第一时间发现的应该是两个色标的数量级都不一样,那肯定有问题了呀。上面这幅图是我根据摸鱼咯大佬的代码出的图,但是画完之后我发现跟另一个班老师的答案不一样啊,但是吧我不相信Eofs的库代码有问题,所以得找原因是吧,所以我去找Eofs库的官方文档,结果发现人家的代码应该就是参考官方的,那怎么办呢?还能怎么办,只好去看源代码了咯。代码看到一半发现问题了,原来是我们平时用的eof和官方例子的eof不一样!官方的例子是用的 - eof1 =
- solver.eofsAsCorrelation(neofs=1)
复制代码即用的是eofsAsCorrelation,从字面意思上看跟相关系数有关系,官方文档对这个函数的解释是这样的Empirical orthogonal functions (EOFs) expressed as the correlationbetween the principal component time series (PCs) and the time series of theEof input dataset at each grid point.也就是说他并不是我们eof空间分解后的V向量,而是某空间点的变化与时间Z向量的相关系数,不知道大家听懂没,没听懂可以回去看看气象统计方法的书。那好不容易找到一个包,难道我们又要重新自己写代码吗?当然不用啦!我们只要把eofsAsCorrelation换成eofs就行啦!我所呈现的这两个eof的关键代码区别就贴在下面啦! - eof = solver.eofsAsCorrelation (neofs=2)
复制代码
- eof = solver.eofs(neofs=2)
复制代码
也不知道是不是国外或者现在学术用相关系数来表示eof的空间场比较多,不然他官网为什么给这样的例子?但是我还是得把课程学的解决了先。 eofs包的安装只需要在我上一篇文章里讲到的Anaconda Prompt里 - conda install -c conda-forge eofs
复制代码或者 实习三四做的是遥相关,以及和我国气温的相关系数,我就简单跳过了,关键代码主要就是求相关系数,直接贴下面了。
from scipy.stats import pearsonr for i in range(37): for j in range(144): r[i,j],p[i,j]=pearsonr(WP,hgt[:,i,j])
如果用过别的绘图软件,然后转Cartopy其实最关心的就是中国地图的绘制,答案肯定是可以的,但是由于我自己的代码也是参考别人的,所以在此就不贴出来了,如果有需要中国九段线shp文件的或者完整脚本的可以在我的公众号后台私信我。最后贴个图 实习五六做的是合成分析,代码部分没啥好整理的。 实习七八是降水预测,做个简单的线性回归,也不讲了。 最后再整理一下字体部分吧,因为发文章一般有字体要求,下面我以Times New Roman为例来讲。从上面的图中大家可以看到,文字的标题以及坐标刻度部分我都改成了Times New Roman,但是色标的字体还是默认的字体。因为我不是直接改默认字体,所以还挺麻烦的,然后当时也没有找到怎么改色标上的字体,所以先讲另外两个吧。 改title是最简单的,直接如下即可 - f1_ax1.set_title('(d) CC bewteen WP and
- temperature JJA',loc='left',fontsize =16)
复制代码如果只有一幅图,那改坐标上刻度的数字也是挺简单的,如下即可 - plt.yticks(fontproperties = 'Times New
- Roman', size = 16)
复制代码但是前段时间我在画有9幅子图的时候,直接用plt就出问题了。怎么改都不管用,九遍上述代码也只有一幅图是改了字体的,所以还得指定具体的ax。但是接下来狗血的就来了,改坐标的字体可以用set_xticklabel方法来改,我们先看一下用法
Parameters ---------- labels : List[str] List of string labels. fontdict : dict, optional Adictionary controlling the appearance of the ticklabels. The default `fontdict` is:: {'fontsize': rcParams['axes.titlesize'], 'fontweight': rcParams['axes.titleweight'], 'verticalalignment': 'baseline', 'horizontalalignment': loc}
这说的是set_xticklabel至少有一个输入是label,也就是说你如果要改fontdict,必须指定label。这就让事情变复杂了,我们需要这样操作:
ax5.set_yticks(np.arange(lowerlat+0.25,upperlat+1,2),crs=ccrs.PlateCarree()) ylabels=ax5.get_yticks() ax5.set_yticklabels(labels=ylabels,fontdict={'family': 'Times New Roman', 'size' : 16})
也就是说,先用set_yticks设置刻度的取法,而set_yticks又不能改字体,我们需要用get_yticks()读取刻度,然后最后用set_yticklabels改字体。确实有点麻烦,但是复制粘贴还是快的哈哈。 至于最后的色标怎么改字体放到下次再说吧,今天先休息了。 没有隐藏内容你们就不回复,那我就只好发金币啦。大家记得关注我的公众号哟!公众号主要是整理一些东西,给自己备忘吧。
公众号
|