[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