#DATE: 06/20/2013                                                                                 

define_crs_from_extent_fun<-function(reg_outline,buffer_dist){

  #Screening values for raster stack

  #input: list_val_range: list of character strings comma separated

  #        e.g.: "mm_12,-15,50","mm_12,-15,50"

  #               variable name, min value, max value

  library(rgeos)

  #Buffer function not in use yet!! need query for specific matching MODIS tile !!! use gIntersection

  if (buffer_dist!=0){

    reg_outline_dissolved <- gUnionCascaded(reg_outline)  #dissolve polygons

    reg_outline <- gBuffer(reg_outline_dissolved,width=buffer_dist*1000)

  }

  #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=WGS84 +datum=WGS84 +units=m +no_defs";

  reg_centroid <- gCentroid(reg_outline)

  reg_centroid_WGS84 <- spTransform(reg_centroid,CRS_locs_WGS84) #get cooddinates of center of region in lat, lon

  reg_outline_WGS84 <- spTransform(reg_outline,CRS_locs_WGS84) #get cooddinates of center of region in lat, lon

  reg_extent <-extent( reg_outline_WGS84) #get boudning box of extent

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

  #Calculate projection parameters

  reg_lat_1 <- ymin(reg_extent)+((ymax(reg_extent)- ymin(reg_extent))/4)

  reg_lat_2 <- ymax(reg_extent)-((ymax(reg_extent)- ymin(reg_extent))/4)

  reg_lon_0 <- coordinates(reg_centroid_WGS84)[1]

  reg_lat_0 <- coordinates(reg_centroid_WGS84)[2]

  reg_x_0 <- 0

  reg_y_0 <- 0

  #Add false northing and false easting calucation for y_0,x_0

  #CRS_interp <- paste("+proj=lcc +lat_1=",43," +lat_2=",45.5," +lat_0=",41.75," +lon_0=",-120.5,

  #                    " +x_0=",0,"+y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs")

  CRS_interp <- paste("+proj=lcc +lat_1=",reg_lat_1," +lat_2=",reg_lat_2," +lat_0=",reg_lat_0," +lon_0=",reg_lon_0,

                      " +x_0=",reg_x_0," +y_0=",reg_y_0," +ellps=WGS84 +datum=WGS84 +units=m +no_defs",sep="")

  reg_outline_interp <- spTransform(reg_outline,CRS(CRS_interp)) #get cooddinates of center of region in lat, lon

  #add part to save projected file??

  #return reg_outline!!!

  reg_outline_obj <-list(reg_outline_interp,CRS_interp)

  names(reg_outline_obj) <-c("reg_outline","CRS_interp")

  return(reg_outline_obj)

} 

  buffer_dist <-list_param$buffer_dist

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

    #now add projection on the fly

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

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

             driver="ESRI Shapefile",overwrite_layer="TRUE")

    reg_outline_obj <- define_crs_from_extent_fun(reg_outline_modis,buffer_dist)

    reg_outline <-reg_outline_obj$reg_outline

    CRS_interp <-reg_outline_obj$CRS_interp

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

    ref_rast <-projectRaster(from=ref_rast_temp,crs=CRS_interp,method="ngb")

    proj4string(ref_rast) <- CRS_interp #Assign given reference system from master script...

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

  #test<-create__m_raster_region(1,list_param_create_region)

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

  #  #test<-create__m_raster_region(1,list_param_create_region)

  #  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

  #}

#DATE: 06/20/2013                                                                                 

source(file.path(script_path,"covariates_production_temperatures_06202013.R"))

infile_reg_outline=""  #input region outline defined by polygon: none for Venezuela

#infile_reg_outline <- "/home/parmentier/Data/IPLANT_project/Oregon_interpolation/Oregon_03142013/OR83M_state_outline.shp"  #input region outline defined by polygon: Oregon

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

#ref_rast_name<-"/home/parmentier/Data/IPLANT_project/Oregon_interpolation/Oregon_03142013/mean_day244_rescaled.rst"  #local raster name defining resolution, exent: oregon

buffer_dist<-0 #not in use yet, must change climatology step to make sure additional tiles are downloaded and LST averages

               #must also be calculated for neighbouring tiles.

CRS_interp<-"+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs";

#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_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";

#list of 18 parameters

                                  buffer_dist,list_val_range,out_suffix,out_suffix_modis,ref_rast_name,hdfdir,covar_names) 

                                      "buffer_dist","list_val_range","out_suffix","out_suffix_modis","ref_rast_name","hdfdir","covar_names")