\documentclass[a4paper]{article}

\usepackage[margin=1.7cm]{geometry}
\usepackage[colorlinks=true,urlcolor=blue]{hyperref}

\author{Mauricio Zambrano-Bigiarini}
%\pdfbookmark[0]{Titlepage}{title} % Sets a PDF bookmark for the title page
\title{Tutorial for Introductory Analysis of \\
       Daily Precipitation Data with hydroTSM}
\date{Last Update: 30-May-2013}
%\date{May 2013, Aug 2011, Apr 2011, Oct 2010}


%\VignetteIndexEntry{Tutorial for introductory analysis of daily precipitation data with hydroTSM}
%\VignetteKeyword{hydrology}
%\VignetteKeyword{hydrological modelling}


\begin{document}

\maketitle

\section{Installation}
\label{sec:Installation}

Installing the latest stable version (from \href{http://cran.r-project.org/web/packages/hydroTSM/}{CRAN}):
<<eval=FALSE>>=
install.packages("hydroTSM")
@

\medskip

\noindent Alternatively, you can also try the under-development version (from  \href{http://www.rforge.net/hydroTSM/}{rforge}):
<<eval=FALSE>>=
install.packages("hydroTSM",, "http://rforge.net/", type="source")
@


\section{Setting Up the Environment}
\label{sec:SettingUptheEnvironment}

\begin{enumerate}

\item Loading the \textit{hydroTSM} library, which contains data and functions used in this analysis.
<<>>=
library(hydroTSM)
@

\item Loading daily precipitation data at the station San Martino di Castrozza, Trento Province, Italy, with data from 01/Jan/1921 to 31/Dec/1990.

<<>>=
data(SanMartinoPPts)
@

\item Selecting only a 6-years time slice for the analysis

<<>>=
x <- window(SanMartinoPPts, start=as.Date("1985-01-01"))
@

\item Monthly values of precipitation
<<>>=
( m <- daily2monthly(x, FUN=sum) )
@

\item Dates of the daily values of 'x'
<<>>=
dates <- time(x)
@

\item Amount of years in 'x' (needed for computations)
<<>>=
( nyears <- yip(from=start(x), to=end(x), out.type="nmbr" ) )
@

\end{enumerate}


\section{Basic Exploratory Data Analysis}
\label{sec:BasicEDA}

\begin{enumerate}

\item Summary statistics
<<>>=
smry(x)
@


\item Using the \textit{hydroplot} function, which (by default) plots 9 different graphs: 3 ts plots, 3 boxplots and 3 histograms summarizing 'x'. For this example, only daily and monthly plots are produced, and only data starting on 01-Jan-1987 are plotted.

\begin{center}
<<fig=TRUE, pdf=TRUE, eps=TRUE, width=10, height=8>>=
hydroplot(x, var.type="Precipitation", main="at San Martino", 
          pfreq = "dm", from="1987-01-01")
@
\end{center}

\item Amount of days with information (not \texttt{NA}) per year
<<>>=
dwi(x)
@

\item Amount of days with information (not \texttt{NA}) per month per year
<<>>=
dwi(x, out.unit="mpy")
@


\item Plotting the monthly precipitation values for each year, useful for identifying dry/wet months.

<<fig=TRUE>>=     
# Daily zoo to monthly zoo
m <- daily2monthly(x, FUN=sum, na.rm=TRUE)
     
# Creating a matrix with monthly values per year in each column
M <- matrix(m, ncol=12, byrow=TRUE)
colnames(M) <- month.abb
rownames(M) <- unique(format(time(m), "%Y"))
     
# Plotting the monthly precipitation values
require(lattice)
print(matrixplot(M, ColorRamp="Precipitation", 
           main="Monthly precipitation at San Martino st., [mm/month]"))
@

\end{enumerate}



\section{Annual Analysis}
\label{sec:AnnualAnalysis}

\begin{enumerate}

\item Annual values of precipitation
<<>>=
daily2annual(x, FUN=sum, na.rm=TRUE)
@

\item Average annual precipitation

Obvious way:

<<>>=
mean( daily2annual(x, FUN=sum, na.rm=TRUE) )
@

Another way (more useful for streamflows, where \texttt{FUN=mean}):

The function \textit{annualfunction} applies \texttt{FUN} twice over \texttt{x}: 
( i) firstly, over all the elements of \texttt{x} belonging to the same year, in order to obtain the corresponding annual values, and
(ii) secondly, over all the annual values of \texttt{x} previously obtained, in order to obtain a single annual value.

<<>>=
annualfunction(x, FUN=sum, na.rm=TRUE) / nyears
@

\end{enumerate}


% % % % % % % % % % % % % % % % % % % % % % % % %
\section{Monthly Analysis}
\label{sec:MonthlyAnalysis}



\begin{enumerate}
\item Median of the monthly values at station 'x'. Not needed, just for looking at these values in the boxplot.

<<>>=
monthlyfunction(m, FUN=median, na.rm=TRUE)
@

\item Vector with the three-letter abbreviations for the month names
<<>>=
cmonth <- format(time(m), "%b")
@

\item Creating ordered monthly factors
<<>>=
months <- factor(cmonth, levels=unique(cmonth), ordered=TRUE)
@

\item Boxplot of the monthly values
<<fig=TRUE, pdf=TRUE, eps=FALSE>>=
boxplot( coredata(m) ~ months, col="lightblue", main="Monthly Precipitation", 
         ylab="Precipitation, [mm]", xlab="Month")
@

\end{enumerate}



\section{Seasonal Analysis}
\label{sec:SeasonalAnalysis}



\begin{enumerate}
\item Average seasonal values of precipitation

<<>>=
seasonalfunction(x, FUN=sum, na.rm=TRUE) / nyears
@

\item Extracting the seasonal values for each year
<<>>=
( DJF <- dm2seasonal(x, season="DJF", FUN=sum) )
( MAM <- dm2seasonal(m, season="MAM", FUN=sum) )
( JJA <- dm2seasonal(m, season="JJA", FUN=sum) )
( SON <- dm2seasonal(m, season="SON", FUN=sum) )
@


\item Plotting the time evolution of the seasonal precipitation values
<<fig=TRUE, pdf=TRUE, eps=FALSE, width=12, height=10>>=
hydroplot(x, pfreq="seasonal", FUN=sum, stype="default")
@


\end{enumerate}




%#########################################
\clearpage
\section{Some Extreme Indices}
\label{sec:SomeExtremeIndices}

Common steps for the analysis of this section:

\begin{enumerate}
\item Loading daily precipitation data at the station San Martino di Castrozza, Trento Province, Italy, with data from 01/Jan/1921 to 31/Dec/1990.
<<>>=
data(SanMartinoPPts)
@
	
\item Selecting only a three-year time slice for the analysis
<<>>=
x <- window(SanMartinoPPts, start=as.Date("1988-01-01"))
@

\item Plotting the selected time series
<<fig=TRUE, pdf=TRUE, eps=FALSE>>=
hydroplot(x,  ptype="ts", pfreq="o", var.unit="mm")
@

\end{enumerate}

%################
\subsection{Heavy Precipitation Days (R10mm)}
\label{subsec:R10mm}

\begin{enumerate}
\item Counting and plotting the number of days in the period where precipitation is > 10 [mm]
<<>>=
( R10mm <- length( x[x>10] ) )
@

\end{enumerate}


%################
\subsection{Very Wet Days (R95p)}
\label{subsec:R95p}

\begin{enumerate}

\item Identifying the wet days (daily precipitation >= 1 mm):
<<>>=
wet.index <- which(x >= 1)
@
	
	
\item Computing the 95th percentile of precipitation on wet days (\textit{PRwn95}):
<<>>=
( PRwn95 <- quantile(x[wet.index], probs=0.95, na.rm=TRUE) )
@

\textbf{Note 1}: this computation was carried out for the three-year time period 1988-1990, not the 30-year period 1961-1990 commonly used. \newline
\textbf{Note 2}: missing values are removed from the computation.


\item Identifying the very wet days (daily precipitation >= \textit{PRwn95})
<<>>=
(very.wet.index <- which(x >= PRwn95))
@


\item Computing the total precipitation on the very wet days:
<<>>=
( R95p <- sum(x[very.wet.index]) )
@
\textbf{Note 3}: this computation was carried out for the three-year time period 1988-1990, not the 30-year period 1961-1990 commonly used \newline

\end{enumerate}



%################
\subsection{5-day Total Precipitation}
\label{subsec:5DayTotalPrecipitation}

\begin{enumerate}


\item Computing the 5-day total (accumulated) precipitation \newline
<<fig=TRUE, pdf=TRUE, eps=FALSE>>=
x.5max <- rollapply(data=x, width=5, FUN=sum, fill=NA, partial= TRUE, 
                    align="center")

hydroplot(x.5max,  ptype="ts+boxplot", pfreq="o", var.unit="mm")
@

\item Maximum annual value of 5-day total precipitation
<<>>=
(x.5max.annual <- daily2annual(x.5max, FUN=max, na.rm=TRUE))
@

\textbf{Note 1}: for this computation, a moving window centred in the current day is used. If the user wants the 5-day total precipitation accumulated in the 4 days before the current day + the precipitation in the current day, the user have to modify the moving window.\newline

\textbf{Note 2}: For the first two and last two values, the width of the window is adapted to ignore values not within the time series


\end{enumerate}


%##################################
\section{Software Details}
\label{sec:SoftwareDetails}
This tutorial was built under: 

<<echo=FALSE>>=
sessionInfo()$platform
sessionInfo()$R.version$version.string 
paste("hydroTSM", sessionInfo()$otherPkgs$hydroTSM$Version)
@



\end{document}
