##Comments and TODO:

#This script is meant to be for general processing tile by tile or region by region.

#We must decide on a local projection. This can best set up from the tile/region extent: for now use

#- produce a distance to ocean layer that is global.

#- do not keep output in memory??

lc_path<-"/home/layers/data/land-cover/lc-consensus-global"

elev_path<-"/home/layers/data/terrain/dem-cgiar-srtm-1km-tif"

infile1<-"worldborder_sinusoidal.shp"

infile2<-"modis_sinusoidal_grid_world.shp"

infile3<-"countries_sinusoidal_world.shp"

infile4<-"srtm_1km.tif"  #this is the global file: replace later with the input produced by the DEM team

list_tiles_modis = c('h11v08','h11v07','h12v07','h12v08','h10v07','h10v08') #tile for Venezuel and surrounding area

infile_reg_outline=""  #input region outline defined by polygon

CRS_interp<-"+proj=lcc +lat_1=43 +lat_2=45.5 +lat_0=41.75 +lon_0=-120.5 +x_0=400000 +y_0=0 +ellps=GRS80 +units=m +no_defs";

CRS_locs_WGS84<-CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +towgs84=0,0,0") #Station coords WGS84

out_suffix<-"_VE_01292013"

ref_rast_name<-""  #local raster name defining resolution, exent, local projection--. set on the fly??

                   #for the processing tile/region? This is a group fo six tiles for now.

  return(reg_outline_poly)

}

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

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

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

  layer_rast<-raster(raster_name)

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

  layer_crop_rast<-crop(layer_rast, region_temp_projected) #crop using the extent from teh 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")

##### STEP 1: Reading region or tile information to set the study or processing region

filename<-sub(".shp","",infile2)       #Removing the extension from file.

modis_grid<-readOGR(".", filename)     #Reading shape file using rgdal library

if (infile_reg_outline!=""){

  filename<-sub(".shp","",infile_reg_outline)   #Removing the extension from file.

if (infile_reg_outline==""){

  reg_outline<-create_modis_tiles_region(modis_grid,list_tiles_modis) #problem...this does not 

                                                                      #align with extent of modis LST!!!

}

#modis_tiles<-create_modis_tiles_region(modis_grid,list_tiles_modis)

##Create covariates for the stuy area: pull everything from the same folder?

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

#1) LST climatology: project, mosaic

pat_str2 <- glob2rx(paste("nobs","*.tif",sep=""))

tmp_str2<- mixedsort(list.files(pattern=pat_str2))

pat_str1 <- glob2rx(paste("mean","*.tif",sep=""))

#list_date_names<-as.character(0:11)

#lsit_date_names<-month.abb

list_date_names<-c("jan","feb","mar","apr","may","jun","jul","aug","sep","oct","nov","dec")

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

  j<-grep(list_date_names[i],mosaic_list,value=FALSE)

  names(mosaic_list)[j]<-list_date_names[i]

}

#reproject and crop if necessary

nobs_m_list<-mosaic_raster_list(mosaic_list,out_rastnames,out_path)

plot(stack(nobs_m_list))

plot(stack(mean_m_list))

#Use this as ref file for now?? Ok for the time being: this will need to change to be a processing tile.

ref_rast<-raster(mean_m_list[[1]]) 

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

##2) Crop and reproject Canopy height data

#Make it a function?

#Use GDAL instead?? system( 'gdalwarp...')

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

#3) Creating elev, aspect, slope from STRM

SRTM_name<-file.path(elev_path,infile4)

#new_proj<-proj4string(SRTM_rast)                  #Assign coordinates reference system in PROJ4 format

#region_temp_projected<-spTransform(modis_tiles,CRS(new_proj))     #Project from WGS84 to new coord. system

###calculate slope and aspect

terrain_rast<-terrain(SRTM_reg, opt=c("slope","aspect"),unit="degrees", neighbors=8) #, filename=\u2019\u2019, ...)

N<-cos(r1)

E<-sin(r1)

Nw<-sin(r2)*cos(r1)   #Adding a variable to the dataframe

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

  lc_name<-lc_list[[i]]

  data_name<-paste("reg_",sub(".tif","",lc_name),"_",sep="") #can add more later...

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

  writeRaster(lc_reg, filename=file.path(out_path,raster_name),overwrite=TRUE)  

setwd(out_path)

lc_reg_list<-mixedsort(list.files(pattern="reg_con.*.tif"))

lc_reg_s<-stack(lc_reg_list)

#Now combine forest classes...in LC1 forest, LC2, LC3, LC4 and LC6-urban...??

#create a local mask for the tile/processing region

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

#rc<-clump(LC12)

#tab_freq<-freq(rc)

#Modify at a later stage:

#raster<-"DISTOC_VE_01292013.rst"  

#Distance calculated in a global layer??

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

#6) X-Y coordinates and LAT-LONG: do not keep in memory?

xy_latlon<-project(xy, CRS_interp, inv=TRUE) # find lat long for projected coordinats (or pixels...)

#coord_s<-stack(x,y,lat,lon)

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

#Create a stack in tif format...

#? output name??

#using bil format more efficient??

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