爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

搜索
查看: 1069|回复: 0

[作图] ncl画北半球降水分布图下载的数据是一维,但程序是三维的如何修改

[复制链接]
发表于 2018-3-31 23:04:14 | 显示全部楼层 |阅读模式

登录后查看更多精彩内容~

您需要 登录 才可以下载或查看,没有帐号?立即注册

x
ncl画北半球降水分布图下载的数据是一维,但程序是三维的如何修改
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"
load "./outnc3d.ncl"
load "./wavede3D.ncl"
load "./month_to_season_DJF.ncl"

begin
;Time Marker
;year marker
ys=1979
ye=2014
ny=ye-ys+1
nt=ny*12
year=ispan(ys,ye,1)
time=ispan(1,nt,1)
month_mark=(/"Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"/)
;Month marker
month=ispan(1,12,1)
times=ispan(1,4,1)
times@unit="season"
timem=ispan(1,12,1)
timem@unit="month"

;Switch
mon=(/2,2,3,3,4,4/)
opt_datas=(/"Z","Z","Z","Z","Z","Z"/)
np=dimsizes(opt_datas)
ilev=(/300,1000,300,1000,300,1000/)
units=(/"m","m","m","m","m","m"/)
opt_addclm=False
opt_wave=False
;land sea mask
nc=addfile("precip.nc","r")
t=short2flt(nc->t(0,:,:))
;Define plot
wks = gsn_open_wks("x11","fig/Fig3-PVO_Z_composite_Mar-May")  
;gsn_define_colormap(wks,"BlAqGrYeOrRe")        
gsn_define_colormap(wks,"BlWhRe")        
plot=new(np,graphic)
;Read data
do ip=0,np-1  ; loop at different levels
opt_data=opt_datas(ip)
if(opt_data.eq."Z")then
  nc=addfile("/media/usb/Reanalysis/ECMWF/Interim/Monthly/Z_1979-2016.nc","r")
  vars=short2flt(nc->z((ys-1979)*12:(ye-1979)*12+11,{ilev(ip)},{0:90},:))
  vars=(/vars/9.8/)
end if
if(opt_data.eq."T")then
  nc=addfile("/media/usb/Reanalysis/ECMWF/Interim/Monthly/T_1979-2016.nc","r")
  vars=short2flt(nc->t((ys-1979)*12:(ye-1979)*12+11,{ilev(ip)},{0:90},:))
end if
if(opt_data.eq."U")then
  nc=addfile("/media/usb/Reanalysis/ECMWF/Interim/Monthly/U_1979-2016.nc","r")
  vars=short2flt(nc->u((ys-1979)*12:(ye-1979)*12+11,{ilev(ip)},{0:90},:))
end if
if(opt_data.eq."SLP")then
  nc=addfile("/media/usb/Reanalysis/ECMWF/Interim/Monthly/SLP_1979-2016.nc","r")
  vars=short2flt(nc->msl((ys-1979)*12:(ye-1979)*12+11,{0:90},:))
end if
if(opt_data.eq."SST")then
  nc=addfile("/media/usb/Satellite/SST/HadleyOI/HadISST_sst_187001-201305.nc","r")
  vars=short2flt(nc->sst((ys-1870)*12:(ye-1870)*12+11,{-30:90},:))
end if
;coordinates
vars!1="lat"
vars!2="lon"
lat1=vars&lat
nlat1=dimsizes(lat1)
lon1=vars&lon
nlon1=dimsizes(lon1)
;climatology
vars0=vars
;============Detrend============
vars=dtrend_n(vars,False,0)
if(opt_wave)then
  varstmp0=wavede3D(vars0,2)
  vars0=varstmp0(:,:,:,0)
  varstmp=wavede3D(vars,2)
  vars=varstmp(:,:,:,0)
  delete([/varstmp0,varstmp/])
end if
;=============Transform=============
var =new((/12,ny,nlat1,nlon1/),float)
varclm=new((/12,ny,nlat1,nlon1/),float)
do m=0,11
  do n=0,ny-1
   var(m,n,:,:)=vars(n*12+m,:,:)
   varclm(m,n,:,:)=vars0(n*12+m,:,:)
  end do
end do
;Composite
opt_mon=(/True,True,True,True,True,True,True,True,True,True,True,True/)
nms=num(opt_mon.eq.True)

;begin composite
do m=mon(ip),mon(ip)           ;!loop over months
  vars_diff     =new((/nlat1,nlon1/),float)
  vars_prob     =new((/nlat1,nlon1/),float)
  vars_clm      =new((/nlat1,nlon1/),float)

  if(opt_mon(m))then  ;!judge at calculated months
   if(month_mark(m).eq."Jan")then
    lw=(/1985,1988,1999,2004,2006,2009,2013/)
    if(ye.ge.2014)then
    hi=(/1996,2000,2011,2014/)
    else
    hi=(/1996,2000,2011/)
    end if
   end if
   if(month_mark(m).eq."Feb")then
    lw=(/1985,1988,1999,2004,2006,2009,2013/)
    if(ye.ge.2014)then
    hi=(/1996,2000,2011,2014/)
    else
    hi=(/1996,2000,2011/)
    end if
   end if
   if(month_mark(m).eq."Mar")then
    lw=(/1985,1988,1999,2004,2006,2009,2013/)
    if(ye.ge.2014)then
    hi=(/1996,2000,2011,2014/)
    else
    hi=(/1996,2000,2011/)
    end if
   end if
   if(month_mark(m).eq."Apr")then
    lw=(/1985,1988,1999,2004,2006,2009,2013/)
    if(ye.ge.2014)then
    hi=(/1996,2000,2011,2014/)
    else
    hi=(/1996,2000,2011/)
    end if
   end if
   if(month_mark(m).eq."May")then
    lw=(/1985,1988,1999,2004,2006,2009,2013/)
    if(ye.ge.2014)then
    hi=(/1996,2000,2011,2014/)
    else
    hi=(/1996,2000,2011/)
    end if
   end if

   lw=lw-1979
   hi=hi-1979
;=================Differences================
   vars_diff(:,:)=dim_avg_n(var(m,lw,:,:),0)-dim_avg_n(var(m,hi,:,:),0)
   vars_clm=dim_avg_n(varclm(m,:,:,:),0)
;T-test significance
   do i=0,nlat1-1
    do j=0,nlon1-1
    if(.not.all(ismissing(var(m,lw,i,j))).and..not.all(ismissing(var(m,hi,i,j))))then
     ave0 = dim_avg(var(m,lw,i,j))
     var0 = dim_variance(var(m,lw,i,j))
     ave1 = dim_avg(var(m,hi,i,j))
     var1 = dim_variance(var(m,hi,i,j))
     ave2 = dim_avg(var(m,:,i,j))
     var2 = dim_variance(var(m,:,i,j))

     vars_prob(i,j) = ttest(ave0,var0,dimsizes(lw), ave1,var1,dimsizes(hi),True, False)
    end if
    end do
   end do
;attribute coordinate info.
   vars_diff!0="lat"
   vars_diff!1="lon"
   vars_diff&lat=lat1
   vars_diff&lon=lon1
   vars_diff=smth9_Wrap(vars_diff,0.5,0.25,False)
   copy_VarCoords(vars_diff,vars_prob)
   copy_VarCoords(vars_diff,vars_clm)
;Plot
   res                       = True     
   res@cnFillOn              = True     
   res@cnLinesOn             = True   
   res@cnLineLabelPlacementMode ="Computed"
   res@cnLineLabelDensityF   = 2
   res@cnLineLabelsOn  = True
   res@gsnSpreadColors       = True   
   res@gsnContourZeroLineThicknessF =2.0
   res@gsnContourNegLineDashPattern = 1
   if(opt_data.ne."SLP")then
    res@tiMainString = month_mark(m)+" Interim "+opt_data+" "+ilev(ip)+"hPa"
   else
    res@tiMainString = month_mark(m)+" Interim "+opt_data
   end if
   res@tmYLMode    = "Explicit"
   FontHeight = 0.03
   res@tmXBLabelFontHeightF = FontHeight
   res@tmYLLabelFontHeightF = FontHeight
   res@tiMainFontHeightF    = FontHeight
   res@lbTitleFontHeightF   = FontHeight
   res@lbLabelFontHeightF   = FontHeight
   res@tiMainOffsetYF =-0.006
   res@lbTopMarginF  = 0.3
   res@lbTitleString = units(ip)
   res@lbTitlePosition ="Right"
   res@lbTitleDirection ="Across"
   res@lbTitleAngleF   = 90
   res@mpMaxLatF =85
   res@mpMinLatF =0
   res@mpMinLonF =90
   res@mpMaxLonF =320
   res@mpCenterLonF=180
   res@mpGeophysicalLineThicknessF = 2.
   res@gsnAddCyclic= False
   res@gsnMaximize = False
   res@tmYRLabelsOn  = False
   res@lbOrientation = "vertical"
   res@lbLabelOffsetF = 0.2
   res@lbLeftMarginF  = 0.33
   res@gsnFrame =False
   res@gsnDraw =False
   res@cnInfoLabelOn       = False      
   res@gsnRightString      = ""
   res@gsnLeftString       = ""
   res@tmXBMode   = "Explicit"
   res@tmXBValues = ispan(0,360,60)
   res@tmXBLabels = res@tmXBValues+"~S~o~N~E"
   res@tmXBMinorValues = ispan(0,360,10)
   res@mpOceanFillColor = -1
   res@cnLevelSelectionMode = "ExplicitLevels"
;   res@cnLevels = ispan(-10,10,2)
;significance shaded
   res0=True
   res0@cnLinesOn            = False
   res0@cnLineLabelsOn       = False
   res0@gsnFrame             = False
   res0@gsnDraw              = False
   res0@cnFillOn             = True
   res0@cnLevels             = (/0.,0.1,.2,.3,.4,.5,.6,.7,.8,.9,1./)
   res0@cnFillColor          = "black"
   res0@cnMonoFillColor      = True
   res0@cnMonoFillPattern    = False            ; want multiple patterns
   res0@cnFillPatterns       = (/17,17 ,-1 ,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1/) ; the patterns
   res0@cnMonoFillScale      = False            ; want different densities
   res0@cnFillScales         = (/0.6,0.6,.3,1.5,2.0,2.5/) ; the densities
   res0@lbLabelBarOn         = False
   res0@cnInfoLabelOn        = False  
;contour
   resc = True
   resc@gsnContourZeroLineThicknessF =2.0
   resc@gsnContourNegLineDashPattern = 1
   resc@gsnDraw  = False
   resc@gsnFrame = False
   resc@cnInfoLabelOn       = False      
   resc@cnLineLabelsOn       = False

   if(opt_addclm)then
    res@cnLinesOn = False
    plotc=gsn_csm_contour(wks,vars_clm,resc)
   end if
   if(opt_data.eq."SST")then
    vars_diff=(/where(abs(vars_diff).ge.2.5,vars_diff@_FillValue,vars_diff)/)
   end if
   if(opt_data.eq."TS")then
    vars_diff=(/mask(vars_diff,vars_diff.eq.0,False)/)
   end if
   if(opt_data.eq."vor")then
    do n=1,3
    vars_diff=smth9_Wrap(vars_diff,0.5,0.25,False)
    end do
   end if
   if(opt_data.eq."prep".or.opt_data.eq."q")then
    vars_diff=(/where(abs(vars_diff).ge.50,vars_diff@_FillValue,vars_diff)/)
    do n=1,4
     vars_diff=smth9_Wrap(vars_diff,0.5,0.25,False)
    end do
   end if
   if(opt_data.eq."prepUSA")then
    do n=1,5
     vars_diff=smth9_Wrap(vars_diff,0.5,0.25,False)
    end do
   end if
   plot(ip)= gsn_csm_contour_map_ce(wks,vars_diff,res)
   plot0= gsn_csm_contour(wks,vars_prob,res0)
   overlay(plot(ip),plot0)
   if(opt_addclm)then
    overlay(plot(ip),plotc)
   end if
   delete([/lw,hi/])
  end if  ;!judge at calculated months
end do  ;loop over different months
delete([/vars,vars0,lat1,lon1,var,varclm,vars_diff,vars_prob,vars_clm/])
end do ;loop at different levels

resp=True
resp@gsnPanelFigureStrings= (/"(a)","(b)","(c)","(d)","(e)","(f)"/)
resp@amJust   = "TopLeft"
resp@amParallelPosF =-0.65
resp@gsnPanelRight  = 0.95
resp@amOrthogonalPosF = -0.8
resp@gsnPanelFigureStringsPerimOn  =False
resp@lbLabelFontHeightF  = 0.015
resp@gsnPanelFigureStringsFontHeightF  =0.016
gsn_panel(wks,plot,(/3,2/),resp)
end


密码修改失败请联系微信:mofangbao
您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表,不代表气象家园立场,禁止在本站发表与国家法律相抵触言论

快速回复 返回顶部 返回列表