登录后查看更多精彩内容~
您需要 登录 才可以下载或查看，没有帐号？立即注册
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改字体。确实有点麻烦，但是复制粘贴还是快的哈哈。 至于最后的色标怎么改字体放到下次再说吧，今天先休息了。 没有隐藏内容你们就不回复，那我就只好发金币啦。大家记得关注我的公众号哟！公众号主要是整理一些东西，给自己备忘吧。
公众号
| 
[复制链接]
窥视卡
雷达卡
发表于 2020-7-3 16:26:25
提升卡
置顶卡
沉默卡
喧嚣卡
变色卡
抢沙发
千斤顶
显身卡
楼主
