通過示例對比Python和NCL，感受兩種語言的區(qū)別。
NCL語言大概是從2018年開始，陪伴我走過了一年的時間，也幫助我發(fā)表了自己的第一篇文章。后來由于自己需要實現(xiàn)的算法過于復雜，自己無力編寫，在朋友的誘惑下不得不轉(zhuǎn)投Python，靠網(wǎng)上各路大神的帖子東拼西湊滿足自己的科研需求。說卸磨殺驢也有點不合適，但在適應了Python之后，我確實基本沒有再打開過NCL了。最近整理文件，找到了自己發(fā)表第一篇學術垃圾時候的NCL腳本，心血來潮的想再用Python實現(xiàn)一遍。也算是對評論里很多讀者的一個回應吧(貌似氣象家園帖子下邊第一個評論就是希望我寫一個兩者對比的文章，被我鴿到現(xiàn)在，咕咕咕)。注：編程水平僅限于能跑就行，warning不管，冗雜語句請忽視。
示例1（EOF）
選擇EOF的原因是圖片內(nèi)容較為豐富，同時方法較為常用
由于原始數(shù)據(jù)過大，只提供處理好的數(shù)據(jù)方便測試(數(shù)據(jù)為每年的寒潮路徑的概率密度分布，為39年×90緯度×360經(jīng)度，由FLEXPART模式追蹤后通過高斯核密度估計算法Gaussian KDE得到。
測試數(shù)據(jù)和代碼見原鏈接
先對比下出圖效果(兩種語言對圖形的渲染有差異，EOF算法可能也有些差異，但是結果是基本相同的，圖只經(jīng)過了基礎的調(diào)整，并不是很好看)。
1 準備工作
NCL：在引入一些特殊函數(shù)(NCL本體不具備的函數(shù))時，通常會加上類似于以下語句：
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
Python：引入各個模塊(庫)：
import xarray as xr
import numpy as np
from eofs.standard import Eof
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt
import cartopy.mpl.ticker as cticker
def contour_map(fig,img_extent,spec):
fig.set_extent(img_extent, crs=ccrs.PlateCarree())
fig.add_feature(cfeature.COASTLINE.with_scale('50m'))
fig.add_feature(cfeature.LAKES, alpha=0.5)
fig.set_xticks(np.arange(leftlon,rightlon+spec,spec), crs=ccrs.PlateCarree())
fig.set_yticks(np.arange(lowerlat,upperlat+spec,spec), crs=ccrs.PlateCarree())
lon_formatter = cticker.LongitudeFormatter()
lat_formatter = cticker.LatitudeFormatter()
fig.xaxis.set_major_formatter(lon_formatter)
fig.yaxis.set_major_formatter(lat_formatter)
2 數(shù)據(jù)讀取
NCL：
f=addfile("out.nc","r");讀取文件
dot=f->cspath(:,:,{0:180});讀入變量
dot:=dot(lat|:,lon|:,year|:);調(diào)整變量維度順序(EOF函數(shù)對維度順序有要求)
x=ispan(1979,2017,1);時間
Python：利用Xarray庫讀取NC文件
f = xr.open_dataset("out.nc" )
dot = np.array(f['cspath'].loc[:,0:90,:])
lat = f['lat'].loc[0:90]
lon = f['lon']
years = range(1979, 2018)
3 EOF
NCL：
eof=eofunc_Wrap(dot,3,False)
pc=eofunc_ts_Wrap(dot,eof,False)
pc=dim_standardize_n(eof_ts,1,1)
copy_VarMeta(dot(:,:,0),eof1);將維度信息重新賦值給新數(shù)組
copy_VarMeta(dot(:,:,0),eof2)
Python：利用EOF庫
solver = Eof(dot)
eof = solver.eofsAsCorrelation(neofs=2)
pc = solver.pcs(npcs=2, pcscaling=1)
var = solver.varianceFraction()
4.1 繪圖：建立畫布
NCL的有些繪圖參數(shù)并不是必要的，由于年代久遠，我記不清起每條語句對應的詳細功能了。
NCL：
wks=gsn_open_wks("pdf","tmp")
res = True
res@gsnMaximize = False
res@gsnDraw = False
res@gsnFrame = False
res@gsnDraw=False
res@gsnFrame=False
respc=res;實際上是設置PC圖形的基礎繪圖參數(shù)
Python：
fig = plt.figure(figsize=(15,15))
4.2 繪圖：繪制EOF
NCL：其繪圖思路是每一條語句指定一個繪圖效果，通過不斷的”疊BUFF“實現(xiàn)全部要素的疊加，先將共同要素疊給res，然后通過res1=res, res2=res賦值后，再對res1和res2分別添加各自的屬性。
res@mpMaxLatF=90;設置經(jīng)緯度邊界
res@mpMinLatF=0
res@mpMaxLonF=160
res@mpMinLonF=0
res@mpFillOn =False;地圖設置
res@mpOutlineOn = True
res@mpDataBaseVersion = "MediumRes"
res@mpDataSetName = "/data/home/Earth..4"
res@cnFillPalette = "MPL_bwr"
res@cnFillOn = True
res@cnFillDrawOrder = "PreDraw"
res@cnLinesOn = False
res@cnLineLabelsOn = False
res@gsnLeftString=""
res@pmTickMarkDisplayMode = "Always"
res@cnLevelSelectionMode="ExplicitLevels"
res@cnLevels =(/-0.05,-0.04,-0.03,-0.02,-0.01,0,0.01,0.02,0.03,0.04,0.05/)
res1=res
res2=res
res1@gsnRightString = "41.84%"
res1@gsnLeftString = "EOF1"
res1@lbLabelBarOn=False
res2@gsnRightString = "14.59%"
res2@gsnLeftString = "EOF2"
res2@lbLabelBarOn=True
map1 = gsn_csm_contour_map(wks,eof1,res1) ;繪圖
map2 = gsn_csm_contour_map(wks,eof2,res2)
Python：與NCL相似的地方在于需要對每個axes添加參數(shù)，不同的地方在于其一條指令可以包含多個效果(比如將所有設置填色繪圖的參數(shù)全部加在f_ax1.contourf里)。
proj = ccrs.PlateCarree(central_longitude=80) #投影
leftlon, rightlon, lowerlat, upperlat = (0,160,10,90) #地圖邊界
img_extent = [leftlon, rightlon, lowerlat, upperlat]
f_ax1 = fig.add_axes([0.1, 0.32, 0.3, 0.4],projection = proj)#EOF1
contour_map(f_ax1,img_extent,20) #利用前邊自定義的繪圖函數(shù)，目的是省去繪制相同圖形時需要再次設置地形，湖泊，軸刻度等參數(shù)
f_ax1.set_title('(a) EOF1',loc='left')#左標題
f_ax1.set_title( '%.2f%%' % (var[0]*100),loc='right')#右標題，解釋方差
f_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)#繪制填色
f_ax2 = fig.add_axes([0.1, 0.1, 0.3, 0.4],projection = proj)#EOF2
contour_map(f_ax2,img_extent,20)
f_ax2.set_title('(c) EOf',loc='left')
f_ax2.set_title( '%.2f%%' % (var[1]*100),loc='right')
c2=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)
#繪制色標
position=fig.add_axes([0.1, 0.17, 0.3, 0.017])
fig.colorbar(c2,cax=position,orientation='horizontal',format='%.1f',)
4.3 繪圖：繪制PC
NCL：
respc@gsnYRefLine=0
respc@trYMinF=-3
respc@trYMaxF=3
res9 = respc
respc@tmYLLabelDeltaF=-1
respc@trXMinF=1979
respc@trXMaxF=2017
respc@tiYAxisOn=False
respc@tmXTOn =False
respc@tmYROn = False
respc@vpHeightF=0.39
respc@vpWidthF=0.6
respc@gsnLeftStringOrthogonalPosF = 0.04
respc1 = respc
respc1@gsnLeftString = "PC1"
respc2 = respc
respc2@gsnLeftString = "PC2"
pc1= gsn_csm_xy(wks,x,eoft1,respc1)
pc2= gsn_csm_xy(wks,x,eoft2,respc2)
pc1_9 = runave_n_Wrap(eoft1, 9, 0, 0);9年滑動平均
pc2_9 = runave_n_Wrap(eoft2, 9, 0, 0)
res9@xyLineColor = "red"
res9@xyLineThicknesses = 4
plotpc9_1 = gsn_csm_xy(wks, x, pc1_9, res9)
plotpc9_2 = gsn_csm_xy(wks, x, pc2_9, res9)
overlay(pc1, plotpc9_1);疊加
overlay(pc2, plotpc9_2)
Python：
f_ax3 = fig.add_axes([0.45, 0.44, 0.3, 0.156])#繪制PC1
f_ax3.set_title('(b) PC1',loc='left')
f_ax3.set_ylim(-3,3)#y軸上下限
f_ax3.axhline(0,linestyle="-",c='k')#水平參考線
f_ax3.plot(years,pc[:,0],c='k')#繪制折線
pc1_9 = np.convolve(pc[:,0], np.repeat(1/9, 9), mode='valid')#計算九年滑動平均
f_ax3.plot(years[4:-4],pc1_9,c='r',lw=2)#繪制滑動平均
f_ax4 = fig.add_axes([0.45, 0.22, 0.3, 0.156])#繪制PC2
f_ax4.set_title('(d) PC2',loc='left')
f_ax4.set_ylim(-3,3)
f_ax4.axhline(0,linestyle="-",c='k')
f_ax4.plot(years,pc[:,1],c='k')
pc2_9 = np.convolve(pc[:,1], np.repeat(1/9, 9), mode='valid')
f_ax4.plot(years[4:-4],pc2_9,c='r',lw=2)
5 收尾工作
NCL：將各個子圖組合起來，并整體添加色標。由于NCL在一開始建立畫布時就指定了輸出文件和格式，因此這里不再需要。
resp=True
resp@gsnPanelRowSpec=True
resp@gsnPanelFigureStrings=(/"a","b","c","d"/)
resp@gsnPanelFigureStringsFontHeightF=0.01
resp@amJust="TopLeft"
gsn_panel(wks,(/map1,pc1,map2,pc2/),(/2,2/),resp)
Python：圖形輸出。
#plt.show()
plt.savefig("eof.pdf")
6 小節(jié)
就個人感覺而言，Python語言還是會更簡潔，思路更清晰一些，尤其是在指定各個繪圖參數(shù)的時候。由于這個示例并不涉及復雜數(shù)據(jù)處理部分，因此沒有很好的體現(xiàn)python的優(yōu)勢。但是就圖形本身而言，NCL畢竟作為專業(yè)的繪圖工具還是有優(yōu)勢的，當然從審美的角度來看各花入個眼，看個人風格喜好吧。本來想的是可以將代碼塊分成兩個縱列，逐條對比，但是首先MD編輯器不允許代碼塊分列，其次兩種語言的差異也沒辦法逐條對比，最終只能妥協(xié)成現(xiàn)在這樣。后邊有時間的話會繼續(xù)補充一些其它的示例，從數(shù)據(jù)處理等其它方面繼續(xù)對比一下兩種語言的差異。
