[ncl-talk] extreme value analysis

[Adam Phillips]

Hi Sri,
I think the best way to diagnose the issue is for you to send your data to ncl-talk. However, before doing so, please closely examine the values that you are passing in to the gsn_csm_xy plotting function. Specifically, print out the xHigh(it), pdf(it) and cdf(it) array values. If you do not see anything amiss with the values send the data file to ncl-talk.
Adam

> On Dec 12, 2019, at 7:13 AM, Sri.durgesh Nandini-Weiss via ncl-talk <ncl-talk at ucar.edu> wrote:
> 
> 
> Hello everyone
> 
> 
> 
> I tried calculating and plotting extreme value analysis on my SSH data, following ncl example extval_6.ncl: 
> 
> My plot was generated with no error, but it doe not look correct.
> 
> See attached, as well as my script. Can someone advice me? I am happy to send my data as well.
> 
> Sri
> 
> ;************************************************
> ; extvalv_6.ncl
> 
>    f     = addfile ("anom_hist_zo_slr.nc", "r")
>    hist_anom    = f->hist_anom
>    printVarSummary(hist_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)
> ;=================================================================
> ; extract data between 95-100percentiles
> ;=================================================================
> 
>   anom1_1d = ndtooned(hist_anom1)                     
>   printVarSummary(anom1_1d)
> 
>   anom1_nomiss=anom1_1d(ind(anom1_1d.ne.-9999))
>   anom1_permsort = dim_pqsort(anom1_nomiss,2)                                                            
>   x95=anom1_nomiss(floattoint(0.95*dimsizes(anom1_nomiss))) 
>   i95  =  ind(anom1_nomiss.ge.x95)
>   xHigh = anom1_nomiss(i95)                 ; values               
>   printVarSummary(xHigh)
> 
> ;***************************************************************
> ;---GEV parameter estimation
> ;***************************************************************
> 
>    vals     = extval_mlegev (xHigh, 0, False)    ; dims=0
>    print(vals) 
> 
>    center   = vals(0)                            ; extract for clarity
>    scale    = vals(1)
>    shape    = vals(2)
>    ;ptype=0
>    pdfcdf   = extval_gev(xHigh, shape, scale, center,0)
>    ;pdfcdf = extval_pareto(xHigh, shape, scale, center, ptype, 0)   ; [/ PDF, CDF /]
>    pdf      = pdfcdf[0]                          ; extract from list for convenience
>    cdf      = pdfcdf[1]
>    delete(pdfcdf)
> 
>    printVarSummary(pdf)
>    printMinMax(pdf,1)
>    print("---")
>    printVarSummary(cdf)
>    printMinMax(cdf,1)
>    print("---")
> 
> ;***************************************************************
> ;--- PLOTS
> ;***************************************************************
> 
>    plotfile="Hist_SSH"
>   wks_type = "pdf"
>   wks = gsn_open_wks(wks_type,plotfile)
> 
>    gsn_define_colormap(wks,"default")        
>    plot = new(2, "graphic")
>    res     = True       
>    res at gsnDraw  = False
>    res at gsnFrame = False
>    res at tiYAxisString = ""
> 
>    it      = dim_pqsort_n(xHigh, 1, 0)        ; for plot reasons make  'x' is asebding 
> 
> ;************************************************
> ; Panel
> ;************************************************
> 
>    resP                     = True                  ; modify the panel plot
>    resP at gsnMaximize         = True                  ; ps, eps, pdf
>    resP at gsnPanelMainString  = "River Flow Rate"     ; use this for NCL V6.4.0 and later
>    resP at txFontHeightF       = 0.020
> 
> ;************************************************
> ;--- Plot PDF and CDF ; create panel
> ;************************************************
> 
>    plot(0) = gsn_csm_xy (wks,xHigh(it),pdf(it),res) ; create plot
> 
>    res at trYMinF  = 0.0   
>    res at trYMaxF  = 1.0  
> 
>    plot(1) = gsn_csm_xy (wks,xHigh(it),cdf(it),res) ; create plot
>  
>    gsn_panel(wks,plot,(/1,2/), resP)             ; now draw as one plot
> 
> 
> 
> 
> -- 
> 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
> <Hist_SSH.pdf>
> _______________________________________________
> ncl-talk mailing list
> ncl-talk at ucar.edu
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20191212/1ffa4365/attachment.html>