#outrastnames <- "reg1_mosaic_weights.tif"

#list_param_mosaic <- list(list_r_weights,out_dir,outrastnames,file_format,NA_flag_val,out_suffix)

#r1_projected <- projectRaster(raster(list_r_weights[[1]]),r)

### First match weights from linear option

lf_files <- unlist(list_linear_r_weights)

list_param_raster_match <- list(lf_files,rast_ref,file_format,out_suffix,out_dir)

names(list_param_raster_match) <- c("lf_files","rast_ref","file_format","out_suffix","out_dir_str")

#debug(raster_match)

#r_test <- raster(raster_match(1,list_param_raster_match))

list_linear_weights_m <- mclapply(1:length(lf_files),FUN=raster_match,list_param=list_param_raster_match,mc.preschedule=FALSE,mc.cores = num_cores)                           

lf_files <- unlist(list_linear_r_weights_prod)

list_param_raster_match <- list(lf_files,rast_ref,file_format,out_suffix,out_dir)

names(list_param_raster_match) <- c("lf_files","rast_ref","file_format","out_suffix","out_dir_str")

num_cores <-11

list_linear_weights_prod_m <- mclapply(1:length(lf_files),FUN=raster_match,list_param=list_param_raster_match,mc.preschedule=FALSE,mc.cores = num_cores)                           

#### Second use use edge (dist) images

lf_files <- unlist(list_edge_r_weights)

list_param_raster_match <- list(lf_files,rast_ref,file_format,out_suffix,out_dir)

names(list_param_raster_match) <- c("lf_files","rast_ref","file_format","out_suffix","out_dir_str")

#debug(raster_match)

#r_test <- raster(raster_match(1,list_param_raster_match))

list_edge_weights_m <- mclapply(1:length(lf_files),FUN=raster_match,list_param=list_param_raster_match,mc.preschedule=FALSE,mc.cores = num_cores)                           

lf_files <- unlist(list_edge_r_weights_prod)

list_param_raster_match <- list(lf_files,rast_ref,file_format,out_suffix,out_dir)

names(list_param_raster_match) <- c("lf_files","rast_ref","file_format","out_suffix","out_dir_str")

num_cores <-11

list_edge_weights_prod_m <- mclapply(1:length(lf_files),FUN=raster_match,list_param=list_param_raster_match,mc.preschedule=FALSE,mc.cores = num_cores)                           

### third match wegihts from sine option

lf_files <- unlist(list_sine_r_weights)

list_param_raster_match <- list(lf_files,rast_ref,file_format,out_suffix,out_dir)

names(list_param_raster_match) <- c("lf_files","rast_ref","file_format","out_suffix","out_dir_str")

#debug(raster_match)

#r_test <- raster(raster_match(1,list_param_raster_match))

list_sine_weights_m <- mclapply(1:length(lf_files),FUN=raster_match,list_param=list_param_raster_match,mc.preschedule=FALSE,mc.cores = num_cores)                           

lf_files <- unlist(list_sine_r_weights_prod)

list_param_raster_match <- list(lf_files,rast_ref,file_format,out_suffix,out_dir)

names(list_param_raster_match) <- c("lf_files","rast_ref","file_format","out_suffix","out_dir_str")

num_cores <-11

list_sine_weights_prod_m <- mclapply(1:length(lf_files),FUN=raster_match,list_param=list_param_raster_match,mc.preschedule=FALSE,mc.cores = num_cores)                           

#get the list of weights into raster objects

list_args_linear_weights <- list_linear_weights_m

#list_args_weights <- (mixedsort(list.files(pattern="r_weights_m_.*.tif")))

list_args_linear_weights <- lapply(1:length(list_args_linear_weights), FUN=function(i,x){raster(x[[i]])},x=list_args_linear_weights)

#get the list of weights product into raster objects

list_args_linear_weights_prod <- list_linear_weights_prod_m

#list_args_weights_prod <- (mixedsort(list.files(pattern="r_weights_prod_m_.*.tif")))

list_args_linear_weights_prod <- lapply(1:length(list_args_linear_weights_prod), FUN=function(i,x){raster(x[[i]])},x=list_args_linear_weights_prod)

###

#get the list of edge weights into raster objects

list_args_edge_weights <- list_linear_weights_m

#list_args_weights <- (mixedsort(list.files(pattern="r_weights_m_.*.tif")))

list_args_edge_weights <- lapply(1:length(list_args_edge_weights), FUN=function(i,x){raster(x[[i]])},x=list_args_linear_weights)

#get the list of weights product into raster objects

list_args_linear_weights_prod <- list_linear_weights_prod_m

#list_args_weights_prod <- (mixedsort(list.files(pattern="r_weights_prod_m_.*.tif")))

list_args_linear_weights_prod <- lapply(1:length(list_args_linear_weights_prod), FUN=function(i,x){raster(x[[i]])},x=list_args_linear_weights_prod)

#get the list of sine weights into raster objects

list_args_sine_weights <- list_sine_weights_m

#list_args_weights <- (mixedsort(list.files(pattern="r_weights_m_.*.tif")))

list_args_sine_weights <- lapply(1:length(list_args_sine_weights), FUN=function(i,x){raster(x[[i]])},x=list_args_sine_weights)

#get the list of weights product into raster objects

list_args_sine_weights_prod <- list_sine_weights_prod_m

#list_args_weights_prod <- (mixedsort(list.files(pattern="r_weights_prod_m_.*.tif")))

list_args_sine_weights_prod <- lapply(1:length(list_args_sine_weights_prod), FUN=function(i,x){raster(x[[i]])},x=list_args_sine_weights_prod)

list_args_linear_weights$fun <- "sum" #we want the sum to compute the weighted mean

list_args_linear_weights$na.rm <- TRUE

list_args_linear_weights_prod$fun <- "sum"

list_args_linear_weights_prod$na.rm <- TRUE

list_args_sine_weights$fun <- "sum" #we want the sum to compute the weighted mean

list_args_sine_weights$na.rm <- TRUE

list_args_sine_weights_prod$fun <- "sum"

list_args_sine_weights_prod$na.rm <- TRUE

r_linear_weights_sum <- do.call(mosaic,list_args_linear_weights) #weights sum image mosaiced

r_linear_prod_sum <- do.call(mosaic,list_args_linear_weights_prod) #weights sum product image mosacied

r_m_linear_weighted_mean <- r_linear_prod_sum/r_linear_weights_sum #this is the mosaic using weighted mean...

r_sine_weights_sum <- do.call(mosaic,list_args_sine_weights) #weights sum image mosaiced

r_sine_prod_sum <- do.call(mosaic,list_args_sine_weights_prod) #weights sum product image mosacied

r_m_sine_weighted_mean <- r_sine_prod_sum/r_sine_weights_sum #this is the mosaic using weighted mean...

raster_name <- file.path(out_dir,paste("r_m_linear_weighted_mean_",out_suffix,".tif",sep=""))

writeRaster(r_m_linear_weighted_mean, NAflag=NA_flag_val,filename=raster_name,overwrite=TRUE)  

raster_name <- file.path(out_dir,paste("r_m_sine_weighted_mean_",out_suffix,".tif",sep=""))

writeRaster(r_m_sine_weighted_mean, NAflag=NA_flag_val,filename=raster_name,overwrite=TRUE)  

r_diff_weighted_mean <- r_m_linear_weighted_mean - r_m_sine_weighted_mean

#r_diff_weighted_mean<-r_diff_weighted_meam

r_diff_mean_linear <- r_m_mean - r_m_linear_weighted_mean 

r_diff_mean_sine <- r_m_mean - r_m_sine_weighted_mean 

r_m_mean_terrain <- terrain(r_m_mean,opt=c("slope","aspect"),unit="degrees")

r_m_sine_weighted_mean_terrain <- terrain(r_m_sine_weighted_mean,opt=c("slope","aspect"),unit="degrees")

r_m_linear_weighted_mean_terrain <- terrain(r_m_linear_weighted_mean,opt=c("slope","aspect"),unit="degrees")

res_pix <- 1200

col_mfrow <- 1 

row_mfrow <- 0.8

png(filename=paste("Figure2_mean_for_region_",region_name,"_",out_suffix,".png",sep=""),

    width=col_mfrow*res_pix,height=row_mfrow*res_pix)

plot(r_m_mean)

dev.off()

res_pix <- 1200

col_mfrow <- 1 

row_mfrow <- 0.8

png(filename=paste("Figure2_linear_weigthed_mean_for_region_",region_name,"_",out_suffix,".png",sep=""),

plot(r_m_linear_weighted_mean )

dev.off()

res_pix <- 1200

col_mfrow <- 1 

row_mfrow <- 0.8

png(filename=paste("Figure2_sine_weigthed_mean_for_region_",region_name,"_",out_suffix,".png",sep=""),

plot(r_m_sine_weighted_mean)

dev.off()

res_pix <- 1200

col_mfrow <- 1 

row_mfrow <- 0.8

png(filename=paste("Figure2_diff_linear_sine_weigthed_mean_for_region_",region_name,"_",out_suffix,".png",sep=""),

plot(r_diff_weighted_mean)

dev.off()

png(filename=paste("Figure2_diff_mean_linear_for_region_",region_name,"_",out_suffix,".png",sep=""),

plot(r_diff_mean_linear)

res_pix <- 1200

col_mfrow <- 1 

row_mfrow <- 0.8

png(filename=paste("Figure2_diff_mean_sine_for_region_",region_name,"_",out_suffix,".png",sep=""),

plot(r_diff_mean_sine)

#### plot terrain to emphasize possible edges..

res_pix <- 1200

col_mfrow <- 1 

row_mfrow <- 0.8

png(filename=paste("Figure2_slope_mean_linear_for_region_",region_name,"_",out_suffix,".png",sep=""),

plot(r_m_linear_weighted_mean_terrain,y=1)

png(filename=paste("Figure2_aspect_mean_linear_for_region_",region_name,"_",out_suffix,".png",sep=""),

    width=col_mfrow*res_pix,height=row_mfrow*res_pix)

plot(r_m_linear_weighted_mean_terrain,y=2)

dev.off()

png(filename=paste("Figure2_slope_mean_sine_for_region_",region_name,"_",out_suffix,".png",sep=""),

plot(r_m_sine_weighted_mean_terrain,y=1)

png(filename=paste("Figure2_aspect_mean_sine_for_region_",region_name,"_",out_suffix,".png",sep=""),

plot(r_m_sine_weighted_mean_terrain,y=2)

png(filename=paste("Figure2_slope_mean_for_region_",region_name,"_",out_suffix,".png",sep=""),

plot(r_m_mean_terrain,y=1)

png(filename=paste("Figure2_aspect_mean_for_region_",region_name,"_",out_suffix,".png",sep=""),

    width=col_mfrow*res_pix,height=row_mfrow*res_pix)

plot(r_m_mean_terrain,y=2)

dev.off()

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

#Ok testing on fake data to experiment and check methods:

##Quick function to generate test dataset

# make_raster_from_lf <- function(i,list_lf,r_ref){

#   vect_val <- list_lf[[i]]

#   r <-  r_ref

#   values(r) <-vect_val

#   #writeRaster...

#   return(r)

# }

# 

# vect_pred1 <- c(9,4,1,3,5,9,9,9,2)

# vect_pred2 <- c(10,3,1,2,4,8,7,8,2)

# vect_pred3 <- c(10,NA,NA,3,5,9,8,9,2)

# vect_pred4 <- c(9,3,2,NA,5,8,9,9,2)

# lf_vect_pred <- list(vect_pred1,vect_pred2,vect_pred3,vect_pred4)

# 

# vect_w1 <- c(0.2,0.5,0.1,0.3,0.4,0.5,0.5,0.3,0.2)

# vect_w2 <- c(0.3,0.4,0.1,0.3,0.4,0.5,0.7,0.1,0.2)

# vect_w3 <- c(0.5,0.3,0.2,0.2,0.3,0.6,0.7,0.3,0.2)

# vect_w4 <- c(0.2,0.5,0.3,0.3,0.4,0.5,0.5,0.2,0.2)

# lf_vect_w <- list(vect_w1,vect_w2,vect_w3,vect_w4)

# df_vect_w <-do.call(cbind,lf_vect_w)

# df_vect_pred <-do.call(cbind,lf_vect_pred)

# 

# tr_ref <- raster(nrows=3,ncols=3)

# 

# r_pred_l <- lapply(1:length(lf_vect_pred),FUN=make_raster_from_lf,list_lf=lf_vect_pred,r_ref=r_ref)

# r_w_l <- lapply(1:length(lf_vect_w),FUN=make_raster_from_lf,list_lf=lf_vect_w,r_ref=r_ref)

# 

# #r_w1<- make_raster_from_lf(2,list_lf=lf_vect_w,r_ref)

# 

# list_args_pred <- r_pred_l

# list_args_pred$fun <- "sum"

# 

# list_args_w <- r_w_l

# list_args_w$fun <- prod

# 

# r_test_val <-do.call(overlay,list_args) #sum

# r_test_w <-do.call(overlay,list_args_w) #prod

# 

# #need to do sumprod

# r1<- r_w_l[[1]]*r_pred_l[[1]]

# r2<- r_w_l[[2]]*r_pred_l[[2]]

# r3<- r_w_l[[3]]*r_pred_l[[3]]

# r4<- r_w_l[[4]]*r_pred_l[[4]]

# 

# r_pred <- stack(r_pred_l)

# r_w <- stack(r_w_l)

# 

# list_args_pred <- r_pred_l

# list_args_pred$fun <- mean

# list_args_pred$na.rm <- TRUE

# #r_sum_pred <-do.call(overlay,list_args_pred) #prod

# 

# #r_sum_pred <-do.call(mean,list_args_pred) #prod

# r_sum_pred <-do.call(mosaic,list_args_pred) #prod

# 

# list_args_pred$na.rm <- FALSE

# r_sum_pred <-do.call(overlay,list_args_pred) #prod

# 

# r_sum_pred <-do.call(overlay,list_args_w) #prod

# 

# list_args_w$fun <- sum

# r_sum_w <-do.call(overlay,list_args_w) #prod

# 

# r_m_w <- ((r1+r2+r3+r4)/(r_sum_w)) #mean weiated

#n33e to check the result!! especially the nubmer of valid pix val

#r_test_val <-do.call(overlay,list_args) #sum

#can do mosaic with sum?? for both weighted sum and val

#

#can use gdal calc...

#r_m <- r1 + r2

#name_method <- paste(interpolation_method,"_",y_var_name,"_",sep="")

##Use python code written by Alberto Guzman

#system("MODULEPATH=$MODULEPATH:/nex/modules/files")

#system("module load /nex/modules/files/pythonkits/gdal_1.10.0_python_2.7.3_nex")

#lf1 <- lf_world_pred_1000x3000

#lf2 <- lf_world_pred_1500x4500

#module_path <- ""

#module_path <- "/nobackupp6/aguzman4/climateLayers/sharedCode/"

#sh /nobackupp6/aguzman4/climateLayers/sharedCode/gdalCalDiff.sh file1.tif file2.tif output.tif

#/nobackupp6/aguzman4/climateLayers/sharedCode/mosaicUsingGdalMerge.py

#l_dates <- paste(day_to_mosaic,collapse=",",sep=" ")

#l_dates <- paste(day_to_mosaic,collapse=",")

## use region 2 first

#lf_out <- paste("diff_world_","1000_3000","by1500_4500_","mod1","_",l_dates,out_suffix,"_",file_format,sep="")

#for (i in 1:length(lf_out)){

#  out_file <- lf_out[i]

#  in_file1 <- lf1[i]

#  in_file2 <- lf2[i]

#    

#  cmd_str <- paste("sh", file.path(module_path,"gdalCalDiff.sh"),

#                 in_file1,

#                 in_file2,

#                 out_file,sep=" ")

#  system(cmd_str)

#

#}