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

Documented in sebkcstack sebkcstack.default

#' @title Stacking, gap-filling, cloud removal of Landsat image ET Model
#' 
#' @description This function prepares landsat 5, 7, and 8 data for  
#' \code{\link{landsat578}}. First, it reads meta data and returns many 
#' input parameters for \code{\link{landsat578}} and \code{\link{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.   
#' 
#' @param folder The folder or directory of landsat files
#' @param meta_pattern The regular expression pattern to identify 
#' metadata file. The default is "MTL.txt". The function works 
#' on the new landsat metadata only.
#' @param 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.
#' @param gap_fill Logical (TRUE or FALSE).It replaces strips in 
#' landsat 7 or NAs in other data set. It uses  \code{\link{raster}{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
#' @param 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. 
#' @param 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 \code{\link{writesebkc}}. 
#' The output file name is "sebkcstack.tif". 
# Whenever the function detects 
# a stacked file ("sebkcstack.tif") in the folder, it skips binding the layers 
# but reads meta data and performs other functions. 
# This is useful especailly after spending time refilling gap and masking 
# clouds. Therefore you can import stacked map from other applications 
# and name it "sebkcstack.tif" and gap-fill strips and/or remove clouds.
# Note that the imported file contains 7 layers. 
# You may delete this file if you want to stack from afresh.       
#' @author 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
#'
#' @return 
#' \itemize{
#' \item{data:} { raster stacked data. If crop is not NULL then  crop Extent object is returned}
#' \item{brescale:} { brescale of the stacked data}
#' \item{grescale:} { grescale of the stacked data}
#' \item{sunelev:} { sun elevation for stacked data}
#' \item{date:} { Acquired date of the raster stacked data}
#' \item{viewangle:} { view angle of landsat 8}
#' \item{azimuth:} { azimuth angle of the image}
#' \item{sensor:} { Type of satellite}
#' \item{K1:} { K1 of landsat 8 thermal band 10}
#' \item{K2:} { K2 of landsat 8 thermal band 10}
#' \item{cc:} { cloud cover percentage}
#' \item{quality:} { Image quality}
#' \item{row:} { WRS_ROW of the image }
#' \item{path:} { WRS_PATH of the image}
#' }
#'
#' @examples
#' \dontrun{
#' folder=system.file("extdata","stack",package="sebkc")
#' stack=sebkcstack(folder=folder)
#' #returns data and grescale
#' stack$data
#' stack$grescale
#' }
#' @export
#' @rdname sebkcstack
sebkcstack<-function(folder,meta_pattern="MTL.txt",remove_cloud=F,
    gap_fill=F,gap_loop=8,clip=NULL) UseMethod ("sebkcstack")
#' @export
#' @rdname sebkcstack
sebkcstack.default<-function(folder,meta_pattern="MTL.txt",remove_cloud=F,
                            gap_fill=F,gap_loop=8,clip=NULL){
  
crop=clip
metafile=list.files(folder,meta_pattern,full.names = TRUE)
metadata=read.delim(metafile,sep="=",stringsAsFactors = F)
SPACECRAFT_ID=(metadata[grep("SPACECRAFT_ID",metadata$GROUP),][[2]])
sensor=as.numeric(strsplit(SPACECRAFT_ID, "_")[[1]][2])
raster1crop=NULL
Mp=NA
Ap=NA



sunelev=as.numeric(metadata[grep("SUN_ELEVATION",metadata$GROUP),][[2]])
viewangle=as.numeric(metadata[grep("ROLL_ANGLE",metadata$GROUP),][[2]])
date=as.Date(metadata[grep("DATE_ACQUIRED",metadata$GROUP),][[2]])
RADIANCE_MULT_BAND_1=as.numeric(metadata[grep("RADIANCE_MULT_BAND_1.",metadata$GROUP),][[2]])[1]
RADIANCE_MULT_BAND_2=as.numeric(metadata[grep("RADIANCE_MULT_BAND_2.",metadata$GROUP),][[2]])[1]

RADIANCE_MULT_BAND_3=as.numeric(metadata[grep("RADIANCE_MULT_BAND_3.",metadata$GROUP),][[2]])[1]
RADIANCE_MULT_BAND_4=as.numeric(metadata[grep("RADIANCE_MULT_BAND_4.",metadata$GROUP),][[2]])[1]

RADIANCE_MULT_BAND_5=as.numeric(metadata[grep("RADIANCE_MULT_BAND_5.",metadata$GROUP),][[2]])[1]
RADIANCE_MULT_BAND_6=as.numeric(metadata[grep("RADIANCE_MULT_BAND_6.",metadata$GROUP),][[2]])[1]

RADIANCE_MULT_BAND_7=as.numeric(metadata[grep("RADIANCE_MULT_BAND_7.",metadata$GROUP),][[2]])[1]
if(sensor==8){
RADIANCE_MULT_BAND_9=as.numeric(metadata[grep("RADIANCE_MULT_BAND_9.",metadata$GROUP),][[2]])[1]
RADIANCE_MULT_BAND_10=as.numeric(metadata[grep("RADIANCE_MULT_BAND_10.",metadata$GROUP),][[2]])[1]
grescale=c(RADIANCE_MULT_BAND_2,RADIANCE_MULT_BAND_3,RADIANCE_MULT_BAND_4,
           RADIANCE_MULT_BAND_5,RADIANCE_MULT_BAND_6,RADIANCE_MULT_BAND_7,RADIANCE_MULT_BAND_10)
}else{
  grescale=c(RADIANCE_MULT_BAND_1,RADIANCE_MULT_BAND_2,RADIANCE_MULT_BAND_3,RADIANCE_MULT_BAND_4,
             RADIANCE_MULT_BAND_5,RADIANCE_MULT_BAND_6,RADIANCE_MULT_BAND_7)
}

RADIANCE_ADD_BAND_1=as.numeric(metadata[grep("RADIANCE_ADD_BAND_1.",metadata$GROUP),][[2]])[1]
RADIANCE_ADD_BAND_2=as.numeric(metadata[grep("RADIANCE_ADD_BAND_2.",metadata$GROUP),][[2]])[1]

RADIANCE_ADD_BAND_3=as.numeric(metadata[grep("RADIANCE_ADD_BAND_3.",metadata$GROUP),][[2]])[1]
RADIANCE_ADD_BAND_4=as.numeric(metadata[grep("RADIANCE_ADD_BAND_4.",metadata$GROUP),][[2]])[1]

RADIANCE_ADD_BAND_5=as.numeric(metadata[grep("RADIANCE_ADD_BAND_5.",metadata$GROUP),][[2]])[1]
RADIANCE_ADD_BAND_6=as.numeric(metadata[grep("RADIANCE_ADD_BAND_6.",metadata$GROUP),][[2]])[1]

RADIANCE_ADD_BAND_7=as.numeric(metadata[grep("RADIANCE_ADD_BAND_7.",metadata$GROUP),][[2]])[1]
if(sensor==8){
  RADIANCE_ADD_BAND_9=as.numeric(metadata[grep("RADIANCE_ADD_BAND_9.",metadata$GROUP),][[2]])[1]
  RADIANCE_ADD_BAND_10=as.numeric(metadata[grep("RADIANCE_ADD_BAND_10.",metadata$GROUP),][[2]])[1]
  brescale=c(RADIANCE_ADD_BAND_2,RADIANCE_ADD_BAND_3,RADIANCE_ADD_BAND_4,
             RADIANCE_ADD_BAND_5,RADIANCE_ADD_BAND_6,RADIANCE_ADD_BAND_7,RADIANCE_ADD_BAND_10)
}else{
  brescale=c(RADIANCE_ADD_BAND_1,RADIANCE_ADD_BAND_2,RADIANCE_ADD_BAND_3,RADIANCE_ADD_BAND_4,
             RADIANCE_ADD_BAND_5,RADIANCE_ADD_BAND_6,RADIANCE_ADD_BAND_7)
}


#if(sensor==8){
REFLECTANCE_MULT_BAND_1=as.numeric(metadata[grep("REFLECTANCE_MULT_BAND_1.",metadata$GROUP),][[2]])[1]
REFLECTANCE_MULT_BAND_2=as.numeric(metadata[grep("REFLECTANCE_MULT_BAND_2.",metadata$GROUP),][[2]])[1]

REFLECTANCE_MULT_BAND_3=as.numeric(metadata[grep("REFLECTANCE_MULT_BAND_3.",metadata$GROUP),][[2]])[1]
REFLECTANCE_MULT_BAND_4=as.numeric(metadata[grep("REFLECTANCE_MULT_BAND_4.",metadata$GROUP),][[2]])[1]

REFLECTANCE_MULT_BAND_5=as.numeric(metadata[grep("REFLECTANCE_MULT_BAND_5.",metadata$GROUP),][[2]])[1]
REFLECTANCE_MULT_BAND_6=as.numeric(metadata[grep("REFLECTANCE_MULT_BAND_6.",metadata$GROUP),][[2]])[1]

REFLECTANCE_MULT_BAND_7=as.numeric(metadata[grep("REFLECTANCE_MULT_BAND_7.",metadata$GROUP),][[2]])[1]
if(sensor==8){
  REFLECTANCE_MULT_BAND_9=as.numeric(metadata[grep("REFLECTANCE_MULT_BAND_9.",metadata$GROUP),][[2]])[1]
  REFLECTANCE_MULT_BAND_10=as.numeric(metadata[grep("REFLECTANCE_MULT_BAND_10.",metadata$GROUP),][[2]])[1]
  Mp=c(REFLECTANCE_MULT_BAND_2,REFLECTANCE_MULT_BAND_3,REFLECTANCE_MULT_BAND_4,
       REFLECTANCE_MULT_BAND_5,REFLECTANCE_MULT_BAND_6,REFLECTANCE_MULT_BAND_7,
       REFLECTANCE_MULT_BAND_10)
}else{
  Mp=c(REFLECTANCE_MULT_BAND_1,REFLECTANCE_MULT_BAND_2,REFLECTANCE_MULT_BAND_3,REFLECTANCE_MULT_BAND_4,
       REFLECTANCE_MULT_BAND_5,REFLECTANCE_MULT_BAND_6,REFLECTANCE_MULT_BAND_7)
}

REFLECTANCE_ADD_BAND_1=as.numeric(metadata[grep("REFLECTANCE_ADD_BAND_1.",metadata$GROUP),][[2]])[1]
REFLECTANCE_ADD_BAND_2=as.numeric(metadata[grep("REFLECTANCE_ADD_BAND_2.",metadata$GROUP),][[2]])[1]

REFLECTANCE_ADD_BAND_3=as.numeric(metadata[grep("REFLECTANCE_ADD_BAND_3.",metadata$GROUP),][[2]])[1]
REFLECTANCE_ADD_BAND_4=as.numeric(metadata[grep("REFLECTANCE_ADD_BAND_4.",metadata$GROUP),][[2]])[1]

REFLECTANCE_ADD_BAND_5=as.numeric(metadata[grep("REFLECTANCE_ADD_BAND_5.",metadata$GROUP),][[2]])[1]
REFLECTANCE_ADD_BAND_6=as.numeric(metadata[grep("REFLECTANCE_ADD_BAND_6.",metadata$GROUP),][[2]])[1]

REFLECTANCE_ADD_BAND_7=as.numeric(metadata[grep("REFLECTANCE_ADD_BAND_7.",metadata$GROUP),][[2]])[1]
if(sensor==8){
  REFLECTANCE_ADD_BAND_9=as.numeric(metadata[grep("REFLECTANCE_ADD_BAND_9.",metadata$GROUP),][[2]])[1]
  REFLECTANCE_ADD_BAND_10=as.numeric(metadata[grep("REFLECTANCE_ADD_BAND_10.",metadata$GROUP),][[2]])[1]
  Ap=c(REFLECTANCE_ADD_BAND_2,REFLECTANCE_ADD_BAND_3,REFLECTANCE_ADD_BAND_4,
       REFLECTANCE_ADD_BAND_5,REFLECTANCE_ADD_BAND_6,REFLECTANCE_ADD_BAND_7,REFLECTANCE_ADD_BAND_10)
}else{
  Ap=c(REFLECTANCE_ADD_BAND_1,REFLECTANCE_ADD_BAND_2,REFLECTANCE_ADD_BAND_3,REFLECTANCE_ADD_BAND_4,
       REFLECTANCE_ADD_BAND_5,REFLECTANCE_ADD_BAND_6,REFLECTANCE_ADD_BAND_7)
}


#}


K1_CONSTANT_BAND_10=as.numeric(metadata[grep("K1_CONSTANT_BAND_10",metadata$GROUP),][[2]])
K2_CONSTANT_BAND_10=as.numeric(metadata[grep("K2_CONSTANT_BAND_10",metadata$GROUP),][[2]])

if(length(K1_CONSTANT_BAND_10)==0){
  K1_CONSTANT_BAND_10=NULL  
}
if(length(K2_CONSTANT_BAND_10)==0){
  K2_CONSTANT_BAND_10=NULL  
}
time=strsplit(
  gsub("Z","",
  (metadata[grep("SCENE_CENTER_TIME",metadata$GROUP),])[2]
  ),":")
#time=round(
#  as.numeric(format(strptime(time, "%I:%M %p"), "%H")) +         # hours component
#    (as.numeric(format(strptime(time, "%I:%M %p"), "%M")) / 60),2  # rounded minutes component
#)
time=as.numeric(time[[1]][1])+(as.numeric(time[[1]][2])/60)
SUN_AZIMUTH=as.numeric(metadata[grep("SUN_AZIMUTH",metadata$GROUP),][[2]])
IMAGE_QUALITY=as.numeric(metadata[grep("IMAGE_QUALITY",metadata$GROUP),][[2]])
CLOUD_COVER=as.numeric(metadata[grep("CLOUD_COVER",metadata$GROUP),][[2]])
row=as.numeric(metadata[grep("WRS_ROW",metadata$GROUP),][[2]])[1]
path=as.numeric(metadata[grep("WRS_PATH",metadata$GROUP),][[2]])[1]
#image files
#chesck first if sebkcstack exists
if(file.exists(paste(folder,"sebkcstack.tif",sep="/"))){
data=brick(paste(folder,"sebkcstack.tif",sep="/"))  
band1=data[[1]]
band2=data[[2]]
band3=data[[3]]
band4=data[[4]]
band5=data[[5]]
band6=data[[6]]
band7=data[[7]]
if(sensor==8){
  band10=data[[10]]
  band11=data[[11]]
}

}else{
imagefiles=list.files(folder,NULL,full.names = FALSE)
band1=raster(paste(folder,grep("_B1.",imagefiles,value = TRUE)[1],sep="/"))
band2=raster(paste(folder,grep("_B2.",imagefiles,value = TRUE)[1],sep="/"))
band3=raster(paste(folder,grep("_B3.",imagefiles,value = TRUE)[1],sep="/"))
band4=raster(paste(folder,grep("_B4.",imagefiles,value = TRUE)[1],sep="/"))
band5=raster(paste(folder,grep("_B5.",imagefiles,value = TRUE)[1],sep="/"))
band6=raster(paste(folder,grep("_B6.",imagefiles,value = TRUE)[1],sep="/"))
band7=raster(paste(folder,grep("_B7.",imagefiles,value = TRUE)[1],sep="/"))
if(sensor==8){
  band9=raster(paste(folder,grep("_B9.",imagefiles,value = TRUE)[1],sep="/"))
  band10=raster(paste(folder,grep("_B10.",imagefiles,value = TRUE)[1],sep="/"))
  data=stack(band2,band3,band4,band5,band6,band7,band10)
  
}else{
  data=stack(band1,band2,band3,band4,band5,band6,band7)
}
}
name1=names(data)

#clousd masking
if((remove_cloud!=FALSE||remove_cloud!=F)&sensor!=8){
  
  if(remove_cloud==2){
  mask<-(data[[1]]>120&data[[1]]<255)|(data[[6]]>=102&data[[6]]<128)
  }else{
  mask<-(data[[1]]>120&data[[1]]<255)&(data[[6]]>=102&data[[6]]<128)
  }
  #band1=data[[1]]
  data[mask==1]<-NA
  names(data)=name1
}

#strips replacing
if(gap_fill==TRUE||gap_fill==T){
  #system.time()
  
  band1[band1==0]=NA
  band2[band2==0]=NA
  band3[band3==0]=NA
  band4[band4==0]=NA
  band5[band5==0]=NA
  band6[band6==0]=NA
  band7[band7==0]=NA
  if(sensor==8){
  band10[band10==0]=NA
  band11[band11==0]=NA
  }
  
  i=1
  while(i<=gap_loop){
    band1<- focal(band1, w=matrix(1,3,3), "mean",na.rm=TRUE,NAonly=TRUE) 
    band2<- focal(band2, w=matrix(1,3,3), "mean",na.rm=TRUE,NAonly=TRUE) 
    band3<- focal(band3, w=matrix(1,3,3), "mean",na.rm=TRUE,NAonly=TRUE) 
    band4<- focal(band4, w=matrix(1,3,3), "mean",na.rm=TRUE,NAonly=TRUE) 
    band5<- focal(band5, w=matrix(1,3,3), "mean",na.rm=TRUE,NAonly=TRUE) 
    band6<- focal(band6, w=matrix(1,3,3), "mean",na.rm=TRUE,NAonly=TRUE) 
    band7<- focal(band7, w=matrix(1,3,3), "mean",na.rm=TRUE,NAonly=TRUE) 
    if(sensor==8){
    band10<- focal(band10, w=matrix(1,3,3), "mean",na.rm=TRUE,NAonly=TRUE) 
    band11<- focal(band11, w=matrix(1,3,3), "mean",na.rm=TRUE,NAonly=TRUE) 
    }
    
    i=i+1
  }
  if(sensor==8){
  data=stack(band1,band2,band3,band4,band5,band6,band7,band10)
  }else{
    data=stack(band1,band2,band3,band4,band5,band6,band7)
  }
  names(data)=name1
  
  
  #system.time()
}
##cropping the image

if(!is.null(crop)){
  data=cropsebkc(x=data,y=crop)$x
  
}
#class(extent(crop))=="Extent"

factor<-list(data=data,brescale=brescale,grescale=grescale,Mp=Mp,Ap=Ap,sunelev=sunelev,
     date=date,viewangle=viewangle,sensor=sensor,path=path,row=row,
     K1=K1_CONSTANT_BAND_10,K2=K2_CONSTANT_BAND_10,cc=CLOUD_COVER,
     quality=IMAGE_QUALITY,azimuth=SUN_AZIMUTH,folder=folder,time=time)


factor$call<-match.call()

class(factor)<-"sebkcstack"
factor
}