# -elevation: SRTM CGIAR 1 km (replace by global DEM?)

# -12 layers from land cover concensus (Jetz lab)

# -distance to coast NASA/NOOA (http://oceancolor.gsfc.nasa.gov/DOCS/DistFromCoast/)

# -24 LST layers: "climatology" produced from MOD11A1, LST (mean and obs) using script in step 1 of workflow

# 3) The output is a multiband file in tif format with projected covariates for the processing region/tile.             

#DATE: 05/27/2013                                                                                 

create__m_raster_region <-function(j,list_param){

  #This functions returns a subset of tiles from the modis grdid.

  #Arguments: raster name of the file,reference file with

  #Output: spatial grid data frame of the subset of tiles

  ## Parse input arguments

  raster_name <- list_param$raster_name[j] #list of raster ot project and crop

  reg_ref_rast <- list_param$reg_ref_rast 

  out_rast_name <- list_param$out_rast_name[j]

  ## Start #

  layer_rast<-raster(raster_name)

  new_proj<-proj4string(layer_rast)                  #Extract coordinates reference system in PROJ4 format

  region_temp_projected<-projectExtent(reg_ref_rast,CRS(new_proj))     #Project from current to region coord. system

  layer_crop_rast<-crop(layer_rast, region_temp_projected) #crop using the extent from the region tile

  #layer_projected_rast<-projectRaster(from=layer_crop_rast,crs=proj4string(reg_outline),method="ngb")

  layer_projected_rast<-projectRaster(from=layer_crop_rast,to=reg_ref_rast,method="ngb")

  writeRaster(layer_projected_rast, filename=out_rast_name,overwrite=TRUE)  

  return(out_rast_name)

}

#Works for now improve later?    

change_names_file_list<-function(list_name,out_suffix,out_prefix,extension,out_path=""){

  #Function to add suffix and prefix to list of file names

  lf_new_names_list<-vector("list",length(list_name)) #this will contain new names for files

  for (i in 1:length(list_name)){

    lf_name<-basename(list_name[[i]])

    lf_out_path<-dirname(list_name[[i]])

    data_name<-paste(out_prefix,sub(extension,"",lf_name),"_",sep="") #can add more later...

    raster_name<-paste(data_name,out_suffix, sep="") #out_suffix must include extension!!!

    if((lf_out_path!="") & (out_path=="")){

      lf_new_names_list[[i]]<-file.path(lf_out_path,raster_name)

    }else{

      lf_new_names_list[[i]]<-file.path(out_path,raster_name)

    }

  }

  return(unlist(lf_new_names_list))

}

  #It requires 16 arguments:

  #3)  lc_path : path to land cover consensus product

  #4)  infile_modis_grid : modis grid litles 

  #7)  infile_distoc : global distance to coast

  #8)  list_tiles_modis : tile for Venezuela and surrounding area

  #9)  infile_reg_outline : input region outline defined by polygon

  #10)  CRS_interp : local projection system

  #11) CRS_locs_WGS84 : CRS_locs_WGS84 #

  #12) out_region_name : generated on the fly?

  #13) out_suffix : added to the covariates stack/brick 

  #14) out_suffix_modis : suffix used in the LST averages used for the covariates (may vary to use older LST averages)

  #15) ref_rast_name: local raster name defining resolution, exent, local projection--. set on the fly??

  #16) hdfdir: directory where the LST averages are stored...

  #17) out_suffix_modis : suffix used in producing LST climatology 

  #18) covar_names : names of covariates

  infile_distoc<-list_param$infile_distoc #Canopy height

  infile_reg_outline <-list_param$infile_reg_outline   #input region outline defined by polygon

    reg_outline<-readOGR(dsn=dirname(infile_reg_outline), filename) # Read in the region outline

  #if no shapefile defining the study/processing area then create one using modis grid tiles

    infile_reg_outline <-paste("outline",out_region_name,"_",out_suffix,".shp",sep="")

    writeOGR(reg_outline,dsn= out_path,layer= sub(".shp","",infile_reg_outline), 

  #Need to clean up this section, later on, this should be written in a function to avoid repetition

  hdfdir_lst_avg<-file.path(hdfdir,"LST_averages")

  tmp_str2 <- mixedsort(list.files(path=hdfdir_lst_avg,pattern=pat_str2,full.names=TRUE)) #note that this assumes the LST averages are stored in hdfdir

  tmp_str1 <- mixedsort(list.files(path=hdfdir_lst_avg,pattern=pat_str1,full.names=TRUE))

  mosaic_list<-split(tmp_str2,list_date_names) #split list of files into groups using lst_date_name_pattern

  mosaic_list_nobs<-new_list

  out_rastnames_nobs<-paste(list_date_names,out_rastnames,sep="")

  mosaic_list_mean <-new_list #list ready for mosaicing

  ### Now mosaic tiles...Note that function could be improved to use less memory

  list_param_mosaic<-list(j,mosaic_list_mean,out_rastnames_mean,out_path)

  names(list_param_mosaic)<-c("j","mosaic_list","out_rastnames","out_path")

  mean_m_list <-mclapply(1:12, list_param=list_param_mosaic, mosaic_m_raster_list,mc.preschedule=FALSE,mc.cores = 12) #This is the end bracket from mclapply(...) statement

  list_param_mosaic<-list(j,mosaic_list_nobs,out_rastnames_nobs,out_path)

  nobs_m_list <-mclapply(1:12, list_param=list_param_mosaic, mosaic_m_raster_list,mc.preschedule=FALSE,mc.cores = 12) #This is the end bracket from mclapply(...) statement

  if (ref_rast_name==""){

  }else{

    ref_rast<-raster(ref_rast_name) #This is the reference image used to define the study/processing area

  ## Project mosaiced tiles if local projection is provided...

  out_suffix_lst <-paste(out_suffix,".tif",sep="")          

  mean_lst_list_outnames<-change_names_file_list(mean_m_list,out_suffix_lst,"reg_",".tif",out_path=out_path)     

  nobs_lst_list_outnames<-change_names_file_list(nobs_m_list,out_suffix_lst,"reg_",".tif",out_path=out_path)    

  if (ref_rast_name!=""){

    #list(mean_m_list)

    list_param_create_region<-list(j,raster_name=mean_m_list,reg_ref_rast=ref_rast,out_rast_name=mean_lst_list_outnames)

    mean_m_list <-mclapply(1:12, list_param=list_param_create_region, create__m_raster_region,mc.preschedule=FALSE,mc.cores = 6) #This is the end bracket from mclapply(...) statement

    list_param_create_region<-list(j,raster_name=nobs_m_list,reg_ref_rast=ref_rast,out_rast_name=nobs_lst_list_outnames)

    nobs_m_list <-mclapply(1:12, list_param=list_param_create_region, create__m_raster_region,mc.preschedule=FALSE,mc.cores = 6) #This is the end bracket from mclapply(...) statement

  }

  # if ref_name!="" need to reproject and clip!!! do this in mclapply for all the list of covar!!!

  ######################################

  #4) LCC land cover

  lc_list<-mixedsort(list.files(path=lc_path,pattern="con_1km_class_.*.tif",full.names=TRUE))

  out_suffix_lc<-paste(out_suffix,".tif",sep="")          

  lc_list_outnames<-change_names_file_list(lc_list,out_suffix_lc,"reg_",".tif",out_path=out_path)    

  list_param_create_region<-list(j,raster_name=lc_list,reg_ref_rast=ref_rast,out_rast_name=lc_list_outnames)

  lc_reg_list <-mclapply(1:length(lc_list), list_param=list_param_create_region, create__m_raster_region,mc.preschedule=FALSE,mc.cores = 9) #This is the end bracket from mclapply(...) statement

  lc_reg_s <- stack(lc_reg_list)

  #Now create mask based on water areas 

  LC12<-raster(lc_reg_s,layer=nlayers(lc_reg_s)) #this is open water

  LC_mask<-LC12

  LC_mask[LC_mask==100]<-NA

  LC_mask <- LC_mask > 100

  ###############################

  #5) DISTOC, distance to coast: Would be useful to have a distance to coast layer ready...

  #Crop and reproject Canopy height data

  lc_list <- c(infile_elev,infile_canheight,infile_distoc) #, lc_list #15 layers to reproject...

  out_suffix_l <-paste(out_suffix,".tif",sep="")          

  lc_list_outnames<-change_names_file_list(lc_list,out_suffix_l,"reg_",".tif",out_path=out_path)    

  list_param_create_region<-list(j,raster_name=lc_list,reg_ref_rast=ref_rast,out_rast_name=lc_list_outnames)

  list_covar_layers <- mclapply(1:3, list_param=list_param_create_region, create__m_raster_region,mc.preschedule=FALSE,mc.cores = 3) #This is the end bracket from mclapply(...) statement

  SRTM_reg <- raster(list_covar_layers[[1]])

  CANHEIGHT <- raster(list_covar_layers[[2]])

  distoc <- raster(list_covar_layers[[3]])

  #3) aspect, slope from STRM

  r<-stack(x,y,lon,lat,N,E,Nw,Ew,SRTM_reg,terrain_rast,CANHEIGHT,distoc)

  ##Write function to screen data values...

  #Screen LST for extreme values?

  #min_val<-(-15+273.16) #if values less than -15C then screen out (note the Kelvin units that will need to be changed later in all datasets)

  #LST[LST < (min_val)]<-NA

  #add screening here...

  #Mask and save layers

  #return reg_outline!!!

  covar_obj <-list(raster_name,infile_reg_outline)

  names(covar_obj) <-c("infile_covariates","infile_reg_outline")

  return(covar_obj)