<html>
<head>
<meta http-equiv="content-type" content="text/html; charset=UTF-8">
</head>
<body>
<p>Hello everyone</p>
<p><br>
</p>
<p>I tried calculating and plotting extreme value analysis on my SSH
data, following ncl example <a
href="https://www.ncl.ucar.edu/Applications/Scripts/extval_6.ncl">extval_6.ncl</a>:
<br>
</p>
<p>My plot was generated with no error, but it doe not look correct.</p>
<p>See attached, as well as my script. Can someone advice me? I am
happy to send my data as well.</p>
<p>Sri</p>
<p>;************************************************<br>
; extvalv_6.ncl<br>
<br>
f = addfile ("anom_hist_zo_slr.nc", "r")<br>
hist_anom = f->hist_anom<br>
printVarSummary(hist_anom)<br>
<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>
; print ("Data read in - start calculations")<br>
; Region domain lat and lon<br>
;==================================================================<br>
<br>
hist_anom1=hist_anom(:,:,{48},{-125}) ;(ntim,
mlon) <br>
printVarSummary(hist_anom1)<br>
;=================================================================<br>
; extract data between 95-100percentiles<br>
;=================================================================<br>
<br>
anom1_1d = ndtooned(hist_anom1) <br>
printVarSummary(anom1_1d)<br>
<br>
anom1_nomiss=anom1_1d(ind(anom1_1d.ne.-9999))<br>
anom1_permsort =
dim_pqsort(anom1_nomiss,2)
<br>
x95=anom1_nomiss(floattoint(0.95*dimsizes(anom1_nomiss))) <br>
i95 = ind(anom1_nomiss.ge.x95)<br>
xHigh = anom1_nomiss(i95) ; values
<br>
printVarSummary(xHigh)<br>
<br>
;***************************************************************<br>
;---GEV parameter estimation<br>
;***************************************************************<br>
<br>
vals = extval_mlegev (xHigh, 0, False) ; dims=0<br>
print(vals) <br>
<br>
center = vals(0) ; extract for
clarity<br>
scale = vals(1)<br>
shape = vals(2)<br>
;ptype=0<br>
pdfcdf = extval_gev(xHigh, shape, scale, center,0)<br>
;pdfcdf = extval_pareto(xHigh, shape, scale, center, ptype,
0) ; [/ PDF, CDF /]<br>
pdf = pdfcdf[0] ; extract from
list for convenience<br>
cdf = pdfcdf[1]<br>
delete(pdfcdf)<br>
<br>
printVarSummary(pdf)<br>
printMinMax(pdf,1)<br>
print("---")<br>
printVarSummary(cdf)<br>
printMinMax(cdf,1)<br>
print("---")<br>
<br>
;***************************************************************<br>
;--- PLOTS<br>
;***************************************************************<br>
<br>
plotfile="Hist_SSH"<br>
wks_type = "pdf"<br>
wks = gsn_open_wks(wks_type,plotfile)<br>
<br>
gsn_define_colormap(wks,"default") <br>
plot = new(2, "graphic")<br>
res = True <br>
res@gsnDraw = False<br>
res@gsnFrame = False<br>
res@tiYAxisString = ""<br>
<br>
it = dim_pqsort_n(xHigh, 1, 0) ; for plot reasons
make 'x' is asebding <br>
<br>
;************************************************<br>
; Panel<br>
;************************************************<br>
<br>
resP = True ; modify the
panel plot<br>
resP@gsnMaximize = True ; ps, eps, pdf<br>
resP@gsnPanelMainString = "River Flow Rate" ; use this for
NCL V6.4.0 and later<br>
resP@txFontHeightF = 0.020<br>
<br>
;************************************************<br>
;--- Plot PDF and CDF ; create panel<br>
;************************************************<br>
<br>
plot(0) = gsn_csm_xy (wks,xHigh(it),pdf(it),res) ; create plot<br>
<br>
res@trYMinF = 0.0 <br>
res@trYMaxF = 1.0 <br>
<br>
plot(1) = gsn_csm_xy (wks,xHigh(it),cdf(it),res) ; create plot<br>
<br>
gsn_panel(wks,plot,(/1,2/), resP) ; now draw as one
plot<br>
<br>
</p>
<p><br>
</p>
<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>