Package 'sebkc'

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

Tday	numeric. Average Day time temperature in degrees Celsius for the growing period. A Remote sensing surface Temperature can also be used
T24	numeric. Average 24-hour (including night hours) temperature in degrees Celsius for the growing period.
Rs	A shortwave radiation data [ca^-2 day^-1]. RasterLayer data can be supplied.
latitude	geographical coordinates in decimal degrees. It should be negative for southern hemisphere
longitude	The longitude of the measurement site i.e. geographical coordinates in decimal degrees for the weather station. It should be negative for West and positive for East.
HI	numeric. Harvest Index. It can be list or spatial
plantdate	character. A date of planting in the form "YYYY-mm-dd"
harvestdate	character. A date of harvest in the form "YYYY-mm-dd"
c30	numeric. A maintenance respiration coefficient (c30) at 30 oC. it takes either a value of 0.0108 for non-legume crop or 0.0283 for legumes
adaptability	numeric. The climate adaptability group number. It takes values of 1,2,3,or 4. it depends on the temperature profile of study area. See references (eg. FAO (1981)) for more explanation.
LAI	Leaf Area Index, dimensionless
crop	character. A crop type. For example "maize"

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

ETo	Numeric. Reference evapotranspiration [mm]. This can be estimated with ETo
P	numeric. Precipitation [mm]
RHmin	numeric. Minimum Relative Humidity [per cent]
soil	character. Name of soil texture as written in FAO56 TABLE 19
crop	character. Name of the crop as written on FAO56 Table 11, TABLE 12 and TABLE 17
I	numeric. Irrigation [mm]
CNII	numeric. The initial Curve Number of the study site. This can be calibrated with cal.kc.
u2	numeric. Wind speed of the nearby weather station [m/s] measured at 2m
FC	numeric. soil water content at field capacity [-][m3 m-3]. This can be calibrated with cal.kc.
WP	numeric. soil water content at wilting point [-][m3 m-3]. This can be calibrated with cal.kc.
h	numeric. Maximum crop height [m]. See details below. This can be calibrated with cal.kc
Zr	numeric. Rooting depth. See details below. This can be calibrated with cal.kc
Ze	numeric. The depth of the surface soil layer that is subject to drying by way of evaporation [0.10-0.15 m]. This can be calibrated with cal.kc
kc	list. The tabulated Kc values for single crop c(kcini,kcmid,kcend) or dual crop model c(kcbini,kcbmid,kcbend). To use with sebal or any of the Surface Energy Balance models, it should be linked as list(model1,model2,model3) or c(model1$EF,model2$EF,model3$EF) but not a mixture of both. This can be calibrated with cal.kc
p	The depletion coefficient. This can be calibrated with cal.kc
lengths	list. The lengths of the crop growth stages. It can take dates in the form of c(start of planting ('YYYY-mm-dd), start of rapid growth, start of mid season, start of maturity, harvesting) or in the form number of days: c(initial days, dev days, mid days, late days). See FAO56 TABLE 11. If performing time series simulation, different dates or one season dates can be set; for instance two seasons simulation can be of the form c(initial days, dev days, mid days, late days, initial days, dev days, mid days, late days) or just c(initial days, dev days, mid days, late days)
fw	list or numeric or character. This can take different fw values throughout the growing season. It can also take one numeric value and this will be transform depending on the type of wetting. If the wetting event is only Precipitation, it is set to 1; if it is only Irrigation it is set to the provided fw value. If it is both Precipitation and Irrigation weights are assigned for calculation of fw. fw can also take a character such as "Drip", "Sprinkler", "Basin", "Border", "alternated furrows", "Trickle", "narrow bed" and "wide bed". In this case fw values from FAO56 are assigned.
fc	numeric. fraction of ground covered by tree canopy
kctype	character. The type of model either "single" or "double". In spatial modelling involving Evaporation Fraction, kctype can be set to "METRIC" or "SSEB" so that Actual Evapotranspiration will be estimated as ETo*EF else it will be estimated with Rn as ETa=(EF*Rn)/2.45. Note that Rn is in MJm-2 day-1.
region	character. The region of the crop as on FAO56 TABLE 11
regexpr	character. Extra term that can be used to search crop names on FAO56 Tables
initgw	numeric. Initial ground water level [m]
LAI	numeric. Leaf Area Index
rc	numeric Runoff coefficient [0-1]. This must be used if there is no information on Curve Number. This can be calibrated with cal.kc
Kb	ground water (base flow) recession parameter. The default is zero where no baseflow occurs
Rn	Net daily radiation [w/m2]. Estimated this with ETo
hmodel	character. The type of crop growth model during the simulation. It takes "gaussian","linear", "exponential"
Zrmodel	character. The type of root growth model during the simulation. It takes "gaussian","linear", "exponential"
xydata	list c(longitude, latitude). The observation points of the model. This should be set only if the model is spatial.
REW	numeric Readily Evaporable Water. It can be calibrated with cal.kc
a1	soil water storage to maximum root depth (Zr) at field capacity
b1	Liu et al. (2006) parameter for computing capillary rise [-0.17]. It can calibrated with cal.kc
a2	Liu et al. (2006) parameter for computing capillary rise [240]. It can be calibrated with cal.kc
b2	Liu et al. (2006) parameter for computing capillary rise [-0.27]. It can be calibrated with cal.kc
a3	Liu et al. (2006) parameter for computing capillary rise [-1.3]. It can be calibrated with cal.kc
b3	Liu et al. (2006) parameter for computing capillary rise [6.6]. It can be calibrated with cal.kc
a4	Liu et al. (2006) parameter for computing capillary rise [4.6]. It can be calibrated with cal.kc
b4	Liu et al. (2006) parameter for computing capillary rise [-0.65]. It can be calibrated with cal.kc
ETa	Measured Actual Evapotranspiration [mm]
report	list. In case of spatial modelling, which days' map should be reported. The default is "stages" where the periods of stages maps are reported as values. In an ascending order it can take values such as 1:12 days or c(1,3,6:12) days.
time	character. Time of the occurrence Irrigation or Precipitation.
HWR	numeric. Default is 1. Height to width ratio of individual plants or groups of plants when viewed from the east or from the west [-].
latitude	geographical coordinates in decimal degrees. It should be negative for southern hemisphere
density	character. density factor for the trees, vines or shrubs on how they are arranged. It can take "full", "row" or "random" or "round".
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
slope	numeric. slope of saturation vapour pressure curve at air temperature T [kPa oC-1].
y	numeric. psychometric constant [kPa oC-1]
r1	leaf resistance FOR STOMATAL CONTROL [s/m]
nongrowing	character. The type of surface during non-growing season. It can take "weeds", "bare", or "mulch"
theta	numeric. The measured soil water content [m3/m3] in the profile. Leave those days without measured data blank.
profile	numeric or character. The profile depth to perform water balance. It takes "variable" when updating the depth with root growth (default). Or taking "fixed" when the maximum root zone is used. It can also take numeric value and this can be longer than maximum root length.
Kcmin	Minimum Kc value
Ky	list or array. yield response factor [-] in the form of c(initial Ky,crop dev't Ky, mid season Ky, late season Ky, total season Ky). For instance Maize Ky is c(0.4,0.4,1.3,0.5,1.25). Ky values for other crops are tabulated in FAO Irrigation and Drainage Paper 33 (Doorenbos and Kassam, 1979) and Paper 56 (Allen et al., 1998). If only one Ky value is supplied it will be repeated for all the seasons.

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 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

Tmin	Numeric. Minimum surface Temperature in map units (eg. Kilvin) that you want computation from. You can keep changing untill you get your desired values
Tmax	Numeric. Maximum surface Temperature in map units (eg. Kilvin) that you want computation from. You can keep changing untill you get your desired values
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.
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
welev	Weather station elevation in meters

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)

---

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

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.
albedo	A RasterLayer data that has albedo values. You can also provide file path or location on your computer
sunangle	sun angle (degree) above the horizon.
DEM	A digital elevation model[m]
cluster	numeric, a value indicating the number of clusters.
extent	A subset of the raster file than include c(xmin, xmax, ymin, ymax). You can set it to "interactive" where you can digitise polygon. You can set to "auto" where it will be automatically generated. You can set it to "full" where all the image dimension will be used. See layout parameter below for more details. see extent
upper	The quantile probability at which NDVI should be included in the selection project. It ranges from 0-1. The default is 0.95
lower	The quantile probability at which Ts should be included in the slection project.The default is 0.2
plot	logical:TRUE or FALSE. To indicate if histograms or triangle should be plotted
layout	character. It takes "portrait" or "landscape". If extent is set to "auto" the middle part of the image will be used so that NAs can be avoided. You indicate whether the automatic slection should indicate west-east (landscape) or north-south direction (portraite)
draw	interactive. The method of digitisation if extent is set "auto". It takes "poly" or "rect"
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
welev	Weather station elevation in meters
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)
iter.max	The number iterations that kmeans cluster should use.
seed	optional single number. Passed to set.seed before the k-means clustering so the pixel selection is reproducible. The default NULL leaves the random-number stream untouched.

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

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.
albedo	A RasterLayer data that has albedo values. You can also provide file path or location on your computer
sunangle	sun angle (degree) above the horizon.
DEM	A digital elevation model[m]
cluster	numeric, a value indicating the number of clusters.
extent	A subset of the raster file than include c(xmin, xmax; serastercond row: ymin, ymax). You can set it to "interactive" where you can digitise polygon. You can set to "auto" where it will be automatically generated. See layout parameter below for more details. see extent
upper	The quantile probability at which NDVI should be included in the selection project. It ranges from 0-1. The default is 0.95
lower	The quantile probability at which Ts should be included in the selection project.The default is 0.2
NDVI.probs	The probability of selecting quantile NDVI that is above the threshold. The value of 0.6 means only NDVI values that is above 60 the analyses.
plot	logical:TRUE or FALSE. To indicate if histograms or triangle should be plotted
layout	character. It takes "portrait" or "landscape". If extent is set to "auto" the middle part of the image will be used so that NAs can be avoided. You indicate whether the automatic selection should indicate west-east (landscape) or north-south direction (portrait)
draw	interactive. The method of digitisation if extent is set "auto". It takes "poly" or "rect"
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
welev	Weather station elevation in meters

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)

---

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

x	object to be extracted
y	extent object or raster object or polygon from which an Extent object can be extracted. A polygon or raster will be reprojected to conform x

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

latitude	geographical coordinates in decimal degrees. It should be negative for southern hemisphere
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
model	character Type of model:"CBM" (Schoolfield, 1982),"BGC" (Running & Coughlan, 1988), "CERES" (Ritchie, 1991) or "FAO56"
Tmax	Numeric. Maximum air Temperature in degree Celsius
Tmin	Numeric. Minimum air Temperature in degree Celsius
p	numeric. CMB parameter

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)

---

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

Tmax	Numeric. Maximum air Temperature in degree Celsius
Tmin	Numeric. Minimum air Temperature in degree Celsius
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
latitude	geographical coordinates in decimal degrees. It should be negative for southern hemisphere
longitude	The longitude of the measurement site i.e. geographical coordinates in decimal degrees for the weather station. It should be negative for West and positive for East.
z	The height above surface where wind speed was measured, in metres.
uz	The wind speed in m/s
altitude	The elevation of the weather station, in metres
RHmax	Maximum Relative Humidity in percent
RHmin	Minimum Relative Humidity in percent
n	Sunshine hours
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
albedo	albedo values. For grass it is 0.23[dimensionless]
as	Regression constant, intercept, expressing the fraction of extra-terrestrial radiation reaching the earth on overcast days (n = 0); either a calibrated value or as = 0.25 can be used.
bs	regression constant, slope, expressing the fraction of extraterrestrial radiation reaching the earth on overcast days (n = 0); either a calibrated value or bs = 0.5 can be used.
Rs	shortwave incoming radiation [MJ/m2/day]
Rn	Net daily radiation in MJ/m2/day. If it is NULL, it will be computed
G	Daily Soil heat flux in W/m2
EF	Spatial. ET fraction [-] from sebal or tseb sebs. It is advisable to link it to any of the two functions.
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]
scale	character The type of scale models to upscale instantaneous ET to daily ET. it takes "Evaporative Fraction" (Bastiaanssen et al 1998); "Sine Function" (Jackson et al 1983); "Solar Radiation" (French et al 2005). Remember the blank spaces between the parameters. You can also respectively use "evaporative" or "solar" or "sine".
Kc	Crop coefficient
surface	character. It is either "grass" or "alfalfa"

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

Tmean	Numeric. Hourly Mean Air Temperature in degree Celsius
RH	Hourly Mean Relative Humidity in percent
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
Lz	longitude of the centre of the local time zone [degrees west of Greenwich]. MUST BE POSITIVE Lz = 0 degrees for Greenwich.
t1	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
Lm	longitude of the measurement site [degrees west of Greenwich]. MUST BE POSITIVE
latitude	geographical coordinates in decimal degrees. It should be negative for southern hemisphere
longitude	The longitude of the measurement site i.e. geographical coordinates in decimal degrees for the weather station. It should be negative for West and positive for East.
n	Sunshine hours
uz	The wind speed in m/s
Rs	solar radiation [MJ/m2/hr]
Rn	net radiation [MJ/m2/hr]
G	Daily Soil heat flux in W/m2
altitude	The elevation of the weather station, in metres
z	The height above surface where wind speed was measured, in metres.
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
Tmax	Numeric. Maximum air Temperature in degree Celsius
Tmin	Numeric. Minimum air Temperature in degree Celsius
albedo	A RasterLayer data that has albedo values. You can also provide file path or location on your computer
as	Regression constant, intercept, expressing the fraction of extra-terrestrial radiation reaching the earth on overcast days (n = 0); either a calibrated value or as = 0.25 can be used.
bs	regression constant, slope, expressing the fraction of extraterrestrial radiation reaching the earth on overcast days (n = 0); either a calibrated value or bs = 0.5 can be used.
ETins	model. It takes ET model from sebal or tseb and retrieve instantaneous ET estimates. This helps to estimate ET24, T24, E24, EF, Kc, EFs, EFc
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]
map	spatial raster data for inverse distance interpolation
surface	character. It is either "grass" or "alfalfa"
period	character. Time of day for the calculation: "daytime" or "nighttime". The default is "daytime".

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

---

Goodness of fit Test

Description

Goodness of fit Test

Usage

fitting(obs, est, P = 1)

Arguments

obs	Observed data
est	Estimated data
P	Number of parameters

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

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.
albedo	A RasterLayer data that has albedo values. You can also provide file path or location on your computer
DEM	A digital elevation model[m]
cluster	numeric, a value indicating the number of clusters.
extent	A subset of the raster file than include c(xmin, xmax; serastercond row: ymin, ymax). You can set it to "interactive" where you can digitise polygon. You can set to "auto" where it will be automatically generated. See layout parameter below for more details. see extent
upper	The quantile probability at which Ts should be included in the selection project. It ranges from 0-1. The default is 0.8
lower	The quantile probability at which NDVI should be included in the selection project. The default is 0.2
plot	logical:TRUE or FALSE. To indicate if histograms or triangle should be plotted
layout	character. It takes "portrait" or "landscape". If extent is set to "auto" the middle part of the image will be used so that NAs can be avoided. You indicate whether the automatic selection should indicate west-east (landscape) or north-south direction (portrait)
draw	interactive. The method of digitisation if extent is set "auto". It takes "poly" or "rect"
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
welev	Weather station elevation in meters
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)
iter.max	maximum iterations of sensible heat calculation.
seed	optional single number. Passed to set.seed before the k-means clustering so the pixel selection is reproducible. The default NULL leaves the random-number stream untouched.

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

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.
albedo	A RasterLayer data that has albedo values. You can also provide file path or location on your computer
sunangle	sun angle (degree) above the horizon.
DEM	A digital elevation model[m]
cluster	numeric, a value indicating the number of clusters.
extent	A subset of the raster file than include c(xmin, xmax; serastercond row: ymin, ymax). You can set it to "interactive" where you can digitise polygon. You can set to "auto" where it will be automatically generated. See layout parameter below for more details. see extent
upper	The quantile probability at which Ts should be included in the selection project. It ranges from 0-1. The default is 0.8
lower	The quantile probability at which NDVI should be included in the slection project.The default is 0.2
NDVI.probs	The probability of selecting quantile NDVI that is above the threshold. The value of 0.6 means only NDVI values that is above 60 the analyses.
plot	logical:TRUE or FALSE. To indicate if histograms or triangle should be plotted
layout	character. It takes "portrait" or "landscape". If extent is set to "auto" the middle part of the image will be used so that NAs can be avoided. You indicate whether the automatic selection should indicate west-east (landscape) or north-south direction (portrait)
draw	interactive. The method of digitisation if extent is set "auto". It takes "poly" or "rect"
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
welev	Weather station elevation in meters

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

longitude	The longitude of the measurement site i.e. geographical coordinates in decimal degrees for the weather station. It should be negative for West and positive for East.
latitude	geographical coordinates in decimal degrees. It should be negative for southern hemisphere
var	numeric variable. The length should be same as the length of latitude and longitude
ext	RasterLayer object.
idp	Neighbourhood parameter

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

ETo	Numeric. Reference evapotranspiration [mm]. This can be estimated with ETo
P	numeric. Precipitation [mm]
RHmin	numeric. Minimum Relative Humidity [per cent]
soil	character. Name of soil texture as written in FAO56 TABLE 19
crop	character. Name of the crop as written on FAO56 Table 11, TABLE 12 and TABLE 17
I	numeric. Irrigation [mm]
CNII	numeric. The initial Curve Number of the study site. This can be calibrated with cal.kc.
u2	numeric. Wind speed of the nearby weather station [m/s] measured at 2m
FC	numeric. soil water content at field capacity [-][m3 m-3]. This can be calibrated with cal.kc.
WP	numeric. soil water content at wilting point [-][m3 m-3]. This can be calibrated with cal.kc.
h	numeric. Maximum crop height [m]. See details below. This can be calibrated with cal.kc
Zr	numeric. Rooting depth. See details below. This can be calibrated with cal.kc
Ze	numeric. The depth of the surface soil layer that is subject to drying by way of evaporation [0.10-0.15 m]. This can be calibrated with cal.kc
kc	list. The tabulated Kc values for single crop c(kcini,kcmid,kcend) or dual crop model c(kcbini,kcbmid,kcbend). To use with sebal or any of the Surface Energy Balance models, it should be linked as list(model1,model2,model3) or c(model1$EF,model2$EF,model3$EF) but not a mixture of both. This can be calibrated with cal.kc
p	The depletion coefficient. This can be calibrated with cal.kc
lengths	list. The lengths of the crop growth stages. It can take dates in the form of c(start of planting ('YYYY-mm-dd), start of rapid growth, start of mid season, start of maturity, harvesting) or in the form number of days: c(initial days, dev days, mid days, late days). See FAO56 TABLE 11. If performing time series simulation, different dates or one season dates can be set; for instance two seasons simulation can be of the form c(initial days, dev days, mid days, late days, initial days, dev days, mid days, late days) or just c(initial days, dev days, mid days, late days)
fw	list or numeric or character. This can take different fw values throughout the growing season. It can also take one numeric value and this will be transform depending on the type of wetting. If the wetting event is only Precipitation, it is set to 1; if it is only Irrigation it is set to the provided fw value. If it is both Precipitation and Irrigation weights are assigned for calculation of fw. fw can also take a character such as "Drip", "Sprinkler", "Basin", "Border", "alternated furrows", "Trickle", "narrow bed" and "wide bed". In this case fw values from FAO56 are assigned.
fc	numeric. fraction of ground covered by tree canopy
kctype	character. The type of model either "single" or "double". In spatial modelling involving Evaporation Fraction, kctype can be set to "METRIC" or "SSEB" so that Actual Evapotranspiration will be estimated as ETo*EF else it will be estimated with Rn as ETa=(EF*Rn)/2.45. Note that Rn is in MJm-2 day-1.
region	character. The region of the crop as on FAO56 TABLE 11
regexpr	character. Extra term that can be used to search crop names on FAO56 Tables
initgw	numeric. Initial ground water level [m]
LAI	numeric. Leaf Area Index
rc	numeric Runoff coefficient [0-1]. This must be used if there is no information on Curve Number. This can be calibrated with cal.kc
Kb	ground water (base flow) recession parameter. The default is zero where no baseflow occurs
Rn	Net daily radiation [w/m2]. Estimated this with ETo
hmodel	character. The type of crop growth model during the simulation. It takes "gaussian","linear", "exponential"
Zrmodel	character. The type of root growth model during the simulation. It takes "gaussian","linear", "exponential"
xydata	list c(longitude, latitude). The observation points of the model. This should be set only if the model is spatial.
REW	numeric Readily Evaporable Water. It can be calibrated with cal.kc
a1	soil water storage to maximum root depth (Zr) at field capacity
b1	Liu et al. (2006) parameter for computing capillary rise [-0.17]. It can calibrated with cal.kc
a2	Liu et al. (2006) parameter for computing capillary rise [240]. It can be calibrated with cal.kc
b2	Liu et al. (2006) parameter for computing capillary rise [-0.27]. It can be calibrated with cal.kc
a3	Liu et al. (2006) parameter for computing capillary rise [-1.3]. It can be calibrated with cal.kc
b3	Liu et al. (2006) parameter for computing capillary rise [6.6]. It can be calibrated with cal.kc
a4	Liu et al. (2006) parameter for computing capillary rise [4.6]. It can be calibrated with cal.kc
b4	Liu et al. (2006) parameter for computing capillary rise [-0.65]. It can be calibrated with cal.kc
report	list. In case of spatial modelling, which days' map should be reported. The default is "stages" where the periods of stages maps are reported as values. In an ascending order it can take values such as 1:12 days or c(1,3,6:12) days.
time	character. Time of the occurrence Irrigation or Precipitation.
HWR	numeric. Default is 1. Height to width ratio of individual plants or groups of plants when viewed from the east or from the west [-].
latitude	geographical coordinates in decimal degrees. It should be negative for southern hemisphere
density	character. density factor for the trees, vines or shrubs on how they are arranged. It can take "full", "row" or "random" or "round".
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
slope	numeric. slope of saturation vapour pressure curve at air temperature T [kPa oC-1].
y	numeric. psychometric constant [kPa oC-1]
r1	leaf resistance FOR STOMATAL CONTROL [s/m]
nongrowing	character. The type of surface during non-growing season. It can take "weeds", "bare", or "mulch"
theta	numeric. The measured soil water content [m3/m3] in the profile. Leave those days without measured data blank.
profile	numeric or character. The profile depth to perform water balance. It takes "variable" when updating the depth with root growth (default). Or taking "fixed" when the maximum root zone is used. It can also take numeric value and this can be longer than maximum root length.
Kcmin	Minimum Kc value
Ky	list or array. yield response factor [-] in the form of c(initial Ky,crop dev't Ky, mid season Ky, late season Ky, total season Ky). For instance Maize Ky is c(0.4,0.4,1.3,0.5,1.25). Ky values for other crops are tabulated in FAO Irrigation and Drainage Paper 33 (Doorenbos and Kassam, 1979) and Paper 56 (Allen et al., 1998). If only one Ky value is supplied it will be repeated for all the seasons.
x	a return object of the function.
output	The kind of plot that should be displayed. It takes 'AW_measured' for plotting Available and modelled soil moisture; 'TS_measured' for Total Available and modelled soil moisture; 'Kc' for displaying Kcb and Kc 'Evaporation' for displaying of Evaporation and Transpiration 'balance' for displaying Irrigation, precipitation, ETa and critical water deficit
main	Title of the plot
cex	the relative font size of the texts. Vlaues: "horiz"; logical: if TRUE, set the legend horizontally rather than vertically
legend	logical. TRUE or FALSE. Whether a legend should be displayed or not.
pos	the position of legend. It takes one of the following 'top', 'topright', 'topleft', 'bottom', 'bottomright' 'bottomleft', 'center'. This must be set so that the legend should not obscure the main plot.
horiz	logical; if TRUE the legend is drawn horizontally rather than vertically.
...	further arguments passed to the plotting functions.

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

data	An original folder containing Landsat data or Raster brick or raster stack. Arrange the bands in ascending order: c(band1,band2,band3,band4,band5,band6,band7) for Landsat 5 and 7. And a similar order for Landsat 8. On folder see sebkcstack See also details.
welev	Weather station elevation in meters
sensor	Numeric. Type of landsat satellite; it takes 5,7, or 8
gain	Band-specific sensor gain. It takes either "high" or "low"
sunelev	Angle of Sun elevation in degrees, as found in the meta data of the satellite image
date	DATE_ACQUIRED of the image, see metadata
weights	albedo bands weights for each band. list as in order c(band1,band2,band3,band4,band5,band6,band7) for landsat 5 and 7. The default is "auto" where sebal weights are used.
grescale	list as in order weights. Band-specific sensor values. The default for landsat 5 and 7 is from Chander et.al (2009). For Landsat 8 use ML Radiance scale in the metadata; i.e. Band specific multiplicative rescaling factor from the metadata (MTL file) (RADIANCE_MULT_BAND_x, where x is the band number).
brescale	list. Band-specific sensor values. The default for landsat 5 and 7 is from Chander et.al (2009). for Landsat 8 use AL Radiance scale in the metadata; i.e. Band specific multiplicative rescaling factor from the metadata (MTL file) (RADIANCE_MULT_BAND_x, where x is the band number).
ESUNx	list. The mean solar exo-atmospheric irradiance for each band. It is used for reflectance correction. The default is taken from SEBAL manual
Mp	Reflectance scale for Landsat 8. Band specific multiplicative rescaling factor from the metadata (MTL file) (REFLECTANCE_MULT_BAND_x, where x is the band number).
Ap	Reflectance scale for Landsat 8. Band specific additive rescaling factor from the metadata (MTL file) (REFLECTANCE_ADD_BAND_x, where x is the band number)
K1	Temperature constants for each Landsat image. Default is taken from landsat handbook
K2	Temperature constants for each Landsat image. Default is taken from landsat handbook
Rp	path radiance in the 10.4-12.5 um band. MODTRAN runs can accurately estimate this. The default is taken from Allen et al(2007) to be 0.91,
Rsky	narrow band downward thermal radiation. MODTRAN runs can accurately estimate this. The default is taken from Allen et al(2007) to be 1.32
tNB	narrow band transmissivity of air (10.4-12.5) um range. MODTRAN runs can accurately estimate this. The default is taken from Allen et al(2007) to be 0.866
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
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

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

x	a longitude or raster or vector object
y	a latitude or raster or vector object
var	numeric variable. The length should be same as the length of latitude and longitude
map	Raster* object that will be used to define the new coordinates.
epsg	coordinate reference number. The default is WGS 84: 4326

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

xydata	x and y data
i	index
map	geographical map

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

condition	logical
trueValue	logical
falseValue	logical

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

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

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

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.
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
latitude	geographical coordinates in decimal degrees. It should be negative for southern hemisphere
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.
iter.max	maximum iterations of sensible heat calculation.
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
welev	Weather station elevation in meters
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
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"
model	character. The type of model. It takes either "SEBAL" or "METRIC"
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

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

expr	an R expression to evaluate

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

folder	The folder or directory of landsat files
meta_pattern	The regular expression pattern to identify metadata file. The default is "MTL.txt". The function works on the new landsat metadata only.
remove_cloud	TRUE FALSE, 1, or 2. To indicate the cloud cover removal from the image. Default is FALSE. If it is set to 1, the overlapping areas of band 1 DN values between 120 and 255 and band 6 DN values of 102 to 128 are used. If it is set to 2 the union of the DN ranges are used; and clouds as well as urban, barren, quarries, rocks, and sand are removed; otherwise only clouds are removed. See Martinuzzi et al(2009) for more details. At the moment only landsat 5 and 7 clouds are removed.
gap_fill	Logical (TRUE or FALSE).It replaces strips in landsat 7 or NAs in other data set. It uses focal neighbourhood approach. It is an automated Version of USGS "Gap-Filling Landsat 7 SLC-off Single Scenes Using Erdas Imagine TM". Like USGS, it also takes approximately one hour to complete. see details at http://landsat.usgs.gov/ERDAS_Approach.php
gap_loop	numeric. The number of loops or repetitions that gap-fill should be run. The default is 8 for one landsat scene. If the landsat image is smaller, it can be reduced.
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)

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)

---