<html>
  <head>
    <meta http-equiv="Content-Type" content="text/html; charset=UTF-8">
  </head>
  <body>
    <p>Hello everyone,</p>
    <p>I wanted to show 2 plots of the distribution of my SSH anomaly
      data.</p>
    <p>The first plot (Hist_RCP45_Vancouver_PDF_zoslr.pdf) shows the 
      PDF of Vancouver, with the 5th and 95th percentiles as vertical
      lines.</p>
    <p>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. <br>
    </p>
    <p>Below is my script for calculating and plotting this.</p>
    <p>Would appreciate any advice.</p>
    <p><br>
         f     = addfile ("anom_hist_zo_slr.nc", "r")<br>
         hist_anom    = f->hist_anom<br>
         printVarSummary(hist_anom)<br>
      <br>
         f1     = addfile ("anom_rcp45_zo_slr.nc", "r")<br>
         rcp45_anom   = f1->rcp45_anom<br>
         printVarSummary(rcp45_anom)<br>
      <br>
         dimx = dimsizes(hist_anom)<br>
         ntim = dimx(0)           ; 240<br>
         nens = dimx(1)           ; 100<br>
         nlat = dimx(2)           ; 45<br>
         nlon = dimx(3)           ; 90<br>
         lat1=nlat<br>
         lon1=nlon<br>
         nmos = 12<br>
         nyrs   = ntim/nmos       ; 20<br>
         printVarSummary(nyrs)<br>
      <br>
      <br>
;==================================================================<br>
         print ("Data read in - start calculations")<br>
      <br>
;==================================================================<br>
      ; Region domain lat and lon<br>
;==================================================================<br>
      <br>
        hist_anom1=hist_anom(:,:,{48},{-125})               ;(ntim,
      mlon) <br>
        printVarSummary(hist_anom1)<br>
        rcp45_anom1=rcp45_anom(:,:,{48},{-125})              ;(ntim,
      mlon)<br>
        printVarSummary(rcp45_anom1)<br>
      <br>
      ;=================================================================<br>
      ; Do dim_avg and dim_stddev on time and ens dimension, only take
      skewness and kurtosis from stat4<br>
      ;=================================================================<br>
      <br>
        anom1_1d = ndtooned(hist_anom1)                     ;to carry
      out the dim_stat4_n on one dimension.<br>
        printVarSummary(anom1_1d)<br>
      <br>
        aveX = avg (anom1_1d)<br>
        printMinMax(aveX,0)                                           <br>
        varstdX=stddev(anom1_1d)<br>
        printVarSummary(varstdX)                          <br>
      <br>
        anom2_1d = ndtooned(rcp45_anom1)<br>
        printVarSummary(anom2_1d)<br>
      <br>
        aveY = avg (anom2_1d)<br>
        printMinMax(aveY,0)<br>
        varstdY=stddev(anom2_1d)<br>
        printVarSummary(varstdY)      <br>
      <br>
;================================================================================<br>
      ; calculate 95-100 percentiles: extract mean+2standard devitations<br>
;===============================================================================<br>
       ; percentiles (requires sorting the array)<br>
      <br>
        x95= aveX+2*varstdX  <br>
        ;no = dim_pqsort(anom1_1d,2)                      ;; sorted
      array is anom1_nomiss!!<br>
        i95  =  ind(anom1_1d.ge.x95)<br>
        xHigh = anom1_1d(i95)                 ; values<br>
        printVarSummary(xHigh)    <br>
        <br>
        y95= aveY+2*varstdY                          <br>
        ;it = dim_pqsort(anom2_1d,2)                      ; sorted array
      is anom1_nomiss!!  <br>
        i95Y  =  ind(anom2_1d.ge.y95)<br>
        YHigh =   anom2_1d(i95Y)              ; values<br>
        printVarSummary(YHigh)   <br>
      <br>
;==================================================================<br>
      <br>
        TTstat1 = dim_stat4(xHigh)   <br>
        printVarSummary(TTstat1)                           <br>
      <br>
        TTstat2 = dim_stat4(YHigh)   <br>
        printVarSummary(TTstat2)                            <br>
        print("calculate PDFs=========================")<br>
          <br>
;================================================================================<br>
      ; calculate empirical PDF<br>
;================================================================================<br>
        opt          = True<br>
        opt@bin_min = -25.<br>
        opt@bin_max =  25.<br>
      <br>
        pdf1 = pdfx(xHigh, 100, opt)<br>
        printVarSummary(pdf1) <br>
        pdf2 = pdfx(YHigh, 100, opt)<br>
        printVarSummary(pdf2) <br>
      <br>
;================================================================================<br>
      ; prepare data for plotting<br>
;================================================================================<br>
      <br>
        nVar    = 2<br>
        nBin    = pdf1@nbins          ; retrieve the number of bins<br>
      <br>
        xx      = new ( (/nVar, nBin/), typeof(pdf1))<br>
        xx(0,:) = pdf1@bin_center     ; assign appropriate "x" axis
      values<br>
        xx(1,:) = pdf2@bin_center<br>
      <br>
        yy      = new ( (/nVar, nBin/), typeof(pdf1))<br>
        yy(0,:) = (/ pdf1 /) * 0.01   ;; convert % to absolut<br>
        yy(1,:) = (/ pdf2 /) * 0.01<br>
      <br>
;================================================================================<br>
      ; Plot PDFs...seperate plots<br>
;================================================================================<br>
      <br>
        plotfile="Hist_RCP45_Vancover_95perc_zoslr"<br>
        wks_type = "pdf"<br>
        wks = gsn_open_wks(wks_type,plotfile)<br>
      <br>
        res                   = True<br>
        res@gsnMaximize       =True<br>
        res@vpHeightF         = 0.7                ; change aspect ratio
      of plot<br>
        res@vpWidthF          = 1.<br>
      <br>
        res@tmYLAutoPrecision    = False          ; no auto precision<br>
        res@tmYLPrecision        = 2              ; set the precision<br>
      <br>
        res@xyLineThicknesses = (/5.,5./)<br>
        res@xyDashPattern = 0               ; dashed/solid(=0)<br>
        res@trXMinF = 0             ; set minimum X-axis value<br>
        res@trXMaxF = 15             ; set maximum X-axis value<br>
        res@trYMaxF = 0.11              ; set maximum Y-axis value<br>
      <br>
        res@tiXAxisFontHeightF = 0.022<br>
        res@tiYAxisFontHeightF = 0.022<br>
        res@tiYAxisString = "probability %"         ; axis string<br>
        res@tiXAxisString = "SSH anomaly [cm]"      ; axis string<br>
      <br>
        res@xyLineColors      = (/"blue","red"/)<br>
        res@gsnXRefLineDashPattern =1<br>
        res@gsnXRefLine           = 0.0<br>
      <br>
        res@gsnDraw           = False              ; don't draw so text
      can be added  <br>
        res@gsnFrame          = False  <br>
        res@gsnStringFontHeightF = 0.026 <br>
        res@gsnCenterString        = "City: "+region +" (Lat:" +lat1 +",
      Lon:"+lon1 +")"<br>
        plot = gsn_csm_xy (wks, xx, yy , res)<br>
      <br>
;================================================================================<br>
      ;  add text<br>
;================================================================================<br>
      <br>
         txres               = True             <br>
         ;txres@txFont        = "duplex_roman" <br>
         txres@txJust        = "CenterLeft"<br>
         txres@txFontHeightF = 0.019<br>
         ;txres@txFuncCode    = "~"<br>
      <br>
      ;; with axis lables<br>
      <br>
         txres@txFontColor = "blue"<br>
         gsn_text_ndc(wks,"1986-2006",0.19,0.645,txres) <br>
         gsn_text_ndc (wks,"Mean     =
      "+decimalPlaces(TTstat1(0),3,True), 0.192, 0.61,txres)<br>
         gsn_text_ndc (wks,"Variance =
      "+decimalPlaces(TTstat1(1),3,True), 0.192, 0.575 ,txres)<br>
         gsn_text_ndc (wks,"Skewness =
      "+decimalPlaces(TTstat1(2),3,True), 0.192, 0.54 ,txres)<br>
         gsn_text_ndc (wks,"Kurtosis =
      "+decimalPlaces(TTstat1(3),3,True), 0.192, 0.505 ,txres)<br>
      <br>
         txres@txFontColor = "red"<br>
         gsn_text_ndc(wks,"RCP4.5: 1980-2100",0.698,0.645,txres) <br>
         gsn_text_ndc (wks,"Mean     =
      "+decimalPlaces(TTstat2(0),3,True), 0.70, 0.61,txres)<br>
         gsn_text_ndc (wks,"Variance =
      "+decimalPlaces(TTstat2(1),3,True), 0.70, 0.575 ,txres)<br>
         gsn_text_ndc (wks,"Skewness =
      "+decimalPlaces(TTstat2(2),3,True), 0.70, 0.54 ,txres)<br>
         gsn_text_ndc (wks,"Kurtosis =
      "+decimalPlaces(TTstat2(3),3,True), 0.70, 0.505 ,txres)<br>
      <br>
         draw (plot)<br>
         frame (wks)<br>
      <br>
      <br>
       <br>
    </p>
    <p><br>
    </p>
    <p><br>
    </p>
    <p><br>
    </p>
    <div class="moz-cite-prefix">On 12/5/19 5:44 PM, Dennis Shea wrote:<br>
    </div>
    <blockquote type="cite"
cite="mid:CAOF1d_5E_a-o-92=-pbn+1CfrqRCgO0x2_FKx5djAyNbib6bRA@mail.gmail.com">
      <meta http-equiv="content-type" content="text/html; charset=UTF-8">
      <div dir="ltr">
        <div>If the source array is 'x[*]' and the distribution in
          normally distributed with xMean and xStd, the.</div>
        <div><br>
        </div>
        <div>   x95 = xMean+2*xStd   ; scalar<br>
        </div>
        <div>   i95  =  ind(x.ge.x95)</div>
        <div>   xHigh = x(i95)              ; values <br>
        </div>
        <div>   print(xHigh)<br>
        </div>
        <div>===</div>
        <div>If 'x' is multidimensional [ 'X' ] then</div>
        <div><br>
        </div>
        <div>   x = ndtooned(X)</div>
        <div><br>
        </div>
        <div><br>
        </div>
        <div><br>
        </div>
      </div>
      <br>
      <div class="gmail_quote">
        <div dir="ltr" class="gmail_attr">On Thu, Dec 5, 2019 at 6:53 AM
          Sri.durgesh Nandini-Weiss via ncl-talk <<a
            href="mailto:ncl-talk@ucar.edu" moz-do-not-send="true">ncl-talk@ucar.edu</a>>
          wrote:<br>
        </div>
        <blockquote class="gmail_quote" style="margin:0px 0px 0px
          0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">
          <div>
            <p>Hello everyone,</p>
            <p>I am trying to extract the upper high end tail (95
              percentiles-100percentiles) (shaded area below) of the
              probability distribution.</p>
            <p>Is one way: Area95-100percentiles= mean+2standard
              deviations? <br>
            </p>
            <p>Does anyone know how to extract (resampled) and plot this
              high end/upper tail from a normal distribution?</p>
            <p><br>
            </p>
            <p><img src="cid:part2.CEEBE360.08F79F80@uni-hamburg.de"
                alt="" class=""></p>
            <p><br>
            </p>
            <pre cols="72">-- 
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</pre>
          </div>
          _______________________________________________<br>
          ncl-talk mailing list<br>
          <a href="mailto:ncl-talk@ucar.edu" target="_blank"
            moz-do-not-send="true">ncl-talk@ucar.edu</a><br>
          List instructions, subscriber options, unsubscribe:<br>
          <a href="http://mailman.ucar.edu/mailman/listinfo/ncl-talk"
            rel="noreferrer" target="_blank" moz-do-not-send="true">http://mailman.ucar.edu/mailman/listinfo/ncl-talk</a></blockquote>
      </div>
    </blockquote>
    <pre class="moz-signature" cols="72">-- 
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</pre>
  </body>
</html>