[ncl-talk] ttest on monthly climatologies
Sri Nandini
snandini at marum.de
Mon Jun 18 03:27:20 MDT 2018
Hello All
I am using the ttest routine on monthly climatology, and plotting only those
values that are significant to the 95% level but the error i get is
warning:ttest: encountered 4692 cases where s1 and/or s2 were less than 2. Output set to missing values in these cases
Upon various searches online i came to understand it means:
This message just indicates that you have spatial locations in your
input arrays which have no non-missing values, and the 'ttest' function
will set the output to missing as well.
Does this mean performing ttest on monthly climatologies or getting the equiv_sample_size is not a good idea?
In the end my plot doesn't contain any stippling dots for ttest values.
my coding looks something
like this:
f= addfile("Total_eva_PI.nc", "r") ;
T41 = f->QFLUX(:,:,:)
T41 at _FillValue = -9.96921e+36 ; (12,lat,lon)
printVarSummary(T41)
aveX = dim_avg_n_Wrap(T41,0) ; (lat,lon)
printVarSummary(aveX)
varX = dim_variance_n_Wrap(T41,0) ; compute variance the 0th dim
printVarSummary(varX) ; (lat,lon)
;=================================================================================
yrStrt = 0950
yrLast = 1000
f1 = addfile("Total_eva_9ka.nc", "r") ;
TIME1 = f1->time
TIME1 = TIME1-31
YYYY1 = cd_calendar(TIME1,-1)/100 ; entire file
iYYYY1 = ind(YYYY1.ge.yrStrt .and. YYYY1.le.yrLast)
air22 = f1->QFLUX(iYYYY1,:,:)
air22 at _FillValue = -9.96921e+36
T2M_4 = clmMonLLT(air22(lat|:,lon|:,time|:)) ; monthly climatology
printVarSummary(T2M_4) ; (12,lat,lon)
aveY = dim_avg_n_Wrap(T2M_4,2)
printVarSummary(aveY) ; (lat,lon)
varY = dim_variance_n_Wrap(T2M_4,2) ; compute variance
printVarSummary(varY) ; (lat,lon)
;================================================================================
;equiv_sample_size: Estimates the number of independent values of a series of correlated observations. Specify a critical significance level to test the lag-one auto-correlation coefficien. Use annual means.
;(Generally, there is no significant year-to-year autocorrelation of monthly data [e.g. successive Januaries].) So here, xX and sY are scalars indicating the number of years used.
;=================================================================================
sigr = 0.05 ; critical sig lvl for r
xEqv = new(dimsizes(aveY),typeof(aveY),aveY at _FillValue)
printVarSummary(xEqv) ; (lat,lon)
xEqv = equiv_sample_size (T41(lat|:,lon|:,time|:), sigr,0)
copy_VarMeta(aveY,xEqv)
printVarSummary(xEqv) ; (lat,lon)
yEqv = new(dimsizes(aveY),typeof(aveY),aveY at _FillValue)
printVarSummary(yEqv)
; (lat,lon)
yEqv = equiv_sample_size (T2M_4(lat|:,lon|:,month|:), sigr,0)
copy_VarMeta(aveY,yEqv)
printVarSummary(yEqv) ; (lat,lon)
;==================================================================================
;Use "ttest" to compute the probabilities
;==================================================================================
probt = new(dimsizes(varY),typeof(varY),varY at _FillValue)
iflag= False ; population variance similar?
probt = 100.*(1. -ttest(aveY,varY,xEqv,aveX,varX,yEqv,iflag,False)) ; all areas whose value > 95 are significant
;printVarSummary(probt) ; (lat,lon)
copy_VarMeta(varY,probt)
printVarSummary(probt) ; (lat,lon)
printMinMax(probt,0) ; min=0.0001486018726404303 max=100
;====================================================================
;calculate the T2M difference
;====================================================================
T2diff = new(dimsizes(aveY),typeof(aveY),aveY at _FillValue)
copy_VarMeta(aveY,T2diff)
T2diff=aveY-aveX
outfile = addfile("E_9_Diff.nc","c") ;
outfile->T1diff = T2diff
;====================================================================
;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","QFLUX_ttest_9ka_PI99%1")
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 = -4 ; set min contour level
res at cnMaxLevelValF = 4 ; set max contour level
res at cnLevelSpacingF = .5 ; set contour spacing
res at gsnDraw = False ; Do not draw plot
res at gsnFrame = False ; Do not advance frome
res at tiMainString = "QFLUX: 9ka-PI"
res at tiMainOffsetYF = 0.02 ; Move title up a little
res at gsnCenterString = " "
res at gsnLeftString = " "
res at gsnRightString = " "
res at cnRasterModeOn = True ; Raster mode shows grid cells
res at lbBoxEndCapStyle = "TriangleBothEnds" ; Added in NCL V6.4.0
res at lbTitlePosition = "bottom" ; title position
res at lbTitleFontHeightF = .01 ; make title smaller
res at lbTitleDirection = "Across" ; title direction
res at lbLabelFontHeightF = 0.01
minlat = 30. ; min lat to mask
maxlat = 80. ; max lat to mask
minlon = -70. ; min lon to mask
maxlon = 63. ; max lon to mask
res at mpProjection = "LambertConformal" ; choose projection
;---masked plot
res at gsnAddCyclic = True ; regional plot
res at mpMinLatF = minlat ; min lat to mask
res at mpMaxLatF = maxlat ; max lat to mask
res at mpMinLonF = minlon ; min lon to mask
res at mpMaxLonF = maxlon ; max lon to mask
res at gsnMaskLambertConformal = True ; turn on lc masking
plot = gsn_csm_contour_map(wks,T2diff,res)
add_lc_labels(wks,plot,minlat,maxlat,minlon,maxlon) ; attach latitude labels
;=====================================================================
;Second plot resources and overlay onto first
;=====================================================================
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 = (/99./) ; only have 1 contour level
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
res2 at gsnRightString = "" ; Turn off subtitles
res2 at gsnLeftString = ""
res at gsnAddCyclic = True ; add cyclic point
plot2 = gsn_csm_contour(wks,probt(:,:),res2)
opt = True ; set up parameters for pattern fill
opt at gsnShadeFillType = "pattern" ; specify pattern fill
opt at gsnShadeHigh = 17 ; stipple pattern
opt at gsnShadeDotSizeF = 1 ; make dots larger
;=====================================Attach the polylines
ADD_SHAPEFILE_OUTLINES = True
sname = "Caspian_Basin.shp"
pres = True
pres at gsLineColor = "black"
pres at gsLineThicknessF = 2.0 ; 3x thickness
shp1 = gsn_add_shapefile_polylines(wks,plot,sname,pres)
printVarSummary(shp1)
;=====================================
plot2 = gsn_contour_shade(plot2,1,99., opt) ; stipple all areas >= 90 contour
overlay (plot, plot2)
draw (plot)
frame(wks)
end
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20180618/60ba67eb/attachment.html>
More information about the ncl-talk
mailing list