[ncl-talk] ttest on monthly climatologies
Adam Phillips
asphilli at ucar.edu
Wed Jun 20 11:21:48 MDT 2018
Hi Sri,
The error messages you are receiving are very likely due to
equiv_sample_size returning sample sizes less than 3 (as the message
states). It is your choice whether to use equiv_sample_size or not. I
personally always use it, as you generally do not want to include points
with very low sample sizes as part of your t-test plots. Did you examine
your probt array to see what the range of values is? Are there any values
greater than 95%? You can use stat_dispersion
<https://www.ncl.ucar.edu/Document/Functions/Contributed/stat_dispersion.shtml>
to
help you with this. You are only attempting to plot those values higher
than 99%, so there's a chance that no point reaches that threshold.
For the record: I don't see anything amiss with your coding.
If you have any further questions please respond to ncl-talk.
Adam
On Mon, Jun 18, 2018 at 3:27 AM Sri Nandini <snandini at marum.de> wrote:
> 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
> _______________________________________________
> ncl-talk mailing list
> ncl-talk at ucar.edu
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>
--
Adam Phillips
Associate Scientist, Climate and Global Dynamics Laboratory, NCAR
www.cgd.ucar.edu/staff/asphilli/ 303-497-1726
<http://www.cgd.ucar.edu/staff/asphilli>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20180620/4727ba25/attachment.html>
More information about the ncl-talk
mailing list