[ncl-talk] ttest values not plotting

Sri Nandini snandini at marum.de
Thu Jan 18 08:28:33 MST 2018


Hello

It appears now the stippling is being added to the map, see pdf attached, however the the latlon seem to be shifted, am i wrong about this?
Here is my modified script:


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  plot: ZeroNegDashLineContour and ShadeLtContour
;           are in shea_util.ncl. 
;=================================================================================
  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) 

 printMinMax(probt,0)                    ;  min=0   max=97.2842
;============================================================
;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_tt")
   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 cnFillMode       = "RasterFill"       ; Raster Mode
   res at cnRasterModeOn              = True              ; Raster mode shows grid cells
   plot = gsn_csm_contour_map(wks,T2diff,res)  
   plot = ZeroNegDashLineContour (plot)            ; a shea_util.ncl function

 aa=stat_dispersion(probt,True)   ; Make sure your data has values less than  or equal to 0.95 that can be plotted
    
  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 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
                                         ; add cyclic point
  ;res at gsnAddCyclic                = True              ; add cyclic point

  plot2   = gsn_csm_contour(wks,gsn_add_cyclic_point(probt(:,:)), res2) 
  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.95, 3., opt) ; stipple all areas < 95 contour                               
  overlay (plot, plot2)

  draw (plot)
  frame(wks)

;=====================================================================
;T2diff contains the temperature difference, probt contains the ttest
;values. 
;=====================================================================






On Jan 17, 2018 9:36:03 PM, Adam Phillips wrote:

> Hi Sri,> I think the issue is that you are letting NCL create the contour levels (for plot2), as you are setting > res2 at cnLevelSelectionMode = "ExplicitLevels" but you are > not setting res2 at cnLevels. 
> > gsn_contour_shade uses the existing contour levels to modify the plot, it does not add them itself. I would suggest doing this:
> 
> 
>   stat_dispersion(probt,True)   ; Make sure your data has values less than  or equal to .05 that can be plotted
> >                                             ; > https://www.ncl.ucar.edu/Document/Functions/Contributed/stat_dispersion.shtml
> 
>   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   = (/.05/)    ; only have 1 contour level
>   ;res2 at cnFillPatterns = (/-1,17/) ; don't fill <0.95, stipple >=0.95
>   res at cnFillOn = True
> 
>   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.05, 3., opt) ; stipple all areas <= 0.05 contour
>   overlay (plot, plot2)


> > or, more simply, do not even use gsn_contour_shade:
> > stat_dispersion(probt,True)   ; Make sure your data has values less than  or equal to .05 that can be plotted
> >                                             ; > https://www.ncl.ucar.edu/Document/Functions/Contributed/stat_dispersion.shtml
> 
> >   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   = (/.05/)    ; only have 1 contour level
>   res2 at cnFillPatterns = (/17,-1/) ; don't fill <0.95, stipple >=0.95
>   res at cnFillOn = True
> 
>   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
>   > overlay (plot, plot2)
> 

> > Hope that helps. If you have any further questions please respond to the ncl-talk email list.
> > Adam
> 
> 
> On Wed, Jan 17, 2018 at 1:18 PM, Sri Nandini via ncl-talk > <> ncl-talk at ucar.edu> >>  wrote:
> > 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)
> > 

> > _______________________________________________
> > 
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 
> 

> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> > 
> 



-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20180118/26db6d69/attachment-0001.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: T2M_tt.pdf
Type: application/pdf
Size: 171332 bytes
Desc: not available
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20180118/26db6d69/attachment-0001.pdf>


More information about the ncl-talk mailing list