[ncl-talk] ttest values not plotting
Sri Nandini
snandini at marum.de
Wed Jan 17 13:18:44 MST 2018
Hello
I am trying to plot ttest values overlay onto the T2m differences, but it is not being plotted nor shaded region according to 95%confidence limits with ttest.
Below is the script, i would be grateful if someone could let me know of my omission or deletion in the resources plotting attributes.
;Use NCL's named dimensions to reorder in time.
;Calculate the temporal means and variances using the dim_avg and dim_variance functions.
;Specify a critical significance level to test the lag-one auto-correlation coefficient and determine the (temporal) number of equivalent sample sizes in each grid point using equiv_sample_size.
;Specify a critical significance level for the ttest and test if the means are different at each grid point.
; This file still has to be loaded manually
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
yrStrt = 0000
yrLast = 0100
; ==============================================================
; Open the file: Read only the user specified period
; ==============================================================
f= addfile("T2M_C.nc", "r") ;Model Control
TIME = f->time
YYYY = cd_calendar(TIME,-1)/100 ; entire file
iYYYY = ind(YYYY.ge.yrStrt .and. YYYY.le.yrLast)
T41 = f->TREFHT(iYYYY,:,:)
printVarSummary(T41) ; (time, lat,lon)
T4 = lonFlip(T41)
printVarSummary(T4) ; (time, lat,lon)
T4 at _FillValue = -9.96921e+36
aveX = dim_avg_n_Wrap(T4, 0) ;; average over the 0th dim
printVarSummary(aveX) ; (lat,lon)
varX = dim_variance_n_Wrap(T4,0) ; compute variance
printVarSummary(varX) ; (lat,lon)
;==============================================================================
f2= addfile("T2M_CC.nc", "r") ;Current Caspian
TIME2 = f2->time
YYYY2 = cd_calendar(TIME2,-1)/100 ; entire file
iYYYY2 = ind(YYYY2.ge.yrStrt .and. YYYY2.le.yrLast)
air22 = f2->TREFHT(iYYYY2,:,:)
air2 = lonFlip(air22)
printVarSummary(air2) ; (time, lat,lon)
air2 at _FillValue = -9.96921e+36
aveY = dim_avg_n_Wrap(air2,0) ; average over the 0th dim
printVarSummary(aveY) ; (lat,lon)
varY = dim_variance_n_Wrap(air2,0) ; compute variance
printVarSummary(varY) ; (lat,lon)
;================================================================================
;Use "ttest" to compute the probabilities
;=================================================================================
sigr = 0.05 ; critical sig lvl for r
xEqv = new(dimsizes(varY),typeof(varY),varY at _FillValue)
printVarSummary(xEqv) ; (lat,lon)
xEqv = equiv_sample_size (T4(lat|:,lon|:,time|:), sigr,0)
copy_VarMeta(varY,xEqv)
printVarSummary(xEqv) ; (lat,lon)
yEqv = new(dimsizes(varY),typeof(varY),varY at _FillValue)
printVarSummary(yEqv) ; (lat,lon)
yEqv = equiv_sample_size (air2(lat|:,lon|:,time|:), sigr,0)
copy_VarMeta(varY,yEqv)
printVarSummary(yEqv) ; (lat,lon)
;==================================================================================
probt = new(dimsizes(varY),typeof(varY),varY at _FillValue)
probt = ttest(aveY,varY,xEqv,aveX,varX,yEqv,True,False) ; values are between 0 and 1
;probt = 100.*(1. -ttest(aveY,varY,xEqv,aveX,varX,yEqv,True,False)) ; all areas whose value > 95 are significant
printVarSummary(probt) ; (lat,lon)
copy_VarMeta(varY,probt)
printVarSummary(probt) ; (lat,lon)
print(probt) ;check to make sure values are between 0 and 1.
;============================================================
;calculate the T2M difference
;============================================================
T2diff = new(dimsizes(aveY),typeof(aveY),aveY at _FillValue)
printVarSummary(T2diff)
copy_VarMeta(aveY,T2diff)
printVarSummary(T2diff)
T2diff=aveY-aveX
printVarSummary(T2diff) ; (lat,lon)
;====================================================================
;Here is a section of my code that draws a color fill plot, and then
;overlays contours of statistical significance.
;====================================================================
wks = gsn_open_wks("pdf","T2M_ttest")
gsn_define_colormap(wks,"BlWhRe")
res = True
res at cnFillOn = True ; turn on color
res at gsnSpreadColors = True ; use full colormap
res at cnLinesOn = False ; turn off contour lines
res at cnLevelSelectionMode = "ManualLevels" ; set manual contour levels
res at cnMinLevelValF = -2. ; set min contour level
res at cnMaxLevelValF = 2. ; set max contour level
res at cnLevelSpacingF = 0.2 ; set contour spacing
res at gsnDraw = False ; Do not draw plot
res at gsnFrame = False ; Do not advance frome
res at tiMainString = "T2M Difference: 50yrs overlay with ttest"
res at gsnCenterString = "5% stippled"
res at gsnLeftString = "K"
res at cnRasterModeOn = True ; Raster mode shows grid cells
plot = gsn_csm_contour_map_ce(wks,T2diff,res)
plot = ZeroNegDashLineContour (plot)
res2 = True ; res2 probability plots
res2 at gsnDraw = False ; Do not draw plot
res2 at gsnFrame = False ; Do not advance frome
res2 at cnLevelSelectionMode = "ExplicitLevels" ; set explicit cnlev
;res2 at cnLevels = (/.95/) ; only have 1 contour level
;res2 at cnFillPatterns = (/-1,17/) ; don't fill <0.95, stipple >=0.95
res at cnFillOn = False
res2 at cnInfoLabelOn = False
res2 at cnLinesOn = False ; do not draw contour lines
res2 at cnLineLabelsOn = False ; do not draw contour labels
res2 at cnFillScaleF = 0.6 ; add extra density
plot2 = gsn_csm_contour(wks,probt,res2) ; add cyclic point
opt = True ; set up parameters for pattern fill
opt at gsnShadeFillType = "pattern" ; specify pattern fill
opt at gsnShadeLow = 17 ; stipple pattern
plot2 = gsn_contour_shade(plot2, 0.091, 3., opt) ; stipple all areas < 0.09 contour
overlay (plot, plot2)
draw (plot)
frame(wks)
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20180117/cae2bb94/attachment.html>
More information about the ncl-talk
mailing list