In ffuentesp7/METRICforR: Actual Evapotranspiration with Energy Balance models
This vignette presents the procedure to estimate the Land Surface Energy Balance (LSEB) using lansat imagery and The
water package.
There are two version of this vignette: Simple and advanced version.
In the first one the procedure is simpler but most of the parameters are selected automatically by the package.
In the second aproach, the procedure is longer because the user has more control on all input parameters and coefficent. Both vignettes follow the METRIC model methodology (Allen et al., 2007) in order to estimate the LSEB using landasat 7 and 8 satellite images.
Also, on the simple procedure we are going to assume flat terrain, but in the advanced procedure we are going to use a digital elevation model (DEM). However, is possible to use DEM in the simple procedure, changing the parameter flat to FALSE and providing a DEM.
Introduction
One of the most cited models to estimate land surface evapotranspiration from satellite-based energy balance is the Mapping EvapoTranspiration at high Resolution with Internalized Calibration (METRIC). This model was developed by Allen et al., (2007) based on the well-known SEBAL model (Bastiaanssen, 1998). It has been widely applied in many countries in the world to estimate crops evapotranspiration (ET) at field scales and large areas using satellite images. The model it has been apply in different vegetation and crops types such us wheat, corn, soybean and alfalfa with good results (3 - 20% of error) and also in recent years over sparse canopies (such us vineyards and olive orchards).
Thus, ET is estimated as a residual of the surface energy equation:
\begin{equation}
\label{eq:EB}
LE = R_n - G - H
\end{equation}
where $LE$ is latent energy consumed by ET ($W \cdot m^{-2}$); $Rn$ is net radiation ($W \cdot m^{-2}$); $G$ is sensible heat flux conducted into the ground ($W \cdot m^{-2}$); and $H$ = sensible heat flux convected to the air ($W \cdot m^{-2}$).
We estimate $Rn$, $G$ and $H$ for each pixel using a Landsat satellite scene and meteorological weather station data.
Then we estimate $LE$ fluxes using the previous equation, and after that, the instantaneous evapotranspiration values were calculated as follows:
\begin{equation}
ET_{inst} = 3600 \cdot \frac{LE}{\lambda \rho_w}
\end{equation}
where $ET_{inst}$ is the instantaneous ET at the time of satellite overpass ($mm \cdot h^{-1}$); 3600 is the convert factor from seconds to hours; $\rho_w$ is density of water = 1000 $kg\cdot m^{-3}$; and $\lambda$ is the water latent heat of vaporization ($J\cdot kg^{-1}$).
Finally the daily ET is computed pixel by pixel (30 x 30 m) as follows:
\begin{equation}
ET_{24} = \frac{ET_{inst}}{ET_r} ET_{r_24}
\label{eq:et24}
\end{equation}
To start the process, we have to install and load the water package:
library(water)
Base data preparation
To compute METRIC crop Evapotranspiration using water package and the simple procedure, we use the following information:
- Raw Landsat 7-8 satellite image (original .TIF data from glovis USGS).
- Weather Station data (.CSV file).
- Area-of-interest (AOI) (if we don´t define an AOI, we will estimate crop evapotranspiration for the entire landsat scene).
First, we create the AOI as a polygon using bottomright and topleft points:
aoi <- createAoi(topleft = c(273110, -3914450),
bottomright = c( 288050, -3926650), EPSG = 32619)
Then, we load weather station data. For that, we use read.WSdata function.
This command transforms our .csv file into a waterWeatherStation object.
Then, we provide a Landsat metadata file (.MTL file).
Doing this, we will be able to estimate the time-specific weather conditions at the time of satellite overpass.
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)
We can visualize our meteorological information using the following command:
print(WeatherStation)
plot(WeatherStation, alldata=FALSE)
After that, we load the Landast satellite image.
Usually, using water, we can use the function called loadImage to load a Landsat image from .TIF files downloaded directly
from Earth Explorer or GLOVIS.
In this vignette we are going to use some Landsat 7 as example data wich, by the way comes with water package as a demo.
image.DN <- L7_Talca[[c(1:5,7)]]
B6 <- L7_Talca[[6]]
Finally we are going to create The Digital Elevation Model (DEM) for the specific satellite image we are processing.
Thus, we needed a image-specifc grid files downloadable from Earth Explorer.
Therefore, check wich grid files we are going to need considering the satellite scene location and the AOI polygon, using:
checkSRTMgrids(image.DN)
You should download all needed grid files, and then you can use the function prepareSRTMdata(extent = image.DN) to create the DEM. In this vignette we are going to load the example data provided with water package.
ARREGLAR EXAMPLE DATASED DEM_Talca
DEM <- DEM_Talca
Net Radiation (R_n) estimation
In order to calculate the Net Radiation from loaded landsat satellite data, first we calculate a surface model (slope + aspect) from the DEM, then we calulate the solar angles (latitude, declination, hour angle and solar incidence angle), and finally we plot the last one. Then we use this function to calculate incoming solar radiation.
surface.model <-METRICtopo(DEM)
solar.angles.r <- solarAngles(surface.model = surface.model,
WeatherStation = WeatherStation, MTL = MTLfile)
plot(solar.angles.r)
Rs.inc <- incSWradiation(surface.model = surface.model,
solar.angles = solar.angles.r,
WeatherStation = WeatherStation)
Then we calculate reflectances at the top-of-atmosphere (TOA), and surface reflectance derived from Landsat image, and use the last one to calculate the broadband albedo 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 = F)
albedo <- albedo(image.SR = image.SR, coeff="Tasumi", sat = "L7")
Later on, we calculate the Leaf Area Index (LAI) using the satellite data, and we plot it. In this step, we can choose diferents methods able in literature to estimate LAI from Landsat data (see package Help for more info.). In this vignette we are going to use the METRIC 2010 method:
LAI <- LAI(method = "metric2010", image = image.TOAr, L=0.1, sat = "L7")
plot(LAI)
Then we estimate land surface temperature (Ts), using computed LAI in order to estimate consequently the surface emissivity, and brightness temperature from landsat`s thermal band (TIR). Then we use this information to compute the incoming and outgoing long-wave radiation as:
Ts <- surfaceTemperature(LAI=LAI, sat = "L7", thermalband = B6,
WeatherStation = WeatherStation)
Rl.out <- outLWradiation(LAI = LAI, Ts=Ts)
Rl.inc <- incLWradiation(WeatherStation,DEM = surface.model$DEM,
solar.angles = solar.angles.r, Ts= Ts)
And finally, we can estimate Net Radiation (R_n) pixel by pixel as follows:
Rn <- netRadiation(LAI, albedo, Rs.inc, Rl.inc, Rl.out)
plot(Rn)
Soil Heat Flux (G) estimation
We estimate the soil heat flux, using as input data the Net radiation, surface reflectance, surface temperature, leaf area index and albedo. in this vignette we are going to use the original METRIC (2007) baded-method, which is:
G <- soilHeatFlux(image = image.SR, Ts=Ts,albedo=albedo,
Rn=Rn, LAI=LAI, sat = "L7")
plot(G)
Sensible Heat Flux (H) estimation
To estimate the sensible heat fluxes derived from the landsat satellite data, first we need to calculate the Momentum momentum roughness length (Zom). In this step, we are able to choose diferents methods able in literature to estimate Zom from Landsat data (see package Help for more information). In this vignette we are going to use the short.crops method.
Then, we are going to use calAnchors function for automatically search end members within the satellite scene (extreme wet and dry conditions). In this vignette we are going to use the CITRA-MCB method.
Finally we estimate the H values at pixel scales using calcH function. This function aplies the CIMEC self-calibration method in order generete an iteration process for the "hot" and "cald" pixels and abdorb all baises in the computation of H. This iteration procces it will be able to see as a dinamic-plot when the function calcH is running. There is a parameter called verbose to control how much information we want to see in the output about the "anchor pixels" and iteration process.
Z.om <- momentumRoughnessLength(LAI=LAI, mountainous = TRUE,
method = "short.crops",
surface.model = surface.model)
hot.and.cold <- calcAnchors(image = image.TOAr, Ts = Ts, LAI = LAI, plots = F,
albedo = albedo, Z.om = Z.om, n = 1,
anchors.method = "CITRA-MCB", sat = "L7",
deltaTemp = 5, verbose = FALSE)
H <- calcH(anchors = hot.and.cold, Ts = Ts, Z.om = Z.om,
WeatherStation = WeatherStation, ETp.coef = 1.05, sat = "L7",
Z.om.ws = 0.0018, DEM = DEM, Rn = Rn, G = G, verbose = FALSE)
Daily Crop Evapotranspiration (ET_24) estimation
To estimate the daily crop evapotranspiration from the Landsat scene we need the daily reference ET (ETr) for our weather station, so we can calculate the daily ETr with:
ET_WS <- dailyET(WeatherStation = WeatherStation)
And finally, we can estimate daily crop ET pixel by pixel:
ET.24 <- ET24h(Rn, G, H$H, Ts, WeatherStation = WeatherStation, ETr.daily=ET_WS)
References
Allen, R. G., Tasumi, M., & Trezza, R. (2007). Satellite-based energy balance for mapping evapotranspiration with internalized calibration (METRIC)—Model. Journal of irrigation and drainage engineering, 133(4), 380-394.
Bastiaanssen, W. G. M., Menenti, M., Feddes, R. A., & Holtslag, A. A. M. (1998). A remote sensing surface energy balance algorithm for land (SEBAL). 1. Formulation. Journal of hydrology, 212, 198-212.
Bastiaanssen, W. G. M., Noordman, E. J. M., Pelgrum, H., Davids, G., Thoreson, B. P., & Allen, R. G. (2005). SEBAL model with remotely sensed data to improve water-resources management under actual field conditions. Journal of irrigation and drainage engineering, 131(1), 85-93.
ffuentesp7/METRICforR documentation built on May 15, 2019, 11:57 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
This vignette presents the procedure to estimate the Land Surface Energy Balance (LSEB) using lansat imagery and The water package.
There are two version of this vignette: Simple and advanced version. In the first one the procedure is simpler but most of the parameters are selected automatically by the package.
In the second aproach, the procedure is longer because the user has more control on all input parameters and coefficent. Both vignettes follow the METRIC model methodology (Allen et al., 2007) in order to estimate the LSEB using landasat 7 and 8 satellite images.
Also, on the simple procedure we are going to assume flat terrain, but in the advanced procedure we are going to use a digital elevation model (DEM). However, is possible to use DEM in the simple procedure, changing the parameter flat to FALSE and providing a DEM.
Introduction
One of the most cited models to estimate land surface evapotranspiration from satellite-based energy balance is the Mapping EvapoTranspiration at high Resolution with Internalized Calibration (METRIC). This model was developed by Allen et al., (2007) based on the well-known SEBAL model (Bastiaanssen, 1998). It has been widely applied in many countries in the world to estimate crops evapotranspiration (ET) at field scales and large areas using satellite images. The model it has been apply in different vegetation and crops types such us wheat, corn, soybean and alfalfa with good results (3 - 20% of error) and also in recent years over sparse canopies (such us vineyards and olive orchards).
Thus, ET is estimated as a residual of the surface energy equation:
\begin{equation} \label{eq:EB} LE = R_n - G - H \end{equation}
where $LE$ is latent energy consumed by ET ($W \cdot m^{-2}$); $Rn$ is net radiation ($W \cdot m^{-2}$); $G$ is sensible heat flux conducted into the ground ($W \cdot m^{-2}$); and $H$ = sensible heat flux convected to the air ($W \cdot m^{-2}$).
We estimate $Rn$, $G$ and $H$ for each pixel using a Landsat satellite scene and meteorological weather station data. Then we estimate $LE$ fluxes using the previous equation, and after that, the instantaneous evapotranspiration values were calculated as follows:
\begin{equation} ET_{inst} = 3600 \cdot \frac{LE}{\lambda \rho_w} \end{equation}
where $ET_{inst}$ is the instantaneous ET at the time of satellite overpass ($mm \cdot h^{-1}$); 3600 is the convert factor from seconds to hours; $\rho_w$ is density of water = 1000 $kg\cdot m^{-3}$; and $\lambda$ is the water latent heat of vaporization ($J\cdot kg^{-1}$).
Finally the daily ET is computed pixel by pixel (30 x 30 m) as follows:
\begin{equation} ET_{24} = \frac{ET_{inst}}{ET_r} ET_{r_24} \label{eq:et24} \end{equation}
To start the process, we have to install and load the water package:
library(water)
Base data preparation
To compute METRIC crop Evapotranspiration using water package and the simple procedure, we use the following information:
- Raw Landsat 7-8 satellite image (original .TIF data from glovis USGS).
- Weather Station data (.CSV file).
- Area-of-interest (AOI) (if we don´t define an AOI, we will estimate crop evapotranspiration for the entire landsat scene).
First, we create the AOI as a polygon using bottomright and topleft points:
aoi <- createAoi(topleft = c(273110, -3914450), bottomright = c( 288050, -3926650), EPSG = 32619)
Then, we load weather station data. For that, we use read.WSdata function.
This command transforms our .csv file into a waterWeatherStation object.
Then, we provide a Landsat metadata file (.MTL file).
Doing this, we will be able to estimate the time-specific weather conditions at the time of satellite overpass.
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)
We can visualize our meteorological information using the following command:
print(WeatherStation) plot(WeatherStation, alldata=FALSE)
After that, we load the Landast satellite image.
Usually, using water, we can use the function called loadImage to load a Landsat image from .TIF files downloaded directly
from Earth Explorer or GLOVIS.
In this vignette we are going to use some Landsat 7 as example data wich, by the way comes with water package as a demo.
image.DN <- L7_Talca[[c(1:5,7)]] B6 <- L7_Talca[[6]]
Finally we are going to create The Digital Elevation Model (DEM) for the specific satellite image we are processing. Thus, we needed a image-specifc grid files downloadable from Earth Explorer. Therefore, check wich grid files we are going to need considering the satellite scene location and the AOI polygon, using:
checkSRTMgrids(image.DN)
You should download all needed grid files, and then you can use the function prepareSRTMdata(extent = image.DN) to create the DEM. In this vignette we are going to load the example data provided with water package.
ARREGLAR EXAMPLE DATASED DEM_Talca
DEM <- DEM_Talca
Net Radiation (R_n) estimation
In order to calculate the Net Radiation from loaded landsat satellite data, first we calculate a surface model (slope + aspect) from the DEM, then we calulate the solar angles (latitude, declination, hour angle and solar incidence angle), and finally we plot the last one. Then we use this function to calculate incoming solar radiation.
surface.model <-METRICtopo(DEM) solar.angles.r <- solarAngles(surface.model = surface.model, WeatherStation = WeatherStation, MTL = MTLfile) plot(solar.angles.r) Rs.inc <- incSWradiation(surface.model = surface.model, solar.angles = solar.angles.r, WeatherStation = WeatherStation)
Then we calculate reflectances at the top-of-atmosphere (TOA), and surface reflectance derived from Landsat image, and use the last one to calculate the broadband albedo 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 = F) albedo <- albedo(image.SR = image.SR, coeff="Tasumi", sat = "L7")
Later on, we calculate the Leaf Area Index (LAI) using the satellite data, and we plot it. In this step, we can choose diferents methods able in literature to estimate LAI from Landsat data (see package Help for more info.). In this vignette we are going to use the METRIC 2010 method:
LAI <- LAI(method = "metric2010", image = image.TOAr, L=0.1, sat = "L7") plot(LAI)
Then we estimate land surface temperature (Ts), using computed LAI in order to estimate consequently the surface emissivity, and brightness temperature from landsat`s thermal band (TIR). Then we use this information to compute the incoming and outgoing long-wave radiation as:
Ts <- surfaceTemperature(LAI=LAI, sat = "L7", thermalband = B6, WeatherStation = WeatherStation) Rl.out <- outLWradiation(LAI = LAI, Ts=Ts) Rl.inc <- incLWradiation(WeatherStation,DEM = surface.model$DEM, solar.angles = solar.angles.r, Ts= Ts)
And finally, we can estimate Net Radiation (R_n) pixel by pixel as follows:
Rn <- netRadiation(LAI, albedo, Rs.inc, Rl.inc, Rl.out) plot(Rn)
Soil Heat Flux (G) estimation
We estimate the soil heat flux, using as input data the Net radiation, surface reflectance, surface temperature, leaf area index and albedo. in this vignette we are going to use the original METRIC (2007) baded-method, which is:
G <- soilHeatFlux(image = image.SR, Ts=Ts,albedo=albedo, Rn=Rn, LAI=LAI, sat = "L7") plot(G)
Sensible Heat Flux (H) estimation
To estimate the sensible heat fluxes derived from the landsat satellite data, first we need to calculate the Momentum momentum roughness length (Zom). In this step, we are able to choose diferents methods able in literature to estimate Zom from Landsat data (see package Help for more information). In this vignette we are going to use the short.crops method.
Then, we are going to use calAnchors function for automatically search end members within the satellite scene (extreme wet and dry conditions). In this vignette we are going to use the CITRA-MCB method.
Finally we estimate the H values at pixel scales using calcH function. This function aplies the CIMEC self-calibration method in order generete an iteration process for the "hot" and "cald" pixels and abdorb all baises in the computation of H. This iteration procces it will be able to see as a dinamic-plot when the function calcH is running. There is a parameter called verbose to control how much information we want to see in the output about the "anchor pixels" and iteration process.
Z.om <- momentumRoughnessLength(LAI=LAI, mountainous = TRUE, method = "short.crops", surface.model = surface.model) hot.and.cold <- calcAnchors(image = image.TOAr, Ts = Ts, LAI = LAI, plots = F, albedo = albedo, Z.om = Z.om, n = 1, anchors.method = "CITRA-MCB", sat = "L7", deltaTemp = 5, verbose = FALSE) H <- calcH(anchors = hot.and.cold, Ts = Ts, Z.om = Z.om, WeatherStation = WeatherStation, ETp.coef = 1.05, sat = "L7", Z.om.ws = 0.0018, DEM = DEM, Rn = Rn, G = G, verbose = FALSE)
Daily Crop Evapotranspiration (ET_24) estimation
To estimate the daily crop evapotranspiration from the Landsat scene we need the daily reference ET (ETr) for our weather station, so we can calculate the daily ETr with:
ET_WS <- dailyET(WeatherStation = WeatherStation)
And finally, we can estimate daily crop ET pixel by pixel:
ET.24 <- ET24h(Rn, G, H$H, Ts, WeatherStation = WeatherStation, ETr.daily=ET_WS)
References
Allen, R. G., Tasumi, M., & Trezza, R. (2007). Satellite-based energy balance for mapping evapotranspiration with internalized calibration (METRIC)—Model. Journal of irrigation and drainage engineering, 133(4), 380-394.
Bastiaanssen, W. G. M., Menenti, M., Feddes, R. A., & Holtslag, A. A. M. (1998). A remote sensing surface energy balance algorithm for land (SEBAL). 1. Formulation. Journal of hydrology, 212, 198-212.
Bastiaanssen, W. G. M., Noordman, E. J. M., Pelgrum, H., Davids, G., Thoreson, B. P., & Allen, R. G. (2005). SEBAL model with remotely sensed data to improve water-resources management under actual field conditions. Journal of irrigation and drainage engineering, 131(1), 85-93.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.