脚本主要内容：
- NetCDF数据读取
- 计算和存储
- 等值线、矢量箭头、图层叠加
脚本略有缺失，完整脚本请购买施宁教授出的《NCL数据处理与绘图实习手册》纸质书籍。版权归施宁教授所有。
效果如下图所示：
另外推荐高端封面配图的书：
由北京大学、中科院化学所等拥有物理化学生物等专业博士学位团队制作，相关文章发表Nat.Commun./Sci. Adv./Adv. Mater.等SCI论文期刊，拥有深厚的科研背景。
代码：
代码语言：javascript复制load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
begin
year=ispan(1979,2013,1) ; 79/80 - 13/14
it_s=197912 ;起始年月
it_e=201411 ;结束年月
refmag = 3 ;参考箭头所表示的风速大小
;;;read data ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; sst
time = f_sst->time ; 读取其日期
; 转换成公历日期
; 截取指定时间段
;
;
;; h300
f_h300 = addfile("./data/h300-197901-201412.nc", "r")
h300 = short2flt(f_h300->hgt(rec_s:rec_e,0,{-90:90},:))
;; u850
f_u850 = addfile("./data/u850-197901-201412.nc", "r")
u850 = short2flt(f_u850->uwnd(rec_s:rec_e,0,{-90:90},:)) ; 850 hPa
;; v850
f_v850 = addfile("./data/v850-197901-201412.nc", "r")
v850 = short2flt(f_v850->vwnd(rec_s:rec_e,0,{-90:90},:)) ; 850 hPa
;; air2m
f_air2m = addfile("./data/air2m-197901-201412.nc", "r")
air2m = short2flt(f_air2m->air(rec_s:rec_e,0,{-90:90},:)) ; T at 2m
;;;DJF 平均 & 异常 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;JFM季节平均，实际是12/1/2月三个月平均，因为从1979年12月开始截取
copy_VarMeta(sst(0,:,:),sst_DJF(0,:,:))
sst_DJF!0 = "year"
sst_DJF&year=year
sst_ano = dim_rmvmean_n_Wrap(sst_DJF,0)
;; h300
h300_DJF = month_to_season(h300, "JFM")
copy_VarMeta(h300(0,:,:),h300_DJF(0,:,:))
h300_DJF!0 = "year"
h300_DJF&year=year
h300_ano = dim_rmvmean_n_Wrap(h300_DJF,0)
;; u850 与h300 同维
u850_DJF = month_to_season(u850, "JFM")
copy_VarMeta(h300_DJF,u850_DJF)
u850_ano = dim_rmvmean_n_Wrap(u850_DJF,0)
;; v850 与h300 同维
v850_DJF = month_to_season(v850, "JFM")
copy_VarMeta(h300_DJF,v850_DJF)
v850_ano = dim_rmvmean_n_Wrap(v850_DJF,0)
;; air2m
air2m_DJF = month_to_season(air2m, "JFM")
copy_VarMeta(air2m(0,:,:),air2m_DJF(0,:,:))
air2m_DJF!0 ="year"
air2m_DJF&year=year
air2m_ano = dim_rmvmean_n_Wrap(air2m_DJF,0)
;;;(3) enso index (5N-5S, 170-120W);;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; 0表示仅用非缺省的数值进行计算
;1 表示标准化时除以的是[N] ; 而0表示除以[N-1]
;; 输出至netcdf文件
path_out = "ENSO-index.nc"
system("rm -f " path_out) ; 若当前路径下有同名文件，则删除
ncdf = addfile(path_out,"c") ; "c" 表示创建 netCDF 文件
ncdf->ensoi = ensoi
;;;(4) composite ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
nnumb = dimsizes(irec_positive)
h300_comp = dim_avg_n_Wrap(h300_ano(irec_positive,:,:),0)
u850_comp = dim_avg_n_Wrap(u850_ano(irec_positive,:,:),0)
v850_comp = dim_avg_n_Wrap(v850_ano(irec_positive,:,:),0)
air2m_comp = dim_avg_n_Wrap(air2m_ano(irec_positive,:,:),0)
;;; (5) t-test ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; sst
;; h300
h300_std = dim_variance_n_Wrap(h300_ano(irec_positive,:,:),0)
h300_std = sqrt(h300_std/nnumb)
t_h300 = h300_comp/h300_std
confi_h300 = h300_comp
confi_h300 = student_t(t_h300, nnumb-1)
;; air2m
air2m_std = dim_variance_n_Wrap(air2m_ano(irec_positive,:,:),0)
air2m_std = sqrt(air2m_std/nnumb)
t_air2m = air2m_comp/air2m_std
confi_air2m = air2m_comp
confi_air2m = student_t(t_air2m, nnumb-1)
;;; (5) plot
wks = gsn_open_wks("eps","plot-comp-enso")
gsn_define_colormap(wks,"rainbow gray") ; 调用rainbow gray色板,，其它色板名称请查阅http://www.ncl.ucar.edu/Document/Graphics/color_table_gallery.shtml
base = new(3,"graphic")
plot = new(3,"graphic")
res = True ; 调整地图及显著性等值线, 每个子图均需该res
res@gsnAddCyclic = True ; 添加循环点，否则会在0度经线左侧出现一根白条
res@gsnDraw = False
res@gsnFrame = False
res@gsnLeftString = ""
res@gsnRightString = ""
resc = res ;拷贝
resv = res ;
rest = res ;
res@mpFillOn = False ; 不填色地图
res@mpCenterLonF = 180 ; 地图的中心经度
res@mpGeophysicalLineThicknessF = 0.5 ; 地图边界的粗细
res@pmTickMarkDisplayMode= "Always" ; 坐标上标签上添加度符号
res@mpGridAndLimbOn = True ; 绘制经纬度线
; 经纬度线间隔
;
; 经纬度线线型取为类型为2的虚线。共17种线型供选择。
; 其粗细
res@cnFillOn = True ; 填色等值线
res@cnLinesOn = True ; 绘制等值线
res@cnLineColor = "white" ; 颜色
res@cnLineThicknessF = 0.3 ; 粗细
res@cnLineLabelsOn = False ; 关闭标签
; 用GMT_gray 进行填色。即调用了第2种色板
; -1 为透明
res@cnInfoLabelOn = False ; 关闭图右下方的等值线信息标签
res@lbLabelBarOn = False ; 关闭labelbar
resc@cnLevelSelectionMode = "ExplicitLevels" ; 指定每根需绘制的等值线
resc@cnLevels = (/-0.75,-0.25,0.25,0.75,1.25/) ;
resc@cnFillOn = False ; 关闭等值线填色
resc@cnLineThicknessF = 2. ; 等值线粗细
resc@gsnContourZeroLineThicknessF = 0. ; 设置0值线粗细。0则不画
resc@cnLineLabelsOn = False ; 关闭标签
resc@cnLineDashPattern = 16 ; 线型为16的虚线
resc@cnInfoLabelOn = True ; 打开图右下方的等值线信息标签
resc@cnInfoLabelOrthogonalPosF = 0.05 ; 移动等值线信息标签的位置
resv@vcPositionMode = "ArrowTail" ;箭头尾部对应着格点的位置
resv@vcGlyphStyle = "Fillarrow" ;其余三种选项为“LineArrow”、“WindBarb” 、“CurlyVector”
resv@vcFillArrowEdgeThicknessF = 2 ; 箭头边界粗细
resv@vcFillArrowEdgeColor = "white" ; 及颜色
resv@vcFillArrowFillColor = "black" ; 箭头内部填充颜色
resv@vcFillArrowWidthF = 0.1 ; 箭头宽度
resv@vcFillArrowHeadXF = 0.6 ; 请参考附录中Fillarrow箭头示意图
resv@vcFillArrowHeadYF = 0.2 ;
resv@vcFillArrowHeadInteriorXF = 0.25 ;
resv@vcMinDistanceF = 0.03 ; 箭头之间的最小距离（在单位平方中）
resv@vcMinMagnitudeF = 1.0 ; 要绘制箭头所表示的最小数值，即小于该数值则不绘制
resv@vcFillArrowMinFracWidthF =1.0
resv@vcFillArrowHeadMinFracXF =1.0
resv@vcFillArrowHeadMinFracYF =1.0
;****设定参考箭头****
resv@vcRefAnnoOn = True
resv@vcRefMagnitudeF = refmag ;标准长度箭头所表示的大小
resv@vcRefLengthF = 0.045 ;标准长度箭头在单位平方中的大小
resv@vcRefAnnoBackgroundColor = "white" ;背景颜色
resv@vcRefAnnoPerimOn = False ;关闭边框
resv@vcRefAnnoFontHeightF = 0.015 ;参考箭头标签字体大小
resv@vcRefAnnoString1On = False ;设定参考箭头上、下的字符
resv@vcRefAnnoString2On = True ; 这里仅设定其下方的字符
resv@vcRefAnnoString2 = refmag " m/s"
resv@vcRefAnnoSide = "Top" ; 参考箭头放至图形上方
resv@vcRefAnnoOrthogonalPosF = -0.12 ; 调整其位置
resv@vcRefAnnoParallelPosF = 0.95
res@gsnCenterString = "sst" ;子图的主标题
res@gsnCenterStringFontHeightF = 0.03 ; 标题字体的大小。由于后面没有修改该值，则每幅图的主标题字体均是此大小
; 只有底图可有地图（map）
; 调用的绘图函数不可带“map”
plot(0) = ColorNegDashZeroPosContour(plot(0),"blue","white","red") ; 负值用蓝色虚线表示，0线用白色实线，正值红色实线
; 带地图的图必须放在最下图层
; 绘制多边形及折线以标明nino 3.4区
plres = True
plres@gsLineColor = "black"
plres@gsLineThicknessF = 1.0
gres = True
gres@gsFillColor = "yellow"
gres@gsFillOpacityF = 0.5
gres@gsLineColor = "black"
latx = (/-5, 5, 5, -5, -5/) ; nino3.4区的坐标位置
lonx = (/190, 190,240, 240, 190/) ;
dum1 = gsn_add_polyline(wks, base(0),lonx,latx,plres)
dum2 = gsn_add_polygon(wks,base(0),lonx,latx,gres)
res@gsnCenterString = "h300&V850"
resc@cnLevelSelectionMode = "AutomaticLevels"
resc@cnLevelSpacingF = 15.
base(1) = gsn_csm_contour_map(wks,confi_h300,res)
plot(1) = gsn_csm_contour(wks,h300_comp,resc)
plot(1) = ColorNegDashZeroPosContour(plot(1),"blue","white","red")
overlay(base(1),plot(1))
plotv = gsn_csm_vector(wks,u850_comp,v850_comp,resv)
overlay(base(1),plotv) ; 也可用gsn_csm_vector_map(wks,h300_comp,u850,v850,res_new)
res@gsnCenterString = "air2m"
resc@cnLevelSelectionMode = "ManualLevels"
resc@cnMaxLevelValF = 2
resc@cnMinLevelValF = -2
resc@cnLevelSpacingF = 0.5
base(2) = gsn_csm_contour_map(wks,confi_air2m,res)
plot(2) = gsn_csm_contour(wks,air2m_comp,resc)
plot(2) = ColorNegDashZeroPosContour(plot(2),"blue","black","red")
overlay(base(2),plot(2))
resP = True ; 绘制panel图
resP@txString = "El nino" ; 添加主标题
resP@txFontHeightF = 0.03 ; 修改其大小
; resP@gsnPanelFigureStrings= (/"a)","b)","c)"/) ;各个子图的标号
resP@gsnPanelFigureStringsFontHeightF = 0.015 ;字体的大小
resP@amJust = "TopLeft" ;摆放的位置，默认是“BottomRight”
; 指定每行绘制的子图的个数
; 第1行绘制1幅，第2行绘制2幅
end
