#MODIFIED ON: 07/02/2014            

function_analyses_paper1 <-"contribution_of_covariates_paper_interpolation_functions_10222013.R" #first interp paper

function_analyses_paper2 <-"multi_timescales_paper_interpolation_functions_05052014.R"

in_dir1 <- "/data/project/layers/commons/NEX_data/output_run3_global_analyses_06192014/output"

out_dir <-"/data/project/layers/commons/NEX_data/output_run3_global_analyses_06192014/"

df_tile_processed <- read.table(file=file.path(out_dir,paste("df_tile_processed_",out_prefix,".txt",sep="")),sep=",")

#in_dir_list <- file.path(in_dir1,read.table(file.path(in_dir1,"processed.txt"))$V1)

y_var_name <- "dailyTmax"

out_prefix<-"run3_global_analyses_06192014"

df_tile_processed <- read.table(file=file.path(out_dir,paste("df_tile_processed_",out_prefix,".txt",sep="")),sep=","),paste("pred_data_month_info_",out_prefix,".txt",sep="")) <- read.table(file=file.path(out_dir,paste("pred_data_month_info_",out_prefix,".txt",sep="")),sep=",")

list_shp_world <- list.files(in_dir_shp,".shp")

l_shp <- unlist(lapply(1:length(list_shp_world),FUN=function(i){paste(strsplit(list_shp_world[i],"_")[[1]][2:3],collapse="_")}))

list_shp_reg_files <- as.character(df_tile_processed$tile_coord)

#matching_index <- match(l_shp,list_shp_reg_files)

matching_index <- match(list_shp_reg_files,l_shp)

df_tile_processed$shp_files <-list_shp_world[matching_index]

list_shp_reg_files <- df_tile_processed$shp_files

list_shp_reg_files<- file.path(in_dir_shp,list_shp_reg_files)

if region_name=="reg1"{

  usa_map <- getData('GADM', country=c('USA'), level=1) #Get US map

  can_map <- getData('GADM', country=c('CAN'), level=1) #Get Canada map

  mex_map <- getData('GADM', country=c('MEX'), level=1) #Get Mexico map

  usa_map <- usa_map[usa_map$NAME_1!="Hawaii",] #remove Hawai 

  #reg_layer <- combine all three?

  proj4string(reg_layer) <- CRS_locs_WGS84

title(paste("Tiles location 10x10 degrees for ", region_name,sep=""))

model_name <- unique(tb$pred_mod)

  title(paste("Average RMSE per tile and by ",model_name[i]))

## Get ID from tile number...

ID_str <- unlist(lapply(1:nrow(df_tile_processed),function(i){unlist(strsplit(as.character(df_tile_processed$tile_id[i]),"_"))[2]}))

assign_FID_spatial_polygons_df <-function(list_spdf,ID_str=NULL){

  list_spdf_tmp <- vector("list",length(list_spdf))

  if(is.null(ID_str)){

    nf <- 0 #number of features

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

    #    shp1 <- list_spdf[[i]]

    #    f <- nrow(shp1)

    #    nf <- nf + f

    #}

    #This assumes that the list has one feature per item list

    nf <- length(list_spdf)

    ID_str <- as.character(1:nf)

  }

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

    #test=spRbind(shps_tiles[[1]],shps_tiles[[2]])

    shp1 <- list_spdf[[i]]

    shp1$FID <- ID_str

    shp1<- spChFIDs(shp1, as.character(shp1$FID)) #assign ID

    list_spdf_tmp[[i]]  <-shp1

  }

  return(list_spdf_tmp)

combine_spatial_polygons_df_fun <- function(list_spdf_tmp,ID_str=NULL){

  if(is.null(ID_str)){

    #call function

    list_spdf_tmp <- assign_FID_spatial_polygons_df

  }

  combined_spdf <- list_spdf_tmp[[1]]

  for(i in 2:length(list_spdf_tmp)){

    combined_spdf <- rbind(combined_spdf,list_spdf_tmp[[i]])

    #sapply(slot(shps_tiles[[2]], "polygons"), function(x) slot(x, "ID"))

    #rownames(as(alaska.tract, "data.frame"))

  }

  return(combined_spdf)

p1 <- spplot(usa_map,"ISO",colorkey=FALSE) #Use ISO instead of NAME_1 to see no color?

p2 <- spplot(can_map,"ISO",colorkey=FALSE) #Use ISO instead of NAME_1 to see no color?

p3 <- spplot(mex_map,"ISO",colorkey=FALSE) #Use ISO instead of NAME_1 to see no color?

p0 +p1+p2+p3

#p1 <- spplot(reg_layer,"ISO",colorkey=FALSE) #Use ISO instead of NAME_1 to see no color?

p0 + p1 + p2 + p3

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

### Figure 10: Plot the number of stations in a processing tile

boxplot(n_mod~tile_id,pred_data_month_info)

boxplot(n_mod~date,pred_data_month_info)