\documentclass[11pt,a4paper]{article}
\usepackage[utf8]{inputenc}
\usepackage{natbib}
\usepackage{authblk}
\usepackage{hyperref}
\usepackage{gensymb}
\usepackage{amsmath}

%\VignetteIndexEntry{User's guide to meteoland}
%\VignettePackage{meteoland}

\title{User's guide to meteoland}
\author[1,2]{Miquel De Cáceres}
\affil[1]{Centre Tecnològic Forestal de Catalunya. Ctra. St. Llorenç de Morunys km 2, 25280, Solsona, Catalonia, Spain}
\affil[2]{CREAF, Cerdanyola del Vallès, 08193, Spain}

\begin{document}
\SweaveOpts{concordance=TRUE}

\maketitle
\tableofcontents

\newpage

<<echo=FALSE>>=
options(width=67)
library(meteoland)
@
\section{Short introduction to the package}
\subsection{Purpose}
Reliable meteorological data are a basic requirement for hydrological and ecological studies at the landscape scale. Given the large spatial variation of meteorology over complex terrains, meteorological records from a single weather station are often not representative of entire landscapes. Studies made on multiple sites over a landscape require different meteorological series for each site; and other studies may require meteorological data series for all grid cells of a landscape, in a continuous way. In these cases, spatial correlation between the meteorology series of different sites or cells must be taken into account. For example, the sequence of days with rain of contiguous cells will normally be the same or very similar, even if precipitation amounts differ. Finally, studies addressing the impacts of climate change on forests and landscapes require downscaling coarse-scale predictions  of global or regional climate models to the landscape scale. When downscaling predictions for several locations in a landscape, spatial correlation of downscaled predictions is also important.

With the aim to assist research of climatic impacts on forests, the R package \texttt{meteoland} provides utilities to estimate daily weather variables at any position in complex terrains:
\begin{enumerate}
\item{Spatial interpolation of daily weather records from meteorological stations.}
\item{Statistical downscaling of coarse-scale meteorological data (coming from regional climate models or re-analyses) to the landscape scale.}
\end{enumerate}


\subsection{Meteorological variables}
Package \texttt{meteoland} assists in the estimation of the following meteorological variables over lanscapes (units in parentheses):
\begin{itemize}
\item{\texttt{DOY}: Day of the year (Julian day).}
\item{\texttt{MeanTemperature}: Mean daily temperature (in degrees Celsius).}
\item{\texttt{MinTemperature}: Minimum daily temperature (in degrees Celsius).}
\item{\texttt{MaxTemperature}: Maximum daily temperature (in degrees Celsius).}
\item{\texttt{Precipitation}: Daily precipitation (in mm of water).}
\item{\texttt{MeanRelativeHumidity}: Mean daily relative humidity (in percent).}
\item{\texttt{MinRelativeHumidity}: Minimum daily relative humidity (in percent).}
\item{\texttt{MaxRelativeHumidity}: Maximum daily relative humidity (in percent).}
\item{\texttt{Radiation}: Incoming radiation (in MJ/m2).}
\item{\texttt{WindSpeed}: Wind speed (in m/s).}
\item{\texttt{WindDirection}: Wind direction (in degrees from North).}
\item{\texttt{PET}: Potential evapo-transpiration (in mm of water).}
\end{itemize}
Although internally it uses specific units, \texttt{meteoland} allows reading and writing meteology data in different units and formats.

\subsection{Spatial variation of meteorology and topography}
Package \texttt{meteoland} deals with two kinds of spatial structures: individual points and grids. The package includes four S4 spatial classes, which are defined as children of classes in package \texttt{sp}. Two classes are defined to represent the variation of topographic features (i.e., elevation, slope and aspect) over space:
\begin{itemize}
\item{\texttt{SpatialPointsTopography} extends \texttt{SpatialPointsDataFrame} and represents the topographic features of a set of points in a landscape.
<<echo=FALSE>>=
showClass("SpatialPointsTopography")
@
}
\item{\texttt{SpatialGridTopography} extends \texttt{SpatialGridDataFrame} and represents the continuous variation of topography over a grid.
<<echo=FALSE>>=
showClass("SpatialGridTopography")
@
}
\end{itemize}
Although they have the same slots as their parent S4 classes, the data frames in \texttt{SpatialPointsTopography} and \texttt{SpatialGridTopography} have only three variables: `\textbf{elevation}' (in meters), `\textbf{slope}' (in degrees) and `\textbf{aspect}' (in degrees from North).

Two analogous spatial classes are used to represent the variation of daily meteorology over space:
\begin{itemize}
\item{\texttt{SpatialPointsMeteorology} extends \texttt{SpatialPoints} and represents daily meteorology series for a set of points in a landscape.
<<echo=FALSE>>=
showClass("SpatialPointsMeteorology")
@
}
\item{\texttt{SpatialGridMeteorology} extends \texttt{SpatialGrid} and represents the continuous variation of daily meteorology over a grid of cells.
<<echo=FALSE>>=
showClass("SpatialGridMeteorology")
@
}
\end{itemize}
In addition to their corresponding inherited slots, \texttt{SpatialPointsMeteorology} and \texttt{SpatialGridMeteorology} have two additional slots: `\texttt{dates}' (a vector of days specifying a time period), and `\texttt{data}' (a vector of data frames with the meteorological data).
Although both have a `\texttt{data}' with data frames, meteorological data is in different form in each class. In objects of \texttt{SpatialPointsMeteorology}, there is one data frame for each point with variables in columns and dates in rows. In objects of \texttt{SpatialGridMeteorology}, each data frames describes the meteorology over the grid for one day, with grid cells in rows and variables in columns.

\subsection{Reading and writing meteorological data}
\subsubsection{Point meteorology}
Objects of class \texttt{SpatialPointsMeteorology} are stored using one \textbf{ascii} data file for each spatial point. Package \texttt{meteoland} provides four input/output functions for point meteorology:
\begin{itemize}
\item{Function \texttt{readmeteorologypoint()} reads the meteorological data stored in one \textbf{ascii} data file and returns a data frame.}
\item{Function \texttt{writemeteorologypoint()} writes the meteorological data of a single point as an \textbf{ascii} file in the file system.}
\item{Function \texttt{readmeteorologypointfiles()} reads several \textbf{ascii} files and returns an object of class \texttt{SpatialPointsMeteorology}.}
\item{Functions \texttt{writemeteorologypointfiles()} writes several \textbf{ascii} files in the disk, one per spatial point. Metadata (i.e. the spatial coordinates of each point and the corresponding file path) is stored in an additional file.}
\end{itemize}

\subsubsection{Grid meteorology}
 Objects of class \texttt{SpatialGridMeteorology} are stored using one \textbf{netCDF} file per day, which also contains the date and spatial projection. The following functions are available for input/output of grid meteorology:
\begin{itemize}
\item{Function \texttt{readmeteorologygrid()} reads the meteorological data stored in one \textbf{netCDF} file and returns a \texttt{SpatialGridDataframe}.}
\item{Function \texttt{writemeteorologygrid()} writes the meteorological data of a single date as a \textbf{netCDF} file.}
\item{Function \texttt{readmeteorologygridfiles()} reads several \textbf{netCDF} files and returns an object of class \texttt{SpatialGridMeteorology}.}
\item{Function \texttt{readmeteorologygridcells()} reads several \textbf{netCDF} files and returns an object of class \texttt{SpatialPointMeteorology} with the meteorological data of a set of specified grid cells.}
\item{Functions \texttt{writemeteorologygridfiles()} writes several \textbf{netCDF} files in the disk, one per date. Metadata (i.e. the dates and their corresponding file path) is stored in an additional file.}
\end{itemize}

\subsection{Meteorology estimation functions}
\subsubsection{Spatial interpolation}
Package \texttt{meteoland} provides two functions for interpolating meteorological data (i.e., one for each data structure):
\begin{itemize}
\item{Function \texttt{interpolationpoints()} interpolates weather for a set of locations given in \texttt{SpatialPointsTopography} and returns an object of class \texttt{SpatialPointsMeteorology}.}
\item{Function \texttt{interpolationgrid()} interpolates weather for a whole grid specified in \texttt{SpatialGridTopography} and returns an object of class \texttt{SpatialGridMeteorology}.}
\end{itemize}
Both functions require an object of class \texttt{`MeteorologyInterpolationData'}, which contains the X-Y coordinates, the meteorological data and topography of a set of weather stations as well as weather interpolation parameters.
<<echo=FALSE>>=
showClass("MeteorologyInterpolationData")
@
When calling functions \texttt{interpolationpoints()} or \texttt{interpolationgrid()}, the user may require interpolation outputs to be written into the file system, instead of being returned in memory. If \texttt{interpolationpoints()} is called with \texttt{export = TRUE}, the function will write the data frame produced for each point into an \textbf{ascii} text file. If \texttt{interpolationgrid()} is called with \texttt{export = TRUE}, the function will write an \textbf{netCDF} file for each day. Metadata files will also be written, so that results can later be loaded in memory.

\subsubsection{Statistical downscaling}
Analogously to interpolation, two functions are available for statistical downscaling of coarse-scale meteorological data (i.e., one function for each data structure):
\begin{itemize}
\item{Function \texttt{downscalingpoints()} performs statistical downscaling of coarse-scale weather data on a set of locations and it returns an object of class \texttt{SpatialPointsMeteorology} containing corrected weather predictions.}
\item{Function \texttt{downscalinggrid()} performs statistical downscaling of coarse-scale weather data over a grid and it returns an object of class \texttt{SpatialPointsMeteorology} containing corrected weather predictions.}
\end{itemize}
Donwscaling functions require an object of class \texttt{`MeteorologyDownscalingData'}, which contains the X-Y coordinates and the coarse-scale meteorological data to be downscaled, which includes a historical (reference) period and future (projected) period:
<<echo=FALSE>>=
showClass("MeteorologyDownscalingData")
@
The historical (reference) period is compared with fine-scale meteorological data of the same period, and the routine uses this information to correct the future (projected) period. Therefore, apart from the \texttt{`MeteorologyDownscalingData'} object, the dowscaling functions require fine-scale historical data (for a set of spatial points or a grid). Normally, these data will be the result of spatial interpolation.

As before, when calling functions \texttt{downscalingpoints()} or \texttt{downscalinggrid()}, the user may require the outputs to be written into the file system, instead of being returned in memory. If \texttt{downscalingpoints()} is called with \texttt{export = TRUE}, the function will write the data frame produced for each point into an \textbf{ascii} text file. If \texttt{downscalinggrid()} is called with \texttt{export = TRUE}, the function will write a \textbf{netCDF} file for each day. Metadata files will also be written, so that results can later be loaded in memory.


\section{Spatial interpolation of weather records}
\subsection{Overview}
Landscape research studies conducted for historical periods can be perfomed using meteorological records obtained from surface weather stations of the area under study. For any target point, minimum temperature, maximum temperature and precipitation are interpolated from weather records using truncated Gaussian filters, while accounting for the relationship between these variables and elevation (Thornton et al. 1997). Relative humidity can be either interpolated (in fact, water vapour pressure is the variable being interpolated) or predicted from temperature estimates, depending on whether it has been measured in weather stations or not. Potential (i.e. top-of-atmosphere) solar radiation is estimated taking into account latitude, seasonality, aspect and slope, following Granier \& Ohmura (1968). Potential solar radiation is then corrected to account for atmosphere transmittance using the predictions of temperature range, relative humidity and precipitation (Thornton \& Running 1999). Finally, the wind vector (wind direction and wind speed) is interpolated by using weather station records and static wind fields. In the following subsections we detail the general algorithm used to obtain interpolation weights and the interpolation/estimation procedure for each variable.

\subsection{Interpolation weights}
Thornton et al. (1997) suggested interpolating meteorological data using a truncated Gaussian filter. Its form with respect to a central point $p$ is:
\begin{equation}
W(r) = e^{-\alpha \cdot (r/R_p)^2} - e^{-\alpha}
\end{equation}
if $r \leq R_p$ and $W(r) = 0$ otherwise. Here $r$ is the radial distance from $p$, $R_p$ is the truncation distance and $\alpha$ is the shape parameter. The spatial convolution of this filter with a set of weather station locations results, for each target point, in a vector of weights associated with observations. The following figure illustrates the Gaussian filter for $R_p = 500$ and either $\alpha = 3.0$ (continuous line) or $\alpha = 6.25$ (dashed line):
\begin{center}
<<fig = TRUE, width=4, height=4, echo = FALSE>>=
r = 0:1000
R_p = 500
gf1 = exp(-3.0*((r/R_p)^2.0)) - exp(-3.0)
gf2 = exp(-6.25*((r/R_p)^2.0)) - exp(-6.25)
gf1[r>R_p] = 0
gf2[r>R_p] = 0
plot(r, gf1, type="l", ylab = "W(r)", xlab ="r")
lines(r, gf2, lty=2)
@
\end{center}
$R_p$ is estimated so that it has lower values in data-rich regions and is increased in data-poor regions. We require the user to specify $N$, the average number of observations to be included for each target point. $R_p$ is then varied as a smooth function of the local density in such a way that this average is achieved over the spatial domain. Estimation of $R_p$ is as follows:
\begin{enumerate}
\item{A user-specified value is used to initialize $R_p$.}
\item{Interpolation weights $W_i$ are calculated for all $i = (1, ..., n)$ stations, and the local station density is calculated as:
\begin{equation}
D_p = \frac{\sum_{i=1}^{n}{(W_i/\hat{W})}}{\pi \cdot R_p^2}
\end{equation}
where $\hat{W}$ is the average weight over the untruncated region of the kernel, calculated as:
\begin{equation}
\hat{W} = \left( \frac{1 - e^{-\alpha}}{\alpha}\right)- e^{-\alpha}
\end{equation}
}
\item{A new $R_p$ value is calculated as a function of $N$ and $D_p$, as:
\begin{equation}
R_p = \sqrt{\frac{N^*}{D_p \cdot \pi}}
\end{equation}
where $N^* = 2N$ for the first $I - 1$ iterations, and $N^* = N$ for the final iteration.
}
\item{ The new $R_p$ is substituted in step (2) and steps (2-4) are iterated a specified number of times $I$. The final $R_p$ value is used to generate interpolation weights $W_i$.}
\end{enumerate}
Thornton et al. (1997) suggested to use this algorithm only once per point (and variable to be estimated), but since missing meteorological values can occur only in some days, we apply the algorithm for each target point and day. The interpolation method for a given set of observations is defined by four parameters $R$, $I$, $N$ and $\alpha$. Following Thornton et al. (1997), we set $R = 140000$ meters and $I = 3$ by default. The other parameters depend on the variable to be interpolated.

\subsection{Temperature}
Predictions for minimum temperature and maximum temperature are done in the same way, so we refer to a general variable $T$. We focus on the prediction of $T_p$, the temperature at a single target point $p$ and for a single day, based on observations $T_i$ and interpolation weights $W_i$ for the $i = (1, ..., n)$ weather stations. Prediction of $T_p$ requires a correction for the effects of elevation differences between observation points $z_1, ..., z_n$ and the prediction point $z_p$. Thornton et al. (1997) established the relationship between elevation and temperature using transformed variables (temporal or spatial moving window averages) for temperature and elevation, instead of the original variables, but we did not implement this feature here. A weighted least-squares regression is used to assess the relationship between temperature and elevation. Instead of regressing $z_i$ on $T_i$, the independent variable is the difference in elevations associated with a pair of stations, and the dependent variable is the corresponding difference in temperatures. This gives a regression of the form:
\begin{equation}
(T_1 - T_2) = \beta_0 + \beta_1 \cdot (z_1 - z_2)
\end{equation}
where subscripts $1$ and $2$ indicate the two stations of a pair and $\beta_0$ and $\beta_1$ are the regression coefficients. Regression is performed using all possible pairs of stations and the regression weight associated with each point is the product of the interpolation weights associated with the stations in a pair. The temperature for the target point, $T_p$ is finally predicted as follows:
\begin{equation}
T_{p} = \frac{\sum_{i=1}^{n}{W_i\cdot (T_i + \beta_0 + \beta_1 \cdot(z_p - z_i))}}{\sum_{i=1}^{n}{W_i}}
\end{equation}
where $z_p$ is the elevation of the target point and $z_i$ is the elevation of the weather station.

\subsection{Relative humidity}
Relative humidity is a parameter not always recorded in weather stations. When input station weather data does not include relative humidity, \texttt{medfate} allows estimating it directly from minimum and maximum temperature (Thornton et al. 1997). Assuming that minimum daily air temperature $T_{min,p}$ at the target point is a good surrogate of dew-point temperature $T_{d,p}$ (i.e. $T_{d,p} = T_{min,p}$; note that this assumption may not be valid in arid climates), one can estimate actual vapor pressure $e_{p}$ (in Pa) as:
\begin{equation}
e_{p} = 610.78 \cdot e^{\left(\frac{17.269\cdot T_{d,p}}{237.3+T_{d,p}}\right)}
\end{equation}
and saturated vapor pressure $e_{s,p}$ (in Pa) as:
\begin{equation}
e_{s,p} = 610.78 \cdot e^{\left(\frac{17.269\cdot T_{a,p}}{237.3+T_{a,p}}\right)}
\end{equation}
where $T_{a,p} = 0.606 \cdot T_{max,p} + 0.394 \cdot T_{min,p}$ is the average daily temperature. Finally, relative humidity $RH_p$ (in percentage) is calculated as:
\begin{equation}
RH_p = 100 \cdot \frac{e_{p}}{e_{s,p}}
\end{equation}

When relative humidity has been measured at weather stations, interpolation should be preferred to estimation from minimum and maximum temperature. However, because relative humidity cannot be corrected for elevation differences between the station and the target climatic grid cell, relative humidity $RH_i$ of each weather station $i$ has to be converted to dew-point temperature $T_{d,i}$ before interpolation (Tymstra et al. 2010). To obtain the dew-point temperature one first needs to calculate vapor pressure:
\begin{equation}
e_i = e_{s,i} \cdot (RH_i / 100)
\end{equation}
where $e_{s,i}$ is the saturated water vapor pressure of station $i$, calculated as indicated above. Then, dew-point temperature of station $i$ is obtained from:
\begin{equation}
T_{d,i} = \frac{237.3\cdot \ln(e_i/610.78)}{17.269 - \ln(e_i/610.78)}
\end{equation}
As for temperature, a weighted least-squares regression is used to assess the relationship between dew-point temperature and elevation. Again, the independent variable is the difference in elevations associated with a pair of stations, and the dependent variable is the corresponding difference in dew-point temperature. This gives a regression of the form:
\begin{equation}
(T_{d,1} - T_{d,2}) = \beta_0 + \beta_1 \cdot (z_1 - z_2)
\end{equation}
where subscripts $1$ and $2$ indicate the two stations of a pair and $\beta_0$ and $\beta_1$ are the regression coefficients. The dew-point temperature for the target point, $T_{d,p}$ is predicted as:
\begin{equation}
T_{d,p} = \frac{\sum_{i=1}^{n}{W_i\cdot (T_{d,i} + \beta_0 + \beta_1 \cdot(z_p - z_i))}}{\sum_{i=1}^{n}{W_i}}
\end{equation}
where $z_p$ is the elevation of the target point and $z_i$ is the elevation of the weather station. From the interpolated dew-point temperature one can obtain actual vapour pressure $e_{p}$ and, together with saturated vapour pressure at point $p$, one calculates relative humidity as indicated above. If saturated vapour pressure is referred to average temperature $T_{a,p}$, then relative humidity is average relative humidity $RH_{a,p}$. If, instead, one refers saturated vapour pressure to minimum and maximum daily temperatures one obtains, respectively, maximum and minimum relative humidity values ($RH_{max,p}$, $RH_{min,p}$). Because of the regression step, one must check that the predicted relative humidity value stays within the physical limits 0\% and 100\%.

\subsection{Precipitation}
Predictions of precipitation are complicated by the need to predict both daily occurrence and, conditioned on this, daily total precipitation. Thornton et al. (1997) define a binomial predictor of spatial precipitation occurrence as a function of the weighted occurrence at surrounding stations. The precipitation occurrence probability $POP_p$ is:
\begin{equation}
POP_p = \frac{\sum_{i=1}^{n}{W_{o,i}\cdot PO_i}}{\sum_{i=1}^{n}{W_{o,i}}}
\end{equation}
where $PO_i$ is the binomial precipitation occurrence in station $i$ (i.e., $PO_i = 0$ if $P_i = 0$ and $PO_i = 1$ if $P_i > 0$) and $W_{o,i}$ is the interpolation weight for precipitation occurrence. Once $POP_p$ is calculated, then precipitation occurs if $POP_p$ is smaller than a critical value (i.e. $PO_p = 1$ if $POP_p < POP_{crit}$ and $PO_p = 0$ otherwise). Conditional on precipitation occurrence we calculate the prediction of daily total precipitation, $P_p$. Like with temperature, Thornton et al. (1997) established the relationship between elevation and precipitation using transformed variables (temporal or spatial moving window averages) for precipitation and elevation. Following their results, we transform precipitation values using a temporal window with side of 5 days. Weighted least-squares, where the weight associated with each point is the product of the interpolation weights associated with the stations in a pair, is used to account for elevation effects on precipitation. Unlike Thornton et al. (1997), who used the same set of interpolation weights (i.e. $W_{o,i}$) for precipitation occurrence and regression, we use a second set of interpolation weights $W_{r,i}$ for the calculation of regression weights. The dependent variable in the regression function is defined as the normalized difference of the precipitation observations $P_i$ for any given pair of stations:
\begin{equation}
\left(\frac{P_1 - P_2}{P_1 + P_2}\right) = \beta_0 + \beta_1 \cdot (z_1 - z_2)
\end{equation}
To obtain the predicted daily total $P_p$ we use the following equation:
\begin{equation}
P_p = \frac{\sum_{i=1}^{n}{W_{o,i}\cdot P_i \cdot PO_i \cdot \left(\frac{1 + f}{1 - f} \right)}}{\sum_{i=1}^{n}{W_{o,i} \cdot PO_i}}
\end{equation}
where $f = \beta_0 + \beta_1 \cdot (z_p - z_i)$. Note the usage of interpolation weight $W_{o,i}$ (and not $W_{r,i}$). The form of prediction requires that $\lvert f \rvert < 1$. A parameter $f_{max}$ (with default $f_{max} = 0.95$ ) is introduced to force $\lvert f \rvert  = f_{max}$ whenever $\lvert f \rvert  > f_{max}$.


\subsection{Radiation}
Incident daily solar radiation is not interpolated, but estimated from topography and measurements of temperature, humidity and precipitation.

Potential solar radiation is estimated from latitude, aspect and slope according to Granier \& Ohmura (1968). In particular, instant potential solar radiation $R_{pot,s}$ is calculated as:
\begin{eqnarray}
R_{pot,s} &=& I_0 [(\sin{\phi}\cos{H})(-\cos{A}\sin{Z_x})  -\sin{H}(\sin{A}\sin{Z_x}) \nonumber \\
 & & + (\cos{\phi}\cos{H})\cos{Z_x}]\cos{\delta} \nonumber \\
 & & + [\cos{\phi}(\cos{A}\sin{Z_x})+ \sin{\phi}\cos{Z_x}]\sin{\delta}
\end{eqnarray}
where $\phi$ is the latitude, $H$ is the hour angle measured from solar noon, positively towards the west, $A$ is the azimuth of the slope (aspect) and $Z_x$ is the zenith angle of the vector normal to the slope (equal to the slope angle) and $\delta$ is the sun's declination, which can be derived from the day of the year $J$ (Julian day):
\begin{equation}
\delta = 0.4093 \cdot \sin(2 \cdot \pi \cdot (J/365) - 1.405)
\end{equation}
Daily potential radiation, $R_{pot}$ is calculated by integrating instant potential radiation over the day between sunrise ($sr$) and sunset ($ss$), using 20 min (i.e. 1200 sec) intervals:
\begin{equation}
R_{pot} = \sum_{s = sr}^{ss}{1200 \cdot R_{pot,s}}
\end{equation}

Improving the method proposed in Thornton et al. (1997), Thornton \& Running (1999) calculate incident daily total solar radiation $R_{g}$ as:
\begin{equation}
R_g = R_{pot} \cdot T_{t,max} \cdot T_{f,max}
\end{equation}
where $T_{t,max}$ is the maximum (cloud-free) daily total transmittance and $T_{f,max}$ is the proportion of $T_{t,max}$ realized on a given day (cloud correction). The maximum daily total transmittance $T_{t,max}$ is estimated as:
\begin{equation}
T_{t,max} = \left[\frac{\sum_{s = sr}^{ss}{R_{pot,s} \cdot \tau^{(P_z/P_0)\cdot m_{\theta}}}}{\sum_{s = sr}^{ss}{R_{pot,s}}}\right] + (\alpha_{e_p} \cdot {e_p})
\end{equation}
where $\tau = 0.87$ is the instantaneous transmittance at sea level, at nadir, for a dry atmosphere; ${e_p}$ is the actual water vapor pressure (in $Pa$), estimated as explained before; $\alpha_{e_p} = -6.1\cdot 10^{-5} Pa^{-1}$ is a parameter describing the effect of vapour pressure on $T_{t,max}$; $m_{\theta} = 1/\cos{\theta}$ is the optical air mass at solar zenith angle $\theta = \sin{\phi}\cdot \sin{\delta}+\cos{\phi}\cdot \cos{\delta} \cdot \cos{H}$; and $P_z/P_0$ is the ratio between air pressure at elevation $z_p$ and air pressure at the sea level, calculated as:
\begin{equation}
(P_z/P_0) = (1.0 -2.2569\cdot 10^{-5}\cdot z_p)^{5.2553}
\end{equation}
In turn, $T_{f,max}$ was empirically related to $\Delta T = T_{\max} - T_{\min}$, the difference between maximum and minimum temperatures for the target point:
\begin{equation}
T_{f,max} = 1.0 - 0.9\cdot e^{-B \cdot {\Delta T}^{C}}
\end{equation}
being $C = 1.5$ and $B$ calculated from:
\begin{equation}
B = b_0 + b_1 \cdot e^{-b_2 \cdot {\hat{\Delta T}}}
\end{equation}
with $b_0 = 0.031$, $b_1 = 0.201$ and $b_2 = 0.185$. In this last equation,  $\hat{\Delta T}$ is a 30-day moving average for the temperature range $\Delta T$. For computational reasons, we do not estimate $\hat{\Delta T}$ from the 30-day moving window average of predicted $\Delta T$ values, but from the interpolation of pre-calculated $\hat{\Delta T}$ values in weather stations. On wet days (i.e. if $P_p > 0$) the estimation of $T_{f,max}$ is multiplied by a factor of $0.75$ to account for clouds.


\subsection{Wind}
\subsubsection{Simple wind interpolation}
Simple wind interpolation for a target point $p$ is as follows. Let $\mathbf{v}_{i}$ be the wind vector in weather station $i$. $\mathbf{v}_{i}$ is initially expressed using polar coordinates. Indeed, we have $u_{i}$ and $\theta_{i}$, the wind speed and wind direction, respectively. If we express $\mathbf{v}_{i}$ in cartesian coordinates we have:
\begin{equation}
x_{i} = u_{i}\cdot \sin(\theta_{i}) \quad y_{i} = u_{i}\cdot \cos(\theta_{i})
\end{equation}
The predicted wind vector $\mathbf{v}_{p}$ is the weighted average of the wind vectors $\{\mathbf{v}_{i}\}$  $i=(1,... , n)$ predicted for point $p$ using the interpolation weights $W_i$ determined from the truncated Gaussian filter:
\begin{equation}
x_{p} = \frac{\sum_{i=1}^{n}{W_i \cdot x_{i}}}{\sum_{i=1}^{n}{W_i}}  \quad y_{p} = \frac{\sum_{i=1}^{n}{W_i \cdot y_{i}}}{\sum_{i=1}^{n}{W_i}}
\end{equation}
The polar coordinates of the predicted wind vector $\mathbf{v}_{p}$  are:
\begin{equation}
u_{p} = \sqrt{x^{2}_{p} + y^{2}_{p}} \quad \theta_{p} = \tan^{-1}(x_{p}/y_{p})
\end{equation}
\subsubsection{Interpolation using wind fields}
More precise wind interpolation requires a set of static wind fields covering the landscape of interest. Each of these wind fields has been calculated assuming a domain-level combination of wind speed and wind direction. The set of domain-level combinations should cover all possible winds in the landscape under study. For example, one could decide to include the combinations of eight different wind directions (i.e., N, NE, E, SE, ...) and three wind speed classes. The wind estimation of a given target point depends on both the wind observations at weather stations and these static wind fields.

In a given day (and before processing target points) we begin by identifying, for each weather station $i = (1, ... ,n)$, the wind field $m_i$ corresponding to a minimum difference between the observed wind vector $\mathbf{v}_i$ and the wind vector of the station in the wind field (i.e., minimum distance between the corresponding cartesian coordinates). The set of wind fields $\{m_i\}$ $i = (1, ... ,n)$ chosen for each weather station conform the information for wind interpolation in a given day.

Actual wind interpolation details for a target point $p$ are as follows. We first draw for each $i=(1,... , n)$ the wind vector $\mathbf{v}_{m_i,p}$ corresponding to the location of the target point $p$ in wind fields $m_i$. Let $u_{m_i,p}$ and $\theta_{m_i,p}$ be the wind speed and wind direction of $\mathbf{v}_{m_i,p}$, respectively. The cartesian coordinates of $\mathbf{v}_{m_i,p}$ are:
\begin{equation}
x_{m_i,p} = u_{m_i,p}\cdot \sin(\theta_{m_i,p}) \quad y_{m_i,p} = u_{m_i,p}\cdot \cos(\theta_{m_i,p})
\end{equation}
The predicted wind vector $\mathbf{v}_{p}$ is the weighted average of the wind vectors $\{\mathbf{v}_{m_i,p}\}$  $i=(1,... , n)$ predicted for point $p$ using the interpolation weights $W_i$ determined from the truncated Gaussian filter:
\begin{equation}
x_{p} = \frac{\sum_{i=1}^{n}{W_i \cdot x_{m_i,p}}}{\sum_{i=1}^{n}{W_i}}  \quad y_{p} = \frac{\sum_{i=1}^{n}{W_i \cdot y_{m_i,p}}}{\sum_{i=1}^{n}{W_i}}
\end{equation}
The polar coordinates of the predicted wind vector $\mathbf{v}_{p}$  are found as before.

\section{Statistical downscaling of coarse-scale weather data}

Statistical downscaling is necessary when meteorological data is available at a spatial scale that is too coarse for landscape-level analysis. This is usually the case when taking predictions from global or regional climate models. The general idea of downscaling to the landscape level is that a fine-scale meteorological series is to be compared to coarse-scale series for the a historical (reference) period. The result of this comparison can be used to correct coarse-scale meteorological series for other periods (normally future projections). In the case of most meteorological variables, the comparison consists in calculating a bias using the reference period and applying this bias to the future period. However, in the case of precipitation a different procedure is needed.

\subsection{Precipitation}
TO DO: Describe precipitation quantile mapping (Gudmundsson et al. 2012).

\subsection{Other variables}
Other meteorological variables (temperature, relative humidity, radiation and wind) are downscaled by simple bias correction applied monthly. Let $x$ be the meteorological variable observed for the target point in the reference period, and $u$ the same variable, but estimated at the coarse-scale for the reference period. Month $m$ bias of the coarse-scale data ($\theta_{m}$) is the average of the difference between $u$ and $x$ over all $n_m$ days of month $m$ in the reference period:
\begin{equation}
\theta_{m} = \sum_{i=1}^{n_m}{(u_i-x_i)}/n_m
\end{equation}
Month biases are used to correct the coarse-scale value of day $i$ in the projected period using:
\begin{equation}
c_i = u_i - \theta_{m(i)}
\end{equation}
where $m(i)$ is the month of day $i$.

Mean temperature, minimum temperature, maximum temperature, radiation and wind speed are downscaled using this general bias correction procedure. In the case of relative humidity, values are transformed to specific humidity, bias correction is applied to this variable and the result is back-transformed to relative humidity. Since historic wind data is often not available, if wind speed data is missing the coarse-scale wind estimate is taken directly without correction.

\section{Miscellaneous functions}
\subsection{Extraction of climatic netCDF data}

\textbf{NetCDF} is a standard data format for meteorological data. In particular, this format is used to store the predictions of global and regional climate models. Function \texttt{extractNetCDF()} parses a set of \textbf{NetCDFs} and extracts the daily meteorological data of a landscape boundary box for use within \texttt{meteoland}. The function first identifies which cells in NetCDF data should be extracted (according to the input boundary box), and the overall period. For each cell to be processed, the function loops over all files (which can describe different variables and time periods) and extracts the corresponding data. The function transforms units to the units used in meteoland. If specific humidity and mean temperature are available, the function also calculates mean relative humidity.

\subsection{Potential evapo-transpiration}
Package \texttt{meteoland} provides a supplementary function to calculate potential evapo-transpiration (PET) using Penman's formulation (Penman 1948; 1956). The code  was taken from package \texttt{`Evapotranspiration'}, which follows McMahon et al. (2013). Penman (1948) proposed an equation to calculate potential evaporation that combined an energy equation based on net incoming radiation with an aerodynamic approach. The Penman or Penman combination equation is:
\begin{equation}
E_{pot} = \frac{\Delta}{\Delta+\gamma}\cdot \frac{R_n}{\lambda}+\frac{\lambda}{\Delta + \lambda}\cdot E_a
\end{equation}
where $PET$ is the daily potential evaporation (in $mm \cdot day^{-1}$) from a saturated surface, $R_{n}$ is the daily radiation to the evaporating surface (in $MJ\cdot m^{-2}\cdot day^{-1}$), $\Delta$ is the slope of the vapour pressure curve ($kPa\cdot \degree C^{-1}$) at air temperature, $\gamma$ is the psychrometric constant ($kPa\cdot \degree C^{-1}$), and $\lambda$ is the latent heat of vaporization (in $MJ\cdot kg^{-1}$). $E_a$  (in $mm \cdot day^{-1}$) is a function of the average daily windspeed ($u$, in $m\cdot s^{-1}$), and vapour pressure deficit ($D$, in $kPa$):
\begin{equation}
E_a = f(u) \cdot D = f(u) \cdot (v_a^*-v_a)
\end{equation}
where $v_a^*$ is the saturation vapour pressure ($kPa$) and $v_a$ the actual vapour pressure ($kPa$) and $f(u)$ is a function of wind speed, for which there are two alternatives (Penman 1948; 1956):
\begin{eqnarray}
f(u) &=& 1.313 + 1.381 \cdot u\\
f(u) &=& 2.626 + 1.381 \cdot u
\end{eqnarray}
If wind speed is not available, an alternative formulation for $E_{pot}$ is used as an approximation (Valiantzas 2006):
\begin{equation}
PET \simeq 0.047\cdot R_s \cdot (T_a+9.5)^{0.5}-2.4\cdot (\frac{R_s}{R_a})^2+0.09\cdot(T_a-20)\cdot(1-\frac{RH_{mean}}{100})
\end{equation}
where $R_s$ is the incoming solar radiation (in $MJ\cdot m^{-2}\cdot day^{-1}$), $T_a$ is the mean daily temperature (in $\degree C$), $R_a$ is the extraterrestrial solar radiation (in $MJ\cdot m^{-2}\cdot day^{-1}$) and $RH_{mean}$ is the mean relative humidity (in percent). Details about the calculation of all these quantities can be found in McMahon et al. (2013).

PET is normally calculated after meteorological data have been interpolated (i.e. within functions \texttt{interpolationpoints()} and \texttt{interpolationgrid()}) or downscaled (i.e. within functions \texttt{downscalingpoints()} and \texttt{downscalinggrid()}), but PET series can also be calculated for a single point using function \texttt{penmanpoint()}. For other formulations of PET, the reader is referred to package \texttt{`Evapotranspiration'}.

\subsection{Obtaining static wind fields}
External software is necessary to calculate the set of wind fields for the study area under different domain-level average situations. For this we recommend using \textbf{\href{http://firelab.org/project/windninja}{WindNinja}}, a computer program that calculates spatially varying wind fields for wildland fire applications. WindNinja allows simulating the spatial variation of wind for one instant in time. It was developed to be used by emergency responders within their typical operational constraints of fast simulation times (seconds), low CPU requirements (single processor laptops), and low technical expertise. WindNinja is typically run on domain sizes up to 50 kilometers by 50 kilometers and at resolutions of around 100 meters. The program is free and can be downloaded from \href{http://www.firemodels.org}{www.firemodels.org}.

The inputs for a basic run of WindNinja are an elevation data file for the study area, a domain-averaged input wind speed and direction and a specification of the dominant vegetation in the area. In order to obtain a set of pre-computed rasters of wind direction and speed, we suggest the following procedure:
\begin{itemize}
\item{Export the elevation raster of the study area in one of the file formats accepted by WindNinja (`.asc', `.tif' or `.img'). In the case of a large study area (e.g. > 100 x 100 km) one should run WindNinja in subsets of the area and then integrate the results (e.g., Sanjuan et al. 2014).}
\item{Run WindNinja, using the elevation of the study area, for all combinations of wind direction and wind speed class (for each wind speed class an mean class value has to be chosen). Several combinations of domain-level wind speed and wind direction can be specified for a single run, and the program can also be run in batch mode.}
\item{Read raster files created by WindNinja (a wind speed file, a wind direction file) for each combination of domain-level wind speed and direction. }
\end{itemize}
Function \texttt{readWindNinjaOutput()} can be used to conduct this last step. The function allows parsing all the ASCII raster files produced by WindNinja for combinations of wind direction (e.g., 0, 45, 90, 135, 180, 225, 270 and 315 degrees) and wind speed (e.g., 2, 7 and 10 m/s). The function returns a list with the following elements:
\begin{itemize}
\item{The vector of domain-level wind directions corresponding to WindNinja Runs}
\item{The vector of domain-level wind speed corresponding to WindNinja Runs}
\item{An object \texttt{SpatialGridDataFrame} containing the wind directions (in degrees from North) for all WindNinja runs.}
\item{An object \texttt{SpatialGridDataFrame} containing the wind speeds (in m/s) for all WindNinja runs.}
\end{itemize}


\section{References}
\begin{itemize}

\item{Garnier, B.J., Ohmura, A., 1968. A method of calculating the direct shortwave radiation income of slopes. J. Appl. Meteorol. 7, 796–800. doi:10.1175/1520-0450(1968)007<0796:AMOCTD>2.0.CO;2.}

\item{Gudmundsson L, Bremnes JB, Haugen JE, Engen-Skaugen T (2012) Technical Note: Downscaling predicted climatic precipitation to the station scale using statistical transformations - A comparison of methods. Hydrology and Earth System Sciences 16, 3383–3390. doi:10.5194/hess-16-3383-2012.}

\item{McMahon, T.A., Peel, M.C., Lowe, L., Srikanthan, R., McVicar, T.R. 2013. Estimating actual, potential, reference crop and pan evaporation using standard meteorological data: a pragmatic synthesis. Hydrology \& Earth System Sciences 17, 1331–1363. doi:10.5194/hess-17-1331-2013.}

\item{Penman, H. L. 1948. Natural evaporation from open water, bare soil and grass. Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences, 193, 120-145.}

\item{Penman, H. L. 1956. Evaporation: An introductory survey. Netherlands Journal of Agricultural Science, 4, 9-29.}

\item{Sanjuan, G., Brun, C., Cort, A., 2014. Wind field uncertainty in forest fire propagation prediction. Procedia Comput. Sci. 29, 1535–1545. doi:10.1016/j.procs.2014.05.139.}

\item{Thornton, P.E., Running, S.W., White, M. a., 1997. Generating surfaces of daily meteorological variables over large regions of complex terrain. J. Hydrol. 190, 214–251. doi:10.1016/S0022-1694(96)03128-9.}

\item{Thornton, P.E., Running, S.W., 1999. An improved algorithm for estimating incident daily solar radiation from measurements of temperature, humidity, and precipitation. Agric. For. Meteorol. 93, 211–228. doi:10.1016/S0168-1923(98)00126-9.}
\item{Tymstra, C., Bryce, R.W., Wotton, B.M., Taylor, S.W., Armitage, O.B., 2010. Development and Structure of Prometheus : the Canadian Wildland Fire Growth Simulation Model, Information report NOR-X-417. doi:ISBN 978-1-100-14674-4}
\item{Valiantzas J.D. 2006. Simplified versions for the Penman evaporation equation using routine weather data. Journal of Hydrology 331, 690–702. doi:10.1016/j.jhydrol.2006.06.012.}
\end{itemize}


\end{document}
