Package {sebkc}

sebkc: Surface Energy Balance and Crop Coefficient Evapotranspiration Estimation

Description

Computes and integrates surface energy balance components of evapotranspiration into the FAO56 water balance model, providing the SEBAL, METRIC, SEBI, SSEB, SEBS and TSEB models.

Details

Generics that the 'raster' package also provides as methods (such as plot, hist, quantile, text, lines and points) are intentionally not imported here, because they are made available through import(raster); importing them again would trigger "replacing previous import" warnings.

Author(s)

Maintainer: George Owusu owusugeorge@ug.edu.gh

Authors:

• George Owusu owusugeorge@ug.edu.gh

See Also

Useful links:

• https://github.com/gowusu/sebkc

• Report bugs at https://github.com/gowusu/sebkc/issues

24 Hour FAO 56 and ASCE-EWRI Reference Evapotranspiration ETo

Description

24 Hour FAO 56 and ASCE-EWRI Reference Evapotranspiration ETo

Usage

ETo(
  Tmax,
  Tmin,
  DOY,
  latitude,
  longitude = NULL,
  z = NULL,
  uz = 2,
  altitude = NULL,
  RHmax = NULL,
  RHmin = NULL,
  n = NULL,
  Krs = 0.16,
  albedo = 0.23,
  as = 0.25,
  bs = 0.5,
  Rs = NULL,
  Rn = NULL,
  G = 0,
  EF = NULL,
  wmo = NULL,
  airport = NULL,
  scale = "sine function",
  Kc = NULL,
  surface = "grass"
)

## Default S3 method:
ETo(
  Tmax,
  Tmin,
  DOY,
  latitude,
  longitude = NULL,
  z = NULL,
  uz = 2,
  altitude = NULL,
  RHmax = NULL,
  RHmin = NULL,
  n = NULL,
  Krs = 0.16,
  albedo = 0.23,
  as = 0.25,
  bs = 0.5,
  Rs = NULL,
  Rn = NULL,
  G = 0,
  EF = NULL,
  wmo = NULL,
  airport = NULL,
  scale = "evaporative fraction",
  Kc = NULL,
  surface = "grass"
)

Arguments

Value

ETo:

24 hour reference ETo [mm/day]

ETa:

24 hour actual ET. It is not NULL if EF has a numeric value

ETc:

24 hour crop ET. It is not NULL if Kc has a numeric value

Rn:

Net daily Radiation [W/m2]

Rn_mj_m2_day:

Net daily Radiation [MJ/m2/day]

y:

psychrometric constant [kPa/degree Celsius]

slope:

slope of saturation vapour pressure curve [kPa/degree Celsius]

vpd:

Vapour pressure deficit[kPa]

Author(s)

George Owusu

References

ALLEN, R. G., PEREIRA, L. S., RAES, D., & SMITH, M. 1996. Crop Evapotranspiration (guidelines for computing crop water requirements) FAO Irrigation and Drainage Paper No. 56: FAO.

See Also

ETohr

Examples

#with minimal data
PET=ETo(Tmax=31,Tmin=28,latitude=5.6,DOY="6-2-2001")
PET$ETo
#with enough data as in FAO56 EXAMPLE 18 on 6 July in Uccle (Brussels, Belgium) 
#located at 50 degrees 48 minutes N and at 100 m above sea level
PET24=ETo(Tmax=21.5,Tmin=12.3,z=10,uz=2.78,altitude=100,
RHmax=84,RHmin=63,n=9.25,DOY=187,latitude=50.8)
PET24$ETo

Hourly FAO-56 and ASCE-EWRI Reference Evapotranspiration

Description

The function calculates hourly or less than hour FAO-56 (grass) and ASCE-EWRI Reference (alfalfa) Evapotranspiration (ETo ETr). The "surface" parameter determines #' the type of ET. The function can also be used to scale hourly or less than hour actual evapotranspiration. The daily ET is implemented in ETo.

Usage

ETohr(
  Tmean,
  RH,
  DOY,
  Lz,
  t1,
  time,
  Lm,
  latitude,
  longitude = NULL,
  n = NULL,
  uz = 2,
  Rs = NULL,
  Rn = NULL,
  G = NULL,
  altitude = NULL,
  z = NULL,
  Krs = 0.16,
  Tmax = NULL,
  Tmin = NULL,
  albedo = 0.23,
  as = 0.25,
  bs = 0.5,
  ETins = NULL,
  ETr24 = NULL,
  wmo = NULL,
  airport = NULL,
  map = NULL,
  surface = "grass",
  period = "daytime"
)

## Default S3 method:
ETohr(
  Tmean,
  RH,
  DOY,
  Lz,
  t1 = 1,
  time,
  Lm,
  latitude,
  longitude = NULL,
  n = NULL,
  uz = 2,
  Rs = NULL,
  Rn = NULL,
  G = NULL,
  altitude = NULL,
  z = NULL,
  Krs = 0.16,
  Tmax = NULL,
  Tmin = NULL,
  albedo = 0.23,
  as = 0.25,
  bs = 0.5,
  ETins = NULL,
  ETr24 = NULL,
  wmo = NULL,
  airport = NULL,
  map = NULL,
  surface = "grass",
  period = "daytime"
)

Arguments

Details

sebal, sebi,sebs,sseb, tseb

Value

ETo:

hour reference ETo [mm/hr]

ETa:

hour actual ET. It is not NULL if EF has a numeric value

ETc:

hour crop ET. It is not NULL if Kc has a numeric value

Rn:

Net Radiation [W/m2]

Rnwm_2:

Net Radiation [MJ/m2/hr]

y:

psychrometric constant [kPa/degree Celsius]

slope:

slope of saturation vapour pressure curve [kPa/degree Celsius]

vpd:

Vapour pressure deficit[kPa]

Author(s)

George Owusu

References

• ALLEN, R. G., PEREIRA, L. S., RAES, D., & SMITH, M. 1998. Crop Evapotranspiration (guidelines for computing crop water requirements) FAO Irrigation and Drainage Paper No. 56: FAO.

• Jackson, R. D., Hatfield, J. L., Reginato, R. J., Idso, S. B., & Pinter Jr, P. J. (1983). Estimation of daily evapotranspiration from one time-of-day measurements. Agricultural Water Management, 7(1-3), 351-362. doi: http://dx.doi.org/10.1016/0378-3774(83)90095-1

• French, A. N., Fitzgerald, G., Hunsaker, D., Barnes, E., Clarke, T., Lesch, S., . . . Pinter, P. (2005). Estimating spatially distributed cotton water use from thermal infrared aerial imagery. Paper presented at the World Water Congress 2005: Impacts of Global Climate Change - Proceedings of the 2005 World Water and Environmental Resources Congress, Reston, Va.

• Colaizzi, P. D., Kustas, W. P., Anderson, M. C., Agam, N., Tolk, J. A., Evett, S. R., . . . O'Shaughnessy, S. A. (2012). Two-source energy balance model estimates of evapotranspiration using component and composite surface temperatures. Advances in Water Resources, 50, 134-151. doi: http://dx.doi.org/10.1016/j.advwatres.2012.06.004

See Also

ETo

Examples

#FAO56 Example 19 14:00-15:00 h using Cd = 0.24 instead 0.34 in FAO56
ET1hr=ETohr(Tmean=38,RH=52,DOY=274,Lz=15,t1=1,time=14.5,Lm=16.25,latitude=16.22,uz=3.3,Rs=2.450)
ET1hr$ETo

Biomass and yield estimation of crops and plants

Description

This function uses Wageningen procedure (De Wit, 1968; Slabbers et al., 1979) to compute the point or spatial possible maximum yield of a standard crop based on a Rs data. This function can also perform internal interpolation of data if longitude and latitude are provided. See datails and examples below.

Usage

biomass(
  Tday = NULL,
  T24,
  Rs,
  latitude,
  longitude = NULL,
  HI = NULL,
  plantdate,
  harvestdate,
  c30 = 0.0108,
  adaptability = 3,
  LAI = 5,
  crop = NULL
)

Arguments

Details

The function can be used in several ways:

Point Estimation

biomass can compute biomass and yield for one location. Each input data must be given only one value

Point Estimation for several locations

The function can compute biomass and yields for more than one location. Input data must be a list or array.

Spatial Estimation

Input data such as Rs, Tday, T24, HI or LAI can be spatial.

Spatial Estimation with internal interpolation

It takes Tday, for example from satellite image, and do interpolation for the rest of the data. Both longitude, latitude, and the variables must be of the same length. If longitude is not provided no interpolation will be done.

@references

• De Wit, C. T. (1968). Plant production Miscellaneous Papers: Landbouw Hogeschool, Wageningen.

• Slabbers, P. J., Herrendorf, V. S., & Stapper, M. (1979). Evaluation of simplified water-crop yield models. Agricultural Water Management, 2(2), 95-129. doi: http://dx.doi.org/10.1016/0378-3774(79)90026-X

• FAO (1981). Report on the Agro-ecological Zones Project: Food and Agriculture Organization of the United Nations.

Value

yield and biomass in t/ha

Author(s)

George Owusu

Examples

## Not run: 
################# Point Estimation ###################
# Identify the climatic input data on 
latitude=11.18 #in decimal degrees
Rs= 452 #Average shortwave radiation over the growing period (cal-2day-1)
Tday = 26.9 #	Average day time difference (oC) over the growing period
T24 =25.3 #Average 24 hour mean temperature over the growing period
# Gather crop Information
crop="maize"
plantdate="2015-05-1" #Start date
harvestdate="2015-08-29" #end date 120 days
HI=0.4 #Harvest index for maize
adaptability=3 #crop climate adaptability group III
LAI = 4 #Leaf Area Index
pointmod=biomass(Tday,T24,Rs,latitude,longitude=NULL,HI,
plantdate,harvestdate,c30=0.0108,adaptability=3,LAI=5,crop=NULL)

################# Point Estimation for several locations#################
latitude=c(10,11.18,12,9)
longitude=c(-1.4,-2,-1,-2.6)
Rs= c(450,452,400,500)
Tday = c(25,26.9,24,20,30)
T24=c(24,25.3,23,28)
HI=c(0.4,0.4,0.5,0.25)
plantdate=c("1-03-2015","1-04-2015","1-05-2015")
harvestdate=c("29-08-2015","29-08-2015","29-08-2015")

# computing several points
pointSmod=biomass(Tday,T24,Rs,latitude,longitude=longitude,HI,
plantdate,harvestdate,c30=0.0108,adaptability=3,LAI=5,crop=NULL)
yield=pointSmod$yield #tabular data
biomass=pointSmod$output #tabular data
#access on biomass
pointSmod$output$yield

## End(Not run)
## Not run: 
############## Spatial Estimation with internal interpolation ###########################
#generate spatial data
folder=system.file("extdata","stack",package="sebkc")
modauto=landsat578(data=folder, welev=362)
Tday=modauto$Ts-273.15 #temperature in degree celcius
LAI=modauto$LAI
latitude=seq(7.544,7.590,0.001) #for interpolation
longitude=seq(-1.211,-1.187,0.001) #for interpolation
latitude=latitude[1:length(longitude)] #to make sure they are of the same length
spatialmodint=biomass(Tday,T24,Rs,latitude,longitude=longitude,HI,
plantdate,harvestdate,c30=0.0108,adaptability=3,LAI=LAI,crop=NULL)
# Rs, T24, and HI are interpolated
plot(spatialmodint$yield)
############## Spatial Estimation  ###########################
Rs= 452
T24 =25.3
HI=0.4 
spatialmod=biomass(Tday,T24,Rs,latitude,longitude=longitude,HI,
plantdate,harvestdate,c30=0.0108,adaptability=3,
LAI=LAI,crop=NULL)
plot(spatialmod$yield)

## End(Not run)

Calibration of crop coefficients

Description

This function can estimate up to 23 parameters of FAO single and dual crop parameters. The parameters include soil parameters: WP,FC, REW; crop parameters: kcini,kcmid,kcend,kcbini,kcbmid,kcbend,h,Zr; the evaporation parameters: p and Ze, runoff coefficients: CNII or rc, as well as groundwater and capillary rise parameters: initgw,a1,b1,a2,b2,a3,b3,a4,b4. For each parameter initial, maximum and range values must be provided if one wants to optimise it. There are also 12 goodness of fit tests for the evaluation of the model with fitting :r-square, index of agreement, F-value, p-value, sse,sst, AIC, BIC, RMSE, NRMSE and nashSutcliffe coefficients. A calibrated model that uses more parameters gets punished by AIC and BIC as well as adjusted R2, RMSE, NRMSE.

Both point and spatial model can be calibrated. The actual model is available at kc

Usage

cal.kc(
  ETo,
  P,
  RHmin,
  soil,
  crop,
  I = 0,
  CNII = c(30, 100, 1),
  u2 = 2,
  FC = NULL,
  WP = NULL,
  h = NULL,
  Zr = NULL,
  Ze = c(0.1, 0.15, 0.01),
  kc = c(c(0.1, 0.3, 0.05), c(0.1, 1, 0.1), c(0.1, 1, 0.1)),
  p = c(0.1, 1, 0.1),
  lengths = NULL,
  fw = 0.8,
  fc = NULL,
  kctype = "dual",
  region = NULL,
  regexpr = "sweet",
  initgw = 50,
  LAI = 3,
  rc = NULL,
  Kb = 0,
  Rn = NULL,
  hmodel = "gaussian",
  Zrmodel = "linear",
  xydata = NULL,
  REW = c(2, 12, 0.5),
  a1 = 360,
  b1 = c(-1, 1, 0.1),
  a2 = 240,
  b2 = c(-1, 1, 0.1),
  a3 = c(-2, 1, 0.1),
  b3 = 6.6,
  a4 = 4.6,
  b4 = c(-1.5, 1, 0.1),
  ETa = NULL,
  report = "stages",
  time = "06:00",
  HWR = 1,
  latitude = NULL,
  density = "full",
  DOY = NULL,
  slope = NULL,
  y = NULL,
  r1 = 100,
  nongrowing = "weeds",
  theta = NULL,
  profile = "variable",
  Kcmin = 0.15,
  Ky = NULL
)

Arguments

Details

There are several ways to use cal.kc:

Measured Data

The kc model can be calibrated against either Actual Evapotranspiration (ETa) or Available Soil Water; the user must specify theta. If the calibration is done on a point data then time series for that location must be provided. If multiple locations are provided with xydata and corresponding multiple locations as input data, then ETa or theta must be provided for all locations. See also kc for more details on locations.

Number of optimised parameters and Goodness of fit

There is a need to limit the number of calibrated parameters. One can compare models by observing the adjusted goodness of fit tests. The more parameters that one uses the less are the values of R-square but RMSE, BIC, AIC increases. Unadjusted and adjusted goodness of fit tests are provided with this function.

Value

The function returns optimised parameters as "parameters", goodness of fit as "fit", and calibrated model as "mod". It also prints candidates of parameters and goodness of fit of the local optimum values. For the explanation of calibrated model output and optimised parameters see the return values of kc. Here, only goodness of fit tests explanations are given.

mod:

The optimised model. See return of kc for more details

parameters:

parameters of the optimised model

fit:

Goodness of fit tests. See fitting for explanation of variables

Author(s)

George Owusu

See Also

kc

Examples

## Not run: 
#fewer iteration
file=system.file("extdata","sys","irrigation.txt",package="sebkc")
data=read.table(file,header=TRUE)  
P=data$P
rc=0
I=data$I
ETo=data$ETo
Zr=data$Zr
p=0.6
FC=0.23
WP=0.10
u2=1.6  
RHmin=35
soil="sandy loam"
crop="Broccoli"
#single kc calibration with fewer iteration
cal=cal.kc(ETo,P,RHmin,soil,crop,I=0,kc=NULL,CNII=67,
u2=2,FC=NULL,WP=NULL,h=NULL,Zr=NULL,lengths=NULL,fc=NULL,
kctype="single",region=NULL,regexpr=NULL,initgw=50,LAI=3,
rc=NULL,Rn=NULL,hmodel="gaussian",Zrmodel="linear",
xydata=NULL,ETa=ETo,a1=360,b1=-0.17,a2=240,b2=-0.27,
a3=-1.3,b3=6.6,a4=4.6,b4=-0.65)
parameters=cal$parameters
output=cal$mod$output
gof=cal$fit
ETc=cal$mod$output$ETc
TAW=cal$mod$output$TAW

#Single kc with different kcs. kc uses default values
#RMSE_adj and NRMSE_adj equal to NAN because the number of parameters more than time steps(days)
cal=cal.kc(ETo,P,RHmin,soil,crop,I=0,CNII=67,u2=2,FC=NULL,WP=NULL,
h=NULL,Zr=NULL,lengths=NULL,fc=NULL,kctype="single",region=NULL,
regexpr=NULL,initgw=50,a1=c(100,400,10),b1=c(-1,-0.1,0.1),
a2=c(100,400,10),b2=c(-1,1,0.1),a3=c(-2,1,0.1),b3=c(1,10,1),
a4=c(1,10,1),b4=c(-1.5,1,0.1),LAI=3,rc=NULL,
Rn=NULL,hmodel="gaussian",Zrmodel="linear",xydata=NULL,ETa=ETo)

#Example with fixed kcs
cal=cal.kc(ETo,P,RHmin,soil,crop,I=0,CNII=67,kc=c(0.1,1.15,1),
u2=2,FC=NULL,WP=NULL,h=NULL,Zr=NULL,lengths=NULL,
fc=NULL,kctype="single",region=NULL,regexpr=NULL,initgw=50,
b1=c(-1,-0.1,0.1),b2=c(-1,1,0.1),LAI=3,rc=NULL,Rn=NULL,
hmodel="gaussian",Zrmodel="linear",xydata=NULL,ETa=ETo)

#Example with spatial kcs
#landsat folder with original files but a subset
folder=system.file("extdata","stack",package="sebkc") 
sebiauto=sebi(folder=folder,welev=317,Tmax=31,Tmin=28) #sebi model
kc2=stack(sebiauto$EF/2,sebiauto$EF,sebiauto$EF/3) #assign EF to kc
xydata=data.frame(cbind(longitude=data$longitude,latitude=data$latitude))
modEF=kc(ETo,P=P,RHmin,soil,crop,I,kc=kc2,Rn=8,lengths = c(4,4,2,2),xydata=xydata)
ETa=modEF$ETcxy[3:14]

calsp=cal.kc(ETo,P,RHmin,soil,crop,I=0,CNII=67,kc=kc2,u2=2,FC=NULL,
 WP=WP,h=NULL,Zr=NULL,lengths = c(4,4,2,2),fc=NULL,kctype="single",
 region=NULL,regexpr=NULL,initgw=50,b1=c(-1,-0.1,0.1),b2=c(-1,1,0.1),
 LAI=3,rc=NULL,Rn=7,hmodel="gaussian",Zrmodel="linear",xydata=xydata,ETa=ETa)
 gof=cal$fit
r2=cal$fit$r2
AIC=cal$fit$AIC
ETcxy=cal$mod$ETcxy
TAW=cal$mod$TAW

## End(Not run)

Automatic computation of cold pixel Radiometric Surface Temperature

Description

This function computes temperature of a cold pixel based on NDVI, Surface Temperature, DEM and albedo. It plots histograms indicating the selected values. The main difference with coldTs2 is that preference is given to temperature selection before considering NDVI and albedo selection.

Usage

coldTs(
  Ts,
  NDVI,
  albedo = NULL,
  sunangle = NULL,
  DEM = NULL,
  cluster = 8,
  extent = "interactive",
  upper = 0.95,
  lower = 0.2,
  plot = TRUE,
  layout = "portrait",
  draw = "poly",
  folder = NULL,
  welev = NULL,
  clip = NULL,
  iter.max = 100,
  seed = NULL
)

## Default S3 method:
coldTs(
  Ts = NULL,
  NDVI,
  albedo = NULL,
  sunangle = NULL,
  DEM = NULL,
  cluster = 8,
  extent = "interactive",
  upper = 0.95,
  lower = 0.2,
  plot = TRUE,
  layout = "portrait",
  draw = "poly",
  folder = NULL,
  welev = NULL,
  clip = NULL,
  iter.max = 100,
  seed = NULL
)

Arguments

Details

The function first divides the selected area into a number of clusters. The mean and standard devaitions are computed. The cluster means with the lowest standard devaitions and qualifying for upper (NDVI) and lower (Ts) limits (see above) are slected. The provided upper and lower probabilities are also used to select the final Ts.

Value

Ts:

The selected or input Ts [K]

Tscold:

The Temperature of cold pixe [K]

TC:

The average temperature of the cold pixels [K]

x:

The x cordinate of the cold pixel

y:

The y cordinate of the cold pixel

xycold:

combined x and y coordinates

candidates:

The similiar candidates' pixels that can be used

Author(s)

George Owusu

References

Allen, R. G., Burnett, B., Kramber, W., Huntington, J., Kjaersgaard, J., Kilic, A., Kelly, C., & Trezza, R. 2013. Automated Calibration of the METRIC-Landsat Evapotranspiration Process. JAWRA Journal of the American Water Resources Association, 49(3): 563-576.

See Also

hotTs

Examples

## Not run: 
#using landsat folder
folder=system.file("extdata","stack",package="sebkc")
modcold=coldTs(folder=folder,welev=170,extent="auto")

#use object of  
folder=system.file("extdata","stack",package="sebkc")
data=landsat578(data=folder, welev=362)
modcold=coldTs(folder=data,welev=170,extent="auto",cluster=3)

modcoldTs=coldTs(data,welev=170,extent="auto",cluster=3)

#using input different input data
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"))
#simulate with default parameters without albedo
modcold=coldTs(Ts=Ts,NDVI=NDVI,sunangle=50.7)
#simulate with albedo but set extent to "auto" and cluster to 7
modcold=coldTs(Ts=Ts,NDVI=NDVI,albedo=albedo,sunangle=50,cluster=7,
extent="auto")

## End(Not run)

Optional Automatic computation of cold pixel Radiometric Surface Temperature

Description

This function computes temperature of a cold pixel based on NDVI, Surface Temperature, DEM and albedo. The main difference with coldTs is that preference is given to temperature selection before considering NDVI and albedo selection. Moreover, this is a simple end member selection where no clustering is done.

Usage

coldTs2(
  Ts,
  NDVI,
  albedo = NULL,
  sunangle = NULL,
  DEM = NULL,
  cluster = 10,
  extent = "interactive",
  upper = 0.2,
  lower = 0.1,
  NDVI.probs = 0.95,
  plot = TRUE,
  layout = "portrait",
  draw = "poly",
  folder = NULL,
  welev = NULL
)

## Default S3 method:
coldTs2(
  Ts = NULL,
  NDVI,
  albedo = NULL,
  sunangle = NULL,
  DEM = NULL,
  cluster = 10,
  extent = "interactive",
  upper = 0.2,
  lower = 0.1,
  NDVI.probs = 0.6,
  plot = TRUE,
  layout = "portrait",
  draw = "poly",
  folder = NULL,
  welev = NULL
)

Arguments

Details

The function first divides the selected area into a number of clusters. The mean and standard deviations are computed. The cluster means with the lowest standard devaitions and qualifying for upper (NDVI) and lower (Ts) limits (see above) are selected. The provided upper and lower probabilities are also used to select the final Ts.

Value

Ts:

The selected or input Ts [K]

Tscold:

The Temperature of cold pixel [K]

TC:

The average temperature of the cold pixels [K]

x:

The x coordinate of the cold pixel

y:

The y coordinate of the cold pixel

xycold:

combined x and y coordinates

candidates:

The similar candidates' pixels that can be used

Author(s)

George Owusu

References

Allen, R. G., Burnett, B., Kramber, W., Huntington, J., Kjaersgaard, J., Kilic, A., Kelly, C., & Trezza, R. 2013. Automated Calibration of the METRIC-Landsat Evapotranspiration Process. JAWRA Journal of the American Water Resources Association, 49(3): 563-576.

See Also

coldTs, hotTs, hotTs2

Examples

## Not run: 
#using landsat folder
folder=system.file("extdata","stack",package="sebkc")
modcold=coldTs2(folder=folder,welev=170,extent="auto")

#use object of  
folder=system.file("extdata","stack",package="sebkc")
data=landsat578(data=folder, welev=362)
modcold=coldTs2(folder=data,welev=170,extent="auto",cluster=3)

modcoldTs2=coldTs2(data,welev=170,extent="auto",cluster=3)

#using input different input data
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"))
#simulate with default parameters without albedo
modcold=coldTs2(Ts=Ts,NDVI=NDVI,sunangle=50.7)
#simulate with albedo but set extent to "auto" and cluster to 7
modcold=coldTs2(Ts=Ts,NDVI=NDVI,albedo=albedo,sunangle=50,cluster=7,
extent="auto")

## End(Not run)

Automatic computation of hot and cold pixels from Radiometric Surface Temperature (Ts)

Description

This function computes temperature coordinates of a hot and cold pixels based on only Surface Temperature. This function uses ranges of temperature.

Usage

coldhot(Tmin, Tmax, Ts = NULL, folder = NULL, welev = NULL)

## Default S3 method:
coldhot(Tmin, Tmax, Ts = NULL, folder = NULL, welev = NULL)

Arguments

Details

The function first divides the selected area into a number of clusters. The mean and standard deviations are computed. The cluster means with the lowest standard deviations and qualifying for upper (Ts) and lower (NDVI) limits (see above) are selected. The provided upper and lower probabilities are also used to select the final Ts.

Value

Ts:

The selected or input Ts

Tshot:

The Temperature of hot pixel

Tscold:

The Temperature of cold pixel

xyhot:

combined x and y coordinates

candidates:

The similar candidates' pixels that can be used

Author(s)

George Owusu

References

Allen, R. G., Burnett, B., Kramber, W., Huntington, J., Kjaersgaard, J., Kilic, A., Kelly, C., & Trezza, R. 2013. Automated Calibration of the METRIC-Landsat Evapotranspiration Process. JAWRA Journal of the American Water Resources Association, 49(3): 563-576.

See Also

coldTs, hotTs, coldTs2

Examples

## Not run: 
 folder=system.file("extdata","stack",package="sebkc")
  model=coldhot(Tmin=300,Tmax=310,folder=folder,welev=170)
  spplot(model$Tshot)
  spplot(model$Tscold)
  modauto=sebal(folder = folder,welev = 380,xycold=model$xycold,xyhot=model$xyhot)
## End(Not run)

Reprojecting and cropping images

Description

This function returns the subset of x object as defined by y extent object. See crop and extent for more details. The difference between crop and cropsebkc is the ability for the latter to reproject y object with x object projection.

Usage

cropsebkc(x, y)

Arguments

Value

cropped x and reprojected y

Author(s)

George Owusu

See Also

crop

Examples

## Not run: 
folder=system.file("extdata","stack",package="sebkc")
stack=sebkcstack(folder=folder)
data=cropsebkc(stack$data,stack$data)

## End(Not run)

Day length calculation

Description

Day length calculation

Usage

dl(latitude, DOY, model = "CBM", Tmax = NULL, Tmin = NULL, p = 0.5)

Arguments

Value

DL

Author(s)

George Owusu

References

• Schoolfield R. (1982). Expressing daylength as a function of latitude and Julian date.

• Running S. W. and Coughlan J. C. (1988). A general model of forest ecosystem processes for regional applications I. Hydrologic balance, canopy gas exchange and primary production processes. Ecological Modelling, 42(2), 125-154 http://dx.doi.org/10.1016/0304-3800(88)90112-3.

• Ritchie J. T. (1991). Wheat Phasic Development. Modeling Plant and Soil Systems. Agronomy Monograph, 31.

• Allen R. G., Pereira L. S., Raes D. and Smith M. (1998). Crop evapotranspiration: Guidelines for computing crop water requirements. FAO Irrigation and Drainage Paper, 56, 300.

Examples

DOY="2001-8-1"
latitude=0
model="CBM"
Tmax=31
Tmin=26
mod=dl(latitude,DOY,model,Tmax,Tmin=Tmin)

Goodness of fit Test

Description

Goodness of fit Test

Usage

fitting(obs, est, P = 1)

Arguments

Value

r2:

R-square [-]

r2_adj:

An adjusted R-square [-]

d:

Index of agreement[-]

p_value:

p value

sse:

Sum of square error

sst:

Sum of square total

AIC:

AIC

BIC:

BIC

RMSE:

Root Mean Square Error

RMSE_adj:

adjusted Root Mean Square Error

NRMSE

Normalised Root Mean Square Error

NRMSE_adj:

adjusted Normalised Root Mean Square Error

E:

Nash Sutcliffe Coefficient

b:

coefficients of regression

AAE:

average absolute error

ARE:

average relative error

MBE:

Mean Bias error

MPE:

Mean Percentage error

Examples

file=system.file("extdata","sys","irrigation.txt",package="sebkc")
data=read.table(file,header=TRUE)  
obs=data$ETo
est=data$ETc
P=4
mod=fitting(obs,est,P)
r2=mod$r2

Automatic computation of hot pixel Radiometric Surface Temperature (Ts)

Description

This function computes temperature of a hot pixel based on NDVI, Surface Temperature, DEM and albedo. The main difference with hotTs2 is that preference is given to temperature selection before considering NDVI and albedo selection. Moreover, this is a simple end member selection where no clustering is done.

Usage

hotTs(
  Ts,
  NDVI,
  albedo = NULL,
  DEM = NULL,
  cluster = 8,
  extent = "interactive",
  upper = 0.8,
  lower = 0.2,
  plot = TRUE,
  layout = "portrait",
  draw = "poly",
  folder = NULL,
  welev = NULL,
  clip = NULL,
  iter.max = 500,
  seed = NULL
)

## Default S3 method:
hotTs(
  Ts = NULL,
  NDVI,
  albedo = NULL,
  DEM = NULL,
  cluster = 8,
  extent = "interactive",
  upper = 0.8,
  lower = 0.2,
  plot = TRUE,
  layout = "portrait",
  draw = "poly",
  folder = NULL,
  welev = NULL,
  clip = NULL,
  iter.max = 500,
  seed = NULL
)

Arguments

Details

The function first divides the selected area into a number of clusters. The mean and standard deviations are computed. The cluster means with the lowest standard deviations and qualifying for upper (Ts) and lower (NDVI) limits (see above) are selected. The provided upper and lower probabilities are also used to select the final Ts.

Value

Ts:

The selected or input Ts

Tshot:

The Temperature of hot pixel

TX:

The average temperature of the hot pixels

x:

The x coordinate of the hot pixel

y:

The y coordinate of the hot pixel

xyhot:

combined x and y coordinates

candidates:

The similar candidates' pixels that can be used

Author(s)

George Owusu

References

Allen, R. G., Burnett, B., Kramber, W., Huntington, J., Kjaersgaard, J., Kilic, A., Kelly, C., & Trezza, R. 2013. Automated Calibration of the METRIC-Landsat Evapotranspiration Process. JAWRA Journal of the American Water Resources Association, 49(3): 563-576.

See Also

coldTs

Examples

## Not run: 
 folder=system.file("extdata","stack",package="sebkc")
  modhot=hotTs(folder=folder,welev=170,extent="auto",cluster=3)
  
   #use object of  landsat578
folder=system.file("extdata","stack",package="sebkc")
data=landsat578(data=folder, welev=170)
modfhot=hotTs(folder=data,welev=170,extent="auto",cluster=3)

#Or Ts = data
modhotTs=hotTs(data,welev=170,extent="auto",cluster=3)
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"))
#simulatte with default parameters without albedo
modhot=hotTs(Ts=Ts,NDVI=NDVI,extent="auto",cluster=3)
#simulate with albedo but set extent to "auto" and cluster to 7
modhot=hotTs(Ts=Ts,NDVI=NDVI,albedo=albedo,cluster=7,
extent="auto")

## End(Not run)

Automatic computation of hot pixel Radiometric Surface Temperature (Ts)

Description

This function computes temperature of a hot pixel based on NDVI, Surface Temperature, DEM and albedo. It plots histograms indicating the selected values

Usage

hotTs2(
  Ts,
  NDVI,
  albedo = NULL,
  sunangle = NULL,
  DEM = NULL,
  cluster = 10,
  extent = "interactive",
  upper = 0.95,
  lower = 0.8,
  NDVI.probs = 0.1,
  plot = TRUE,
  layout = "portrait",
  draw = "poly",
  folder = NULL,
  welev = NULL
)

## Default S3 method:
hotTs2(
  Ts = NULL,
  NDVI,
  albedo = NULL,
  sunangle = NULL,
  DEM = NULL,
  cluster = 10,
  extent = "interactive",
  upper = 0.95,
  lower = 0.8,
  NDVI.probs = 0.2,
  plot = TRUE,
  layout = "portrait",
  draw = "poly",
  folder = NULL,
  welev = NULL
)

Arguments

Details

The function first divides the selected area into a number of clusters. The mean and standard deviations are computed. The cluster means with the lowest standard deviations and qualifying for upper (Ts) and lower (NDVI) limits (see above) are selected. The provided upper and lower probabilities are also used to select the final Ts.

Value

Ts:

The selected or input Ts

Tshot:

The Temperature of hot pixel

TX:

The average temperature of the hot pixels

x:

The x cordinate of the hot pixel

y:

The y cordinate of the hot pixel

xyhot:

combined x and y coordinates

candidates:

The similar candidates' pixels that can be used

Author(s)

George Owusu

References

Allen, R. G., Burnett, B., Kramber, W., Huntington, J., Kjaersgaard, J., Kilic, A., Kelly, C., & Trezza, R. 2013. Automated Calibration of the METRIC-Landsat Evapotranspiration Process. JAWRA Journal of the American Water Resources Association, 49(3): 563-576.

See Also

coldTs, hotTs, coldTs2

Examples

## Not run: 
 folder=system.file("extdata","stack",package="sebkc")
  modhot=hotTs2(folder=folder,welev=170,extent="auto",cluster=3)
  
   #use object of  landsat578
folder=system.file("extdata","stack",package="sebkc")
data=landsat578(data=folder, welev=170)
modfhot=hotTs2(folder=data,welev=170,extent="auto",cluster=3)

#Or Ts = data
modhotTs2=hotTs2(data,welev=170,extent="auto",cluster=3)
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"))
#simulate  with default parameters without albedo
modhot=hotTs2(Ts=Ts,NDVI=NDVI,extent="auto",cluster=3)
#simulate with albedo but set extent to "auto" and cluster to 7
modhot=hotTs2(Ts=Ts,NDVI=NDVI,albedo=albedo,cluster=7,
extent="auto")

## End(Not run)

Inverse distance interpolation

Description

This function performs inverse distance interpolation based on projection and geometry of input spatial data 'ext'.

Usage

invdist(longitude, latitude, var, ext, idp = 0.5)

Arguments

Value

var: interpolated values

Author(s)

George Owusu

Examples

 ## Not run: 
#Get ext data
folder=system.file("extdata","stack",package="sebkc")
data=landsat578(data=folder, welev=362)
temp=data$Ts
latitude=seq(7.544,7.590,0.001)
longitude=seq(-1.211,-1.187,0.001)
var=seq(20,30,0.2)
latitude=latitude[1:length(longitude)]
var=var[1:length(longitude)]
map=invdist(longitude,latitude,var,ext=temp)

## End(Not run)

Estimating FAO56 Single and Dual Crop coefficients with a water balance model in R

Description

This function computes FAO56 crop coefficients and water balance. The output of the model includes crop evapotranspiration, water stress coefficient, irrigation requirements, soil available water, etc. In addition to point estimation, spatial modelling is integrated using Evaporation fractions from sebal, sebi, ETo, sseb, sebs. cal.kc can be used to calibrate the model.

Usage

kc(
  ETo,
  P,
  RHmin,
  soil = "Sandy Loam",
  crop,
  I = 0,
  CNII = 67,
  u2 = 2,
  FC = NULL,
  WP = NULL,
  h = NULL,
  Zr = NULL,
  Ze = 0.1,
  kc = NULL,
  p = NULL,
  lengths = NULL,
  fw = 0.8,
  fc = NULL,
  kctype = "dual",
  region = NULL,
  regexpr = "sweet",
  initgw = 50,
  LAI = 3,
  rc = NULL,
  Kb = 0,
  Rn = NULL,
  hmodel = "gaussian",
  Zrmodel = "linear",
  xydata = NULL,
  REW = NULL,
  a1 = 360,
  b1 = -0.17,
  a2 = 240,
  b2 = -0.27,
  a3 = -1.3,
  b3 = 6.6,
  a4 = 4.6,
  b4 = -0.65,
  report = "stages",
  time = "06:00",
  HWR = 1,
  latitude = NULL,
  density = "full",
  DOY = NULL,
  slope = NULL,
  y = NULL,
  r1 = 100,
  nongrowing = "weeds",
  theta = NULL,
  profile = "variable",
  Kcmin = 0.15,
  Ky = NULL
)

## Default S3 method:
kc(
  ETo,
  P,
  RHmin,
  soil = "Sandy Loam",
  crop,
  I = 0,
  CNII = 67,
  u2 = 2,
  FC = NULL,
  WP = NULL,
  h = NULL,
  Zr = NULL,
  Ze = 0.1,
  kc = NULL,
  p = NULL,
  lengths = NULL,
  fw = 1,
  fc = NULL,
  kctype = "dual",
  region = NULL,
  regexpr = "sweet",
  initgw = 50,
  LAI = 3,
  rc = NULL,
  Kb = 0,
  Rn = NULL,
  hmodel = "gaussian",
  Zrmodel = "linear",
  xydata = NULL,
  REW = NULL,
  a1 = 360,
  b1 = -0.17,
  a2 = 240,
  b2 = -0.27,
  a3 = -1.3,
  b3 = 6.6,
  a4 = 4.6,
  b4 = -0.65,
  report = "stages",
  time = "06:00",
  HWR = 1,
  latitude = NULL,
  density = "full",
  DOY = NULL,
  slope = NULL,
  y = NULL,
  r1 = 100,
  nongrowing = "weeds",
  theta = NULL,
  profile = "variable",
  Kcmin = 0.15,
  Ky = NULL
)

## S3 method for class 'kc'
plot(
  x,
  output = "AW_measured",
  main = NULL,
  cex = 1,
  legend = TRUE,
  pos = "top",
  horiz = FALSE,
  ...
)

fit.kc(x, output = "AW_measured")

Arguments

Details

Using Existing Tables

The function include backend database of FAO Tables. The user can use the parameter values such as Zr, h, p, FC, WP, etc from FAO56 tables. Most of the parameters with NULL values can take default values from FAO56 Tables. When a crop name results in 3 crops, the user can differentiate that with regexpr parameter. For example, a crop parameter with a name "maize" can match two crops on FAO Table 12. Including regexpr parameter with "sweet" with match only one.

Modelling crop height and Rooting Depth

At each time step a crop height or root depth is estimated. The crop or rooting depth growth model can take exponential model of the form hday=0+0.75*max(h)*(1-exp(-day/0.75*max(days))) or Zrday=0+0.75*Zr*(1-exp(-day/0.75*max(days))) The gaussian model is of the form hday=0+0.75*max(h)*(1-exp(-day^2/(0.75*max(days)^2))) or Zrday=0+0.75*(1-exp(-day^2/(0.75*max(day))^2)) The linear model is of the form hday=(day*max(h))/max(days), Zrday=(day*Zr)/max(days)

Point modeling with Single and dual crop coefficients

The function can estimate kc curve for the growing period by using the traditional kcini, kcmid and kcend and crop lengths growth stages, in case of single crop coefficients. Just like single crop coefficients, the dual crop coefficients kcbini, kcbmid, kcbend are embedded in kc parameter. The "kc" parameter must take 3 values as c(initial, mid, end). The values of "lengths" parameter must also be specified as c(initial, dev, mid, late). See example section for more details.

Spatial modeling With Evaporation Fraction

Instead of supplying kcs values as described above, Evaporation Fraction (EF) from can be used. The corresponding EF at initital, mid amd end seasons can be supplied. The EF curve will be developed for each location at each time step and ETa derived from (EF*Rn)/2.45.

Interpolation of Input data

Once EF is used instead of Kcs, one can provide longitudes and latitudes in xydata parameter data frame. A 2D column data of P, I,ETo,RHmin, u2, Rn for each time step for each location can be provided. At each time step an inverse distance interpolation will be perform to get continuous surface. In case of initial and static input column data such as CNII, FC, WP, initgw, Zr, h can also be interpolated.

Calibration

A calibration of unmeasured variables can be done with cal.kc

Value

The output of the function depends of the kctype and kc data type. For single and dual crop coefficient, an "output" dataframe contains variables such as TAW,RAW,Runoff,Ks,Kcb,h,Kc,ETc,DP,CR,Drend,Kr,De,AW,E,Ke,CN. In case kc is set to EF the function returns Kc,Kcxy,ETc,ETcxy,Drend,Drendxy,TAW,TAWxy,RAW,RAWxy,AW,Awcxy,Ineed,Ineedxy

output:

dataframe of the return values

ETc_adj:

The soil water stress adjusted ETc [mm]

Ya_by_Ym:

actual yield by potential yield

yield_decrease:

The relative yield decrease: 1-Ya/Ym

ETdecrease:

The relative evapotranspiration deficit :1- ETc_adj/ETc

yield_by_ET_decrease:

he relative yield decrease: 1-Ya/Ym and the relative evapotranspiration deficit in one dataframe

Awcxy:

Available Soil water for each location at each time step [mm]

CR:

capillary rise [mm/day]

GW:

Updated ground water depth

CN:

Curve Number

De:

cumulative depth of evaporation (depletion) from the soil surface layer [mm] at the end of the day [mm/day]

DP:

Deep Percolation [mm/day]

Drend:

cumulative depth of evapotranspiration (depletion) from the root zone at the end of the day [mm]. This is equal to the amount of water needed to bring soil water to field capacity

Drendxy:

cumulative depth of evapotranspiration (depletion) from the root zone at the end of the day [mm] for each xydata location at each time step. This is equal to the amount of water needed to bring soil water to field capacity

E:

evaporation [mm/day]

ETc:

crop evapotranspiration [mm/day]. If water stress coefficient (Ks) is less than 1, ETc is crop evapotranspiration under non-standard conditions else ETc is crop evapotranspiration under standard conditions

ETcxy:

crop evapotranspiration [mm/day] for each location at each time step. crop evapotranspiration under standard conditions

h:

crop height [m]

Ineed:

Amount of soil water that is below maximum allowed depletion [mm]

Ineedxy:

Amount of soil water that is below maximum allowed depletion for each xydata location at each time step [mm]

Kc:

standard crop coefficient [-]

Kc_adj:

adjusted crop coefficient [-]

EF:

Evaporation Fraction [-]

EFxy:

Evaporation Fraction for each monitored location [-]

kcb:

basal crop coefficient for the growing season [-]

Kcxy:

crop coefficient for each location at each time step [-]

Ke:

soil evaporation coefficient

Kr:

soil evaporation reduction coefficient [-]

Ks:

Water stress coefficient [-]

RAW:

Readily Available Water [mm]

RAwcxy:

Readily Available Soil water for each xydata location at each time step [mm]

Runoff:

Runoff [mm/day]

TAW:

Total Available Water [mm]

TAwcxy:

Total Available Soil water for each location at each time step [mm]

Author(s)

George Owusu

References

• Allen, R. G., PEREIRA, L. S., RAES, D., & SMITH, M. (1998). Crop Evapotranspiration (guidelines for computing crop water requirements) FAO Irrigation and Drainage Paper No. 56: FAO.

• Allen, R.G., Wright, J.L., Pruitt, W.O., Pereira, L.S., Jensen, M.E., 2007. Water requirements. In: Hoffman, G.J., Evans, R.G., Jensen, M.E., Martin, D.L., Elliot, R.L. (Eds.), Design and Operation of Farm Irrigation Systems. , 2nd edition. ASABE, St. Joseph, MI, pp. 208-288.

• Liu, Y., Pereira, L.S., Fernando, R.M., 2006. Fluxes through the bottom boundary of the root zone in silty soils: parametric approaches to estimate groundwater contribution and percolation. Agric. Water Manage. 84, 27-40.

See Also

cal.kc

Examples

## Not run: 
#FAO56 Example 32: Calculation of the crop coefficient (Kcb + Ke) under sprinkler irrigation
FAO56Example32=kc(ETo=7,P=0,RHmin=20,soil="sandy loam",crop="Broccoli",I=10,
kctype = "dual",u2=3,kc=c(0.9,0,0),h=1,Zr=0.1)
FAO56Example32$output$fc
FAO56Example32$output$Kc_adj
#FAO56 Example 33: Calculation of the crop coefficient (Kcb + Ke) under furrow irrigation
FAO56Example33=kc(ETo=7,P=0,RHmin=20,soil="sandy loam",crop="Broccoli",I=10,
         kctype = "dual",u2=3,kc=c(0.9,0,0),h=1,Zr=0.1,fw=0.3)
#FAO56 Example 34: Calculation of the crop coefficient (Kcb + Ke) under drip irrigation
FAO56Example34=kc(ETo=7,P=0,RHmin=20,soil="sandy loam",crop="Broccoli",I=10,
           kctype = "dual",u2=3,kc=c(0.9,0,0),h=1,Zr=0.1,fw="drip")
         
file=system.file("extdata","sys","irrigation.txt",package="sebkc")
data=read.table(file,header=TRUE)  
P=data$P
rc=0
I=data$I
ETo=data$ETo
Zr=data$Zr
p=0.6
FC=0.23
WP=0.10
u2=1.6  
RHmin=35
soil="sandy loam"
crop="Broccoli"
############Single crop coefficient model###########
modsingle=kc(ETo,P=P,RHmin,soil,crop,I)
#dataframe of the return values
modsingle$output 
modsingle$output$TAW #Total Avalaible water [mm]
#Single with runoff coefficient
############Dual crop coefficient model###########
moddual=kc(ETo,P=P,RHmin,soil,crop,I,kctype = "dual",FC=FC,p=p,WP=WP,u2=u2,kc=c(0.3,0,0),Zr=Zr)
#Reproducing FAO56 Example 36
fc=(0.00667*1:12)+0.07333
h=1:12/1:12*0.3
 moddual=kc(ETo,P=P,RHmin,soil,crop,I,kctype = "dual",FC=FC,p=p,WP=WP,u2=u2,rc=0,
 kc=data$Kcb,Zr=Zr,lengths = c(4,4,2,2),h=h,hmodel="linear",fc=fc,fw=0.8)
 plot(moddual$output$Day,moddual$output$Kc_adj,type="l",ylim=c(0,1.2))
 lines(moddual$output$Day,moddual$output$Kcb,type="l", lty=2)
###########Spatial model I: Evaporation Fraction (EF) Model###################
#landsat folder with original files but a subset
folder=system.file("extdata","stack",package="sebkc") 

sebiauto=sebi(folder=folder,welev=317,Tmax=31,Tmin=28) #sebi model
kc2=stack(sebiauto$EF/2,sebiauto$EF,sebiauto$EF/3) #assign EF to kc
modEF=kc(ETo,P=P,RHmin,soil,crop,I,kc=kc2,Rn=8,lengths = c(4,4,2,2))
spplot(modEF$EF)

#############Spatial model II: interpolating precipitation ######################
h=1:12
Zr=1:12/12
xydata=data.frame(cbind(longitude=data$longitude,latitude=data$latitude))
print(xydata) #take a look at xydata
P2=data.frame(matrix(rep(P,12),ncol=12))
print(P2) #take a look at the format of precipitation data
modEFh=kc(ETo,P=P2,RHmin,soil,crop,I,kc=kc2,Rn=8,lengths = c(4,4,2,2),h=h,Zr=Zr,xydata=xydata)
#Amount of irrigation water needed 
spplot(modEFh$Drend)
#Amount of water Irrigation that is critically needed to reach maximum allowed depletion level
modEFh$Ineedxy
ETcxy=modEFh$ETcxy ##get ETc for each location at each time step
TAW=modEFh$TAW ## 


############FAO56 Example 41: Estimation of mid-season crop coefficient for Trees or weeds ########
FAOExample41 =kc(fc=0.3,DOY=200,latitude=40,u2=1.5,h=2,density="row",lengths=c(0,0,1,0),
ETo=0,P=0,RHmin=55,soil="sand",crop="Broccoli",kctype="dual")
FAOExample41$output$Kcb

#FAO56 Example 43: Estimation of Kcb mid from ground cover with reduction for stomatal control
FAOExample43 =kc(fc=0.196,DOY=180,latitude=30,u2=2,h=5,density="round",y=0.0676,slope=0.189,
r1=420,RHmin=25,lengths=c(0,0,1,0),ETo=0,P=0,soil="sand",crop="Broccoli",kctype="dual")
FAOExample43$output$Kcb

## End(Not run)

Historical Weather data of Kumasi from 1960 to 2015

Description

A dataset containing Rain, maximum and minimum temperature, wind speed, sunshine hours, and minimum and maximum Relative Humidity.

Format

A data frame with 20454 rows and 8 variables:

DOY

Day of the Year in dates

Tmax

Maximum Temperature [oC]

Tmin

Minimum Temperature [oC]

RHmin

Minimum Relative Humidity [per cent]

RHmax

Maximum Relative Humidity [per cent]

n

sunshine hours

uz

wind speed [m/s]

P

Daily precipitation [mm]

Author(s)

George Owusu

A 2.6 RCP downscaled Weather data of Kumasi from 2016 to 2050

Description

A dataset containing Rain, maximum and minimum temperature, wind speed, sunshine hours, and minimum and maximum Relative Humidity. The RCPs 2.6 assumes that green house gas will peak in between 2010-2020 and declines throughout 21st Century.

Format

A data frame with 12784 rows and 8 variables:

DOY

Day of the Year in dates

Tmax

Maximum Temperature [oC]

Tmin

Minimum Temperature [oC]

RHmin

Minimum Relative Humidity [per cent]

RHmax

Maximum Relative Humidity [per cent]

n

sunshine hours

uz

wind speed [m/s]

P

Daily precipitation [mm]

Author(s)

George Owusu

A 4.5 RCP downscaled Weather data of Kumasi from 2016 to 2050

Description

A dataset containing Rain, maximum and minimum temperature, wind speed, sunshine hours, and minimum and maximum Relative Humidity. The RCP 4.5 assumes that green house gas will peak in 2040 and declines throughout 21st Century.

Format

A data frame with 12784 rows and 8 variables:

DOY

Day of the Year in dates

Tmax

Maximum Temperature [oC]

Tmin

Minimum Temperature [oC]

RHmin

Minimum Relative Humidity [per cent]

RHmax

Maximum Relative Humidity [per cent]

n

sunshine hours

uz

wind speed [m/s]

P

Daily precipitation [mm]

Author(s)

George Owusu

A 8.5 RCP downscaled Weather data of Kumasi from 2016 to 2050

Description

A dataset containing Rain, maximum and minimum temperature, wind speed, sunshine hours, and minimum and maximum Relative Humidity. The RCP 8.5 assumes that green house gas will rise throughout 21st Century.

Format

A data frame with 12784 rows and 8 variables:

DOY

Day of the Year in dates

Tmax

Maximum Temperature [oC]

Tmin

Minimum Temperature [oC]

RHmin

Minimum Relative Humidity [per cent]

RHmax

Maximum Relative Humidity [per cent]

n

sunshine hours

uz

wind speed [m/s]

P

Daily precipitation [mm]

Author(s)

George Owusu

Radiometric correction and processing of Landsat Images

Description

landsat578 performs radiometric collection of landsat 5,7,and 8 images. It returns albedo, NDVI, fractional vegetation cover, surface temperature SAVI,radiance and reflectance of the image

Usage

landsat578(
  data,
  welev,
  sensor = NULL,
  gain = "low",
  sunelev = NULL,
  date = NULL,
  weights = "auto",
  grescale = "auto",
  brescale = "auto",
  ESUNx = "auto",
  Mp = NULL,
  Ap = NULL,
  K1 = NULL,
  K2 = NULL,
  Rp = 0.91,
  Rsky = 1.32,
  tNB = 0.886,
  DOY = NULL,
  clip = NULL,
  folder = NULL
)

## Default S3 method:
landsat578(
  data = NULL,
  welev,
  sensor = NULL,
  gain = "low",
  sunelev = NULL,
  date = NULL,
  weights = "auto",
  grescale = "auto",
  brescale = "auto",
  ESUNx = "auto",
  Mp = NULL,
  Ap = NULL,
  K1 = NULL,
  K2 = NULL,
  Rp = 0.91,
  Rsky = 1.32,
  tNB = 0.886,
  DOY = NULL,
  clip = NULL,
  folder = NULL
)

Arguments

Details

If you set the data to landsat folder that contains landsat band and meta data sebkcstack is used to compute the parameters. The user must specify welev and may use default parameters. The data must always be a raster brick or stack. To use the default values please order bind the bands c(band1,band2,band3,band4,band5,band6,band7). The defaults values are for landsat 5 and 7. For landsat 8, use the meta data.

Value

albedo:

albedo values

Ts:

Surface Temperature

NDVI:

NDVI values

SAVI:

SAVI values

fc:

Fractional vegetation cover

radiance:

radiance values

reflectance:

TOA reflectance

Author(s)

George Owusu

References

• Chander, G., Markham, B. L., & Helder, D. L. 2009. Summary of current radiometric calibration coefficients for Landsat MSS, TM, ETM+, and EO-1 ALI sensors. Remote Sensing of Environment, 113(5): 893-903.

• Waters, R. 2002. SEBAL: Surface Energy Balance Algorithms for Land. Idaho Implementation, Advanced Training and Users Manual. Kimberly, Idaho: University of Idaho.

• Survey, U. S. G. 2015. Landsat 8 (L8) data users handbook: 97p. from https://landsat.usgs.gov/documents/Landsat8DataUsersHandbook.pdf

Examples

## Not run: 
#Automatic detection of parameters
folder=system.file("extdata","stack",package="sebkc")
modauto=landsat578(data=folder, welev=362)

#Manual estimation with input parameters
landsat7_6feb2004=brick(system.file("extdata","landsat7_6feb2004.grd",
package="sebkc"))
mod=landsat578(data=landsat7_6feb2004,welev=317.1,sensor=7,gain="high",
sunelev=50.71154048,weights="auto",grescale="auto",brescale="auto",
ESUN="auto",K1=666.09,K2=1282.71,date="2002-2-6")
#acces albedo value
albedo=mod$albedo

## End(Not run)

Reproject WGS 84 longitude and latitude into a new raster coordinates and vice versa

Description

Reproject WGS 84 longitude and latitude into a new raster coordinates and vice versa

Usage

lonlatreproject(x, y = NULL, var = NA, map, epsg = 4326)

Arguments

Details

If x is a list of longitudes then new horizontal coordinates based on map coordinates. If x is a raster object with horizontal coordinates then it will be projected to WGS 84.

Value

New x, y and map

Author(s)

George Owusu

Examples

folder=system.file("extdata","stack",package="sebkc")
data=landsat578(data=folder, welev=362)
temp=data$Ts
latitude=seq(7.544,7.590,0.001)
longitude=seq(-1.211,-1.187,0.001)
var=seq(20,30,0.2)
latitude=latitude[1:length(longitude)]
var=var[1:length(longitude)]
newlonlat=lonlatreproject(longitude,latitude,var=var,map=temp)
longitude=newlonlat$longitude
map=newlonlat$map

Internal function for doing interpolation

Description

Internal function for doing interpolation

Usage

mapmatrix(xydata, i, map)

Arguments

Value

map

Author(s)

George Owusu

Examples

## Not run: 
folder=system.file("extdata","stack",package="sebkc")
sebiauto=sebi(folder=folder,welev=317,Tmax=31,Tmin=28)
points=sampleRandom(sebiauto$EF,100,sp=TRUE)
pt=cbind(points@coords,points@data)
longitude=pt$x
latitude=pt$y
value=pt$layer
mapmatrix=mapmatrix(xydata,2)
map=invdist(longitude,latitude,var=value,ext=sebiauto$EF)

## End(Not run)

Spatial Conditional function on raster object

Description

It computes equivalent of logical if function in R

Usage

rastercon(condition, trueValue, falseValue)

Arguments

Value

Return True or False

Author(s)

George Owusu

Examples

## Not run: 
NDVI=raster(system.file("extdata","NDVI.grd",package="sebkc"))
albedo=raster(system.file("extdata","albedo.grd",package="sebkc"))
LAI=raster(system.file("extdata","LAI.grd",package="sebkc"))
  eNB=rastercon( NDVI<0 & albedo<0.47,0.99,rastercon(LAI>=3,0.98,0.97+(LAI*0.0033)))
  
## End(Not run)

R computation of Surface Energy Balance Algorithm for Land and METRIC

Description

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

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

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

Author(s)

George Owusu

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.

Examples

## Not run: 
#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)

## End(Not run)

SEBI Evapotranspiration Fraction Estimation in R

Description

sebi computes Evapotranspiration Fraction Estimation based on Menenti& Choudhury (1993) parameterization of land surface evaporation by means of location dependent potential evaporation and surface temperature range.

Usage

sebi(
  albedo,
  Ts,
  Tmax,
  Tmin,
  RHmax = NULL,
  RHmin = NULL,
  latitude = NULL,
  NDVI,
  SAVI,
  iter.max = 7,
  xyhot = "auto",
  xycold = "auto",
  DOY,
  sunelev,
  welev,
  zx = 10,
  u = 2,
  zomw = 2,
  zom = NULL,
  LAI = NULL,
  DEM = NULL,
  lapse = 0.0065,
  Rn24 = NULL,
  wmo = NULL,
  airport = NULL,
  Krs = 0.16,
  surface = "grass",
  model = "sebi",
  clip = NULL,
  folder = NULL
)

## Default S3 method:
sebi(
  albedo = NULL,
  Ts,
  Tmax,
  Tmin,
  RHmax = NULL,
  RHmin = NULL,
  latitude = NULL,
  NDVI,
  SAVI,
  iter.max = 7,
  xyhot = "auto",
  xycold = "auto",
  DOY,
  sunelev,
  welev,
  zx = 10,
  u = 2,
  zomw = 2,
  zom = NULL,
  LAI = NULL,
  DEM = NULL,
  lapse = 0.0065,
  Rn24 = NULL,
  wmo = NULL,
  airport = NULL,
  Krs = 0.16,
  surface = "grass",
  model = model,
  clip = NULL,
  folder = NULL
)

Arguments

Details

Function ETo will be automatically used to estimate parameters y (psychrometric constant), slope (slope of saturation vapour pressure curve) and vpd (Vapour pressure deficit).

Value

EF Evapotranspiration Fraction

Author(s)

George Owusu

References

Menenti, M., & Choudhury, B. 1993. Parameterization of land surface evaporation by means of location dependent potential evaporation and surface temperature range. Paper presented at the Proceedings of IAHS conference on Land Surface Processes.

See Also

sebal , ETo, sebs

Examples

## Not run: 
#using folder parameter 
folder=system.file("extdata","stack",package="sebkc")
sebiauto=sebi(folder=folder,welev=317,Tmax=31,Tmin=28)

#manual input of parameters
#data
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"))
#model
modsebi=sebi(albedo=albedo,Ts=Ts,Tmax=31,Tmin=28,RHmax=84,
RHmin=63,NDVI=NDVI,SAVI=NULL,iter.max=7,xyhot="full",
xycold="full",DOY=37,sunelev=50.71154048,welev=317.1,zx=10,
u=2,zomw=2,zom=NULL,LAI=LAI,Rn24=NULL,model="sebi")


#Another interactive example
#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="auto",cluster=2)
#determine the cold. Digitize polygon on the Ts map
modcold=coldTs(data,welev=275,extent="auto",cluster=2)
xyhot=modhot$xyhot
xycold=modcold$xycold
#use object of \code{\link{coldTs}} or \code{\link{hotTs}} in different ways
modsebi=sebi(data,xyhot=modhot,xycold=xycold,welev=welev,
Tmax=Tmax,Tmin=Tmin,RHmax=RHmax,RHmin=RHmin)

## End(Not run)

Internal function to Catch errors and warnings

Description

Internal function to Catch errors and warnings

Usage

sebkc.tryCatch(expr)

Arguments

Value

value a list with 'value' and 'warning', where value may be an error caught

Examples

thiserror=str( sebkc.tryCatch( log( "a" ) ) )$value

Stacking, gap-filling, cloud removal of Landsat image ET Model

Description

This function prepares landsat 5, 7, and 8 data for landsat578. First, it reads meta data and returns many input parameters for landsat578 and sebal. Second, it stacks 7 landsat layers: it stacks bands 1-7 for landsat 5 and 7; and bands 2-7,10 for landsat 8.. Third, it removes landsat clouds based on Martinuzzi et al(2009) method. Fourth, it gap-fills landsat 7 strips due to SLC-off as well as gap-filling NAs in other landsat sensors. It can also process any selected bands.

Usage

sebkcstack(
  folder,
  meta_pattern = "MTL.txt",
  remove_cloud = FALSE,
  gap_fill = FALSE,
  gap_loop = 8,
  clip = NULL
)

## Default S3 method:
sebkcstack(
  folder,
  meta_pattern = "MTL.txt",
  remove_cloud = FALSE,
  gap_fill = FALSE,
  gap_loop = 8,
  clip = NULL
)

Arguments

Details

The function can simulate any of the activities described in the description section. It takes several minutes to replace strips (1 hour) and remove clouds (30 minutes) on one landsat scene. But stacking of bands and reading of metadata take few seconds. One can save the stacked and processed data in the same "folder" with writesebkc. The output file name is "sebkcstack.tif".

Value

data:

raster stacked data. If crop is not NULL then crop Extent object is returned

brescale:

brescale of the stacked data

grescale:

grescale of the stacked data

sunelev:

sun elevation for stacked data

date:

Acquired date of the raster stacked data

viewangle:

view angle of landsat 8

azimuth:

azimuth angle of the image

sensor:

Type of satellite

K1:

K1 of landsat 8 thermal band 10

K2:

K2 of landsat 8 thermal band 10

cc:

cloud cover percentage

quality:

Image quality

row:

WRS_ROW of the image

path:

WRS_PATH of the image

Author(s)

George Owusu

References

Martinuzzi, S., Gould, W.A., Ramos Gonzales, O.M. 2007. Creating Cloud-Free Landsat ETM+ Data Sets in Tropical Landscapes: Cloud and Cloud-Shadow Removal. USDA Forest Service General Technical Report IITF-GTR-32.

USGS Using Landsat 7 Data http://landsat.usgs.gov/ERDAS_Approach.php

Examples

## Not run: 
folder=system.file("extdata","stack",package="sebkc")
stack=sebkcstack(folder=folder)
#returns data and grescale
stack$data
stack$grescale

## End(Not run)

SEBS Estimation of Evapotranspiration Fraction

Description

sebs computes Su (2002) Surface Energy Balance System (SEBS) for estimation of turbulent heat fluxes

Usage

sebs(
  albedo,
  Ts,
  Tmax,
  Tmin,
  RHmax = NULL,
  RHmin = NULL,
  n = NULL,
  NDVI,
  SAVI,
  iter.max = 7,
  xyhot,
  xycold,
  DOY,
  sunelev,
  welev,
  zx,
  u,
  zomw,
  zom = NULL,
  LAI = NULL,
  DEM = NULL,
  lapse = 0.0065,
  Rn24 = NULL,
  wmo = NULL,
  airport = NULL,
  Krs = 0.16,
  surface = "grass",
  latitude = NULL,
  model = "SEBS",
  clip = NULL,
  folder = NULL
)

## Default S3 method:
sebs(
  albedo = NULL,
  Ts,
  Tmax,
  Tmin,
  RHmax = NULL,
  RHmin = NULL,
  n = NULL,
  NDVI,
  SAVI,
  iter.max = 7,
  xyhot = "auto",
  xycold = "auto",
  DOY = NULL,
  sunelev,
  welev,
  zx = 10,
  u = 2,
  zomw = 2,
  zom = NULL,
  LAI = NULL,
  DEM = NULL,
  lapse = 0.0065,
  Rn24 = NULL,
  wmo = NULL,
  airport = NULL,
  Krs = 0.16,
  surface = "grass",
  latitude = NULL,
  model = "SEBS",
  clip = NULL,
  folder = NULL
)

Arguments

Details

Function ETo will be automatically used to estimate parameters y (psychrometric constant), slope (slope of saturation vapour pressure curve) and vpd (Vapour pressure deficit). For information on data input procedure see 'details' at sebal

Value

EF Evapotranspiration Fraction

Author(s)

George Owusu

References

Su, Z. 2002. The Surface Energy Balance System (SEBS) for estimation of turbulent heat fluxes. Hydrol. Earth Syst. Sci., 6(1): 85-100.

See Also

sebal , ETo, sebi

Examples

## Not run: 
#define Input data
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"))
Tmax=31
Tmin=28
RHmax=84
RHmin=63
#perform simulation
modsebs=sebs(albedo=albedo,Ts=Ts,Tmax=Tmax,Tmin=Tmin,RHmax=RHmax,
RHmin=RHmin,NDVI=NDVI,SAVI=NULL,iter.max=7,xyhot="full",
xycold="full",DOY=37,sunelev=50.71154048,welev=317.1,zx=10,
u=2,zomw=2,zom=NULL,LAI=LAI,Rn24=NULL,model="SEBS")


#use original landsat 7 data by specifying folder path
folder=system.file("extdata","stack",package="sebkc")
modauto=sebs(folder = folder,welev = 380,Tmax=Tmax,
Tmin=Tmin,RHmax=RHmax,RHmin=RHmin)
#plot ET fraction
plot(modauto$EF)

#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="auto",cluster=2)
#determine the cold. Digitize polygon on the Ts map
modcold=coldTs(data,welev=275,extent="auto",cluster=2)
xyhot=modhot$xyhot
xycold=modcold$xycold
#use object of \code{\link{coldTs}} or \code{\link{hotTs}} in different ways
modsebs=sebs(data,xyhot=modhot,xycold=xycold,welev=welev,
Tmax=Tmax,Tmin=Tmin,RHmax=RHmax,RHmin=RHmin)
 
## End(Not run)

Simplified Surface Energy Balance

Description

Senay etal (2011) Simplified Surface Energy Balance (SSEB) approach for estimating landscape ET

Usage

sseb(
  Ts,
  welev = NULL,
  TH = "auto",
  TC = "auto",
  sunelev = NULL,
  NDVI = NULL,
  DEM = NULL,
  albedo = NULL,
  NDVImax = 0.7,
  cc = 0.65,
  KL = 0.0065,
  ETo = NULL,
  wmo = NULL,
  airport = NULL,
  Krs = 0.16,
  surface = "grass",
  x = 1,
  Tmax = NULL,
  Tmin = NULL,
  zx = NULL,
  u = 2,
  DOY = NULL,
  latitude = NULL,
  n = NULL,
  RHmax = NULL,
  RHmin = NULL,
  clip = NULL,
  folder = NULL
)

## Default S3 method:
sseb(
  Ts = NULL,
  welev = NULL,
  TH = "auto",
  TC = "auto",
  sunelev = NULL,
  NDVI = NULL,
  DEM = NULL,
  albedo = NULL,
  NDVImax = 0.7,
  cc = 0.65,
  KL = 0.0065,
  ETo = NULL,
  wmo = NULL,
  airport = NULL,
  Krs = 0.16,
  surface = "grass",
  x = 1.2,
  Tmax = NULL,
  Tmin = NULL,
  zx = NULL,
  u = 2,
  DOY = NULL,
  latitude = NULL,
  n = NULL,
  RHmax = NULL,
  RHmin = NULL,
  clip = NULL,
  folder = NULL
)

Arguments

Value

ET24:

Actual 24 hour ET

EF:

ET fraction

ETo:

24 hour reference ETo

Author(s)

George Owusu

References

Senay, G. B., Budde, M. E., & Verdin, J. P. 2011. Enhancing the Simplified Surface Energy Balance (SSEB) approach for estimating landscape ET: Validation with the METRIC model. Agricultural Water Management, 98(4): 606-618.

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"))
#minimal data
modsseb=sseb(Ts=Ts,TH="full",TC="full",sunelev=50,NDVI=NDVI,DEM=NULL,
albedo=albedo,NDVImax=0.7,cc=0.65,KL=0.0065,
ETo="auto",x=1.2,Tmax=31,Tmin=28,zx=10,
u=2,DOY=37,latitude=5.6,n=NULL,RHmax=NULL,RHmin=NULL)

#uaing the folder parameters
folder=system.file("extdata","stack",package="sebkc")
ssebauto=sseb(folder=folder,welev=317,Tmax=31,Tmin=28, latitude=5.6)

#Another interactive example
#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="auto",cluster=2)
#determine the cold. Digitize polygon on the Ts map
modcold=coldTs(data,welev=275,extent="auto",cluster=2)
xyhot=modhot$Tshot
xycold=modcold$Tscold
#use object of \code{\link{coldTs}} or \code{\link{hotTs}} in different ways
modsseb=sseb(Ts=Ts,TH=xyhot,TC=modcold,sunelev=50,NDVI=NDVI,DEM=NULL,
albedo=albedo,NDVImax=0.7,cc=0.65,KL=0.0065,
ETo="auto",x=1.2,Tmax=31,Tmin=28,zx=10,
u=2,DOY=37,latitude=5.6,n=NULL,
RHmax=NULL,RHmin=NULL)

## End(Not run)

S-SEBI Evapotranspiration Fraction

Description

ssebi computes Roerink et. al (2000) S-SEBI: A simple remote sensing algorithm to estimate the surface energy balance.

Usage

ssebi(
  Ts,
  albedo = NULL,
  threshold = 0.2,
  plot = TRUE,
  Rn24 = NULL,
  folder = NULL,
  welev = NULL
)

## Default S3 method:
ssebi(
  Ts = NULL,
  albedo = NULL,
  threshold = 0.2,
  plot = TRUE,
  Rn24 = NULL,
  folder = NULL,
  welev = NULL
)

Arguments

Value

EF, amax,bmax,amin,bmin

Author(s)

George Owusu

References

Roerink, G. J., Su, Z., & Menenti, M. 2000. S-SEBI: A simple remote sensing algorithm to estimate the surface energy balance. Physics and Chemistry of the Earth, Part B: Hydrology, Oceans and Atmosphere , 25(2): 147-157.

Examples

## Not run: 
#Manual data specification
albedo=raster(system.file("extdata","albedo.grd",package="sebkc"))
Ts=raster(system.file("extdata","Ts.grd",package="sebkc"))
mod=ssebi(Ts=Ts,albedo=albedo,threshold=0,plot=TRUE)

#Using landsat folder
folder=system.file("extdata","stack",package="sebkc")
ssebiauto=ssebi(folder=folder,welev=278)

## End(Not run)

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

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="auto",cluster=2)
#determine the cold. Digitize polygon on the Ts map
modcold=coldTs(data,welev=275,extent="auto",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)

Water Deficit Index (WDI)

Description

Moran et. al (1994) Estimation of crop water deficit using surface-air temperature and spectral vegetation index.

Usage

wdi(
  Ts = NULL,
  Ta,
  NDVI,
  albedo = NULL,
  Tsmax = "auto",
  Tsmin = "auto",
  sunelev = NULL,
  folder = NULL,
  welev = NULL,
  Rn24 = NULL
)

wdi(
  Ts = NULL,
  Ta,
  NDVI,
  albedo = NULL,
  Tsmax = "auto",
  Tsmin = "auto",
  sunelev = NULL,
  folder = NULL,
  welev = NULL,
  Rn24 = NULL
)

Arguments

Details

Note that surface temperature (Tsmin and Tsmax) are different from air temperature (Tmin and Tmax) in ETo, sseb, sebs

You can get Ta from sebal Ta

Value

EF:

standardised wdi excluding negatives and one plus

index:

raw wdi including negatives and more than one

Author(s)

George Owusu

References

Moran, M. S., Clarke, T. R., Inoue, Y., & Vidal, A. 1994. Estimating crop water deficit using the relation between surface-air temperature and spectral vegetation index. Remote Sensing of Environment, 49(3): 246-263.

See Also

sseb and sebi

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"))
modWDI=wdi(Ts=Ts,Ta=299,NDVI=NDVI,albedo=NULL,Tsmax="auto",
Tsmin="auto",sunelev=50)
plot(modWDI$EF)

#Using folder parameter
folder=system.file("extdata","stack",package="sebkc")
wdiauto=wdi(folder=folder,welev=317,Tsmax=31,Tsmin=28,Ta=290)

#' #Another interactive example
#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="auto",cluster=2)
#determine the cold. Digitize polygon on the Ts map
modcold=coldTs(data,welev=275,extent="auto",cluster=2)
Tsmax=modhot$Tshot
Tsmin=modcold$Tscold
#use object of \code{\link{coldTs}} or \code{\link{hotTs}} in different ways
modWDI2=wdi(data,Ta=299,Tsmax=Tsmax,Tsmin=modcold,sunelev=50)
plot(modWDI2$EF)

## End(Not run)

Retrieval of Weather Data from WMO and NASA Surface meteorology and Solar Energy

Description

The function retrieves weather data from the World Meteorological Organization and NASA SSE. The retrieval is based on airports IATA or ICAO codes or WMO weather station codes. The actual WMO data is from Weather Underground . Internet connection is therefore needed before data can be accessed

Usage

weather(
  data = NULL,
  wmo = NULL,
  airport = NULL,
  date = "YYYY-m-d",
  time = NULL,
  latitude = NULL,
  longitude = NULL,
  NASA.SSE = list(from = "YYYY-m-d", to = "YYYY-m-d"),
  folder = NULL
)

Arguments

Value

WMO and NASE.SSE. The WMO contains hour and day data. The NASE.SSE contains meta and data. See example below.

Author(s)

George Owusu

References

• Worldwide Station List.https://www.wetterzentrale.de/ [Accessed on 2016-5-6]

• List of airports by IATA code. https://en.wikipedia.org/wiki/List_of_airports_by_IATA_code [Accessed on 2016-5-6]

• Weather Underground https://www.wunderground.com/ [Accessed on 2016-5-6]

Examples

## Not run: 
#all hours data at Lincoln Airport, Nabraska
lincoln=weather(airport = "KLNK", date="2011-08-04")
lincoln$WMO$hour
lincoln$WMO$day
#current hour data at London
london.now=weather(airport = "LON", date=format(Sys.time(), "%Y-%m-%d"),
time=format(Sys.time(), "%H:%M"))
london.now$WMO$hour
#NASA SSE data for Kumasi, Ghana
kumasi=weather(longitude = -1.5,latitude = 6.7,NASA.SSE = list(from="2001-04-01",
to="2002-04-01"))

#Hour, daily and SSE data at JFK 
JFK=weather(airport = "JFK", date="2001-02-02",time="17:20", latitude=40.6413, 
longitude=73.7781)
JFK$WMO$hour
JFK$WMO$day
JFK$NASA.SEE$data
JFK$NASA.SEE$meta

#' #Hour, daily and SSE slice data at Kotoka International Airport, ACCRA 
accra=weather(airport = "ACC", date="2001-02-02",time="17:20", latitude=40.6413, longitude=73.7781, 
NASA.SSE=list(from="2001-01-01",to="2004-01-01"))
accra$WMO$hour
accra$WMO$hour
accra$NASA.SEE$data
accra$NASA.SEE$meta

#Using wmo code
kumasi.now=weather(wmo="65442", date=format(Sys.time(), "%Y-%m-%d"),time=format(Sys.time(), 
"%H:%M"))
kumasi.now$WMO$hour

#retrieve daily data for a specific period i.e 1 month
lincoln=weather(airport = "KLNK", date=c("2011-08-04","2011-09-04"))
#kumasi and accra data for a time slice. May not work for different 
#geographical regions; write the results to a file
accrakumasi=weather(airport ="ACC", date=c("2011-08-04","2011-10-04"),
folder="C:/Users/george/Documents/")

## End(Not run)

Writing sebkc object to a folder

Description

This function writes spatial output of the following functions sebal, sebi,ssebi, wdi,sebs,sseb, tseb and landsat578

Usage

writesebkc(object, folder = NULL, xy = NULL, overwrite = TRUE)

Arguments

Value

Writes output to a folder

Examples

## Not run: 
folder=system.file("extdata","stack",package="sebkc")
stack=landsat578(data=folder, welev=362)
writesebkc(stack)
sebaloutput=sebal(folder = folder,welev = 317)
writesebkc(sebaloutput)


## End(Not run)

Reprojecting and extracting Raster values

Description

The function reprojects longitudes and latitudes into map coordinates and extract raster values

Usage

xyextract(map, longitude, latitude)

Arguments

Value

data

Author(s)

George Owusu

Examples

## Not run: 
lonlat=read.table(system.file("extdata","sys","xyvalues.txt",package="sebkc"),header=TRUE)
longitude=lonlat$longitude
latitude=lonlat$latitude
folder=system.file("extdata","stack",package="sebkc")
data=landsat578(data=folder, welev=362)
Ts=data$Ts
Tsdata=xyextract(Ts,longitude,latitude)


## End(Not run)