- 积分
- 20
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2018-3-31
- 最后登录
- 1970-1-1
|
登录后查看更多精彩内容~
您需要 登录 才可以下载或查看,没有帐号?立即注册
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
|
|