[ncl-talk] The 99th percentile analysis method with GEV distribution and return value calculation
Sri.durgesh Nandini-Weiss
sri.durgesh.nandini-weiss at uni-hamburg.de
Wed Nov 27 04:10:19 MST 2019
Hello everyone,
I was hoping someone would be able to help my with my script for
extracting the 99th percentiles (peak-over-threshold method) from my
data, fit to a GEV distribution using the maximum likelihood method. and
than from this distribution i want to obtain the uncertainty estimates
of the extreme value distribution (i.e the 5% and the 95% confidence
bounds).
Can someone please advice me?
Really appreciate this!
Sri
On 11/26/19 3:40 PM, Sri.durgesh Nandini-Weiss via ncl-talk wrote:
>
> Hello everyone,
>
>
> Im trying to calculate extreme value analysis and plot the CDF of the
> long term changes in extreme sea level by 2100 using the 99^th
> percentiles and the GEV distribution method.
>
> According to previous studies (Vogel et al., 1993
> <https://www.frontiersin.org/articles/10.3389/feart.2018.00216/full#B33>;
> Huang et al., 2008
> <https://www.frontiersin.org/articles/10.3389/feart.2018.00216/full#B16>;
> Feng and Jiang, 2015
> <https://www.frontiersin.org/articles/10.3389/feart.2018.00216/full#B10>)
> the GEV distribution was used to get the return levels of the extreme
> sea levels. Percentile analysis method is widely used to assess the
> extreme sea level changes (Menéndez and Woodworth, 2010
> <https://www.frontiersin.org/articles/10.3389/feart.2018.00216/full#B28>;
> Feng et al., 2015
> <https://www.frontiersin.org/articles/10.3389/feart.2018.00216/full#B11>;
> Marcos and Woodworth, 2017
> <https://www.frontiersin.org/articles/10.3389/feart.2018.00216/full#B26>).
>
> I modified NCL example extvalv_6 but it uses the Annual maxima series
> and I have used the highest percentile method (threshold method): see
> my script attached below. My steps are:
>
> 1) Calculate and extract the 99th percentile value of the anomaly data
> using stat_dispersion2) Basic statistics of the original sample
>
> 3) Using extval_mlegev and extval_gev
>
> 4) Estimate GEV distribution parameters: Can I use the GEV
> distribution to get the return levels of the extreme sea levels or use
> this?
>
> Pr = 0.99 ; user specified
> Nr = extval_return_period(anom1_1d, Pr) ; Nr=?? years at 0.99 probability
> print(Nr)
> The risks associated with extreme sea levels can be assessed from the
> estimates of return levels and return periods. The return period of
> extreme sea level is defined as the sea level statistically expected
> to be equaled or exceeded every specific year.
> 5) Plot the CDF and the return value with uncertainty estimates
>
> At the moment i have errors when calculating my return value: "(23998)
> extval_return_period: Tr must be > 0: Tr=4.20953" and so on.
>
> Can someone please help me with this, or use my script on a sample
> temperature data?
>
> f = addfile ("anom_hist_zo_slr.nc", "r")
> hist_anom = f->hist_anom ;(time|240,ens|100,lat|90,lon|45)
> printVarSummary(hist_anom)
>
> f1 = addfile ("anom_rcp85_zo_slr.nc", "r")
> rcp45_anom = f1->rcp85_anom ;(time|240,ens|100,lat|90,lon|45)
> printVarSummary(rcp45_anom)
>
>
> ;==================================================================
> print ("Data read in - start Extreme value calculations")
> ;==================================================================
> ; Region domain lat and lon
> ;==================================================================
>
> hist_anom1=hist_anom(:,:,{18},{178}) ;(ntim, ens)
> printVarSummary(hist_anom1)
>
> rcp45_anom1=rcp45_anom(:,:,{18},{178}) ;(ntim, ens)
> printVarSummary(rcp45_anom1)
>
> ;=================================================================
>
> anom1_1d = ndtooned(hist_anom1)
> printVarSummary(anom1_1d)
> anom2_1d = ndtooned(rcp45_anom1)
> printVarSummary(anom2_1d)
>
> ;=================================================================
> ; Use stat_dispersion over each region: to calculate 99 percentiles
> ; (a) For each time period
> ; (b) For each statistic,
> ; (c) For each gridpoint? or for all gridpoints? [45*90] sort into
> ascending order and pick the highest 1%
> ; STATX(22 to 29) added for v6.1.0
> ; STATX(22 to 27) require more
> than 1000 observations
> ; ; Otherwise, they are set to
> _FillValue
> ;
> ; Lower 0.1%="+STATX(22) ; lower 0.1% => 0.1th percentile
> ; Lower 1.0%="+STATX(23) ; lower 1.0% => 1th
> ; Lower 5.0%="+STATX(24) ; lower 5.0% => 5th
> ; Upper 5.0%="+STATX(25) ; upper 5.0% => 95th percentile
> ; Upper 1.0%="+STATX(26) ; upper 1.0% => 99th
> ; Upper 0.1%="+STATX(27) ; Upper 0.1% => 99.9th
>
> ;=================================================================
>
> hist_perc = stat_dispersion(anom1_1d, True)
> printVarSummary(hist_perc)
> hist_99= hist_perc(26)
>
> rcp45_perc = stat_dispersion(anom2_1d, True)
> printVarSummary(rcp45_perc)
> rcp45_99= hist_perc(26)
>
> ; In extreme analysis, a decision has to be made regarding the
> definition of the extreme events,either in the form of the number of
> extremes (when a GEV distribution is used) the method of maximum
> likelihood is used to estimate GEV parameters following Kharin and
> Zwiers [2005].
>
>
> vals = extval_mlegev ( anom1_1d, 0, False) ; dims=0
> print(vals) ; vals(6)
>
> center = vals(0) ; extract for
> clarity
> scale = vals(1)
> shape = vals(2)
>
>
> ;================================================================================
> ; extval_return_period:Calculates the period of an event (eg, flood,
> heat wave, drought) occurring given an average event recurrence
> interval and specified probability level.
> ; The risks associated with extreme sea levels can be assessed from
> the estimates of return levels and return periods. The return period
> of extreme sea level is defined as the sea level statistically
> expected to be equaled or exceeded every specific year.
> ;================================================================================
>
> Pr = 0.99 ; user specified
> Nr = extval_return_period(anom1_1d, Pr) ; Nr=?? years at 0.99
> probability
> print(Nr)
>
> ;================================================================================
> pdfcdf = extval_gev( anom1_1d, shape, scale, center,0)
> pdf = pdfcdf[0] ; extract from list
> for convenience
> cdf = pdfcdf[1]
> delete(pdfcdf)
>
> printVarSummary(pdf)
> printMinMax(pdf,1)
>
> printVarSummary(cdf)
> printMinMax(cdf,1)
>
> ;================================================================================
> ; Plot the PDF and CDF
> ;================================================================================
>
> nscl = dimsizes( scale )
> labels = new( nscl , "string")
> plot = new( 2 , "graphic")
>
> wks = gsn_open_wks ("png","extval")
>
> res = True
> res at gsnDraw = False
> res at gsnFrame = False
>
> res at xyLineThicknessF = 2.0
> res at xyMonoDashPattern = True
> res at xyLineColors = (/"green"/)
> res at trYMinF =-0.02 ; min value on y-axis
> res at trYMaxF = 0.5 ; max value on y-axis
>
> res at pmLegendDisplayMode = "Always" ; turn on legend
> res at pmLegendSide = "Top" ; Change location of
> res at pmLegendParallelPosF = 0.675 ; move units right
> res at pmLegendOrthogonalPosF = -0.30 ; more neg = down
> res at pmLegendWidthF = 0.12 ; Change width and
> res at pmLegendHeightF = 0.150 ; height of legend.
> res at lgLabelFontHeightF = 0.0175 ; change font height
> res at lgPerimOn = True
> res at lgItemOrder = ispan(nscl-1,0,1)
> res at xyExplicitLegendLabels = " shp="+sprintf("%3.1f", shape)+ \
> "; scl="+sprintf("%3.1f", scale)+";
> ctr="+sprintf("%3.1f", center)
> it = dim_pqsort_n(anom1_1d, 1, 0) ; for plot reasons
> make 'x' is ascending
> printVarSummary(it)
>
> res at tiYAxisString = "PDF: extval_gev"
> plot(0)= gsn_csm_xy (wks,anom1_1d(it),pdf(it),res) ; create plot
>
> res at trYMinF = 0.0
> res at trYMaxF = 1.0
> res at pmLegendOrthogonalPosF = -0.95 ; more neg = down
> res at tiYAxisString = "CDF: extval_gev"
> plot(1)= gsn_csm_xy (wks,anom1_1d(it),cdf(it),res) ; create plot
>
> resP = True ; modify the
> panel plot
> resP at gsnMaximize = True ; ps, eps, pdf
> resP at gsnPanelMainString = "GEV: PDF and CDF"
> gsn_panel(wks,plot,(/1,2/),resP)
>
>
>
> 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
> 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/20191127/76e9ca58/attachment.html>
More information about the ncl-talk
mailing list