tseb: Computing Two Source Surface Energy Balance (TSEB) in R

In gowusu/sebkc: The Surface Energy Balance and Crop Coefficient Estimation with R

Computing Two Source Surface Energy Balance (TSEB) in R

Description

This function estimates evaporation (E), transpiration (T), and evapotranspiration(ET) of land surface. The computation is based on Norman et al (1995) and Colaizzi et al(2012) models.

Usage

tseb(albedo, Ts, LAI, DOY, xyhot = "auto", xycold = "auto", TA = "sebal",
  Tmax, Tmin, RHmax = NULL, RHmin = NULL, n = NULL, latitude = NULL,
  zom = NULL, NDVI, SAVI, hc = 2, DEM = NULL, viewangle = 0, sunelev,
  welev, u = 2, s = 2, C = 90, lapse = 0.0065, Rn24 = NULL,
  ETr = NULL, ETr24 = NULL, wmo = NULL, airport = NULL, Krs = 0.16,
  surface = "grass", time = NULL, Lz = NULL, Lm = NULL, t1 = 1, zx,
  zomw, xPT = 1.3, rc = c(70, 1000, 100), series = "xPT", clump = 1,
  fg = 1, fc = "auto", network = "parallel", clip = NULL,
  folder = NULL, iter.max = 7)

## Default S3 method:
tseb(albedo = NULL, Ts = NULL, LAI = NULL, DOY = NULL,
  xyhot = "auto", xycold = "auto", TA = "sebal", Tmax, Tmin,
  RHmax = NULL, RHmin = NULL, n = NULL, latitude = NULL, zom = NULL,
  NDVI, SAVI, hc = 2, DEM = NULL, viewangle = 0, sunelev, welev, u = 2,
  s = 2, C = 90, lapse = 0.0065, Rn24 = NULL, ETr = NULL,
  ETr24 = NULL, wmo = NULL, airport = NULL, Krs = 0.16,
  surface = "grass", time = NULL, Lz = NULL, Lm = NULL, t1 = 1,
  zx = 200, zomw = 2, xPT = c(1.26, 0, 0.1), rc = c(50, 1000, 100),
  series = "PM", clump = 1, fg = 1, fc = "auto", network = "series",
  clip = NULL, folder = NULL, iter.max = 7)

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.
LAI	Leaf Area Index, dimensionless
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
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.
TA	Air temperature[K]. It can take "sebal" or 'Tmax' or 'Tmin' or numeric value. When it is set to sebal, TA = Ts-dT. The 'Tmax' and 'Tmin' correspond to maximum and minimum temperature.
Tmax	Numeric. Maximum air Temperature in degree Celsius
Tmin	Numeric. Minimum air Temperature in degree Celsius
RHmax	Maximum Relative Humidity in percent
RHmin	Minimum Relative Humidity in percent
n	Sunshine hours
latitude	geographical coordinates in decimal degrees. It should be negative for southern hemisphere
zom	the momentum roughness length for each pixel [m]. You can also point to a 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.
hc	numeric. it is canopy height[m]
DEM	A digital elevation model[m]
viewangle	angle of view of the sensor in degrees
sunelev	Angle of Sun elevation in degrees, as found in the meta data of the satellite image
welev	Weather station elevation in meters
u	The satellites overpass time wind speed at the weather station [m/s]
s	numeric. Mean leaf size; it is four times the leaf area divided by the perimeter. The default is 2 (Grace et al., 1981).
C	coefficient. Weighting a coefficient in the equation for leaf boundary layer. A value of 90 s^0.5/m (Norman et al, 1995).
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/about/faq/international_cities.asp [Accessed on 2016-5-6] or http://www.wetterzentrale.de/klima/stnlst.html [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"
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° for Greenwich.
Lm	longitude of the measurement site [degrees west of Greenwich]. MUST BE POSITIVE
t1	numeric. The length of the calculation period in hour; 1 for hour, 0.5 for 30 minutes 0.25 for 15 minutes
zx	The height above the weather station where the wind speed is measured [m]
zomw	the roughness length for the weather station surface [m]
xPT	coefficient(s). The Priestley-Taylor equation coefficient (Norman et al, 2005). It can take a single coefficient or in the form c(initial, minimum, increment). Example is xPT=c(1.26,1,0.1). It must be noted this is a decreasing increment from the initial value .i.e 1.26. The subsequent increment is only meant for a cell where LEs is negative.
rc	coefficient. Penman-Monteith equation coefficient (Colaizzi et al,2012). It can take a single coefficient or in the form c(initial, maximum, increment). Example is rc=c(0,1000,100). Unlike xPT, this is positively incremented.
series	character.The type of series. It takes "xPT" for Priestley-Taylor equation(Norman, 1995) or "PM" for Penman-Monteith equation (Colaizzi et al, 2012)
clump	coefficient. A canopy clumping coefficient used for calculation of directional fc (fractional cover of vegetation)
fg	numeric. The fraction of LA1 that is green. If no information is available on fg, then it is assumed to be unity (1).
fc	fraction (0-1).Directional vegetation cover. The default uses Norman et al(1995) fc equation
network	character. It takes either "parallel" (Norman et al, 1995) or "series" (Colaizzi et al, 2012; Norman et al, 1995).
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
iter.max	maximum iterations of sensible heat calculation.

Details

Theory

This function separately computes LEs, LEs using Net radiation (Rn), sensible heat flux (H) and soil heat flux (G). The Rn was first computed with sebal. The soil net radiation computed as Rns=(sebal$Rn*exp(0.9*log(1-fc))) . The canopy radiation was computed as Rnc=sebal$Rn-Rns . The G was computed from sebal. The canopy sensible heat (Hc) and soil sensible heat (Hs) were computed based on corresponding canopy Temperature (Tc) and soil temperature (Tsoil). Two different approaches were employed to sensible heat fluxes. The first is Norman et al(1995) parallel network. The second is Norman (1995) secant series but uses Priestley-Taylor equation (xPT).

Default Parameter

These parameters have been set based on literature

Value

• LEc: canopy latent heat flux [W/m2]

• LEs: Soil latent heat flux[W/m2]

• T24: 24 hour Transpiration [mm/day]

• E24: 24 hour Evaporation [mm/day]

• ET24: 24 hour Evapotranspiration [mm/day]

• EF: ET24 (daily Evapotranspiration) Evaporative Fraction. EF corresponds to the type of scaling method used: see ETohr for details

• EFs: E(Evaporation) Evaporative Fraction

• EFc: T (Transpiration) Evaporative Fraction

Author(s)

George Owusu

References

• Norman JM, Kustas WP, Humes KS (1995). A two-source approach for estimating soil and vegetation energy fluxes from observations of directional radiometric surface temperature. Agric For Meteorol 1995;77:263-93.

• Colaizzi, Paul; Kustas, William P.; Anderson, Martha C.; Agam, Nurit; Tolk, Judy A.; Evett, Steven R.; Howell, Terry A.; Gowda, Prasanna H.; and O'Shaughnessy, Susan A., "Two-source energy balance model estimates of evapotranspiration using component and composite surface temperatures" (2012). Publications from USDA-ARS / UNL Faculty. Paper 1147. http://digitalcommons.unl.edu/usdaarsfacpub/1147

• Grace, J., 1981. Some effects of wind on plants. In: J. Grace, E.D. Ford and P.G. Jarvis (Editors), Plants and their Atmospheric Environment. Blackwell Scientific, London, pp. 31-56.

• Ronglin Tang, Zhao-Liang Li, Yuanyuan Jia, Chuanrong Li, Kun-Shan Chen, Xiaomin Sun & Jinyong Lou (2012): Evaluating one- and two-source energy balance models in estimating surface evapotranspiration from Landsat-derived surface temperature and field measurements, International Journal of Remote Sensing. http://dx.doi.org/10.1080/01431161.2012.716529

• KUSTAS, W. & ANDERSON, M. 2009. Advances in thermal infrared remote sensing for land surface modeling. Agricultural and Forest Meteorology, 149, 2071-2081.

Examples

## Not run: 
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"))
Parallel=tseb(Ts=Ts,LAI=LAI,DOY=37,xyhot="full",
xycold="full",albedo=albedo,Tmax=31,
Tmin=28,RHmax=NULL,RHmin=NULL,zom=NULL,NDVI=NDVI,SAVI=NULL,hc=20,
DEM=NULL,viewangle=0,sunelev=50,welev=345,u=2,
s=2,C=90,lapse=0.0065,Rn24=NULL,zx=200,zomw=2,
xPT=1.3,clump=1,fg=1,fc="auto",network="parallel",
latitude=5.6,n=6.5)

series=tseb(Ts=Ts,LAI=LAI,DOY=37,xyhot="full",
xycold="full",albedo=albedo,Tmax=31,
Tmin=28,RHmax=NULL,RHmin=NULL,zom=NULL,NDVI=NDVI,SAVI=NULL,hc=20,
DEM=NULL,viewangle=0,sunelev=50,welev=345,u=2,
s=2,C=90,lapse=0.0065,Rn24=NULL,zx=200,zomw=2,
xPT=1.3,clump=1,fg=1,fc="auto",network="series",
latitude=5.6,n=6.5)

#using landsat folder parameters
folder=system.file("extdata","stack",package="sebkc")
tsebautoseries=tseb(folder=folder,welev=317,Tmax=31,Tmin=27,
latitude=5.6,n=6.5)
tsebautoparallel=tseb(folder=folder,welev=317,Tmax=31,Tmin=27,
latitude=5.6,n=6.5,network="parallel")

#3 step simulation using landsat 7 with SLC OFF
##step one: prepare the input data 
### point to landsat 7 data with strips 

#####get the data stacked and remove strips
rawdata=sebkcstack(folder)
###### 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 
# Determine xyhot. Digitize polygon on the Ts map
modhot=hotTs(data,welev=300,extent="digitize",cluster=2)
#determine the cold. Digitize polygon on the Ts map
modcold=coldTs(data,welev=275,extent="digitize",cluster=2)
xyhot=modhot$xyhot
xycold=modcold$xycold
modTSEB2=tseb(data,Tmax=31,Tmin=28,
latitude=5.6,n=6.5)
#the same
modTSEB2=tseb(data,Tmax=33,Tmin=28,welev=345,xyhot=xyhot,xycold=xycold,
latitude=5.6,n=6.5)
#the same
modTSEB2=tseb(data,Tmax=33,Tmin=28,welev=345,xyhot=modhot,xycold=modcold,
latitude=5.6,n=6.5)

## End(Not run)

gowusu/sebkc documentation built on July 28, 2023, 11:44 a.m.