[ncl-talk] ttest values not plotting
Adam Phillips
asphilli at ucar.edu
Thu Jan 18 10:08:16 MST 2018
Hi Sri,
To clarify my previous response, res2 at cnMonoFillPattern = False also needs
to be set when not using gsn_contour_shade:
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 cnMonoFillPattern = False
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)
Looking at the new modifications to your script, you are doing this for
your t-test calculation:
probt = 100.*(1. -ttest(aveY,varY,xEqv,aveX,varX,yEqv,True,False)) ; all
areas whose value > 95 are significant
As you are multiplying by 100 after subtracting the values from one, your
comment is correct, all values greater than 95 would be considered
significant. However, you are also setting your contour levels like this:
res2 at cnLevels = (/.95/) ; only have 1 contour level
Thus, you are setting your single contour level to be at .95, and not 95.
So you should change that resource to 95.
Also, you are calling gsn_contour_shade like this:
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
Those settings will result in all values less than or equal to .95 (matches
your original contour level of .95) being shaded, which is not what you
want. Those setting should be switched to this:
res2 at cnLevels = (/95./)
opt = True ; set up parameters for pattern fill
opt at gsnShadeFillType = "pattern" ; specify pattern fill
opt at gsnShadeHigh = 17 ; stipple pattern
plot2 = gsn_contour_shade(plot2, 1, 95., opt) ; stipple all areas >= 95
contour
Also note that when you turn stippling on that not every grid box will
necessarily have one stipple in it, and some (depending on your grid) might
have more.
If you have any further questions please respond to ncl-talk.
Adam
On Thu, Jan 18, 2018 at 8:28 AM, Sri Nandini via ncl-talk <ncl-talk at ucar.edu
> wrote:
> 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 <(303)%20497-1726>
>
> <http://www.cgd.ucar.edu/staff/asphilli>
>
>
>
> _______________________________________________
> 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/20180118/495f24d3/attachment.html>
More information about the ncl-talk
mailing list