Skip to contents

R computation of Surface Energy Balance Algorithm for Land and METRIC

Source: R/sebal.R
sebal.Rd

This R function computes Bastiaanssen (1998) surface energy balance components of evapotranspiration (ET), Sensible heat(H), soil heat flux(G) and Net Radiation(Rn).

Usage

sebal(
  albedo,
  Ts,
  NDVI,
  SAVI,
  welev,
  xyhot = "auto",
  xycold = "auto",
  DOY = NULL,
  sunelev,
  zx = 10,
  u = 2,
  zomw = 2,
  zom = NULL,
  LAI = NULL,
  DEM = NULL,
  lapse = 0.0065,
  Rn24 = NULL,
  ETr = NULL,
  ETr24 = NULL,
  wmo = NULL,
  airport = NULL,
  Krs = 0.16,
  surface = "grass",
  latitude = NULL,
  t1 = 1,
  time = NULL,
  Lz = NULL,
  Lm = NULL,
  model = "SEBAL",
  iter.max = 7,
  clip = NULL,
  folder = NULL
)

# Default S3 method
sebal(
  albedo = NULL,
  Ts = NULL,
  NDVI = NULL,
  SAVI = NULL,
  welev,
  xyhot = "auto",
  xycold = "auto",
  DOY = NULL,
  sunelev = NULL,
  zx = 10,
  u = 2,
  zomw = 2,
  zom = NULL,
  LAI = NULL,
  DEM = NULL,
  lapse = 0.0065,
  Rn24 = NULL,
  ETr = NULL,
  ETr24 = NULL,
  wmo = NULL,
  airport = NULL,
  Krs = 0.16,
  surface = "grass",
  latitude = NULL,
  t1 = 1,
  time = NULL,
  Lz = NULL,
  Lm = NULL,
  model = "SEBAL",
  iter.max = 7,
  clip = NULL,
  folder = NULL
)

Arguments

albedo

A RasterLayer data that has albedo values. You can also provide file path or location on your computer

Ts

A RasterLayer data that indicates radiometric surface temperature values from remote sensing image, preferably in Kelvin (K). You can also point to raster file on your computer.

NDVI

A RasterLayer data that indicates NDVI (Normalized Difference Vegetation Index) values. You can also point to a file on your computer.

SAVI

A RasterLayer data that indicates SAVI(Soil-adjusted Vegetation Index) values. You can also point to file on your computer.

welev

Weather station elevation in meters

xyhot

numeric or "auto". A list of x and y coordinates of a hot pixel in the form c(x,y). If it is set to "auto", hotTs will be used to compute it.

xycold

numeric or "auto". A list of x and y coordinates of a cold pixel in the form of c(x,y). If it set to "auto", coldTs will be used to compute it.

DOY

Numeric or Date [YYYY-mm-dd]. Day of the Year. If you give data in the form of date [YYYY-mm-dd], it will be converted to DOY

sunelev

Angle of Sun elevation in degrees, as found in the meta data of the satellite image

zx

The height above the weather station where the wind speed is measured [m]

u

The satellites overpass time wind speed at the weather station [m/s]

zomw

the roughness length for the weather station surface [m]

zom

the momentum roughness length for each pixel [m]. You can also point to a raster file on your computer

LAI

Leaf Area Index, dimensionless

DEM

A digital elevation model[m]

lapse

A local lapse rate that is applied to correct DEM [K/m]. Default value is 0.0065

Rn24

A 24 hour net radiation [W/m2]. It is needed to estimate daily ET. It can be computed with ETo with even a minimal data of Tmax and Tmin. You can also point to raster file on your computer

ETr

an alfalfa or grass reference ET at the time of satellites overpass. It only needed if you set model to "metric" [mm/hr]

ETr24

daily grass or alfalfa reference Evapotranspiration (mm/day)

wmo

numeric. World Meteorological Organization weather station code. You can use the Worldwide Station List at https://www.wunderground.com/ [Accessed on 2016-5-6] or https://www.wetterzentrale.de/ [Accessed on 2016-5-6]

airport

numeric. IATA or ICAO code. They are alphabetically listed at https://en.wikipedia.org/wiki/List_of_airports_by_IATA_code:_A [Accessed on 2016-5-6]

Krs

Relationship between the fraction of extraterrestrial radiation that reaches the earth's surface, Rs/Ra, and the air temperature difference Tmax - Tmin for interior (Krs = 0.16) and coastal (Krs = 0.19) regions

surface

character. It is either "grass" or "alfalfa"

latitude

geographical coordinates in decimal degrees. It should be negative for southern hemisphere

t1

numeric. The length of the calculation period in hour; 1 for hour, 0.5 for 30 minutes 0.25 for 15 minutes

time

The midpoint of the time of measurement[hour]; for example time is 12.5 for a period between 12:00 and 13:00

Lz

longitude of the centre of the local time zone [degrees west of Greenwich]. MUST BE POSITIVE Lz = 0 degrees for Greenwich.

Lm

longitude of the measurement site [degrees west of Greenwich]. MUST BE POSITIVE

model

character. The type of model. It takes either "SEBAL" or "METRIC"

iter.max

maximum iterations of sensible heat calculation.

clip

extent object or raster object or polygon from which an Extent object can be extracted. A polygon or raster will be reprojected to conform the data to be cropped. for example clip can take the form of c(xmin, xmax, ymin, ymax)

folder

An original directory of the images that contains all landsat 5,7 or 8 bands and metadata. At the moment only Landsat folder is supported

Value

LE:

latent heat flux [W/m2]

Rn:

Net Radiation [W/m2]

H:

Sensible Heat Flux [W/m2]

G:

Soil heat flux [W/m2]

ETins:

instantaneous ET if model = "METRIC" [mm/hr]

ET24:

24 hour Evapotranspiration only if Rn24 is not NULL [mm/day]

ETrF:

METRIC ET fraction [dimensionless]

EF:

ET fraction [dimensionless]

Ta:

Air temperature = Ts-dT [K]

iter.coef:

Parameters from sensible heat iteration

Details

Computing SEBAL or METRIC

sebal can compute either SEBAL(Bastiaanssen et al,1998) or METRIC (Allen et al 2007). In the case of METRIC, input data such as DEM, local lapse rate, and ETr must be provided. See Allen et al (2007) for more details.

Using Landsat Folder as input data

Though this sebkc function is not sensor dependent, one can point the folder parameter to a Landsat TM(5), ETM(7) or 8 folder that contains landsat bands and metadata. One can also point the folder parameter to object of landsat578

Automatic computation

When xycold and xyhot are set to "auto", about 80 per cent of the input image pixels are used to compute cold and hot pixels. This avoids using NAs in classification of pixels.

Taking long time to Simulate?

If you point the folder to landsat files and set xycold and xyhot to "auto", it will take about 2 hours for the function to run. Without these "autos", it will take about 45 minutes. It is advisable to perform 3 step analyses of using three functions as follows: sebkcstack, landsat578, hotTs, coldTs. You need not to go through this if you have your albedo, Ts, NDVI, SAVI from MODIS, ASTER etc

References

  • Bastiaanssen, W. G. M., Menenti, M., Feddes, R. A., & Holtslag, A. A. M. 1998. A remote sensing surface energy balance algorithm for land (SEBAL), part 1: formulation. Journal of Hydrology: 212-213: 198-212.

  • Bastiaanssen, W. G. M., Pelgrum, H., Wang, J., Ma, Y., Moreno, J. F., Roerink, G. J., & van der Wal, T. 1998. A remote sensing surface energy balance algorithm for land (SEBAL).: Part 2: Validation. Journal of Hydrology, 212-213(0): 213-229.

  • Allen, R., 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.

  • Allen, R., Tasumi, M., Morse, A., Trezza, R., Wright, J., Bastiaanssen, W., Kramber, W., Lorite, I., & Robison, C. 2007. Satellite-Based Energy Balance for Mapping Evapotranspiration with Internalized Calibration (METRIC)-Applications. Journal of Irrigation and Drainage Engineering, 133(4): 395-406.

Author

George Owusu

Examples

if (FALSE) { # \dontrun{
#use original landsat 8 data by specifying folder path
folder=system.file("extdata","stack",package="sebkc")
modauto=sebal(folder = folder,welev = 380,xycold="full",xyhot="full")
#plot ET fraction
plot(modauto$EF)

#Semi-auto
#load individual input files, for example from MODIS, ASTER etc
albedo=raster(system.file("extdata","albedo.grd",package="sebkc"))
Ts=raster(system.file("extdata","Ts.grd",package="sebkc"))
NDVI=raster(system.file("extdata","NDVI.grd",package="sebkc"))
LAI=raster(system.file("extdata","LAI.grd",package="sebkc"))
mod=sebal(albedo=albedo,Ts=Ts,NDVI=NDVI,SAVI=NULL,
iter.max=7,xyhot="full",xycold="full",
DOY=37,sunelev=50.71154048,welev=317.1,zx=10,
u=2,zomw=3,zom=NULL,LAI=LAI,DEM=NULL,
lapse=0.0065,Rn24=NULL,ETr=NULL,model="SEBAL")
#get latent heat flux
mod$LE
#evaporation fraction
mod$EF

#3 step simulation using landsat 7 with SLC OFF
##step one: prepare the input data 
### point to landsat 7 data with strips 
folder=system.file("extdata","slc",package="sebkc")
#####get the data stacked and remove strips
rawdata=sebkcstack(folder,remove_cloud = 2,gap_fill = TRUE)
###### check if strips are removed
plot(rawdata$data)
#get Ts, Albedo, NDVI, SAVI, sunelev, DOY from landsat data 
welev=278
 data=landsat578(rawdata,welev=welev)
 #perform semi-auto simulation mod=sebal(data,welev=welev)
# Determine xyhot. Digitize polygon on the Ts map
modhot=hotTs(data,welev=welev,extent="auto",cluster=2)
#determine the cold. Digitize polygon on the Ts map
modcold=coldTs(data,welev=welev,extent="auto",cluster=2)
xyhot=modhot$xyhot
xycold=modcold$xycold
modsebal=sebal(folder=data,xyhot=xyhot,xycold=xycold,welev=welev)
} # }