Date: Sun May 06 2012 - 14:25:30 MDT
Dear all,
i used wavelet analysis to perform olr in intraseasonal scale.
I just wanted to eliminate some colour (the blue).can some one help me?
This is the script i used and plot result.
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 "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/time_axis_labels.ncl"
;load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/diagnostics_cam.ncl"
;--------------------------------------------------------------------
; CREATE SAME WAVELET FIGURE AS TORRENCE & COMPO
;--------------------------------------------------------------------
begin
latS = -5.
latN = 10.
lonL = 5.
lonR = 27.5
; latS = 3.5
; latN = 7.
; lonL = -5.
; lonR = 5.
minprd = 1.
maxprd = 150.
;ZONE = (/"ZONE3","ZONE5"/)
ZONE = (/"ZONE1"/)
NZONE = dimsizes(ZONE)
VARNAME = (/"olr","pr_wtr"/) ;
; VARNAME = (/"olr","data","Tb","precip"/) ;
NVAR = dimsizes(VARNAME)
runmeannts = 5
SEASON = "MAMJ"
csig = "ok"
;********************************************************************************
; initial resource settings
;********************************************************************************
if (csig.eq."ok") then
wks = gsn_open_wks("ps","AIC.WAVELET_zone_AnCy_PLOTfull_"+SEASON)
; open ps file
else
wks = gsn_open_wks("ps","WAVELET_zone_AnCy_PLOTfull") ; open ps file
end if
gsn_define_colormap(wks,"BlAqGrYeOrReVi200")
ncolor=192
do ivar=0,NVAR-1 ; #################### DO VAR
VAR = VARNAME(ivar)
if (VAR.eq."olr") then
diri = "/media/VERBATIM/CLIMSERV_DATA/"
fili = "1bpf_10_60dayolr.79.10.nc"
f = addfile (diri+fili, "r")
timeUnits = f->time@units
printVarSummary(timeUnits)
startDate = ut_inv_calendar( 1979, 01, 01, 00, 0, 0, timeUnits, 0 )
endDate = ut_inv_calendar( 2010, 12, 31, 00, 0, 0, timeUnits, 0 )
byear=2006
eyear=2009
in_datapier = f->$VAR$(:,:,:)
printVarSummary(in_datapier)
delete(in_datapier)
in_data0 = short2flt(f->$VAR$(:,:,:))
in_data0 = lonFlip(in_data0)
end if
if (VAR.eq."pr_wtr") then
diri = "/media/VERBATIM/CLIMSERV_DATA/"
fili = "premean.nc"
f = addfile (diri+fili, "r")
timeUnits = f->time@units
printVarSummary(timeUnits)
startDate = ut_inv_calendar( 1979, 01, 01, 00, 0, 0, timeUnits, 0 )
endDate = ut_inv_calendar( 2010, 12, 31, 00, 0, 0, timeUnits, 0 )
byear=2006
eyear=2009
in_data0 = short2flt(f->$VAR$(:,:,:))
end if
in_dataf = in_data0({startDate:endDate},{latS:latN},:)
Gtime = in_data0&time
Gtime!0="time"
Gtime&time= Gtime
Gtunits= timeUnits
delete(in_data0)
olr = in_dataf(:,:,{lonL:lonR})
N = dimsizes(olr&time)
delete(in_dataf)
TALON = olr(0,:,:) *0.0 + 3.
X1D0 = ndtooned(TALON)*0.0 + 3.
idxp = ind_resolve(ind(X1D0.eq.3.), dimsizes(TALON))
delete(TALON)
Nidxp = dimsizes(idxp(:,0))
X1 = new((/N,Nidxp/),"float")
X = new(N,"float")
if (Nidxp.gt.0)
do iNidxp = 0,Nidxp-1
print("******** grid's point number "+iNidxp+" / "+Nidxp+" ********")
X = olr(:,idxp(iNidxp,0),idxp(iNidxp,1))
X1(:,iNidxp) = X
;************************************
; compute wavelet
;************************************
mother = 0 ;An integer giving the mother wavelet to use
param = 6.0
dt = 1. ;timesteps in units of years
s0 = dt
dj = 0.25
jtot = 1+floattointeger(((log10(N*dt/s0))/dj)/log10(2.))
npad = N
nadof = 0
noise = 1
siglvl = .05
isigtest = 0
w = wavelet(X,mother,dt,param,s0,dj,jtot,npad,noise,isigtest,siglvl,nadof)
power = onedtond( w@power, (/jtot,N/) ) ; 2-D ARRAY CONTAINING
THE SQUARED SUM OF THE REAL IMAGINARY PART OF W
power!0 = "period" ; Y axis
power&period = w@period ; 1-D ARRAY CONTAINING THE FOURIER
PERIOD CORRESPONDING TO W@SCALE
power!1 = "time" ; X axis
power&time = olr&time
power@long_name = "Power Spectrum"
power@units = "C^2"
SIG = power ; transfer meta data
SIG = power/conform (power,w@signif,0)
SIG@long_name = "Significance"
SIG@units = " "
SCALE = w@scale
do it=0,N-1
power(:,it) = power(:,it)/SCALE
end do
if (iNidxp.eq.0)
POWERzone = new((/jtot,N,Nidxp/),"float")
SIGzone = new((/jtot,N,Nidxp/),"float")
end if
POWERzone(:,:,iNidxp) = power
SIGzone(:,:,iNidxp) = SIG
minprd = 1.1; min(w@period)
maxprd = max(w@period)
delete(SCALE)
if (iNidxp.ne.Nidxp-1)
delete(w)
end if
delete(SIG)
delete(power)
delete(X)
end do
power = dim_avg_n(POWERzone,2)
power!0 = "period" ; Y axis
power&period = w@period ; 1-D ARRAY CONTAINING THE FOURIER
PERIOD CORRESPONDING TO W@SCALE
power!1 = "time" ; X axis
power&time = olr&time
power@long_name = "Power Spectrum"
power@units = "C^2"
SIG = dim_avg_n(SIGzone,2)
SIG!0 = "period" ; Y axis
SIG&period = w@period ; 1-D ARRAY CONTAINING THE FOURIER
PERIOD CORRESPONDING TO W@SCALE
SIG!1 = "time" ; X axis
SIG&time = olr&time
SIG@long_name = "Significance"
delete(POWERzone)
delete(SIGzone)
X = dim_avg_n(X1,1)
delete(X1)
X!0 = "time"
X&time = olr&time
w = wavelet(X,mother,dt,param,s0,dj,jtot,npad,noise,isigtest,siglvl,nadof)
WSTDEV = w@stdev
sqrWSTDEV = WSTDEV*WSTDEV
SDOF = w@dof(0)
TRESHOLD = 1. ;WSTDEV*WSTDEV ;chiinv(pHigh,SDOF) / (1. * SDOF)
;sqrWSTDEV * chiinv(pHigh,SDOF) / (1. * SDOF)
SDOF@xlag1 = w@r1
idxprd = ind_resolve(ind(w@period.le.maxprd.and.w@period.ge.minprd),
dimsizes(w@period))
WFRQ = w@period(idxprd(:,0))
fft_theor = w@fft_theor(idxprd(:,0))
signif = w@signif(idxprd(:,0))
delete(idxprd)
NWFRQ = dimsizes(WFRQ)
FRQ = new((/NWFRQ/),"float")
do ijtot=0,NWFRQ-1
FRQ (ijtot) = 1./WFRQ(ijtot)
end do
SDOF@frq = FRQ
delete(FRQ)
delete(WFRQ)
powersig = power
powersig@long_name = "Global Power Spectrum"
Nyear= eyear-byear+1
Wplot = new(Nyear,graphic)
vl11 = new(Nyear,graphic)
vl12 = new(Nyear,graphic)
vl21 = new(Nyear,graphic)
vl22 = new(Nyear,graphic)
do iyear=byear,eyear ; #################### DO YEAR
print("****** "+VAR+" "+iyear+" *****" )
power = powersig
sig = SIG
;################ SERIE TEMPORELLE DU CYCLE REGIONAL MOYEN
##########################################
GTunits = "days since 1800-01-01 00:00:00"
NstartDate = ut_inv_calendar( iyear, 01, 01, 00, 0, 0, timeUnits, 0 )
NendDate = ut_inv_calendar( iyear, 12, 31, 00, 0, 0, timeUnits, 0 )
xver11 = ut_inv_calendar( iyear, 06, 01, 00, 0, 0, timeUnits, 0 )
xver11 = ut_convert( xver11, GTunits )
xver12 = ut_inv_calendar( iyear, 09, 30, 00, 0, 0, timeUnits, 0 )
xver12 = ut_convert( xver12, GTunits )
if (SEASON.eq."annual") then
NstartDate = ut_inv_calendar( iyear, 01, 01, 00, 0, 0, timeUnits, 0 )
NendDate = ut_inv_calendar( iyear, 12, 31, 00, 0, 0, timeUnits, 0 )
else
NstartDate = ut_inv_calendar( iyear, 06, 01, 00, 0, 0, timeUnits, 0 )
NendDate = ut_inv_calendar( iyear, 09, 30, 00, 0, 0, timeUnits, 0 )
end if
POWERprdmax = power({minprd:maxprd},{NstartDate:NendDate})
SIGprdmax = sig({minprd:maxprd},{NstartDate:NendDate})
if (VAR.eq."olr") then
olrTS = POWERprdmax(time|:,period|:)
olrTS&time = ut_convert( olrTS&time, GTunits )
olrTSsig = SIGprdmax(time|:,period|:)
olrTSsig&time = ut_convert( olrTSsig&time, GTunits )
else
end if
delete(POWERprdmax)
delete(power)
delete(SIGprdmax)
delete(sig)
res = True ; plot mods desired
res@gsnDraw = False ; Do not draw plot
res@gsnFrame = False ; Do not advance frome
res@cnFillOn = True ; turn on color
res@cnFillMode = "RasterFill" ; turn on raster mode
res@cnRasterSmoothingOn = True ; turn on raster smoothing
res@cnLinesOn = False ; turn off contour lines
res@cnLineLabelsOn = False
res@cnInfoLabelOn = False
res@lbLabelBarOn = False ; turn off individual lb's
res@gsnSpreadColors = True ; use full colormap
res@gsnSpreadColorStart = 1
res@gsnSpreadColorEnd = -1
res@tmLabelAutoStride = False
res@trYReverse = True ; reverse y-axis
res@vpHeightF = .5 ;
res@vpWidthF = 1.
res@gsnRightString = "Wavelet Power"
res@tiYAxisString = "Period (Days)"
; =================== Set special resources for the time axis
===================
resTick = True
resTick@ttmFormat = "%c"
resTick@ttmAxis = "XB"
resTick@ttmMajorStride = 31
if (VAR.eq."olr") then
res@gsnLeftString = "Olr-NOAA"
res@gsnCenterString = " ";zone
res@gsnRightString = iyear
power = olrTS(period|:,time|:) ;/1000.
power!0 = "period"
power&period = olrTS&period
power!1 = "time"
power&time = olrTS&time
sig = power
sig = olrTSsig(period|:,time|:)
delete(olrTS)
delete(olrTSsig)
printVarSummary(power)
end if
res@tmYLMode = "Explicit"
; res@tmYLValues = (/2,5,10,25,40,60,90/)
res@tmYLValues = (/10,30,60,90/) ;power&period
res@tmYLLabels = "" + res@tmYLValues
res@cnLevelSelectionMode = "ManualLevels"; manual set levels so lb consistent
res@cnMinLevelValF = min(power) ; min level
res@cnMaxLevelValF = max(power) ; max level
res@cnLevelSpacingF =
(res@cnMaxLevelValF-res@cnMinLevelValF)/(ncolor+1) ;
contour interval
;================================ RESOURCES SIGNIFICANCE
res2 = True ; res2 probability plots
res2@gsnDraw = False ; Do not draw plot
res2@gsnFrame = False ; Do not advance frome
res2@cnLevelSelectionMode = "ManualLevels" ; set manual contour levels
res2@cnMinLevelValF = 1. ;TRESHOLD ; set min contour level
res2@cnMaxLevelValF = 6. ;max(sig) ; set max contour level
res2@cnLevelSpacingF = 2. ; set contour spacing
res2@cnLinesOn = True ; draw contour lines
res2@cnLineLabelsOn = False ; draw contour labels
res2@cnInfoLabelOn = False
;res2@cnFillScaleF = . ; add extra density
res2@gsnLeftString = ""
res2@gsnRightString = ""
time_axis_labels( power&time, res, resTick )
Wplot(iyear-byear) = gsn_csm_contour(wks,power,res)
iplot = gsn_csm_contour(wks,sig,res2)
overlay(Wplot(iyear-byear),iplot) ; overlay probability plot onto power plot
;===================Resources marking the rainy season
resvl = True
resvl@gsLineColor = 155
resvl@gsLineThicknessF = 2.0 ; twice as thick
XVER11 = flt2dble(power&period)
XVER11 = (/xver11/)
XVER12 = flt2dble(power&period)
XVER12 = (/xver12/)
vl11(iyear-byear) =
gsn_add_polyline(wks,Wplot(iyear-byear),XVER11,power&period,resvl)
vl12(iyear-byear) =
gsn_add_polyline(wks,Wplot(iyear-byear),XVER12,power&period,resvl)
delete(XVER11)
delete(XVER12)
gws = dim_sum_n_Wrap(power,1)
sdof = SDOF
sdof@spcx = gws
Nperiod = dimsizes(power&period)
;************************************************
; calculate confidence interval [here 5 and 95%]
; return 4 curves to be plotted
;************************************************
;printVarSummary(sdof)
;splt = specx_ci (sdof, 0.1, 0.90)
xHigh = chiinv (1.-siglvl, SDOF)/SDOF
splt = new((/3,Nperiod/),typeof(gws)) ; 4 spec curves
splt(0,:) = gws /power&period ; input spectrum
sclfct = sum(gws)/sum(fft_theor)
splt(1,:) = fft_theor * sclfct / power&period ; red noise spectrum
splt(2,:) = splt(1,:)*xHigh ; high ci for red noise spectrum
resl = True
resl@gsnFrame = False
resl@gsnDraw = False
resl@trYAxisType = "LogAxis"
resl@trYReverse = True ; reverse y-axis
resl@tmYLMode = "Explicit"
;resl@tmYLValues = (/2,5,10,15,20,40,60,90/)
resl@tmYLValues = (/10,30,60,90/)
;res@tmYLValues = (/2,5,10,25,40,60,90/)
resl@tmYLLabels = "" + resl@tmYLValues
resl@xyLineThicknesses = (/2.,1.,1./) ; Define line thicknesses
resl@xyDashPatterns = (/0,0,1/) ; Dash patterns
resl@xyLineColors = (/"foreground","green","blue"/)
resl@tiXAxisString = "GWS*frequency"
;plotg = gsn_csm_xy(wks,gws,power&period,resl)
plotg = gsn_csm_xy(wks,splt,power&period,resl)
plotc = gsn_attach_plots(Wplot(iyear-byear),plotg,res,resl)
Nperiod = dimsizes(power&period)
Ntime = dimsizes(power&time)
if (iyear.eq.byear)
powerALL = new((/Nyear,Nperiod,Ntime/),typeof(power))
;sigALL = new((/Nyear,Nperiod,Ntime/),typeof(power))
gwsALL = new((/Nyear,Nperiod/),typeof(gws))
POWERALL = power
;sigALL = sig
GWSALL = gws
SPLTALL = splt
r = res
r2 = res2
rl = resl
rTick = resTick
;rTick@ttmFormat = "%Y%c"
end if
powerALL(iyear-byear,:,:) = power(:,:)
gwsALL(iyear-byear,:) = gws
delete(res)
delete(res2)
delete(resl)
delete(resvl)
delete(resTick)
delete(power)
delete(sig)
delete(gws)
delete(sdof)
delete(splt)
end do ; #################### DO YEAR
pres = True
pres@gsnMaximize = True
;pres@gsnPaperOrientation = "portrait"
;pres@gsnPanelLabelBar = True ; add common colorbar
;pres@lbLabelAutoStride = True ; auto stride on labels
;pres@txString = zone+" lat="+latS+"/"+latN+" lon="+lonL+"/"+lonR
if (VAR.eq."olr") then
gsn_panel(wks,Wplot,(/Nyear,1/),pres)
else
if (VAR.eq."Tb") then
gsn_panel(wks,Wplot,(/Nyear,1/),pres)
else
gsn_panel(wks,Wplot,(/Nyear,1/),pres)
end if
end if
POWERALL = dim_avg_n(powerALL,0)
GWSALL = dim_avg_n(gwsALL,0)
SIGALL = POWERALL
SIGALL = POWERALL/conform (POWERALL,signif,0)
SIGALL@long_name = "Significance"
SIGALL@units = " "
;printVarSummary(SIGALL)
;print(SIGALL(0,:))
SDOFALL = SDOF
SDOFALL@spcx = GWSALL
;************************************************
; calculate confidence interval [here 5 and 95%]
; return 4 curves to be plotted
;************************************************
;printVarSummary(sdof)
;SPLTALL = specx_ci (SDOFALL, 0.1, 0.90)
xHigh = chiinv (1.-siglvl, SDOF)/SDOF
SPLTALL(0,:) = GWSALL / POWERALL&period ; input spectrum
sclfct = sum(GWSALL)/sum(fft_theor)
SPLTALL(1,:) = fft_theor * sclfct / POWERALL&period ; red
noise spectrum
SPLTALL(2,:) = SPLTALL(1,:)*xHigh ; high ci for red noise spectrum
;################ set ressources for zonal climatological plot
#########################
r@gsnRightString = byear+"-"+eyear
r@cnMinLevelValF = min(POWERALL) ; min level
r@cnMaxLevelValF = max(POWERALL) ; max level
r@cnLevelSpacingF =
(r@cnMaxLevelValF-r@cnMinLevelValF)/(ncolor+1) ; contour
interval
time_axis_labels( POWERALL&time, r, rTick )
PWRALLplt = gsn_csm_contour(wks,POWERALL,r)
iplot = gsn_csm_contour(wks,SIGALL,r2)
overlay(PWRALLplt,iplot)
;plotG = gsn_csm_xy(wks,GWSALL,POWERALL&period,rl)
plotG = gsn_csm_xy(wks,SPLTALL,POWERALL&period,rl)
;=================== Resources marking the rainy season for
CLIM WAVELET PLOT
resvl = True
resvl@gsLineColor = 155
resvl@gsLineThicknessF = 2.0 ; twice as thick
xver11 = ut_inv_calendar( byear, 06, 01, 00, 0, 0, timeUnits, 0 )
xver11 = ut_convert( xver11, GTunits )
xver12 = ut_inv_calendar( byear, 09, 30, 00, 0, 0, timeUnits, 0 )
xver12 = ut_convert( xver12, GTunits )
XVER11 = flt2dble(POWERALL&period)
XVER11 = (/xver11/)
XVER12 = flt2dble(POWERALL&period)
XVER12 = (/xver12/)
vl11ALL = gsn_add_polyline(wks,PWRALLplt,XVER11,POWERALL&period,resvl)
vl12ALL = gsn_add_polyline(wks,PWRALLplt,XVER12,POWERALL&period,resvl)
delete(XVER11)
delete(XVER12)
plotC = gsn_attach_plots(PWRALLplt,plotG,r,rl)
delete(r)
delete(rl)
delete(resvl)
delete(rTick)
;if (zone.eq."ZONE3".or.zone.eq."ZONE4".or.zone.eq."ZONE5") then
; gsn_panel(wks,PLTPWRALL(izone-2:izone-2+1),(/1,2/),pres)
;else
gsn_panel(wks,PWRALLplt,(/1,1/),pres)
;end if
delete(pres)
delete(Wplot)
delete(vl11)
delete(vl12)
delete(vl21)
delete(vl22)
delete(X)
delete(w)
delete(fft_theor)
delete(signif)
delete(idxp)
delete(Nidxp)
delete(powersig)
delete(SIG)
delete(powerALL)
delete(POWERALL)
delete(SIGALL)
delete(gwsALL)
delete(GWSALL)
delete(SDOF)
delete(SDOFALL)
delete(SPLTALL)
end if
;delete(X0class)
delete(X1D0)
;delete(olr0)
PPWR = True
PPWR@gsnMaximize = True
gsn_panel(wks,PWRALLplt,(/3,2/),PPWR)
delete(PPWR)
delete(PWRALLplt)
delete(vl11ALL)
delete(vl12ALL)
;delete(vl21ALL)
;delete(vl22ALL)
delete(Gtime)
delete(olr)
;delete(lat)
;delete(lon)
;delete(X1D)
;delete(Xclass)
end do
; #################### DO VAR
end
-- Sinclair ZEBAZE Laboratory for Environmental modelling & Atmospheric Physics Department of Physics, Faculty of science - University of YAOUNDE I P.O.BOX 812 YAOUNDE - CAMEROON Phone: (237) 77834713
_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
