[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
Tue Nov 26 07:40:03 MST 2019
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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20191126/5b60c615/attachment.html>
More information about the ncl-talk
mailing list