Python vs NCL （之前数据有误，已更换）

气海无涯 2021-09-27

1687

1、前言

NCL语言大概是从2018年开始，陪伴我走过了一年的时间，也帮助我发表了自己的第一篇文章。后来由于自己需要实现的算法过于复杂，自己无力编写，在朋友的诱惑下不得不转投Python，靠网上各路大神的帖子东拼西凑满足自己的科研需求。说卸磨杀驴也有点不合适，但在适应了Python之后，我确实基本没有再打开过NCL了。最近整理文件，找到了自己发表第一篇学术垃圾时候的NCL脚本，心血来潮的想再用Python实现一遍。也算是对评论里很多读者的一个回应吧(貌似气象家园帖子下边第一个评论就是希望我写一个两者对比的文章，被我鸽到现在，咕咕咕)。注：编程水平仅限于能跑就行，warning不管，冗杂语句请忽视。

2、示例对比

下面通过示例对比Python和NCL，感受两种语言的区别。

示例1（EOF）

选择EOF的原因是图片内容较为丰富，同时方法较为常用

由于原始数据过大，只提供处理好的数据方便测试(数据为每年的寒潮路径的概率密度分布，为39年×90纬度×360经度，由FLEXPART模式追踪后通过高斯核密度估计算法Gaussian KDE得到。

测试数据和代码如下（请回复关键字“NCL”获取）：

1.测试数据 2. NCL代码 3. Python代码

先对比下出图效果(两种语言对图形的渲染有差异，EOF算法可能也有些差异，但是结果是基本相同的，图只经过了基础的调整，并不是很好看)。

1 准备工作

NCL：在引入一些特殊函数(NCL本体不具备的函数)时，通常会加上类似于以下语句：

1load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"

Python：引入各个模块(库)：

 1import xarray as xr
 2import numpy as np
 3from eofs.standard import Eof
 4import cartopy.crs as ccrs
 5import cartopy.feature as cfeature
 6import matplotlib.pyplot as plt
 7import cartopy.mpl.ticker as cticker
 8def contour_map(fig,img_extent,spec):
 9    fig.set_extent(img_extent, crs=ccrs.PlateCarree())
10    fig.add_feature(cfeature.COASTLINE.with_scale('50m')) 
11    fig.add_feature(cfeature.LAKES, alpha=0.5)
12    fig.set_xticks(np.arange(leftlon,rightlon+spec,spec), crs=ccrs.PlateCarree())
13    fig.set_yticks(np.arange(lowerlat,upperlat+spec,spec), crs=ccrs.PlateCarree())
14    lon_formatter = cticker.LongitudeFormatter()
15    lat_formatter = cticker.LatitudeFormatter()
16    fig.xaxis.set_major_formatter(lon_formatter)
17    fig.yaxis.set_major_formatter(lat_formatter)

2 数据读取

NCL：

1f=addfile("out.nc","r");读取文件
2dot=f->cspath(:,:,{0:180});读入变量
3dot:=dot(lat|:,lon|:,year|:);调整变量维度顺序(EOF函数对维度顺序有要求)
4x=ispan(1979,2017,1);时间

Python：利用Xarray库读取NC文件

1f = xr.open_dataset("out.nc" )
2dot = np.array(f['cspath'].loc[:,0:90,:])
3lat = f['lat'].loc[0:90]
4lon = f['lon']
5years = range(1979, 2018)

3 EOF

NCL：

1eof=eofunc_Wrap(dot,3,False)
2pc=eofunc_ts_Wrap(dot,eof,False)
3pc=dim_standardize_n(eof_ts,1,1)
4copy_VarMeta(dot(:,:,0),eof1);将维度信息重新赋值给新数组
5copy_VarMeta(dot(:,:,0),eof2)

Python：利用EOF库

1solver = Eof(dot)
2eof = solver.eofsAsCorrelation(neofs=2)
3pc = solver.pcs(npcs=2, pcscaling=1)
4var = solver.varianceFraction()

4.1 绘图：建立画布

NCL的有些绘图参数并不是必要的，由于年代久远，我记不起每条语句对应的详细功能了。

NCL：

1wks=gsn_open_wks("pdf","tmp")
2res                            = True            
3res@gsnMaximize                = False
4res@gsnDraw                    = False
5res@gsnFrame                   = False
6res@gsnDraw=False
7res@gsnFrame=False
8respc=res;实际上是设置PC图形的基础绘图参数

Python：

1fig = plt.figure(figsize=(15,15))

4.2 绘图：绘制EOF

NCL：其绘图思路是每一条语句指定一个绘图效果，通过不断的”叠BUFF“实现全部要素的叠加，先将共同要素叠给res，然后通过res1=res, res2=res赋值后，再对res1和res2分别添加各自的属性。

 1res@mpMaxLatF=90;设置经纬度边界
 2res@mpMinLatF=0
 3res@mpMaxLonF=160
 4res@mpMinLonF=0
 5res@mpFillOn =False;地图设置
 6res@mpOutlineOn = True
 7res@mpDataBaseVersion         = "MediumRes"             
 8res@mpDataSetName              = "/data/home/Earth..4"
 9res@cnFillPalette               = "MPL_bwr"
10res@cnFillOn                   = True
11res@cnFillDrawOrder            = "PreDraw"
12res@cnLinesOn                  = False
13res@cnLineLabelsOn             = False
14res@gsnLeftString=""
15res@pmTickMarkDisplayMode      = "Always"  
16res@cnLevelSelectionMode="ExplicitLevels"
17res@cnLevels            =(/-0.05,-0.04,-0.03,-0.02,-0.01,0,0.01,0.02,0.03,0.04,0.05/)
18res1=res
19res2=res
20res1@gsnRightString              = "41.84%"
21res1@gsnLeftString              = "EOF1"
22res1@lbLabelBarOn=False 
23res2@gsnRightString              = "14.59%"
24res2@gsnLeftString              = "EOF2"
25res2@lbLabelBarOn=True
26map1 = gsn_csm_contour_map(wks,eof1,res1) ;绘图
27map2 = gsn_csm_contour_map(wks,eof2,res2)

Python：与NCL相似的地方在于需要对每个axes添加参数，不同的地方在于其一条指令可以包含多个效果(比如将所有设置填色绘图的参数全部加在f_ax1.contourf里)。

 1proj = ccrs.PlateCarree(central_longitude=80) #投影
 2leftlon, rightlon, lowerlat, upperlat = (0,160,10,90) #地图边界
 3img_extent = [leftlon, rightlon, lowerlat, upperlat]
 4f_ax1 = fig.add_axes([0.1, 0.32, 0.3, 0.4],projection = proj)#EOF1
 5contour_map(f_ax1,img_extent,20) #利用前边自定义的绘图函数，目的是省去绘制相同图形时需要再次设置地形，湖泊，轴刻度等参数
 6f_ax1.set_title('(a) EOF1',loc='left')#左标题
 7f_ax1.set_title( '%.2f%%' % (var[0]*100),loc='right')#右标题，解释方差
 8f_ax1.contourf(lon,lat, eof[0,:,:], levels=np.arange(-0.9,1.0,0.1), zorder=0, extend='both',transform=ccrs.PlateCarree(), cmap=plt.cm.RdBu_r)#绘制填色
 9f_ax2 = fig.add_axes([0.1, 0.1, 0.3, 0.4],projection = proj)#EOF2
10contour_map(f_ax2,img_extent,20)
11f_ax2.set_title('(c) EOf',loc='left')
12f_ax2.set_title( '%.2f%%' % (var[1]*100),loc='right')
13c2=f_ax2.contourf(lon,lat, eof[1,:,:], levels=np.arange(-0.9,1.0,0.1), zorder=0, extend='both', transform=ccrs.PlateCarree(), cmap=plt.cm.RdBu_r)
14#绘制色标
15position=fig.add_axes([0.1, 0.17, 0.3, 0.017])
16fig.colorbar(c2,cax=position,orientation='horizontal',format='%.1f',)

4.3 绘图：绘制PC

NCL：

 1respc@gsnYRefLine=0
 2respc@trYMinF=-3
 3respc@trYMaxF=3
 4res9 = respc
 5respc@tmYLLabelDeltaF=-1
 6respc@trXMinF=1979
 7respc@trXMaxF=2017
 8respc@tiYAxisOn=False
 9respc@tmXTOn =False
10respc@tmYROn = False
11respc@vpHeightF=0.39
12respc@vpWidthF=0.6
13respc@gsnLeftStringOrthogonalPosF =  0.04
14respc1 = respc
15respc1@gsnLeftString              = "PC1"
16respc2 = respc
17respc2@gsnLeftString              = "PC2"
18pc1= gsn_csm_xy(wks,x,eoft1,respc1)
19pc2= gsn_csm_xy(wks,x,eoft2,respc2)
20pc1_9 = runave_n_Wrap(eoft1, 9, 0, 0);9年滑动平均
21pc2_9 = runave_n_Wrap(eoft2, 9, 0, 0)
22res9@xyLineColor = "red"
23res9@xyLineThicknesses = 4
24plotpc9_1 = gsn_csm_xy(wks, x, pc1_9, res9)
25plotpc9_2 = gsn_csm_xy(wks, x, pc2_9, res9)
26overlay(pc1, plotpc9_1);叠加
27overlay(pc2, plotpc9_2)

Python：

 1f_ax3 = fig.add_axes([0.45, 0.44, 0.3, 0.156])#绘制PC1
 2f_ax3.set_title('(b) PC1',loc='left')
 3f_ax3.set_ylim(-3,3)#y轴上下限
 4f_ax3.axhline(0,linestyle="-",c='k')#水平参考线
 5f_ax3.plot(years,pc[:,0],c='k')#绘制折线
 6pc1_9 = np.convolve(pc[:,0], np.repeat(1/9, 9), mode='valid')#计算九年滑动平均
 7f_ax3.plot(years[4:-4],pc1_9,c='r',lw=2)#绘制滑动平均
 8f_ax4 = fig.add_axes([0.45, 0.22, 0.3, 0.156])#绘制PC2
 9f_ax4.set_title('(d) PC2',loc='left')
10f_ax4.set_ylim(-3,3)
11f_ax4.axhline(0,linestyle="-",c='k')
12f_ax4.plot(years,pc[:,1],c='k')
13pc2_9 = np.convolve(pc[:,1], np.repeat(1/9, 9), mode='valid')
14f_ax4.plot(years[4:-4],pc2_9,c='r',lw=2)

5 收尾工作

NCL：将各个子图组合起来，并整体添加色标。由于NCL在一开始建立画布时就指定了输出文件和格式，因此这里不再需要。

1resp=True
2resp@gsnPanelRowSpec=True
3resp@gsnPanelFigureStrings=(/"a","b","c","d"/)
4resp@gsnPanelFigureStringsFontHeightF=0.01
5resp@amJust="TopLeft"
6gsn_panel(wks,(/map1,pc1,map2,pc2/),(/2,2/),resp)

Python：图形输出。

1#plt.show()
2plt.savefig("eof.pdf")

6 小节

就个人感觉而言，Python语言还是会更简洁，思路更清晰一些，尤其是在指定各个绘图参数的时候。由于这个示例并不涉及复杂数据处理部分，因此没有很好的体现python的优势。但是就图形本身而言，NCL毕竟作为专业的绘图工具还是有优势的，当然从审美的角度来看各花入个眼，看个人风格喜好吧。本来想的是可以将代码块分成两个纵列，逐条对比，但是首先MD编辑器不允许代码块分列，其次两种语言的差异也没办法逐条对比，最终只能妥协成现在这样。后边有时间的话会继续补充一些其它的示例，从数据处理等其它方面继续对比一下两种语言的差异。

如需获取测试数据和代码请回复关键字“NCL”，系统将向您推送下载链接。感谢几个群友的测试并提醒数据出错，目前已经更换。

有问题可以到QQ群里进行讨论，我们在那边等大家。

最后，如果文章对您有帮助的话，欢迎大家花几秒时间点一下赞，在看、收藏和文末的广告、并转发给有需要的同学们，谢谢！

python

文章转载自气海无涯，如果涉嫌侵权，请发送邮件至：contact@modb.pro进行举报，并提供相关证据，一经查实，墨天轮将立刻删除相关内容。