<html><head><meta http-equiv="Content-Type" content="text/html; charset=UTF-8"><title></title><style type="text/css">.felamimail-body-blockquote {margin: 5px 10px 0 3px;padding-left: 10px;border-left: 2px solid #000088;} </style></head><body>Hello<br><br>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.<br>Below is the script, i would be grateful if someone could let me know of my omission or deletion in the resources plotting attributes.<br><br>;Use NCL's named dimensions to reorder in time.<br>;Calculate the temporal means and variances using the dim_avg and dim_variance functions.<br>;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.<br>;Specify a critical significance level for the ttest and test if the means are different at each grid point. <br>; This file still has to be loaded manually<br>load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"<br><br> yrStrt = 0000<br> yrLast = 0100<br>; ==============================================================<br>; Open the file: Read only the user specified period <br>; ==============================================================<br>f= addfile("T2M_C.nc", "r") ;Model Control<br>TIME = f->time<br> YYYY = cd_calendar(TIME,-1)/100 ; entire file<br> iYYYY = ind(YYYY.ge.yrStrt .and. YYYY.le.yrLast)<br>T41 = f->TREFHT(iYYYY,:,:)<br>printVarSummary(T41) ; (time, lat,lon)<br>T4 = lonFlip(T41)<br>printVarSummary(T4) ; (time, lat,lon)<br>T4@_FillValue = -9.96921e+36 <br>aveX = dim_avg_n_Wrap(T4, 0) ;; average over the 0th dim<br>printVarSummary(aveX) ; (lat,lon) <br>varX = dim_variance_n_Wrap(T4,0) ; compute variance<br>printVarSummary(varX) ; (lat,lon) <br>;==============================================================================<br><br>f2= addfile("T2M_CC.nc", "r") ;Current Caspian<br>TIME2 = f2->time<br> YYYY2 = cd_calendar(TIME2,-1)/100 ; entire file<br> iYYYY2 = ind(YYYY2.ge.yrStrt .and. YYYY2.le.yrLast)<br>air22 = f2->TREFHT(iYYYY2,:,:)<br>air2 = lonFlip(air22)<br>printVarSummary(air2) ; (time, lat,lon)<br>air2@_FillValue = -9.96921e+36 <br>aveY = dim_avg_n_Wrap(air2,0) ; average over the 0th dim<br>printVarSummary(aveY) ; (lat,lon) <br>varY = dim_variance_n_Wrap(air2,0) ; compute variance<br>printVarSummary(varY) ; (lat,lon) <br>;================================================================================<br>;Use "ttest" to compute the probabilities<br>;=================================================================================<br><br> sigr = 0.05 ; critical sig lvl for r<br><br> xEqv = new(dimsizes(varY),typeof(varY),varY@_FillValue)<br>printVarSummary(xEqv) ; (lat,lon) <br> xEqv = equiv_sample_size (T4(lat|:,lon|:,time|:), sigr,0) <br> copy_VarMeta(varY,xEqv)<br>printVarSummary(xEqv) ; (lat,lon) <br><br> yEqv = new(dimsizes(varY),typeof(varY),varY@_FillValue)<br>printVarSummary(yEqv) ; (lat,lon) <br> yEqv = equiv_sample_size (air2(lat|:,lon|:,time|:), sigr,0)<br>copy_VarMeta(varY,yEqv)<br>printVarSummary(yEqv) ; (lat,lon) <br> <br>;==================================================================================<br> probt = new(dimsizes(varY),typeof(varY),varY@_FillValue)<br><br> probt = ttest(aveY,varY,xEqv,aveX,varX,yEqv,True,False) ; values are between 0 and 1<br><br>;probt = 100.*(1. -ttest(aveY,varY,xEqv,aveX,varX,yEqv,True,False)) ; all areas whose value > 95 are significant<br>printVarSummary(probt) ; (lat,lon) <br>copy_VarMeta(varY,probt)<br>printVarSummary(probt) ; (lat,lon) <br>print(probt) ;check to make sure values are between 0 and 1.<br>;============================================================<br>;calculate the T2M difference <br>;============================================================ <br> T2diff = new(dimsizes(aveY),typeof(aveY),aveY@_FillValue)<br>printVarSummary(T2diff) <br>copy_VarMeta(aveY,T2diff)<br>printVarSummary(T2diff) <br>T2diff=aveY-aveX<br> printVarSummary(T2diff) ; (lat,lon) <br>;====================================================================<br>;Here is a section of my code that draws a color fill plot, and then<br>;overlays contours of statistical significance.<br>;====================================================================<br> wks = gsn_open_wks("pdf","T2M_ttest")<br> gsn_define_colormap(wks,"BlWhRe")<br><br> res = True<br> res@cnFillOn = True ; turn on color<br> res@gsnSpreadColors = True ; use full colormap<br><br> res@cnLinesOn = False ; turn off contour lines<br> res@cnLevelSelectionMode = "ManualLevels" ; set manual contour levels<br> res@cnMinLevelValF = -2. ; set min contour level<br> res@cnMaxLevelValF = 2. ; set max contour level<br> res@cnLevelSpacingF = 0.2 ; set contour spacing<br><br> res@gsnDraw = False ; Do not draw plot<br> res@gsnFrame = False ; Do not advance frome<br><br> res@tiMainString = "T2M Difference: 50yrs overlay with ttest"<br> res@gsnCenterString = "5% stippled"<br> res@gsnLeftString = "K"<br><br> res@cnRasterModeOn = True ; Raster mode shows grid cells<br> plot = gsn_csm_contour_map_ce(wks,T2diff,res) <br> plot = ZeroNegDashLineContour (plot)<br><br> res2 = True ; res2 probability plots<br> res2@gsnDraw = False ; Do not draw plot<br> res2@gsnFrame = False ; Do not advance frome<br> res2@cnLevelSelectionMode = "ExplicitLevels" ; set explicit cnlev<br> ;res2@cnLevels = (/.95/) ; only have 1 contour level<br> ;res2@cnFillPatterns = (/-1,17/) ; don't fill <0.95, stipple >=0.95<br> res@cnFillOn = False<br><br> res2@cnInfoLabelOn = False<br> res2@cnLinesOn = False ; do not draw contour lines<br> res2@cnLineLabelsOn = False ; do not draw contour labels<br> res2@cnFillScaleF = 0.6 ; add extra density<br><br> plot2 = gsn_csm_contour(wks,probt,res2) ; add cyclic point<br><br> opt = True ; set up parameters for pattern fill<br> opt@gsnShadeFillType = "pattern" ; specify pattern fill<br> opt@gsnShadeLow = 17 ; stipple pattern<br> plot2 = gsn_contour_shade(plot2, 0.091, 3., opt) ; stipple all areas < 0.09 contour<br> overlay (plot, plot2)<br><br> draw (plot)<br> frame(wks)<br></body></html>