#### FUNCTION USED IN SCRIPT

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

#### Parameters and constants  

infile_reg_outline <- "/data/project/layers/commons/data_workflow/inputs/region_outlines_ref_files/OR83M_state_outline.shp"  #input region outline defined by polygon: Oregon

met_stations_outfiles_obj_file<-"/data/project/layers/commons/data_workflow/output_data_365d_gam_fus_lst_test_run_07172013/met_stations_outfiles_obj_gam_fusion__365d_gam_fus_lst_test_run_07172013.RData"

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

y_var_name <- "dailyTmax"

out_prefix<-"_07182014"

out_dir<-"/home/parmentier/Data/IPLANT_project/paper_contribution_covar_analyses_tables_fig"

#out_dir<-"/home/parmentier/Data/IPLANT_project/paper_multitime_scale__analyses_tables"

create_out_dir_param = TRUE

#figure 2: methodological worklfow (generated outside R)

#figure 3: LST daily and monthly climatology

#figure 4: MAE/RMSE and distance to closest fitting station.

#Figure 7a: Spatial pattern of prediction for one day

#figure 7b: Spatial autocorrelation profile : Moran's vs lag distance

### Figure 1: Oregon study area

ghcn_dat <- readOGR(dsn=dirname(met_stations_obj$monthly_covar_ghcn_data),

        sub(".shp","",basename(met_stations_obj$monthly_covar_ghcn_data)))

interp_area <- readOGR(dsn=dirname(infile_reg_outline),sub(".shp","",basename(infile_reg_outline)))

interp_area_WGS84 <-spTransform(interp_area,CRS_locs_WGS84)         # Project from WGS84 to new coord. system

usa_map_2 <- usa_map[usa_map$NAME_1!="Alaska",] #remove Alaska

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

usa_map_OR <- usa_map_2[usa_map_2$NAME_1=="Oregon",] #get OR

title("Mean for month of January")

dev.off()

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

################### Figure 4. MAE/RMSE and distance to closest fitting station.

res_pix<-480

col_mfrow<-2

row_mfrow<-1

png(filename=file.path(out_dir,png_file_name),

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

par(mfrow=c(row_mfrow,col_mfrow))

list_param_plot<-list(list_dist_obj,col_t,pch_t,legend_text,mod_name,x_tick_labels,metric_name,

                      title_plot,y_lab_text,add_CI)

names(list_param_plot)<-c("list_dist_obj","col_t","pch_t","legend_text","mod_name","x_tick_labels","metric_name",

plot_dst_MAE(list_param_plot)

title(xlab="Distance to closest fitting station (km)")

barplot(l1$n_tb$res_mod1,names.arg=limit_val,

        ylab="Number of observations",

        xlab="Distance to closest fitting station (km)")

dev.off()

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

#Using baseline 2: lat,lon and elev

mod_compk1 <-kruskal.test(prop_tb$rmse ~ as.factor(prop_tb$method_interp)+as.factor(prop_tb$prop)) #Kruskal Wallis test

mod_prop <-lm(prop_tb$rmse ~ as.factor(prop_tb$method_interp)+as.factor(prop_tb$prop)) #This is significant!!

col_t<-c("red","blue","black")

pch_t<- 1:length(col_t)

CI_bar <- c(TRUE,TRUE,TRUE)

#add_CI <- c(TRUE,FALSE,FALSE)

####

res_pix<-480

col_mfrow<-2

row_mfrow<-1

png(filename=file.path(out_dir,png_file_name),

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

par(mfrow=c(row_mfrow,col_mfrow))

      xlab="Hold out validation/testing proportion",

      ylab="MAE (°C)")

metric_name<-"rmse"

list_param_plot<-list(list_prop_obj,col_t,pch_t,legend_text,mod_name,metric_name,add_CI,CI_bar)

names(list_param_plot)<-c("list_prop_obj","col_t","pch_t","legend_text","mod_name","metric_name","add_CI","CI_bar")

dev.off()

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

#read in relevant data:

## Calculate average difference for RMSE for all three methods

plot(prop_obj_gam_s$avg_tb$rmse ~ prop_obj_gam_s$avg_tb$prop, type="b",)

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

#set up the output file to plot

res_pix<-480

col_mfrow<-2

row_mfrow<-1

png(filename=file.path(out_dir,png_file_name),

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

par(mfrow=c(row_mfrow,col_mfrow))

pch_t<- 1:length(col_t)

##Make this a function???

y_range<-range(prop_obj_kriging$avg_tb$rmse,prop_obj_kriging_s$avg_tb$rmse)

plot_ac_holdout_prop<- function(l_prop,l_col_t,l_pch_t,add_CI,y_range){

dev.off()

mae_tmp<- data.frame(gam=tb1[tb1$pred_mod=="mod1",c("mae")],

                     kriging=tb3[tb3$pred_mod=="mod1",c("mae")],

                      gwr=tb4[tb4$pred_mod=="mod1",c(metric_name,"month")])

y_range<-range(unlist(month_data_list))

res_pix<-480

col_mfrow<-2

#row_mfrow<-2

row_mfrow<-1

png(filename=file.path(out_dir,png_file_name),

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

par(mfrow=c(row_mfrow,col_mfrow))

#Second plot

ylab_text<-"MAE (°C)"

xlab_text<-"Month"

boxplot(mae~month,data=month_data_list$gam,main="Monthly MAE boxplot", xlab=xlab_text,ylab=ylab_text,outline=FALSE)

dev.off()

#Now generate table 5

min_val <- 0

layout_m<-c(1,3) #one row two columns

#    height=480*layout_m[1],width=480*layout_m[2])

p<-levelplot(pred_temp_s,main="Interpolated Surfaces Method Comparison", ylab=NULL,xlab=NULL,

#col.regions=temp.colors(25))

#dev.off()

lf2 <- raster_prediction_obj_2$method_mod_obj[[index]][[y_var_name]]

lf2 #contains the models for gam

#min_val <- 0

layout_m<-c(4,3) #one row two columns

    height=480*layout_m[1],width=480*layout_m[2])

p<- levelplot(pred_temp_s,main="Interpolated Surfaces Model Comparison baseline 1", ylab=NULL,xlab=NULL,

###baseline 2: s(lat,lon) + s(elev)

l_obj<-vector("list",length=2)

l_obj[[2]]<-raster_prediction_obj_2

l_table <- vector("list",length=length(l_obj))   

l_s_table_term_tb <- vector("list",length(l_obj))

for (k in 1:length(l_obj)){

  raster_prediction_obj<- l_obj[[k]]

  #extract models for every day

  s.table_term_tb <- do.call(rbind.fill,s.table_term_list) 

  #Adding month to df

  #s.table_term_tb <- add_month_tag(s.table_term_tb)

  threshold_val<-c(0.01,0.05,0.1)

  s.table_term_tb$p_val_rec1 <- s.table_term_tb[["p-value"]] < threshold_val[1]

  s.table_term_tb$p_val_rec2 <- s.table_term_tb[["p-value"]] < threshold_val[2]

tb_mod4_LST_rec3 <- aggregate(s_table_LST_mod4$p_val_rec3~s_table_term_mod4$month,FUN=mean)

plot(tb_mod4_LST_rec3,type="l",ylim=c(0.2,1))

s_table_elev_mod4 <- subset(s.table_term_tb,mod_name=="mod4" & term_name=="s(elev_s)")

plot(tb_mod4_elev_rec3)

lines(tb_mod4_elev_rec3)

test1 <- subset(s.table_term_tb,mod_name=="mod1" & term_name=="s(elev_s)")

#### Now plot figure for paper:

names(tb_mod4_LST_rec3) <- c("month","p_val_rec3")

res_pix<-480

layout_m <- c(1,3) # works if set to this?? ok set the resolution...      

#date_selected <- c("20100101","20100901")

    height=res_pix*layout_m[1],width=res_pix*layout_m[2])

    #height=480*6,width=480*4)

          col=c("black"),

          type="b",

          ylab=list(label="\u0394RMSE between mod1 and mod4",cex=1.5),

          xlab=list(label="Month",cex=1.5),

          main=list(label="Proportion of significant LST term in mod4",cex=1.8))

p_dif <- xyplot(rmse_dif ~ 1:12,

          col=c("black"),

          type="b",

          ylab=list(label="\u0394RMSE between mod1 and mod4",cex=1.5),

          xlab=list(label="Month",cex=1.5),

p_dif <- update(p_dif, panel = function(...) {

            panel.abline(h = 0, lty = 2, col = "gray")

            panel.xyplot(...)

            panel.abline(h = 0, lty = 2, col = "gray")

            panel.xyplot(...)

        })

#grid.arrange(p1,p2,p3,ncol=1)

grid.arrange(p_prop,p_dif,p_cor,ncol=3)

dev.off()

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

### Prepare table 7: correlation matrix between covariates      

#Plot a sample variogram, and possibly a fitted model

#model 1 lat,lon and elev

layout_m <- c(1,2) # works if set to this?? ok set the resolution...

date_selected <- c("20100101","20100901")

png(paste("Figure9_paper_","_variogram_",date_selected[1],"_",date_selected[2],"_",out_prefix,".png", sep=""),

    #height=480*6,width=480*4)

#p3 <- list_plots_spt[[3]]

         ylim=c(0,9),

         ylab=list(label="Semivariance",cex=1.5),

         xlab=list(label="Distance (meter)",cex=1.5),

#plot(p1)

p241<-plot(raster_prediction_obj_3$method_mod_obj[[241]]$mod[[1]]$exp_var,raster_prediction_obj_3$method_mod_obj[[241]]$mod[[1]]$var_model,

         ylim=c(0,9),

         ylab=list(label="Semivariance",cex=1.5),

         xlab=list(label="Distance (meter)",cex=1.5),

grid.arrange(p1,p241,ncol=2)

dev.off()

#Combine both plot?     + plot info on sill, nugget and range? and most frequent model selected

list_kg_var_model <- lapply(1:365,function(i){raster_prediction_obj_3$method_mod_obj[[i]]$mod[[1]]$var_model})

list_kg_var_model <- lapply(1:365,function(i){raster_prediction_obj_3$method_mod_obj[[i]]$mod[[1]]$var_model})

list_kg_var_model

test <- do.call(rbind,list_kg_var_model)

tt2 <- subset(test,model=="Nug")      

tb_variogram["Nug"] <- (tt2$psill)

dates<-(extract_from_list_obj(raster_prediction_obj_1$method_mod_obj,"sampling_dat"))$date #get vector of dates

tb_variogram["date"] <- dates

tb_variogram$month <- add_month_tag(tb_variogram)

layout_m <- c(2,2) # works if set to this?? ok set the resolution...      

#date_selected <- c("20100101","20100901")

    height=480*layout_m[1],width=480*layout_m[2])

    #height=480*6,width=480*4)

dev.off()

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

index <- 244

names_mod <- names(raster_prediction_obj_2$method_mod_obj[[index]][[y_var_name]]) #names of models to plot

#in_dir2 <-"/home/parmentier/Data/IPLANT_project/Oregon_interpolation/Oregon_03142013/output_data_365d_gam_day_lst_comb4_08152013"

  #p1 <- levelplot(elev,scales = list(draw = FALSE), colorkey = FALSE,col.regions=rev(terrain.colors(255)),contour=T)

  #add legend..par.settings = GrTheme

  cx <- ((data_v[[res_model_name]])*2)

  #p2 <- spplot(data_v,res_model_name, cex=1 * cx,main=paste("Residuals for ",res_model_name," ",datelabel,sep=""),

  #             col.regions=matlab.like(25))

  p3 <- p2 + p1 + p2 #to force legend...

  #p1 <- spplot(interp_area)

  #p3 <- p1+p2 #to force legend...

  list_p[[k]] <- p3

}

    height=480*layout_m[1],width=480*layout_m[2])

dev.off()   

l_formulas<-(extract_from_list_obj(raster_prediction_obj_2$method_mod_obj,"formulas")) #get vector of dates

#                            y_var ~ s(lat,lon)

#                y_var ~ s(lat,lon) + s(elev_s)

}

layout_m<-c(1,3) # works if set to this?? ok set the resolution...

    height=480*layout_m[1],width=480*layout_m[2])

grid.arrange(list_p_mae[[1]],list_p_mae[[2]],list_p_mae[[3]],ncol=3)

dev.off()   

list_data_v <- lapply(1:365,function(i){raster_prediction_obj_2$validation_mod_obj[[i]]$data_v}) #testing with residuals

l_formulas<-(extract_from_list_obj(raster_prediction_obj_2$method_mod_obj,"formulas")) #get vector of dates

test <- do.call(rbind,list_data_v)        

data_v_ag <-test                

cor(test$res_mod1,test$elev)

cor(test$res_mod2,test$elev,use="complete.obs") #decrease in corellation when using elev

cor(test$res_mod1,test$LST,,use="complete.obs")

cor(test$res_mod5,test$LST,use="complete.obs") #decrease in correlation when using LST

brks<-c(0,500,1000,1500,2000)

lab_brks<-1:4

elev_rcstat<-cut(data_v_ag$elev,breaks=brks,labels=lab_brks,right=F)

#Now set up plotting device

layout_m <- c(1,3) # works if set to this?? ok set the resolution...      

    height=480*layout_m[1],width=480*layout_m[2])

    #height=480*6,width=480*4)

p_bw1<-bwplot(data_v_ag$res_mod1~elev_rcstat,do.out=F,ylim=c(-15,15),

         ylab=list(label="Residuals (deg C)",cex=1.5),

         xlab=list(label="Elevation classes (meter)",cex=1.5),

p_bw2 <- bwplot(data_v_ag$res_mod5~elev_rcstat,do.out=F,ylim=c(-15,15),

         ylab=list(label="Residuals (deg C)",cex=1.5),

         xlab=list(label="Elevation classes (meter)",cex=1.5),

         scales = list(x = list(at = c(1, 2, 3, 4), 

         ylab=list(label="Residuals (deg C)",cex=1.5),

         xlab=list(label="Elevation classes (meter)",cex=1.5),

         scales = list(x = list(at = c(1, 2, 3, 4), 

grid.arrange(p_bw1,p_bw2,p_bw3,ncol=3)

dev.off()

###################### END OF SCRIPT #######################