water: Tools and Functions to Estimate Actual Evapotranspiration Using Land Surface Energy Balance Models in R

1 Introduction and motivation

The crop water requirement is a key factor in the agricultural process. It is usually estimated throughout actual evapotranspiration ($E{T}_a$). An accurate quantification of $E{T}_a$ helps to develop irrigation strategies, improve the efficiency of water use and increase the irrigated area and the production (Millar 1993), (Baruch and M. Fisher 1991), (Ferreyra, G. Sellés, and J. Tosso 1985).

Traditional methods to estimate $E{T}_a$ are based on (a) direct measurements by sophisticated instruments, such as lysimeters (Payero and S. Irmak 2008), (López-Urrea, a. Montoro, P. López-Fuster, and E. Fereres 2009), Eddy covariance systems (Paço, M. Ferreira, and N. Conceição 2006), (Parent and F. Anctil 2012),(Poblete-Echeverría and S. Ortega-Farias 2013) or Bowen ratios (Cragoa and W. Brutsaert 1996), (Ortega-Farı́as, R. Cuenca, and M. English 1995), (Twine, W. Kustas, J. Norman, D. Cook, P. Houser, T. Meyers, J. Prueger, P. Starks, and M. Wesely 2000), or on (b) empirical methods, such as the FAO-56 approach (Allen, L. S. Pereira, D. Raes, and M. Smith 1998). This method uses a reference evapotranspiration ($E{T}_r$) from an automatic weather station multiplied by crop coefficients ($K_c$) from literature (Allen, I. A. Walter, E. R., T. Howell, D. Itenfisu, and M. Jensen 2005). Although all these methods can be accurate enough, they are restrictive to be extrapolated to a farm or a regional level since they do not take into account the effect that the spatial and temporal variation of the soil, the climate and the crop have over the $E{T}_a$ (Allen, A. Irmak, R. Trezza, J. M. H. Hendrickx, W. Bastiaanssen, and J. Kjaersgaard 2011).

However, new physical methods to estimate $E{T}_a$ have been developed using remote sensing data, considering the land spatial and temporal patterns. A major restriction for the estimation of $E{T}_a$ using remote sensing is the need of an absolute surface temperature for calibration. One of the first methods to be developed and applied worldwide was the Surface Energy Balance Algorithms for Land (SEBAL) model (Bastiaanssen, M. Menenti, R. Feddes, and A. Holtslag 1998),(Bastiaanssen, M. Menenti, R. a. Feddes, and a. a. M. Holtslag 1998). This model calculates $E{T}_a$ from satellite-based land surface energy balance (LSEB) equation and uses a near-surface temperature gradient ($dT$) for calibration. $dT$ is computed by taking two pixels with extreme water condition (anchor pixels) selected in the scene to generate a linear relationship between surface temperature and the difference between surface and air temperatures. Based on this, (Allen, M. Tasumi, and R. Trezza 2007) developed the Mapping EvapoTranspiration at High Resolution with Internalized Calibration (METRIC) model. The main difference between SEBAL and METRIC is that the latter uses the $E{T}_r$ from a weather station, incorporating climatic conditions, while SEBAL uses the potential evaporation from a water body in the scene considering that sensible heat and soil heat fluxes are zero.

METRIC has been widely applied to estimate $E{T}_a$ at field and regional scale over different crops such as wheat, corn, soybean and alfalfa, with errors ranging between $3$ and $20\mathrm{\%}$ (Allen, M. Tasumi, A. Morse, R. Trezza, J. L. Wright, W. Bastiaanssen, W. Kramber, I. Lorite, and C. W. Robison 2007), (Choi, W. P. Kustas, M. C. Anderson, R. G. Allen, F. Li, and J. H. Kjaersgaard 2009), (Mkhwanazi and J. L. Chávez 2012). In recent years METRIC has been used to compute $E{T}_a$ over sparse woody canopies such as vineyards and olive orchards (Carrasco-Benavides, S. Ortega-Farı́as, L. O. Lagos, J. Kleissl, L. Morales, C. Poblete-Echeverrı́a, and R. G. Allen 2012), (Carrasco-Benavides, S. Ortega-Farı́as, L. Lagos, J. Kleissl, L. Morales-Salinas, and A. Kilic 2014), (Santos, I. J. Lorite, R. G. Allen, and M. Tasumi 2012), (Pôças, T. a. Paço, M. Cunha, J. a. Andrade, J. Silvestre, A. Sousa, F. L. Santos, L. S. Pereira, and R. G. Allen 2014) in both flat and mountainous terrains (Allen, R. Trezza, A. Kilic, M. Tasumi, and H. Li 2013).

The current implementations of SEBAL and METRIC imply the need to use more than one software to run the model (Allen, T. M., R. Trezza, and J. Kjaersgaard 2010) and they involve multiple steps. There are many different sub-models published for the estimation of some parameters (e.g. leaf area index, momentum roughness length, land surface temperature) that are not integrated into the current implementation of METRIC. (Allen, B. Burnett, W. Kramber, J. Huntington, J. Kjaersgaard, A. Kilic, C. Kelly, and R. Trezza 2013) proposed a methodology for an automation procedure by using statistical conditions and expert knowledge. This technique reduced the effect of human criteria helping to increase the model robustness. However, a software tool for the automatic selection of anchor pixels has not been published yet.

About land surface energy balance models and crop evapotranspiration

As mentioned above, METRIC estimates $E{T}_a$ as the residual from the surface energy balance equation considering information from satellite images and weather stations located near to the study site. Bellow, the key equations are detailed, beginning with the estimation of $E{T}_a$ as the residual from the surface energy balance equation:

$$\begin{array}{}\text{(1)}& LE=R_n-G-H\end{array}$$

where $LE$ is latent heat flux consumed by $E{T}_a$ ($W\cdot m^{-2}$); $R_n$ is net radiation ($W\cdot m^{-2}$); $G$ is soil heat flux ($W\cdot m^{-2}$); and $H$ is the sensible heat flux convected to the air ($W\cdot m^{-2}$).

$R_n$ is calculated considering information obtained at the time of satellite overpass. Some correction processes are necessary, such as radiometric and atmospheric corrections. $G$ is estimated using an empirical equation that considers mainly $R_n$, surface temperature, normalized difference vegetation index (NDVI), soil-adjusted vegetation index (SAVI) and albedo. More detailed information concerning the equations and models used in METRIC can be found in (Allen, M. Tasumi, and R. Trezza 2007).

$H$ is the general equation of heat transport and is estimated using an approach called “calibration using inverse modeling at extreme conditions” (CIMEC) (Allen, B. Burnett, W. Kramber, J. Huntington, J. Kjaersgaard, A. Kilic, C. Kelly, and R. Trezza 2013). This method involves the selection of pixels with near extreme conditions (hot and cold anchor pixels) from which the $E{T}_a$ can be estimated and assigned. $H$ is computed as follows:

$$H=\frac{\rho \cdot c_p\cdot dT}{r_{ah}}$$

where $dT$ is the difference between land surface and near-surface air temperatures, $r_{ah}$ is the aerodynamic resistance to heat transport ($s\cdot m^{-1}$), $\rho$ is the air density ($kg\cdot m^{-3}$) and $c_p$ is the specific heat of air ($1004\cdot J\cdot k{g}^{-1}{\cdot }^{\circ }{K}^{-1}$). $dT$ is solved by using a linear relationship between air temperature and the estimated surface temperature of the the anchor pixels (Bastiaanssen, M. Menenti, R. Feddes, and A. Holtslag 1998). To calculate $r_{ah}$, wind speed is extrapolated to a height at which forces of buoyancy and mechanical mix are equal (about 200 meters), using an iterative correction process based on the Monin-Obhukov equations (Allen 1996), (Bastiaanssen, M. Menenti, R. Feddes, and A. Holtslag 1998).

After $LE$ from Equation (1), it is possible to compute the instantaneous evapotranspiration values:

$$\begin{array}{}\text{(2)}& E{T}_{inst}=3600\cdot \frac{LE}{\lambda {\rho }_w}\end{array}$$

where $E{T}_{inst}$ is the instantaneous $E{T}_a$ at the satellite overpass ($mm\cdot h^{-1}$); $3600$ is the conversion factor from seconds to hours; ${\rho }_w$ is the density of water ($1000kg\cdot m^{-3}$); and $\lambda$ is the water latent heat of vaporization ($J\cdot k{g}^{-1}$).

Finally, the daily ET is computed pixel by pixel as: $$\begin{array}{}\text{(3)}& E{T}_{24}=\frac{E{T}_{inst}}{E{T}_r}E{T}_{r\mathrm{\_}24}\end{array}$$

where $E{T}_{inst}$ in the instantaneous $E{T}_a$ estimated on equation (2); $E{T}_r$ is the standardized 0.5 m tall alfalfa reference evapotranspiration at the image time and $E{T}_{r\mathrm{\_}24}$ is the cumulative 24 h $E{T}_r$ for the image day. The relationship between $E{T}_{inst}$ and $E{T}_r$ is the reference ET fraction and is the same as the alfalfa based coefficient, $K_c$, and is used to extrapolate $E{T}_a$ from the image time to periods of 24 hours or longer (Allen, M. Tasumi, and R. Trezza 2007).

Sensible heat flux

As (Allen, M. Tasumi, and R. Trezza 2007) mentioned, the computation of latent heat flux ($LE$) is only as accurate as the summed estimates for $R_n$, $G$, and $H$. Table 1 shows some errors reported by METRIC models for different crops. It can be seen that net radiation ($R_n$) and soil heat flux ($G$) present the lowest estimation errors, while sensible heat flux ($H$), is the hardest component of the surface energy balance to estimate.

One of the weaknesses of METRIC model reported in the literature is the selection of the anchor pixels. (Long and V. P. Singh 2013) and (Morton 2013) indicated that the selection of anchor pixels is subjective and depends on the ability of the operator to search and isolate the most appropriates hot and cold pixels. This process produces important biases in the estimation of $H$. Also, (Choragudi 2011), (Wang, T. Sammis, V. Gutschick, M. Gebremichael, and D. Miller 2009) mentioned that METRIC was very sensitive to the selection of the hot pixel. A group of possible candidates could have minimal differences in some attributes, but these can generate a big bias in the estimations. It means that the estimation of $H$ in METRIC is very sensitive to the selection of anchor pixels.

Table 1: Root mean square error (RMSE) of energy balance components estimated using METRIC for different crops

Crop	Validation tool	$R_n$ ($W{m}^{-2}$(%))	$G$ ($W{m}^{-2}$(%))	$H$ ($W{m}^{-2}$(%))	$LE$ ($W{m}^{-2}$(%))	$ET$ ($mm{h}^{-1}$(%))	Source
grass	lysimeter	nr	nr	nr	nr	($4$-$22\mathrm{\%}$) (mean $=4\mathrm{\%}$)	$1$
sugar beet	lysimeter	nr	nr	nr	nr	($6$-$137\mathrm{\%}$) (mean $=1\mathrm{\%}$)	$1$
soybean	lysimeter	$22.1$ ($4.1\mathrm{\%}$)	$14.2$ ($27.6\mathrm{\%}$)	nr	nr	$0.14$ ($17.6\mathrm{\%}$)	$2$
corn, soybean	EC	nr	nr	$39$-$48$	$34$-$44$	$0.58$-$0.89$	$3$
olive	EC	nr	nr	nr	nr	$0.14$-$1.2$	$4$
vineyard	EC	$24$ ($3.8\mathrm{\%}$)	$16$ ($9.4\mathrm{\%}$)	$39$-$59$ ($10$-$26.0\mathrm{\%}$)	$33$-$54$ ($14$-$27.2\mathrm{\%}$)	$0.9$ ($9\mathrm{\%}$)	$5$

nr: not reported; EC: Eddy covariance

Sources: 1: (Allen, M. Tasumi, A. Morse, R. Trezza, J. L. Wright, W. Bastiaanssen, W. Kramber, I. Lorite, and C. W. Robison 2007); 2: (Mkhwanazi and J. L. Chávez 2012),(Mkhwanazi and J. L. Chávez 2013); 3: (Gonzalez-Dugo, C. Neale, L. Mateos, W. Kustas, J. Prueger, M. Anderson, and F. Li 2009); 4: (Santos, I. J. Lorite, R. G. Allen, and M. Tasumi 2012),(Pôças, T. a. Paço, M. Cunha, J. a. Andrade, J. Silvestre, A. Sousa, F. L. Santos, L. S. Pereira, and R. G. Allen 2014); 5: (Poblete-Echeverría and S. Ortega-Farias 2012), (Carrasco-Benavides, S. Ortega-Farı́as, L. O. Lagos, J. Kleissl, L. Morales, C. Poblete-Echeverrı́a, and R. G. Allen 2012), (Carrasco-Benavides, S. Ortega-Farı́as, L. Lagos, J. Kleissl, L. Morales-Salinas, and A. Kilic 2014)

The objective of this article is to propose an open implementation of a land surface energy balance model (LSEB) as an R package that integrates most of the METRIC sub-models and allows automatic selection of the anchor pixels. In this version of the package, specific functions for METRIC model (Allen, M. Tasumi, and R. Trezza 2007) are provided. Apart from the previous features, this package: i) is written in R, one of the most used scientific programming languages; ii) can be automated for batch processing of many satellite images; iii) provides functions for loading and processing satellite images; iv) provides functions and a new class object to manage weather station data; and v) is a free and open software.

2 About the water package

Package organization

The water package is developed in R to estimate actual evapotranspiration ($E{T}_a$) from Landsat satellite scenes using Land Surface Energy Balance (LSEB) models, such as METRIC (Allen, M. Tasumi, and R. Trezza 2007).

Functions in water package are arranged in three groups: i) general functions to estimate sub-components of LSEB (e.g. leaf area index, albedo, land surface temperature, momentum roughness length); ii) specific functions to estimate the components of LSEB; iii) internal functions and methods to handle data from a weather station using a new proposed S3 class and functions to control global options such as saving results to disk, overwrite files, etc. 

The first group of functions consist of models and equations to estimate the sub-components, which generally are controlled by the argument method. Most of the models available here are presented in the Appendix. In the second group, there are three functions: i) METRIC.Rn(); ii) METRIC.G() and iii) METRIC.EB(). The first one estimates net radiation using METRIC model. The second one estimates soil heat flux, and the third one estimate all the components of energy balance: $R_n$, $G$, $H$ and $LE$ according to METRIC model.

Three example datasets are provided with the water package: i) a subset of a Landsat 7 scene path 223, row 85 from 15th February 2013, bands 1 to 7; ii) data from a weather station from the same Landsat subset in CSV file format (comma separated values); iii) a subset of NASA SRTM digital elevation model, with the same spatial extent as the example image. These datasets are used in examples in Sections 3.2 and 3.3, and also in the vignettes included with the package.

The water package is available on The Comprehensive R Archive Network at
https://CRAN.R-project.org/package=water. This software is made freely available under the terms and conditions of the GNU General Public License.

Key functions in the water package

The key functions included in the package are:

read.WSdata()

This function allows to import weather station data from a table in comma-separate values (CSV) format. The result is a new object of class "waterWeatherStation". The main input arguments are the CSV file and a vector with the order of the needed variables (radiation, temperature, wind speed and relative humidity) called columns. An optional argument is a Landsat metadata file (MTL). When this function is used with the CSV and the MTL files, it will interpolate the weather conditions to the moment of the satellite overpass.

METRIC.EB()

This is the main function of water package. It runs each of the sub-models needed to get all the components of the LSEB (equation (1)), from satellite and weather station data. The input arguments are: a satellite scene and a "waterWeatherStation" object. The arguments alb.coeff, LAI.method, Zom.method and anchors.method allow choosing between the different sub-models or coefficients used. More information about the sub-models available in water is presented in the Appendix. An optional logical argument is plain, which allows to use a digital elevation model or to consider that the surface is flat. When using a digital elevation model, the net radiation estimation and the land surface temperature are corrected using the elevation, slope and aspect of the surface. The output of this function is a raster layer object (from raster package) with 4 different layers: net radition, soil heat flux, sensible heat flux and surface temperature.

ET24h()

This function estimates the $E{T}_r$ for a 24-hour period. The input arguments are the components of the LSEB ($R_n$, $G$, $H$) and a "waterWeatherStation" object. The argument ET allows to select between two $E{T}_r$ methods. ET="ETo", is the method for short crops, similar to clipped, cool-season grass and ET="ETr", is the method for tall crops, similar to 0.5 m tall full-cover alfalfa. By default, it will use ET="ETr", but the user should choose according to the conditions of the weather station. The output of this function is a raster layer object (from raster package) with the 24-hour $E{T}_a$ in $mm\cdot da{y}^{-1}$.

METRIC.EB() uses many different steps to estimate the parameters needed to calculate the components of the LSEB. These steps are available as individual functions like:

albedo()

This function is used to calculate the broadband albedo from narrowband satellite data. This process involves applying a weighting function with empirical coefficients. water includes models and coefficients described by (Tasumi, R. G. Allen, and R. Trezza 2008) and (Liang 2001), as well as new coefficients for Landsat 8 estimated by The Simple Model of the Atmospheric Radiative Transfer of Sunshine (SMARTS2) version 2.9.5 (Gueymard 1995). coeff="Olmedo" computes clear sky spectral irradiances for the spectral range of each Landsat 8 OLI band (see Appendix, Section 6.1). The output of this function is a raster layer object (from raster package) with the albedo.

LAI()

This function estimates the leaf area index (LAI) from Landsat data. water includes many different models to estimate LAI (see Appendix, Section 6.2). In the Examples Section of this article, the METRIC 2010 method (Allen, T. M., R. Trezza, and J. Kjaersgaard 2010) will be used. This method estimates LAI from SAVI (Huete 1988), as follows:

$$SAVI=(1+L)({\rho }_{NIR}-{\rho }_R)/(L+{\rho }_{NIR}+{\rho }_R)$$

where $\rho$ is the reflectance at the top-of-atmosphere from $R$, the red band, and from $NIR$, the near infrared band; $L$ is a soil correction factor. By default water uses L=0.5, although (Allen, M. Tasumi, and R. Trezza 2007) suggested a value of L=0.1 in METRIC applications for western USA. This value varies by the amount or coverage of green vegetation: in very high coverage vegetation regions, L=0 and SAVI = NDVI; in areas with no green vegetation, L=1. Then SAVI is used to estimate LAI as follows:

$$LAI=11\cdot SAV{I}^3$$

The output of this function is a raster layer object (from raster package) with the estimated $LAI$.

As it was mentioned before, the critical part in the estimation of the LSEB is the sensible heat flux estimation. The key functions used to estimate this are: momentumRoughnessLength(); calcAnchors() and calcH().

momentumRoughnessLength()

This function estimates the Momentum Roughness Length ($Z_{om}$) from the average vegetation height around the weather station. water includes several methods to estimate $Z_{om}$ from Landsat data (see Appendix, Section 6.4). For example, when method="short-crops" (Allen, M. Tasumi, and R. Trezza 2007) $Z_{om}$ is estimated as:

$$Z_{om}=0.0018\text{LAI}$$

And if mountainous=TRUE, this value of $Z_{om}$ is corrected as:

$$Z_{om\mathrm{\_}mtn}=Z_{om}\cdot (1+\frac{(180/\pi )\cdot \text{slope}-5}{20})$$

The output of this function is a raster layer object with the estimated $Z_{om}$

calcAnchors()

This function automatically select anchor pixels that represent the dry and wet ends of the ET spectrum within the satellite scene using information of land surface characteristics (LAI, albedo, $Z_{om}$, $T_s$). When the anchor pixels are found, it assigns the $E{T}_a$ estimates.

The criteria to select “hot” and “cold” pixels were adapted from (Allen, A. Irmak, R. Trezza, J. M. H. Hendrickx, W. Bastiaanssen, and J. Kjaersgaard 2011), (Tasumi 2003) and shown in Table 2. The methodology searches for image-specific pixels with the lowest and highest temperature values that match with these criteria. The output of this function is a data frame with the coordinates of the selected anchor pixels.

Table 2: Ranges of variables used for selection of anchor (cold and hot) pixels in calcAnchors(method="CITRA-MCB")

Variable	Cold pixel	Hot pixel
albedo${}^{\ast }$	0.18 - 0.25	0.13 - 0.15
NDVI${}^{\ast }$	0.76 - 0.84	0.10 - 0.28
LAI ($m^2\cdot m^{-2}$)	3 - 6	-
$Z_{om}$ ($m$)	0.03 - 0.08	$\le$ 0.005

* dimensionless

calcH()

This function applies the CIMEC self-calibration method in order to generate an iterative process for the “hot” and “cold” pixels and absorb all biases in the computation of $H$. A near surface temperature difference ($dT$) is used in place of a surface air temperature difference to drive the determination of sensible heat flux, in an iterative correction process. The convergence of the function is reached when the change in the estimated aerodynamic resistance is less than $1\mathrm{\%}$ for the cold and hot condition. There is an argument called verbose to control how much information about the iterative process is shown in the output. The output of this function is a raster layer object with the estimated soil heat flux.

Input data requirements

METRIC uses two sources to estimate the LSEB: a satellite image and a weather station with hourly data.

METRIC can be run using different satellite sensors. Currently, the coefficients needed for Landsat 7 and 5 are available in Allen, M. Tasumi, and R. Trezza (2007). Those coefficients can also be applied to Landsat 8 data. In the Appendix Section 6.1, we propose specific coefficients for the estimation of albedo using this satellite sensor. Other coefficients to run METRIC using MODIS images are available in (Allen, M. Tasumi, and R. Trezza 2007), and will be included in future versions of water package.

The weather station data should include near-surface air temperature, wind speed, relative humidity and solar radiation. The water package uses the function read.WSdata() to convert the weather station data from a comma-separate values table to a special R object. The input data should be on an hourly or shorter time basis. The expected units are ${}^{\circ }C$ for temperature; $m\cdot s^{-1}$ for wind speed; $\mathrm{\%}$ for relative humidity and $W\cdot m^{-2}$ for solar radiation. However, this function uses a parameter cf when conversion factors are needed to convert the variables from different units (e.g. if wind speed is in $km\cdot h^{-1}$ the conversion factor for this variable should be $0.278$). To estimate $E{T}_r$ from the weather station data using (Allen, I. A. Walter, E. R., T. Howell, D. Itenfisu, and M. Jensen 2005) equation, more information is needed: position in latitude and longitude, and the wind sensor height in meters. When the weather station data is being imported, a satellite metadata file can be included as a function argument. This allows interpolating the weather conditions at the exact moment of the satellite overpass.

The water package only uses one weather station to estimate $E{T}_r$ for an entire Landsat scene of $180\times 180$ $km$. The model uses $E{T}_r$ to derive the ET reference fraction ($E{T}_{r}F$) at image time (using equation (3)). This assumes that $E{T}_a$ in the entire area changes in proportion to the change in $E{T}_r$ at the weather station (Allen, M. Tasumi, and R. Trezza 2007). This means that $E{T}_r$ is only used as an index of the relative change and this is retained through the $E{T}_{r}F$. Any biases caused by variation in weather conditions should be canceled by using the same $E{T}_{r}F$ for both instantaneous and 24 h period. Nevertheless, it is recommended to use a spatial mask when the weather conditions are heterogeneous, for example in irrigated areas surrounded by deserts.

Performance and memory use

The water package uses large temporal memory in order to obtain the results and sub products. Most of the results are "RasterLayer" or "RasterStack" objects from raster package (Hijmans 2015). Processing an entire Landsat scene could need more than 2 gigabytes of memory to store the temporal data. One approach to solve this is to write products and sub-products to the disk. There is an option writeResults=TRUE in function waterOptions() to force water to store the results on the disk, instead of on temporal memory.

If water runs out of memory while processing data, it will usually stop working without a warning message. We suggest processing only a portion of a Landsat scene using an area-of-interest (aoi) polygon or storing results to disk.

Example code and datasets

Two different approaches to estimate land surface energy balance demonstrate the features and procedures in water. The first example in Section 3.2 is a simple procedure, and the second one in Section 3.3 refers to an advanced procedure. Finally, the estimation of $E{T}_a$ from the output of any of the previous procedures is demonstrated in Section 3.4.

In Section 3 functions createAoi() and read.WSdata() are summarized. Then, in Section 3.2 function METRIC.EB() is shown. In Section 3.3 the LSEB is estimated step by step. And later, in Section 3.4 daily $E{T}_r$ is estimated using functions dailyET() and ET24h().

3 Estimating ETa using METRIC model and water package

Base data preparation

Figure 1: Plot of the hourly weather station data. The conditions at the time of the sattelite overpass are marked by a gray bar. This is the default plot for waterWeatherStation objects.

To calculate METRIC Actual Evapotranspiration using water package, three sources are needed:

1. A raw Landsat 7-8 satellite image.

2. A Weather Station data (.csv file).

3. A polygon with our Area-of-interest (AOI) Spatial-Polygon object, to run the model using only a portion of the satellite scene.

First, AOI is created as a polygon using bottomright and topleft coordinates:

aoi <- createAoi(topleft = c(272955, 6085705), 
                 bottomright = c(288195, 6073195), EPSG = 32719)    

Then, the weather station data is loaded using the function read.WSdata(). This function converts the CSV file into a "waterWeatherStation" object. Then, if a Landsat metadata file (MTL file) is provided, the time-specific weather conditions at the time of satellite overpass will be calculated. This is shown on Figure 1 as a gray bar. Files apples.csv and L7.MTL.txt are included in the package as raw data. In R, system.file() is used to call this files.

csvfile <- system.file("extdata", "apples.csv", package = "water")
MTLfile <- system.file("extdata", "L7.MTL.txt", package = "water")
        
WeatherStation <- read.WSdata(WSdata=csvfile, 
                              date.format = "%d/%m/%Y",
                              lat = -35.42222, long = -71.38639, 
                              elev = 201, height = 2.2, 
                              MTL = MTLfile)

Next, the Landsat satellite image is loaded. water provides a function to load a Landsat image (loadImage()) from TIFF files. Landsat images can be downloaded directly from USGS archives in Earth Explorer (http://earthexplorer.usgs.gov/). In this article, an example dataset will be used which comes with water package as demonstration data.

image.DN <- L7_Talca

Finally, Digital Elevation Model (DEM) will be created for the area being processed. water provides two functions to do this: checkSRTMgrids() will search for the downloadable grid files in http://earthexplorer.usgs.gov/. However, this function will only print the links to the files. The downloading process has to be done manually. After this, prepareSRTMdata() can be used to mosaic and clip those files using the same extent of the image. In this article, the example data, provided with water package, will be loaded.

DEM <- DEM_Talca

Simple procedure

The simple procedure is summarized on Figure 2. The function METRIC.EB() will be used to estimate the land surface energy balance. This function has many parameters to choose from the different METRIC model equations. e.g. changes can be made in:

Figure 2: Schematic diagram with the functions, data inputs and outputs used in the simple procedure when running the METRIC model with water package. The green rounded boxes represent data, the blue boxes represent the functions and the yellow box is the final result.

• Coefficients used to estimate broadband albedo from narrowband data.

• Model to estimate Leaf Area Index (LAI) from satellite data.

• Model to estimate momentum roughness length ($Z_{om}$)

• Automatic method for the selection of anchor pixels

• Reference ET coefficient and momentum roughness length estimated for the weather station

When this function is run, the energy balance and the surface temperature ($T_s$) used are assigned to the Energy.Balance object. This function prints the position and characteristics of the anchor pixels to the console. Also, a plot with the values of the aerodynamic resistance during the iterative process is generated after every iteration. Here, the logical argument verbose controls how much information is shown in the output, and the plotting of the diagnostic graph.

Energy.Balance <- METRIC.EB(image.DN = image.DN, 
                            plain = FALSE, DEM = DEM, 
                            WeatherStation = WeatherStation, ETp.coef = 1.2,
                            MTL = MTLfile, sat = "L7", 
                            thermalband = image.DN$thermal.low)

The results of the energy balance estimated using this function are shown in Figure 3. The console output with information related to the anchor pixels goes like this:

pixel      X        Y       Ts  LAI type
1 139253 282420 -3922830 323.1587 0.13  hot
2 121566 274710 -3921780 310.0151 4.40 cold

Figure 3: Land surface energy balance ($W\cdot m^{-2}$) estimated using METRIC with water package

Advanced procedure

The advanced procedure involves running many different functions one-by-one, this is summarized in figures 4 and 5. These functions were run inside the code of METRIC.EB() in the previous example. Running water with this procedure allows to have more control in the different arguments used.

Figure 4: Schematic diagram with the functions, data inputs and outputs used in the advanced procedure when estimating Net Radiation and Soil Heat Flux running the METRIC model with water package. The green rounded boxes represent data and the blue boxes represent the functions. The functions marked with “*” have multiple methods available.

Net Radiation estimation

In order to calculate the $R_n$ for the loaded Landsat satellite data, a surface model (slope + aspect) from the DEM is calculated, then the solar angles (latitude, declination, hour angle and solar incidence angle) are calculated. Then incSWradiation() is used to calculate incoming solar radiation.

surface.model <-METRICtopo(DEM)

solar.angles.r <- solarAngles(surface.model = surface.model, 
                              WeatherStation = WeatherStation, MTL = MTLfile)

Rs.inc <- incSWradiation(surface.model = surface.model, 
                         solar.angles = solar.angles.r, 
                         WeatherStation = WeatherStation)

After this, reflectances are calculated at the top-of-atmosphere (TOA), and surface reflectance derived from the Landsat image as:

image.TOAr <- calcTOAr(image.DN = image.DN, sat = "L7", MTL = MTLfile, 
                       incidence.rel = solar.angles.r\$incidence.rel)
        
image.SR <- calcSR(image.TOAr = image.TOAr, sat = "L7", 
                   surface.model = surface.model, 
                   incidence.hor = solar.angles.r\$incidence.hor, 
                   WeatherStation = WeatherStation, ESPA = FALSE)

Following this, broadband albedo is calculated as the sum of visible to near infrared narrowband satellite bands and coefficients related to atmospheric transmittance of global solar beam radiation. In this example coeff="Tasumi" was used.

albedo <- albedo(image.SR=image.SR,  coeff="Tasumi")

Later on, Leaf Area Index (LAI) is calculated using the satellite data. In this example method=metric2010 is used:

LAI <- LAI(method="metric2010", image=image.TOAr, L=0.1)

Land surface temperature ($T_s$) is estimated using computed LAI values in order to estimate consequently the surface emissivity and brightness temperature from Landsat’s thermal band (TIR). Then this information is used to compute the incoming and outgoing long-wave radiation as:

Ts <- surfaceTemperature(LAI = LAI, sat = "L7",
                         thermalband = image.DN\$thermal.low, 
                         WeatherStation = WeatherStation)
        
Rl.out <- outLWradiation(LAI = LAI, Ts = Ts)
        
Rl.inc <- incLWradiation(WeatherStation, DEM = surface.model\$DEM, 
                         solar.angles = solar.angles.r, Ts = Ts)

Finally, Net Radiation ($R_n$) can be estimated pixel by pixel as follows:

Rn <- netRadiation(LAI, albedo, Rs.inc, Rl.inc, Rl.out)

Soil Heat Flux estimation

Soil heat flux is estimated $G$, using as input data the $R_n$, surface reflectance, $T_s$, LAI and albedo. In this example the original METRIC (2007) based-method will be used, which is:

G <- soilHeatFlux(image = image.SR, Ts = Ts, albedo = albedo, 
                  Rn = Rn, LAI = LAI)

Sensible Heat Flux estimation

To estimate the sensible heat fluxes derived from the Landsat satellite data, first, the calculation of the momentum roughness length ($Z_{om}$) is needed.

Figure 5: Schematic diagram with the functions, data inputs and outputs used in the advanced procedure when estimating the sensible heat flux running the METRIC model with water package. The green rounded boxes represent data, the blue boxes represent the functions and the yellow box is the final result. The functions marked with * have multiple methods available.

Z.om <- momentumRoughnessLength(LAI = LAI, mountainous = TRUE, 
                                method = "short.crops", 
                                surface.model = surface.model)      
    

Then, calcAnchors() is used to search for the anchor pixels within the Landsat scene. And finally, calcH() is used to run the CIMEC process and estimate the sensible heat flux:

hot.and.cold <- calcAnchors(image = image.TOAr, Ts = Ts,
                            LAI = LAI, plots = FALSE, albedo = albedo, 
                            Z.om = Z.om, n = 1, 
                            anchors.method = "CITRA-MCB", 
                            deltaTemp = 5, verbose = FALSE)
        
H <- calcH(anchors = hot.and.cold, Ts = Ts, Z.om = Z.om, 
           WeatherStation = WeatherStation, ETp.coef = 1.05, 
           Z.om.ws = 0.0018, DEM = DEM, Rn = Rn, G = G, verbose = TRUE)
    

When the function calcH() is runnig, and verbose=TRUE, the output shows the intermediate values of the CIMEC process parameters. Also, the value for the aerodynamic resistance and its change for every iteration is plotted iteratively by this function. This plot is shown on Figure 6. In this example, the change in the value of the aerodynamic resistance goes down after iteration #9, and in iteration #14 is less than $1\mathrm{\%}$.

Figure 6: Leaf area index $m^2{m}^{-2}$ and the position of the hot pixel (red cross), cold pixel (blue X) and the weather station (circle with a X) (left). Change on aerodynamic resistance convergence in the iterative process for hot (red) and cold (blue) conditions (right).

Daily crop evapotranspiration estimation

To estimate the daily actual evapotranspiration from the Landsat scene, the daily reference ET ($E{T}_r$) is needed. The daily $E{T}_r$ can be calculated with dailyET(). This function calculates the cumulative 24h standardized reference evapotranspiration for the day of the image using the ASCE standardized Penman-Monteith method (Allen, I. A. Walter, E. R., T. Howell, D. Itenfisu, and M. Jensen 2005). Finally, 24h crop ET can be estimated for every pixel of the Landsat scene using the function ET24h() (Figure 7):

ET_WS <- dailyET(WeatherStation = WeatherStation, ET = "ETr")

ET.24 <- ET24h(Rn = Rn, G = G, H = H\$H, 
               Ts = Ts, WeatherStation = WeatherStation, 
               ETr.daily = ET_WS)
    

Figure 7: Crop evapotranspiration in mm/day, estimated using METRIC and water package

4 Conclusions

The water package offers a fast and reliable platform for estimating actual evapotranspiration using the land surface energy balance. It includes different methods for many sub-models. The simple procedure showed in this article allows to estimate the $E{T}_a$ in a simple and fast way. The advanced procedure allows to have more control in the different available methods. Further versions of this package will implement other LSEB models such as SEBAL. Also, other satellite sensors such as MODIS will be included. Because water is written in R, a language very popular in the scientific community and published as free software, further developments could come from a wide community of users.

5 Acknowledgments

This work was supported by the Argentinian government through the projects INTA 1133043 and 1133042, and Universidad de Talca, Chile through the research program “Adaptation of Agriculture to Climate Change (A2C2)”. The authors wish to thank Danlu Guo and and an anonymous reviewer for their careful reading of the manuscript and their enriching comments. Finally, we thank Marcos Angelini, Rosana Vallone and Marcos Carrasco-Benavides for their valuable support and suggestions to improve this work.

6 Appendix

Included models for albedo

1. coeff="Tasumi" (Tasumi, R. G. Allen, and R. Trezza 2008)

   $$\begin{array}{rl}\text{albedo}=& {\rho }_{s,B}\cdot 0.254+{\rho }_{s,G}\cdot 0.149+{\rho }_{s,R}\cdot 0.147+{\rho }_{s,NIR}\cdot 0.311\\ & +{\rho }_{s,SWIR1}\cdot 0.103+{\rho }_{s,SWIR2}\cdot 0.036\end{array}$$

2. coeff="Liang" (Liang 2001)

   $$\begin{array}{rl}\text{albedo}=& {\rho }_{s,B}\cdot 0.356+{\rho }_{s,R}\cdot 0.130+{\rho }_{s,NIR}\cdot 0.373+{\rho }_{s,SWIR1}\cdot 0.085\\ & +{\rho }_{s,SWIR2}\cdot 0.072-0.0018\end{array}$$

3. coeff="Olmedo"

   $$\begin{array}{rl}\text{albedo}=& {\rho }_{s,B}\cdot 0.246+{\rho }_{s,G}\cdot 0.146+{\rho }_{s,R}\cdot 0.191+{\rho }_{s,NIR}\cdot 0.304\\ & +{\rho }_{s,SWIR1}\cdot 0.105+{\rho }_{s,SWIR2}\cdot 0.008\end{array}$$

   where ${\rho }_{s,b}$ is the surface reflectance for band $b$.

Included models for Leaf Area Index

1. method="metric" (Allen, M. Tasumi, and R. Trezza 2007)

   $${\text{SAVI}}_{ID}=(1+L)({\rho }_{t,NIR}-{\rho }_{t,R})/(L+{\rho }_{t,NIR}+{\rho }_{t,R})$$

   where $\rho$ is the reflectance at top-of-atmosphere, and the subindex refers to bands $R$:red or $NIR$:near infrared; and L is a soil correction factor. The default value used for L is 0.5, and in METRIC applications in western us, (Allen, M. Tasumi, and R. Trezza 2007) suggested a value of L=0.1. And Leaf Area Index is:

   $$\text{LAI}=-\frac{ln((0.69-{\text{SAVI}}_{ID})/0.59)}{0.91}$$

2. method="metric2010" (Pôças, T. a. Paço, M. Cunha, J. a. Andrade, J. Silvestre, A. Sousa, F. L. Santos, L. S. Pereira, and R. G. Allen 2014)

   $$\text{LAI}=11\cdot {\text{SAVI}}_{ID}^3$$

3. method="vineyard" (Johnson, D. E. Roczen, S. K. Youkhana, R. R. Nemani, and D. F. Bosch 2003)

   $$\text{NDVI}=({\rho }_{t,NIR}-{\rho }_{t,R})/({\rho }_{t,NIR}+{\rho }_{t,R})$$

   where $\rho$ is the reflectance at top-of-atmosphere, and the subindex refers to bands $R$:red or $NIR$:near infrared. And Leaf Area Index is:

   $$\text{LAI}=4.9\cdot \text{NDVI}-0.46$$

4. method="MCB" (Carrasco-Benavides, S. Ortega-Farı́as, L. Lagos, J. Kleissl, L. Morales-Salinas, and A. Kilic 2014)

   $$\text{LAI}=1.2-3.08\cdot \exp (-2013.35\cdot {\text{NDVI}}^6.41)$$

5. method="turner" (Turner, W. B. Cohen, R. E. Kennedy, K. S. Fassnacht, and J. M. Briggs 1999)

   $$\text{NDVI}=({\rho }_{s,NIR}-{\rho }_{s,R})/({\rho }_{s,NIR}+{\rho }_{s,R})$$

   where $\rho$ is the reflectance at surface level, and the subindex refer to bands $R$:red or $NIR$:near infrared. And Leaf Area Index is:

   $$\text{LAI}=0.5724+0.0989\cdot \text{NDVI}-0.0114\cdot {\text{NDVI}}^2+0.0004\cdot {\text{NDVI}}^3$$

Included models for Land Surface Temperature

1. method="metric" (Allen, M. Tasumi, and R. Trezza 2007) (Landsat 7)

   $${ϵ}_{NB}=0.97+0.0033LAI$$

   and ${ϵ}_{NB}=0.98$ when $LAI>3$; where ${ϵ}_0$ is the broadband surface emissivity (dimesionless); and $LAI$ is the leaf area index ($m^2\cdot m^{-2}$).

   $$T_s=\frac{K_2}{\ln [({ϵ}_{NB}{K}_1/R_c)+1]}$$

   where $T_s$ is the land surface temperature ($K$); $K_1$ and $K_2$ are specific constants for Landsat 7 ($K_1=666.1$ and $K_2=1283W\cdot m^{-2}\cdot s{r}^{-1}\cdot \mu m$); and $R_c$ is the corrected thermal radiance from the surface using the spectral radiance from band 6 of Landsat following (Wukelic, D. E. Gibbons, L. M. Martucci, and H. P. Foote 1989).

2. method="SC" (Jimenez-Munoz, J. Cristobal, J. A. Sobrino, G. S??ria, M. Ninyerola, and X. Pons 2009), (Jimenez-Munoz, J. A. Sobrino, D. Skokovic, C. Mattar, and J. Cristobal 2014) (Landsat 8)

   $$T_s=\gamma [\frac{1}{ϵ}({\upsilon }_{1}L+{\upsilon }_2)+{\upsilon }_3]+\delta$$

   where $ϵ$ is the broadband surface emissivity; and $\gamma$ and $\delta$ are two parameters given by

   $$\gamma =\frac{T^2}{b_{\gamma }\cdot L}$$

   $$\delta =T-\frac{T^2}{b_{\gamma }}$$

   where $T$ is the at-sensor brightness temperature; $b_{\gamma }$ is a coefficient equal to $1324K$ for L8 band 10, or $1199K$ for band 11. And ${\upsilon }_1$, ${\upsilon }_2$ and ${\upsilon }_3$ are the atmosferic functions, given by

   $${\upsilon }_1=\frac{1}{{\tau }_{NB}};{\upsilon }_2=-R_{sky}-\frac{R_p}{{\tau }_{NB}};{\upsilon }_3=R_{sky}$$

   where ${\tau }_{NB}$ is the narrow band transmissivity of air; $R_{sky}$ is the narrow band downward thermal radiation from a clear sky $(W{m}^{-2}s{r}^{-1}\mu m^{-1})$; and $R_p$ is path radiance in the 10.4–12.54 $\mu m$ band $(W{m}^{-2}s{r}^{-1}\mu m^{-1})$.

3. method="SW" (Jimenez-Munoz, J. A. Sobrino, D. Skokovic, C. Mattar, and J. Cristobal 2014) (Landsat 8)

   $$

   $$\begin{array}{rlr}{T}_s& =& T_i+1.378(T_{10}-T_{11})+0.183(T_{10}-T_{11}{)}^2-0.268+(54.30-2.238w)(1-ϵ)\\ & & +(-129.20+16.40w)\mathrm{\Delta }ϵ\end{array}$$

   $$

   where $T_s$ is the land surface temperature ($K$); $T_{10}$ and $T_{11}$ are the at-sensor brightness temperatures for bands 10 and 11 of Landsat 8 ($K$); $ϵ$ is the mean emissivity; $w$ is the total atmospheric water vapor content (in $g\cdot c{m}^{-2}$) and $\mathrm{\Delta }ϵ$ is the emissivity difference.

Included models for Momentum Roughness Length

1. method="short-crops" (Allen, M. Tasumi, and R. Trezza 2007)

   $$Z_{om}=0.018\ast \text{LAI}$$

2. method="custom" (Allen, M. Tasumi, and R. Trezza 2007)

   $$Z_{om}=\exp ((a\ast \text{NDVI}/\text{albedo})+b)$$

   where $a$ and $b$ are the regression coefficients derived by adjusting a lineal model between $\log Z_{om}\sim \text{NDVI}/\text{albedo}$ for points inside the Landsat scene representing specific vegetation types.

3. method="Perrier" (Santos, I. J. Lorite, R. G. Allen, and M. Tasumi 2012), (Pôças, T. a. Paço, M. Cunha, J. a. Andrade, J. Silvestre, A. Sousa, F. L. Santos, L. S. Pereira, and R. G. Allen 2014)

   $$Z_{om}=((1-\exp (-a\ast \text{LAI}/2))\ast exp(-a\ast \text{LAI}/2){)}^h$$

   where $h$ is the crop height in meters, and $a$ is:

   $$a<-(2\ast (1-fLAI){)}^{-1}$$

   when $fLAI>0.5$, or

   $$a<-2\ast fLAI$$

   when $fLAI<0.5$. And $fLAI$ is the proportion of LAI lying above $h/2$.