[ncl-talk] Extract upper high tail of probability distribution
Sri.durgesh Nandini-Weiss
sri.durgesh.nandini-weiss at uni-hamburg.de
Mon Dec 9 07:05:34 MST 2019
Hello everyone,
I wanted to show 2 plots of the distribution of my SSH anomaly data.
The first plot (Hist_RCP45_Vancouver_PDF_zoslr.pdf) shows the PDF of
Vancouver, with the 5th and 95th percentiles as vertical lines.
I have problems with the second
plot(Hist_RCP45_Vancouver_95perc_zoslr.pdf), where i just want to
extract the high end/upper tail and plot this (between 95-100%), but it
doesn't show the correct data.
Below is my script for calculating and plotting this.
Would appreciate any advice.
f = addfile ("anom_hist_zo_slr.nc", "r")
hist_anom = f->hist_anom
printVarSummary(hist_anom)
f1 = addfile ("anom_rcp45_zo_slr.nc", "r")
rcp45_anom = f1->rcp45_anom
printVarSummary(rcp45_anom)
dimx = dimsizes(hist_anom)
ntim = dimx(0) ; 240
nens = dimx(1) ; 100
nlat = dimx(2) ; 45
nlon = dimx(3) ; 90
lat1=nlat
lon1=nlon
nmos = 12
nyrs = ntim/nmos ; 20
printVarSummary(nyrs)
;==================================================================
print ("Data read in - start calculations")
;==================================================================
; Region domain lat and lon
;==================================================================
hist_anom1=hist_anom(:,:,{48},{-125}) ;(ntim, mlon)
printVarSummary(hist_anom1)
rcp45_anom1=rcp45_anom(:,:,{48},{-125}) ;(ntim, mlon)
printVarSummary(rcp45_anom1)
;=================================================================
; Do dim_avg and dim_stddev on time and ens dimension, only take
skewness and kurtosis from stat4
;=================================================================
anom1_1d = ndtooned(hist_anom1) ;to carry out the
dim_stat4_n on one dimension.
printVarSummary(anom1_1d)
aveX = avg (anom1_1d)
printMinMax(aveX,0)
varstdX=stddev(anom1_1d)
printVarSummary(varstdX)
anom2_1d = ndtooned(rcp45_anom1)
printVarSummary(anom2_1d)
aveY = avg (anom2_1d)
printMinMax(aveY,0)
varstdY=stddev(anom2_1d)
printVarSummary(varstdY)
;================================================================================
; calculate 95-100 percentiles: extract mean+2standard devitations
;===============================================================================
; percentiles (requires sorting the array)
x95= aveX+2*varstdX
;no = dim_pqsort(anom1_1d,2) ;; sorted array is
anom1_nomiss!!
i95 = ind(anom1_1d.ge.x95)
xHigh = anom1_1d(i95) ; values
printVarSummary(xHigh)
y95= aveY+2*varstdY
;it = dim_pqsort(anom2_1d,2) ; sorted array is
anom1_nomiss!!
i95Y = ind(anom2_1d.ge.y95)
YHigh = anom2_1d(i95Y) ; values
printVarSummary(YHigh)
;==================================================================
TTstat1 = dim_stat4(xHigh)
printVarSummary(TTstat1)
TTstat2 = dim_stat4(YHigh)
printVarSummary(TTstat2)
print("calculate PDFs=========================")
;================================================================================
; calculate empirical PDF
;================================================================================
opt = True
opt at bin_min = -25.
opt at bin_max = 25.
pdf1 = pdfx(xHigh, 100, opt)
printVarSummary(pdf1)
pdf2 = pdfx(YHigh, 100, opt)
printVarSummary(pdf2)
;================================================================================
; prepare data for plotting
;================================================================================
nVar = 2
nBin = pdf1 at nbins ; retrieve the number of bins
xx = new ( (/nVar, nBin/), typeof(pdf1))
xx(0,:) = pdf1 at bin_center ; assign appropriate "x" axis values
xx(1,:) = pdf2 at bin_center
yy = new ( (/nVar, nBin/), typeof(pdf1))
yy(0,:) = (/ pdf1 /) * 0.01 ;; convert % to absolut
yy(1,:) = (/ pdf2 /) * 0.01
;================================================================================
; Plot PDFs...seperate plots
;================================================================================
plotfile="Hist_RCP45_Vancover_95perc_zoslr"
wks_type = "pdf"
wks = gsn_open_wks(wks_type,plotfile)
res = True
res at gsnMaximize =True
res at vpHeightF = 0.7 ; change aspect ratio of plot
res at vpWidthF = 1.
res at tmYLAutoPrecision = False ; no auto precision
res at tmYLPrecision = 2 ; set the precision
res at xyLineThicknesses = (/5.,5./)
res at xyDashPattern = 0 ; dashed/solid(=0)
res at trXMinF = 0 ; set minimum X-axis value
res at trXMaxF = 15 ; set maximum X-axis value
res at trYMaxF = 0.11 ; set maximum Y-axis value
res at tiXAxisFontHeightF = 0.022
res at tiYAxisFontHeightF = 0.022
res at tiYAxisString = "probability %" ; axis string
res at tiXAxisString = "SSH anomaly [cm]" ; axis string
res at xyLineColors = (/"blue","red"/)
res at gsnXRefLineDashPattern =1
res at gsnXRefLine = 0.0
res at gsnDraw = False ; don't draw so text can
be added
res at gsnFrame = False
res at gsnStringFontHeightF = 0.026
res at gsnCenterString = "City: "+region +" (Lat:" +lat1 +",
Lon:"+lon1 +")"
plot = gsn_csm_xy (wks, xx, yy , res)
;================================================================================
; add text
;================================================================================
txres = True
;txres at txFont = "duplex_roman"
txres at txJust = "CenterLeft"
txres at txFontHeightF = 0.019
;txres at txFuncCode = "~"
;; with axis lables
txres at txFontColor = "blue"
gsn_text_ndc(wks,"1986-2006",0.19,0.645,txres)
gsn_text_ndc (wks,"Mean = "+decimalPlaces(TTstat1(0),3,True),
0.192, 0.61,txres)
gsn_text_ndc (wks,"Variance = "+decimalPlaces(TTstat1(1),3,True),
0.192, 0.575 ,txres)
gsn_text_ndc (wks,"Skewness = "+decimalPlaces(TTstat1(2),3,True),
0.192, 0.54 ,txres)
gsn_text_ndc (wks,"Kurtosis = "+decimalPlaces(TTstat1(3),3,True),
0.192, 0.505 ,txres)
txres at txFontColor = "red"
gsn_text_ndc(wks,"RCP4.5: 1980-2100",0.698,0.645,txres)
gsn_text_ndc (wks,"Mean = "+decimalPlaces(TTstat2(0),3,True),
0.70, 0.61,txres)
gsn_text_ndc (wks,"Variance = "+decimalPlaces(TTstat2(1),3,True),
0.70, 0.575 ,txres)
gsn_text_ndc (wks,"Skewness = "+decimalPlaces(TTstat2(2),3,True),
0.70, 0.54 ,txres)
gsn_text_ndc (wks,"Kurtosis = "+decimalPlaces(TTstat2(3),3,True),
0.70, 0.505 ,txres)
draw (plot)
frame (wks)
On 12/5/19 5:44 PM, Dennis Shea wrote:
> If the source array is 'x[*]' and the distribution in normally
> distributed with xMean and xStd, the.
>
> x95 = xMean+2*xStd ; scalar
> i95 = ind(x.ge.x95)
> xHigh = x(i95) ; values
> print(xHigh)
> ===
> If 'x' is multidimensional [ 'X' ] then
>
> x = ndtooned(X)
>
>
>
>
> On Thu, Dec 5, 2019 at 6:53 AM Sri.durgesh Nandini-Weiss via ncl-talk
> <ncl-talk at ucar.edu <mailto:ncl-talk at ucar.edu>> wrote:
>
> Hello everyone,
>
> I am trying to extract the upper high end tail (95
> percentiles-100percentiles) (shaded area below) of the probability
> distribution.
>
> Is one way: Area95-100percentiles= mean+2standard deviations?
>
> Does anyone know how to extract (resampled) and plot this high
> end/upper tail from a normal distribution?
>
>
>
> --
> Dr. Sri Nandini-Weiß
> Research Scientist
>
> Universität Hamburg
> Center for Earth System Research and Sustainability (CEN)
> Cluster of Excellence 'Climate, Climatic Change, and Society' (CLICCS)
>
> Bundesstrasse 53, 20146 Hamburg
> Tel: +49 (0) 40 42838 7472
>
> _______________________________________________
> ncl-talk mailing list
> ncl-talk at ucar.edu <mailto:ncl-talk at ucar.edu>
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>
--
Dr. Sri Nandini-Weiß
Research Scientist
Universität Hamburg
Center for Earth System Research and Sustainability (CEN)
Cluster of Excellence 'Climate, Climatic Change, and Society' (CLICCS)
Bundesstrasse 53, 20146 Hamburg
Tel: +49 (0) 40 42838 7472
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20191209/73ea87f8/attachment-0001.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: igmpcfpjcpojomll.jpeg
Type: image/jpeg
Size: 43896 bytes
Desc: not available
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20191209/73ea87f8/attachment-0001.jpeg>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: Hist_RCP45_Vancouver_PDF_zoslr.pdf
Type: application/pdf
Size: 91281 bytes
Desc: not available
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20191209/73ea87f8/attachment-0002.pdf>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: Hist_RCP45_Vancover_95perc_zoslr.pdf
Type: application/pdf
Size: 83519 bytes
Desc: not available
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20191209/73ea87f8/attachment-0003.pdf>
More information about the ncl-talk
mailing list