Confirm both this file ("2016AppliedHydrologyAssignment.nb") and "AppliedHydrology2016dailyRainfall_1.xls" or “AppliedHydrology2016dailyRainfall_2.xls” are located in the same folder in your PC.

How to use this file :
   to execute command  :  Shift + Enter key
      to quit evaluation cells  :  Command + .  key

You can execute each cell one by one and also at once. In order to eveluate or execute several cells, first, select cells and then press Shift + Enter key.

With next command, you can set working directory to the folder which include this notebook.

SetDirectory[NotebookDirectory[]]

FileNames[]

Read data from xls file

targetData=Import["AppliedHydrology2016dailyRainfall_1.xls"][[1]];

(*targetData=Import["AppliedHydrology2016dailyRainfall_2.xls"][[1]];*)

targetData[[;;2]](*check the first 2 lines*)

Usefull Modules

grids[min_,max_]:=Join[Range[Ceiling[min,10],Floor[max,10],10],Table[{j-5,Dashed},{j,Round[min,10],Round[max,10],10}]]

marks

square={{-0.5,-0.5},{0.5,-0.5},{0.5,0.5},{-0.5,0.5},{-0.5,-0.5}};triangl={{0,1},{Cos[7Pi/6],Sin[7Pi/6]},{Cos[11Pi/6],Sin[11Pi/6]},{0,1}};
colmark=(({{m1, m2, m3, m4, m5, m6}, {m1r, m2r, m3r, m4r, m5r, m6r}, {m1b, m2b, m3b, m4b, m5b, m6b}, {m1g, m2g, m3g, m4g, m5g, m6g}, {m1p, m2p, m3p, m4p, m5p, m6p}, {m1c, m2c, m3c, m4c, m5c, m6c}, {m1m, m2m, m3m, m4m, m5m, m6m}, {m1br, m2br, m3br, m4br, m5br, m6br}, {m1o, m2o, m3o, m4o, m5o, m6o}}))=Graphics/@{{#,Circle[{0,0},1]},{#,Disk[{0,0},1]},{#,Line[square]},{#,Polygon[square]},{#,Line[triangl]},{#,Polygon[triangl]}}&/@{Black,Red,Blue,Green,Pink,Cyan,Magenta,Brown,Orange};

trend test

RunLengthEncode[x_List]:={First[#],Length[#]}&/@Split[x];
(*MannKendall[x_List,ansfull_:False]:=Module[{n=Length[x],Vs2,dd,ddorder,tr2,ti,Sx,Z2,tievec=Table[0,{Length[x]}]},dd=RunLengthEncode[Sort[x]];ddorder=dd[[Ordering[dd[[All,2]]]]];tr2=Transpose[ddorder];ti=Length/@Split[tr2[[2]]];tievec=Join[ti,Drop[tievec,Length[ti]]];Vs2=(1/18.)(n(n-1)(2n+5)-(∑_(i=1)^n)tievec[[i]]i(i-1)(2i+5));Sx=(∑_(k=1)^(n-1))(∑_(j=k+1)^n)Sign[x[[j]]-x[[k]]];Z2=((Sx-Sign[Sx])/(Sqrt[Vs2]));If[ansfull,{Z2,Sx,(Sx/(n(n-1)/2))},Z2]];*)(*20140128にチェックの上、新しいモジュールはtau=S/Dとしている。comment out*)
MannKendall[x_List,ansfull_:False]:=Module[{n=Length[x],SRunLen,Vs2,S,Z2,D,tt2,ttti},SRunLen=SortBy[RunLengthEncode[Sort[x]],Last][[All,2]];ttti=#(#-1)(2#+5)&/@SRunLen//Total;tt2=#(#-1)&/@SRunLen//Total;
Vs2=(1/18)(n(n-1)(2n+5)-ttti);S=(∑_(k=1)^(n-1))(∑_(j=k+1)^n)Sign[x[[j]]-x[[k]]];Z2=((S-Sign[S])/(Sqrt[Vs2]));D=(Sqrt[(1/2)n(n-1)-(1/2)tt2])(Sqrt[(1/2)n(n-1)]);If[ansfull,{Z2,S,(S/D),Vs2}//N,N[Z2]]];

threshold selection

SampleMEFatu[dataX_,u_]:=If[Sort[dataX,Greater][[1]]≥u,Mean[Select[dataX-u,#≥0&]],Null];(*for Real u*)
SampleVarEFatu[dataX_,u_]:=If[Sort[dataX,Greater][[2]]≥u,Variance[Select[dataX-u,#≥0&]],Null];(*for Real u*)

L - moment

Lmoment[data_List]:=Module[{X,n=Length[data],B0,B1,B2,B3,λ1,λ2,λ3,λ4},
(* L moment method from HandBook of Hydrology Chap.18 FREQUENCY ANALYSIS OF EXTREME EVENTS *)
(* unbiased PWM estimators *)
	X=Sort[data,Greater];
B0:=Mean[X];
B1:=(∑_(j=1)^(n-1))(((n-j)X[[j]])/(n(n-1)));
B2:=(∑_(j=1)^(n-2))(((n-j)(n-j-1)X[[j]])/(n(n-1)(n-2)));                                         (* 18.1.14 *)
B3:=(∑_(j=1)^(n-3))(((n-j)(n-j-1)(n-j-2)X[[j]])/(n(n-1)(n-2)(n-3)));
(* L-moment estimators *)
λ1=B0;
λ2=2B1-B0;
λ3=6B2-6B1+B0;
λ4=20B3-30B2+12B1-B0;                                        (* 18.1.16 *)
{λ1,λ2,λ3,λ4,λ2/λ1,λ3/λ2,λ4/λ2}];

Plotting position

PlottingPosition[i_,n_,α_]:=((i-α)/(n+1-2 α));(*Weibull:0, Cunnane:0.4,Hazen:0.5*)
ipos[α_,n_,Ny_]:=α-((Ny )/n)(1+n-2 α) Log[((1-α)/(1+Ny-2 α))]

AMS analysis

AMSDistFit[data_,locName_,sz_:3.5 72,xrange_:{0,400},mark_:{m1,0.04},Tsw_:False,Rotatelbl_:False,imgpos_:{0.18,0.86},GEVcolor_:Purple,LGumbelcolor_:Orange]:=Show[GumbelProbPlot[data,locName,0.4,"Rainfall (mm)",mark,xrange,Tsw,Rotatelbl,sz,imgpos],GEVplot[data,GEVcolor],LGumbelplot[data,LGumbelcolor]];

GumbelProbPlot[x_List,title_:"",α_:0.4,vUnit_:"Rainfall (mm)",mark_:{m1,0.04},xrange_:{0,0},Tsw_:False,Rotatelbl_:False,imgsz_:4 72,imgpos_:{0.18,0.86}]:=Module[{n=Length[x],ynorm,Tnen,grf,orderx,ytick,lab1,lab2,yp},
ynorm[data_]:=-Log[-Log[data]];
		orderx=Sort[x];ytick=0.01{99.5,99,98,95,90,80,70,50,30,20,10};
lab1=Transpose[{ynorm[ytick],ytick}];
Tnen={200,100,50,20,10,5};
lab2=Transpose[{Take[ynorm[ytick],6],Tnen}];
yp=Table[N[(i-α)/(n+1-2α)],{i,n}];      grf=ListPlot[Transpose[{orderx,ynorm[yp]}],PlotRange->{If[xrange=={0,0},{Max[0,Min[x]-0.1(Max[x]-Min[x])],Max[x]+0.3(Max[x]-Min[x])},xrange],{-1.54,6}},PlotMarkers->mark,AspectRatio->1.3,GridLines->{Automatic,ynorm[ytick]},Frame->True,Axes->None,FrameTicks->{{lab1,If[Tsw,lab2,None]},Automatic},FrameLabel->{{ "F",If[Tsw,"T",None]}, {vUnit,None}},PlotLabel->"Gumbel Probability Paper",ImageSize->imgsz,RotateLabel->Rotatelbl,PlotRangePadding->None,Epilog->{Black,Text[Style[title,FontSize->12],ImageScaled[imgpos],{Left,Top},Background->White]}]
	];

GEV  Generalized Extreme Value Distribution

GEVprm[X_List]:=Module[{κ,α,ξ,c,λ},
			 λ=Lmoment[X];
   		c=((2λ[[2]])/(λ[[3]]+3λ[[2]]))-((Log[2])/(Log[3]));
        κ=7.859c + 2.9554c^2;                                                   (* 18.2.22a *)
        α=((κ  λ[[2]])/(Gamma[1+κ](1-(2^(-κ)))));
        ξ=λ[[1]]+((α( Gamma[1+κ]-1))/κ);
		{α,κ,ξ}];
GEVXq[X_List,ReturnPeriod_]:=Module[{κ,α,ξ,prm,Xq},
		prm=GEVprm[X];α=prm[[1]];κ=prm[[2]];ξ=prm[[3]];
		Xq=ξ+(α/κ)(1-(((-Log[1-(1/ReturnPeriod)]))^κ))] ;
GEVplot[X_List,col_:Purple,prmsw_:False]:=Module[{κ,α,ξ,prm,xstrt,GEV,graph,gumgraph,x},
		GEV[x_]:=Exp[-(((1-((κ(x-ξ))/α)))^(1/κ))];    (* κ≠0 18.2.18 *)
		prm=GEVprm[X];α=prm[[1]];κ=prm[[2]];ξ=prm[[3]];
				If[prmsw,Print["GEV   parameters   : α: ",α," κ: ",κ," ξ: ",ξ]];

			xstrt=Max[Min[X]/4,ξ+(α/κ)(1-(((-Log[.004]))^κ))];
		gumgraph=Plot[ -Log[-Log[ GEV[x]]],{x,xstrt,ξ+(α/κ)(1-(((-Log[.9975]))^κ))},PlotStyle->col]];

Gumbel Distribution

LGumbelprm[X_List]:=Module[{α,ξ,λ},
				 λ=Lmoment[X];
          α=λ[[2]]/Log[2.];
          ξ=λ[[1]]-0.5772α ;
     	{α,ξ}];
LGumbelXq[X_List,ReturnPeriod_]:=Module[{α,ξ,prm,Xq},
		prm=LGumbelprm[X];α=prm[[1]];ξ=prm[[2]];
		Xq=ξ-α  Log[-Log[1-(1/ReturnPeriod)]]] ;
LGumbelplot[X_List,col_:Orange]:=Module[{α,ξ,prm,xstrt,Fgumbel,sp,graph,gumgraph,x},
	Fgumbel[x_]:=Exp[-Exp[-((x-ξ)/α)]];
			sp[p_]:=-Log[-Log[p]];
		prm=LGumbelprm[X];α=prm[[1]];ξ=prm[[2]];
			xstrt=Max[Min[X]/4,ξ+α sp[.004]];
		gumgraph=Plot[ -Log[-Log[ Fgumbel[x]]],{x,xstrt,ξ+α sp[.9975]},PlotStyle->col]];

POT analysis

GumbelpaperPOTtoAME[x_List,ny_,α_:0.4,unit_:"Rainfall [mm]",mark_:{m2r,0.03},xrange_:{0,0},imgsz_:4 72,title_:"",tpos_:{0.3,0.83}]:=Module[{n=Length[x],ynorm,Tnen,grf,orderx,ytick,lab1,lab2,yp},
ynorm[data_]:=-Log[-Log[data]];
		orderx=Sort[x];ytick=0.01{99.5,99,98,95,90,80,70,50,30,20,10,5,1};
lab1=Transpose[{ynorm[ytick],ytick}];
Tnen={200,100,50,20,10,5};
lab2=Transpose[{Take[ynorm[ytick],6],Tnen}];
yp= PlottingPosition[Range[n],n,α];    
				yp=Exp[-(n/ny)(1-yp)];  grf=ListPlot[Transpose[{orderx,ynorm[yp]}],PlotMarkers->mark,PlotRange->{If[xrange=={0,0},{Max[0,Min[x]-0.1(Max[x]-Min[x])],Max[x]+0.2(Max[x]-Min[x])},xrange],{-1.54,6}},AspectRatio->1.3,GridLines->{Automatic,ynorm[ytick]},Frame->True,Axes->None,FrameTicks->{{lab1,lab2},Automatic},FrameLabel->{ { "F","T"},{unit,None}},RotateLabel->False,ImageSize->imgsz,PlotLabel->"Gumbel Probability Paper",Epilog->{Black,Text[title,ImageScaled[tpos],{-1,0},Background->White]}]
	];

Exponential Distribution

Expprm[X_List,threshold_:-9999]:=Module[{α,ξ},
  ξ=If[threshold==-9999, Min[X],threshold]; α=Mean[X]-ξ;
		{α,ξ}];(*L moment and Product moment is same *)

(* Using ny information, draw Exponential distribution on Gumbel Probability Papaer *)ExpPlotonGumPaper[POT_List,ny_,threshold_:-9999,col_:{Orange,Dashed}]:=Module[{n=Length[POT],α,ξ,prm,xstrt,Gexp,sp,x},
	Gexp[x_]:=1-Exp[-(((x-ξ)/α))];
		sp[p_]:=-Log[-Log[p]];
	{α,ξ}=Expprm[POT,threshold];
	Plot[ -Log[(n/ny)(1-Gexp[x])],{x,ξ,ξ+α sp[0.9999]},PlotStyle->col]];

Generalized Pareto Distribution

GPprm[X_List,threshold_:-9999]:=Module[{κ,α,ξ,λ},  ξ=If[threshold==-9999, Min[X],threshold];
			 λ=Lmoment[X];
   		κ=((λ[[1]]-ξ)/(λ[[2]]))-2.;
        α=(λ[[1]]-ξ)(1+κ);
		{α,κ,ξ}];

(*Using ny information, draw GP distribution on Gumbel Probability Papaer with appropriate range *)GPPlotonGumPaper[POT_List,ny_,threshold_:-9999,col_:{Purple,Dotted}]:=Module[{n=Length[POT],α,κ,ξ,Ggp,sp,x},
	Ggp[x_]:=1-(((1-((κ(x-ξ))/α)))^(1/κ));
		sp[p_]:=-Log[-Log[p]];
{ α,κ,ξ}=GPprm[POT,threshold];
Plot[ -Log[(n/ny)(1-Ggp[x])],{x,ξ,Which[κ>0,ξ+(α/κ) ,κ≤0,3Max[POT]]},PlotStyle->col]];

Practice for each assignment

Example: Suppose data of your assignment is in column No.4 (Sapporo) in “AppliedHydrology2016dailyRainfall_ 1.xls”

Pick up all of your concerning data

col=4(*put columun No. of your data *)

MyData=targetData[[2;;,{1,col}]];(*processing col-4*)

Dimensions[MyData]

yearly=SplitBy[MyData,DateList[#[[1,1,1]]]&];

Dimensions[yearly]

yearly[[1]]

yearly[[-1]]

This means 25 years, 365 days each year.

2) Pick up AMS

MyAMS=Max/@yearly[[All,All,2]]

Length[MyAMS]

Min[MyAMS]

Move white button to find suitable bin size.

Manipulate[Histogram[MyAMS,{bin},"PDF",PlotLabel->targetData[[1,col]]<>"  AMS",PlotRange->All],{bin,3,20,1}]

3) Assess with AMS

3 - 1) probability plot on a Gumbel probability paper using plotting position Cunnane

plotting position  : ((i-α)/(n+1-2 α))      in this case, n=25, i=1,2,3,,,25, α=0.4

First, Sort AMS data,

gams=Sort[MyAMS,Greater]

Next, prepare {{x1, y1}, {x2, y2}, , , , {xn, yn}}

gamsxy=Table[{gams[[i]],1-((i-α)/(Length[gams]+1-2 α))/.α->0.4},{i,Length[gams]}]

ListPlot[gamsxy,Frame->True,GridLines->Automatic]

Change Y axis into Gumbel scale; -Log(-Log(y))

gamsxy2=Table[{gams[[i]],-Log[-Log[1-((i-α)/(Length[gams]+1-2 α))/.α->0.4]]},{i,Length[gams]}]

ListPlot[gamsxy2,Frame->True,GridLines->Automatic]

ListPlot[gamsxy2,PlotMarkers->{m1,0.04},PlotRange->{{0,300},{-1.6,6}},Frame->True,GridLines->Automatic,AspectRatio->1.3]

Change frame ticks

AMSPP=GumbelProbPlot[MyAMS,targetData[[1,col]]]

3 - 2) estimate parameters of distributions

gumprm=LGumbelprm[MyAMS]

gevprm=GEVprm[MyAMS]

3 - 3) plot the distributions onto  3 - 1)

gumplt=Plot[-Log[-Log[Exp[-Exp[-((x-ξ)/α)]]]]/.{α->gumprm[[1]],ξ->gumprm[[2]]},{x,0,500},PlotStyle->Orange]

gevplt=Plot[-Log[-Log[Exp[-(((1-((κ(x-ξ))/α)))^(1/κ))]]]/.{α->gevprm[[1]],κ->gevprm[[2]],ξ->gevprm[[3]]},{x,0,500}]

Show[AMSPP,gumplt,gevplt]

3 - 4) estimate quantile

Missing