Project

General

Profile

« Previous | Next » 

Revision f7349270

Added by Benoit Parmentier almost 11 years ago

modifications figures for paper and adding many functions including for Moran's lags profiles

View differences:

climate/research/oregon/interpolation/multi_timescale_paper_interpolation.R
5 5
#Figures, tables and data for the  paper are also produced in the script.
6 6
#AUTHOR: Benoit Parmentier 
7 7
#CREATED ON: 10/31/2013  
8
#MODIFIED ON: 12/18/2013            
9
#Version: 1
8
#MODIFIED ON: 12/23/2013            
9
#Version: 2
10 10
#PROJECT: Environmental Layers project                                     
11 11
#################################################################################################
12 12

  
......
39 39
#### FUNCTION USED IN SCRIPT
40 40

  
41 41
function_analyses_paper1 <-"contribution_of_covariates_paper_interpolation_functions_10222013.R"
42
function_analyses_paper2 <-"multi_timescales_paper_interpolation_functions_12182013.R"
42
function_analyses_paper2 <-"multi_timescales_paper_interpolation_functions_12232013.R"
43 43

  
44 44
##############################
45 45
#### Parameters and constants  
......
88 88
raster_obj_file_8 <- "raster_prediction_obj_kriging_fusion_dailyTmax_365d_kriging_fss_lst_comb5_11052013.RData"
89 89
raster_obj_file_9 <- "raster_prediction_obj_gwr_fusion_dailyTmax_365d_gwr_fss_lst_comb5_11052013.RData"
90 90

  
91
## holdout
91
## Results from monthly holdout 0 to 70%
92 92

  
93 93
raster_obj_file_10 <- "raster_prediction_obj_gam_fusion_dailyTmax_365d_gam_fss_lst_mults_0_70_comb5_11082013.RData"
94 94
raster_obj_file_11 <- "raster_prediction_obj_kriging_fusion_dailyTmax_365d_kriging_fss_lst_mults_0_70_comb5_11132013.RData"
......
97 97
raster_obj_file_14 <- "raster_prediction_obj_kriging_CAI_dailyTmax_365d_kriging_cai_lst_mults_0_70_comb5_11272013.RData"
98 98
raster_obj_file_15 <- "raster_prediction_obj_gwr_CAI_dailyTmax_365d_gwr_cai_lst_mults_0_70_comb5_11222013.RData"
99 99

  
100
out_dir<-"/home/parmentier/Data/IPLANT_project/paper_multitime_scale__analyses_tables_fig_09032013"
100
#List of object with 0% monthly hold out
101
list_raster_obj_files  <- list(file.path(in_dir1,raster_obj_file_1),file.path(in_dir2,raster_obj_file_2),
102
                               file.path(in_dir3a,raster_obj_file_3a),file.path(in_dir3b,raster_obj_file_3b),
103
                               file.path(in_dir4,raster_obj_file_4),file.path(in_dir5,raster_obj_file_5),
104
                               file.path(in_dir6,raster_obj_file_6),file.path(in_dir7,raster_obj_file_7),
105
                               file.path(in_dir8,raster_obj_file_8),file.path(in_dir9,raster_obj_file_9))
106
 
107
names(list_raster_obj_files)<- c("gam_daily","kriging_daily","gwr_daily_a","gwr_daily_b",
108
                                 "gam_CAI","kriging_CAI","gwr_CAI",
109
                                 "gam_fss","kriging_fss","gwr_fss")
110

  
111
y_var_name <- "dailyTmax"
112
out_prefix<-"analyses_12232013"
113
out_dir<-"/home/parmentier/Data/IPLANT_project/paper_multitime_scale__analyses_tables"
114
out_dir <-paste(out_dir,"_",out_prefix,sep="")
115

  
116
if (!file.exists(out_dir)){
117
  dir.create(out_dir)
118
  #} else{
119
  #  out_path <-paste(out_path..)
120
}
101 121
setwd(out_dir)
102 122

  
103 123
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
104 124
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"
105 125
CRS_locs_WGS84<-CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +towgs84=0,0,0") #Station coords WGS84
106
y_var_name <- "dailyTmax"
107
out_prefix<-"analyses_12022013"
108 126
ref_rast_name<- "/data/project/layers/commons/data_workflow/inputs/region_outlines_ref_files/mean_day244_rescaled.rst"                     #This is the shape file of outline of the study area. #local raster name defining resolution, exent, local projection--. set on the fly??
109 127
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
110 128
ref_rast_name <-"/data/project/layers/commons/data_workflow/inputs/region_outlines_ref_files/mean_day244_rescaled.rst"  #local raster name defining resolution, exent: oregon
......
119 137
met_obj_file_1 <- "met_stations_outfiles_obj_gam_daily__365d_gam_daily_lst_comb5_11012013.RData"
120 138

  
121 139
#Load objects containing training, testing, models objects 
122

  
123 140
met_stations_obj <- load_obj(file.path(in_dir1,met_obj_file_1))
124 141
covar_obj <-load_obj(file.path(in_dir1,covar_obj_file_1)) #Reading covariates object for GAM daily method
125 142
infile_covariates <- covar_obj$infile_covariates
......
129 146
s_raster <- brick(infile_covariates)
130 147
names(s_raster)<-covar_names
131 148

  
132
#raster_prediction_obj_1 <-load_obj(file.path(in_dir1,raster_obj_file_1)) #comb5 gam_daily
149
###################################################################
150
############### BEGIN SCRIPT ##################################
133 151

  
134
############### BEGIN SCRIPT #################
152
#############################PART I: Generate tables for paper ##########
153
#Table 1: Covariates used (not produced in R)
154
#Table 2: Covariates models and functional form (not produced in R)
155
#Table 3: Combinations of models tested (not produced in R)
156
#Table 4. Average RMSE for all models acroos year 2010 wiht 30% daily hold out
157
#Table 5: Average monhtly RMSE for CAI model 7, climatology surfaces,  CAI and single time scale methods 
135 158

  
136
############
137
##### 1) Generate: Table 4. Contribution of covariates using validation accuracy metrics
159
##### Table 4: Contribution of covariates using validation accuracy metrics
138 160
## This is a table of accuracy  
139 161

  
140
list_raster_obj_files  <- list(file.path(in_dir1,raster_obj_file_1),file.path(in_dir2,raster_obj_file_2),
141
                               file.path(in_dir3a,raster_obj_file_3a),file.path(in_dir3b,raster_obj_file_3b),
142
                               file.path(in_dir4,raster_obj_file_4),file.path(in_dir5,raster_obj_file_5),
143
                               file.path(in_dir6,raster_obj_file_6),file.path(in_dir7,raster_obj_file_7),
144
                               file.path(in_dir8,raster_obj_file_8),file.path(in_dir9,raster_obj_file_9))
145
 
146
names(list_raster_obj_files)<- c("gam_daily","kriging_daily","gwr_daily","gwr_daily",
147
                                 "gam_CAI","kriging_CAI","gwr_CAI",
148
                                 "gam_fss","kriging_fss","gwr_fss")
149

  
150 162
summary_metrics_v_list<-lapply(list_raster_obj_files,FUN=function(x){x<-load_obj(x);x[["summary_metrics_v"]]$avg$rmse})                           
151 163
names(summary_metrics_v_list)
152 164

  
......
160 172
table_kriging <- table_kriging[1:7,]                         
161 173

  
162 174
#for gwr models
163
table_gwr <- summary_metrics_v_list[grep("gwr",names(summary_metrics_v_list))]
175
table_gwr <- summary_metrics_v_list[grep("gwr",names(summary_metrics_v_list))] #select models related to gwr
176
table_gwr$gwr_daily <- c(table_gwr[[1]],table_gwr[[2]]) #combine both gwr object results
177
table_gwr <- table_gwr[c(5,3,4)] #reorder list and drop first two objects
164 178
table_gwr <- do.call(cbind,table_gwr)
165 179
table_gwr <- table_gwr[1:7,]                         
166 180

  
......
176 190
#strsplit(list_formulas,"~")
177 191
                             
178 192
table4_paper$Formulas<-rep(list_formulas,3)                             
179
table4_paper<-table4_paper[(c(5,4,1,2,3))]                             
193
table4_paper<-table4_paper[(c(5,4,1,2,3))]  #reordering columns                           
180 194

  
181 195
#Testing siginificance between models
182 196
                             
......
195 209
file_name<-paste("table4_multi_timescale_paper","_",out_prefix,".txt",sep="")
196 210
write.table(table4_paper,file=file_name,sep=",")
197 211

  
198
####################################
199
####### Now create figures #############
212
################ Table 5 ################
213
#### Monthly RMSE for GAM CAI Climatology, GAM CAI daily and Single  time scale
214
## This is a table of accuracy  
215

  
216
#raster_prediction_obj_4 <-load_obj(file.path(in_dir4,raster_obj_file_4)) #comb5 gam_dir
217

  
218
tb_m_gam_CAI <- load_obj(file.path(in_dir13,raster_obj_file_13))$tb_month_diagnostic_v #getting objet
219
tb_gam_CAI <- load_obj(file.path(in_dir13,raster_obj_file_13))$tb_diagnostic_v #getting objet
220
tb_gam_daily <- load_obj(file.path(in_dir1,raster_obj_file_1))$tb_diagnostic_v #getting objet
221

  
222
names_mod <- paste("mod",1:7,sep="")
223
list_stat_tb_gam_CAI <-calc_stat_by_month_tb(names_mod,tb_gam_CAI,month_holdout=T)
224
list_stat_tb_gam_daily <-calc_stat_by_month_tb(names_mod,tb_gam_daily,month_holdout=F)
225

  
226
list_stat_tb_gam_CAI$avg
227

  
228
clim_rmse <- subset(tb_m_gam_CAI,prop==30 & pred_mod== "mod7",select="rmse")
229
CAI_rmse <- subset(list_stat_tb_gam_CAI$avg,prop_month==30 & pred_mod== "mod7",select="rmse")
230
daily_rmse <- subset(list_stat_tb_gam_daily$avg, pred_mod== "mod7",select="rmse")
231

  
232
table5 <- data.frame(clim=clim_rmse,CAI_month=CAI_rmse,daily_month=daily_rmse)
233
table5 <- round(table5,digit=3) #roundto three digits teh differences
234

  
235
table5$month <- month.abb
236
table5 <- table5[,c(4,1,2,3)]
237

  
238
###Table 5, writeout
239
names_table_col <-c("Month","Long term Climatology RMSE",
240
                    "Average Monthly RMSE from Daily pred (CAI)"
241
                    ,"Average Monthly RMSE from Daily pred (Single time scale)") # 
242

  
243
names(table5)<- names_table_col
244
print(table5) #show resulting table
245

  
246
#Now write out table 5
247

  
248
file_name<-paste("table5_multi_timescale_paper","_",out_prefix,".txt",sep="")
249
write.table(table5,file=file_name,sep=",")
250

  
251
#######################################################
252
####### PART 2: generate figures for paper #############
200 253

  
201 254
#figure 1: Study area OR
202 255
#Figure 2: LST climatology production: daily mean compared to monthly mean
......
208 261
#Figure 8: Transect profiles for one day 
209 262
#Figure 9: Image differencing and land cover: spatial patterns   
210 263
#Figure 10: LST and Tmax at stations, elevation and land cover.
211
#Figure 11: Spatial lag profiles and stations data  
264
#Figure 11: Spatial lag profiles for predicted surfaces 
212 265
#Figure 12: Daily deviation
213 266
#Figure 13: Spatial correlogram at stations, LST and elevation every 5 lags
214 267

  
......
264 317
dev.off()
265 318

  
266 319
################################################
267
######### Figure 2:  Method comparison workflow 
320
######### Figure 2: LST averaging: daily mean compared to monthly mean
268 321

  
269
# Workflow figure is not generated in R
270

  
271
################################################
272
######### Figure 3: LST averaging: daily mean compared to monthly mean
273

  
274
lst_md<-raster(ref_rast_name)
322
lst_md <- raster(ref_rast_name)
275 323
projection(lst_md)<-projection(s_raster)
276 324
lst_mm_09<-subset(s_raster,"mm_09")
277 325

  
......
279 327
lst_md<- lst_md - 273.16
280 328
lst_mm_01<-subset(s_raster,"mm_01")
281 329

  
282
png(filename=paste("Comparison_daily_monthly_mean_lst",out_prefix,".png",sep=""),width=960,height=480)
330
png(filename=paste("Figure2_paper_Comparison_daily_monthly_mean_lst",out_prefix,".png",sep=""),width=960,height=480)
283 331
par(mfrow=c(1,2))
284 332
plot(lst_md)
285 333
plot(interp_area,add=TRUE)
286
title("Mean January 1")
334
title("Mean for January 1")
287 335
plot(lst_mm_01)
288 336
plot(interp_area,add=TRUE)
289 337
title("Mean for month of January")
290 338
dev.off()
291 339

  
340
################################################
341
######### Figure 3: Prediction procedures: direct, CAI and FSS 
342
#(figure created outside R)
343

  
344

  
292 345
################################################
293 346
######### Figure 4. RMSE multi-timescale mulitple hold out for FSS and CAI
294 347

  
......
317 370

  
318 371
list_tb <- c(tb_s_list,tb_v_list,tb_ms_list,tb_mv_list) #combined in one list
319 372
ac_metric <- "rmse"
373
#plot_names <- c("RMSE for gam_FSS","RMSE for kriging_FSS","RMSE for gwr_FSS",
374
#               "RMSE for gam_CAI","RMSE for kriging_CAI","RMSE for gwr_CAI")
375
names_mod <- paste("mod",1:7,sep="")
376
plot_names <- names(list_tb) #this is the default names for the plots
377
#now replace names for relevant figure used later on.
378
plot_names[7:12] <- c("RMSE for gam_FSS","RMSE for kriging_FSS","RMSE for gwr_FSS",
379
               "RMSE for gam_CAI","RMSE for kriging_CAI","RMSE for gwr_CAI")
320 380
#Quick function to explore accuracy make this a function to create solo figure...and run only a subset...
321 381

  
322
list_plots <- plot_accuracy_by_holdout_fun(list_tb,ac_metric)
323
names(list_tb)
324
tb_v_list 
325
#For paper...
326
#Combine figures... tb_v for GWR, Kriging and GAM for both FSS and CAI
327
#grid.arrange(p1,p2, ncol=2)
382
list_plots <- plot_accuracy_by_holdout_fun(list_tb,ac_metric,plot_names,names_mod)
383

  
384
##For paper...combine figures... tb_v for GWR, Kriging and GAM for both FSS and CAI
328 385

  
329 386
layout_m<-c(2,3) #one row two columns
330 387
#par(mfrow=layout_m)
331 388
    
332 389
##add option for plot title? 
333
png(paste("Figure4__accuracy_",ac_metric,"_prop_month","_",out_prefix,".png", sep=""),
390
png(paste("Figure4_paper_accuracy_",ac_metric,"_prop_month","_",out_prefix,".png", sep=""),
334 391
    height=480*layout_m[1],width=480*layout_m[2])
335 392
p1<-list_plots[[7]]
336 393
p2<-list_plots[[8]]
......
380 437
## Now create boxplots...
381 438
layout_m<-c(2,2) #one row two columns
382 439

  
383
png(paste("Figure_5_boxplot_overtraining_",out_prefix,".png", sep=""),
440
png(paste("Figure5_paper_boxplot_overtraining_",out_prefix,".png", sep=""),
384 441
    height=480*layout_m[1],width=480*layout_m[2])
385 442
#boxplot(diff_kriging_CAI$rmse,diff_gam_CAI$rmse,diff_gwr_CAI$rmse,names=c("kriging_CAI","gam_CAI","gwr_CAI"),
386 443
#        main="Difference between training and testing daily rmse")
......
388 445

  
389 446
#monthly CAI
390 447
boxplot(list_m_diff$kriging_CAI$rmse,list_m_diff$gam_CAI$rmse,list_m_diff$gwr_CAI$rmse,
391
        names=c("kriging_CAI","gam_CAI","gwr_CAI"),
448
        names=c("kriging_CAI","gam_CAI","gwr_CAI"),ylab="RMSE (°C)",xlab="Interpolation Method",
392 449
        main="Difference between training and testing for monhtly rmse")
450
legend("topleft",legend=c("a"),bty="n") #bty="n", don't put box around legend
393 451

  
394 452
#daily CAI
395 453
boxplot(list_diff$kriging_CAI$rmse,list_diff$gam_CAI$rmse,list_diff$gwr_CAI$rmse,
396
        names=c("kriging_CAI","gam_CAI","gwr_CAI"),
454
        names=c("kriging_CAI","gam_CAI","gwr_CAI"),ylab="RMSE (°C)",xlab="Interpolation Method",
397 455
        main="Difference between training and testing daily rmse")
456
legend("topleft",legend=c("b"),bty="n")
457

  
398 458
#monthly fss
399 459
boxplot(list_m_diff$kriging_fss$rmse,list_m_diff$gam_fss$rmse,list_diff$gwr_fss$rmse,
400
        names=c("kriging_fss","gam_fss","gwr_fss"),
460
        names=c("kriging_FSS","gam_FSS","gwr_FSS"),ylab="RMSE (°C)",xlab="Interpolation Method",
401 461
        main="Difference between training and testing for monhtly rmse")
462
legend("topleft",legend=c("c"),bty="n")
402 463

  
403 464
#daily fss
404 465
boxplot(list_diff$kriging_fss$rmse,list_diff$gam_fss$rmse,list_diff$gwr_fss$rmse,
405
        names=c("kriging_fss","gam_fss","gwr_fss"),
466
        names=c("kriging_FSS","gam_FSS","gwr_FSS"),ylab="RMSE (°C)",xlab="Interpolation Method",
406 467
        main="Difference between training and testing daily rmse")
468
legend("topleft",legend=c("d"),bty="n")
407 469

  
408 470
dev.off()
409 471

  
......
427 489

  
428 490
names_layers <-c("mod1 = lat*long","mod2 = lat*long + LST","mod3 = lat*long + elev","mod4 = lat*long + N_w*E_w",
429 491
                 "mod5 = lat*long + elev + DISTOC","mod6 = lat*long + elev + LST","mod7 = lat*long + elev + LST*FOREST")
430
nb_fig<- c("7a","7b","7c")
492
nb_fig<- c("6a","6b","6c")
431 493
list_plots_spt <- vector("list",length=length(lf))
432 494
for (i in 1:length(lf)){
433 495
  pred_temp_s <-stack(lf[[i]])
......
462 524
layout_m<-c(2,4) # works if set to this?? ok set the resolution...
463 525
#layout_m<-c(2*3,4) # works if set to this?? ok set the resolution...
464 526

  
465
png(paste("Figure7_","_spatial_pattern_tmax_prediction_models_gam_levelplot_",date_selected,out_prefix,".png", sep=""),
466
    height=480*layout_m[1],width=480*layout_m[2])
527
png(paste("Figure6_paper_","_spatial_pattern_tmax_prediction_models_gam_levelplot_",date_selected,out_prefix,".png", sep=""),
528
    height=960*layout_m[1],width=960*layout_m[2])
467 529
    #height=480*6,width=480*4)
468 530

  
469 531
p1 <- list_plots_spt[[1]]
......
479 541
#######Figure 7a: Map of transects
480 542

  
481 543
## Transects image location in OR             
482
png(paste("Fig7_elevation_transect_paths_",date_selected,out_prefix,".png", sep=""),
544
png(paste("Figrue7_paper_elevation_transect_paths_",date_selected,out_prefix,".png", sep=""),
483 545
    height=480*layout_m[1],width=480*layout_m[2])
484 546

  
485 547
plot(elev)
......
512 574
layers_names3 <- names(rast_pred3)<-c("mod1","mod3","mod6","elev")
513 575
#pos<-c(1,2) # postions in the layer prection
514 576
#transect_list
515
list_transect2[[1]]<-c(transect_list[1],paste("figure_3_tmax_elevation_transect1_OR_",date_selected,
577
list_transect2[[1]]<-c(transect_list[1],paste("Figure8a_paper_tmax_elevation_transect1_OR_",date_selected,
516 578
                                           paste("mod1_2_6",collapse="_"),out_prefix,sep="_"))
517
list_transect2[[2]]<-c(transect_list[2],paste("figure_3_tmax_elevation_transect2_OR_",date_selected,
579
list_transect2[[2]]<-c(transect_list[2],paste("Figure8b_tmax_elevation_transect2_OR_",date_selected,
518 580
                                           paste("mod1_2_6",collapse="_"),out_prefix,sep="_"))
519
list_transect2[[3]]<-c(transect_list[3],paste("figure_3_tmax_elevation_transect3_OR_",date_selected,
581
list_transect2[[3]]<-c(transect_list[3],paste("Figure8c_paper_tmax_elevation_transect3_OR_",date_selected,
520 582
                                           paste("mod1_2_6",collapse="_"),out_prefix,sep="_"))
521 583

  
522
list_transect3[[1]]<-c(transect_list[1],paste("figure_3_tmax_elevation_transect1_OR_",date_selected,
584
list_transect3[[1]]<-c(transect_list[1],paste("Figure8a_tmax_elevation_transect1_OR_",date_selected,
523 585
                                           paste("mod1_3_6",collapse="_"),out_prefix,sep="_"))
524
list_transect3[[2]]<-c(transect_list[2],paste("figure_3_tmax_elevation_transect2_OR_",date_selected,
586
list_transect3[[2]]<-c(transect_list[2],paste("Figure8b_tmax_elevation_transect2_OR_",date_selected,
525 587
                                           paste("mod1_3_6",collapse="_"),out_prefix,sep="_"))
526
list_transect3[[3]]<-c(transect_list[3],paste("figure_3_tmax_elevation_transect3_OR_",date_selected,
588
list_transect3[[3]]<-c(transect_list[3],paste("Figure8c_tmax_elevation_transect3_OR_",date_selected,
527 589
                                           paste("mod1_3_6",collapse="_"),out_prefix,sep="_"))
528 590

  
529
names(list_transect2)<-c("transect_OR1","transect_OR2","transect_OR3")
530
names(list_transect3)<-c("transect_OR1","transect_OR2","transect_OR3")
591
names(list_transect2)<-c("Oregon Transect 1","Oregon Transect 2","Oregon Transect 3")
592
names(list_transect3)<-c("Oregon Transect 1","Oregon Transect 2","Oregon Transect 3")
531 593

  
532 594
names(rast_pred2)<-layers_names2
533 595
names(rast_pred3)<-layers_names3
534 596

  
535
title_plot2<-paste(names(list_transect2),date_selected,sep=" ")
536
title_plot2<-paste(rep("Oregon transect on ",3), date_selected,sep="")
537
title_plot3<-paste(names(list_transect3),date_selected,sep=" ")
538
title_plot3<-paste(rep("Oregon transect on ",3), date_selected,sep="")
597
title_plot2<-paste(names(list_transect2),"on",date_selected,sep=" ")
598
#title_plot2<-paste(names(list_transect2),"on",date_selected,sep=" ")
599

  
600
#title_plot2<-paste(rep("Oregon transect on ",3), date_selected,sep="")
601
#title_plot3<-paste(names(list_transect3),date_selected,sep=" ")
602
#title_plot3<-paste(rep("Oregon transect on ",3), date_selected,sep="")
539 603

  
540 604
#r_stack<-rast_pred
541 605
#m_layers_sc<-c(3) #elevation in the third layer
......
545 609
#rast_pred2
546 610
#debug(plot_transect_m2)
547 611
trans_data2 <-plot_transect_m2(list_transect2,rast_pred2,title_plot2,disp=FALSE,m_layers_sc)
548
trans_data3 <-plot_transect_m2(list_transect3,rast_pred3,title_plot3,disp=FALSE,m_layers_sc)
612
trans_data3 <-plot_transect_m2(list_transect3,rast_pred3,title_plot2,disp=FALSE,m_layers_sc)
549 613

  
550 614
################################################
551
#Figure 8: Spatial pattern: Image differencing and land cover  
615
#Figure 9: Spatial pattern: Image differencing  
552 616
#Do for january and September...?
553 617

  
554 618
#names_layers <-c("mod1 = lat*long","mod2 = lat*long + LST","mod3 = lat*long + elev","mod4 = lat*long + N_w*E_w",
......
574 638
temp.colors <- colorRampPalette(c('blue', 'white', 'red'))
575 639

  
576 640
layout_m<-c(1,1) #one row two columns
577
png(paste("Figure_9_difference_image_",out_prefix,".png", sep=""),
641
png(paste("Figure9_paper_difference_image_",out_prefix,".png", sep=""),
578 642
    height=480*layout_m[1],width=480*layout_m[2])
579 643
plot(r_stack_diff,col=temp.colors(25))
580 644
#levelplot(r_stack_diff)
581 645
dev.off()
582 646
###
583 647

  
648

  
649
####################################################################
650
#Figure 10: LST and Tmax at stations, elevation and land cover.
651

  
584 652
LC_subset <- c("LC1","LC5","LC6","LC7","LC9","LC11")  
585 653
LC_names <- c("LC1_forest", "LC5_shrub", "LC6_grass", "LC7_crop", "LC9_urban","LC11_barren")
586 654
LC_s <- subset(s_raster,LC_subset)
......
594 662
#r_subset_name <- "elev_s"
595 663
#r_names <- c("elev_s")
596 664
#stat_list_elev <- extract_diff_by_landcover(r_stack_diff,s_raster,LC_subset,LC_names,avl)
597

  
598 665
#write_out_raster_fun(s_raster,out_suffix=out_prefix,out_dir=out_dir,NA_flag_val=-9999,file_format=".rst")
599 666

  
600 667
#show correlation with LST by day over the year, ok writeout s_raster of coveriate??
......
605 672
layout_m<-c(2,3) #one row two columns
606 673
#savePlot(paste("fig6_diff_prediction_tmax_difference_land cover",mf_selected,mc_selected,date_selected,out_prefix,".png", sep="_"), type="png")
607 674

  
608
png(paste("Figure_9_diff_prediction_tmax_difference_land cover,ac_metric","_",out_prefix,".png", sep=""),
675
png(paste("Figure10_paper_diff_prediction_tmax_difference_land cover,ac_metric","_",out_prefix,".png", sep=""),
609 676
      height=480*layout_m[1],width=480*layout_m[2])
610 677
par(mfrow=layout_m)    
611 678
#funciton plot
......
630 697
  points(zones_stat$zones,zones_stat[,6],col="purple",lty="dashed",pch=6)
631 698

  
632 699
  breaks_lab<-zones_stat$zones
700
  #make it slanted...
633 701
  tick_lab<-c("0","1-10","","20-30","","40-50","","60-70","","80-90","90-100") #Not enough space for  
634 702
  #tick_lab<-c("0","10-20","30-40","60-70","80-90","90-100")
635 703
  axis(side=1,las=1,tick=TRUE,
636 704
       at=breaks_lab,labels=tick_lab, cex.axis=1.2,font=2) #reduce number of labels to Jan and June
637 705
  #text(tick_lab, par(\u201cusr\u201d)[3], labels = tick_lab, srt = 45, adj = c(1.1,1.1), xpd = TRUE, cex=.9)
638
  axis(2,cex.axis=1.2,font=2)
706
  axis(2,cex.axis=1.4,font=2)
639 707
  box()
640 708
  legend("topleft",legend=names(zones_stat)[-7], 
641
        cex=1, col=c("black","red","green","blue","darkgreen","purple"),bty="n",
709
        cex=1.4, col=c("black","red","green","blue","darkgreen","purple"),bty="n",
642 710
        lty=1,pch=1:7)
643
  title(paste(title_plots_list[i],sep=""),cex=1.4, font=2)
711
  title(paste(title_plots_list[i],sep=""),cex=1.6, font=2)
644 712
  #title(paste("Prediction tmax difference (",mf_selected,"-",mc_selected,") and land cover ",sep=""),cex=1.4,font=2) 
645 713
}
646 714
dev.off()
647 715

  
648
#### Now elev?
649
#LC1<-mask(LC1,mask_ELEV_SRTM)
650
#  cellStats(LC1,"countNA")        #Check that NA have been assigned to water and areas below 0 m
651
  
652
#LC1_50_m<- LC1>50
653
#LC1_100_m<- LC1>=100
654
#LC1_50_m[LC1_50_m==0]<-NA
655
#LC1_100_m[LC1_100_m==0]<-NA
656
#LC1_50<-LC1_50_m*LC1
657
#LC1_100<-LC1_100_m*LC1
658
#avl<-c(0,500,1,500,1000,2,1000,1500,3,1500,2000,4,2000,4000,5)
659
#rclmat<-matrix(avl,ncol=3,byrow=TRUE)
660
#elev_rec<-reclass(ELEV_SRTM,rclmat)  #Loss of layer names when using reclass
661
  
662
#elev_rec_forest<-elev_rec*LC1_100_m
663
#avg_elev_rec<-zonal(rast_diff,zones=elev_rec,stat="mean",na.rm=TRUE)
664
#std_elev_rec<-zonal(rast_diff,zones=elev_rec,stat="sd",na.rm=TRUE)
665
#avg_elev_rec_forest<-zonal(rast_diff,zones=elev_rec_forest,stat="mean",na.rm=TRUE)
666
#std_elev_rec_forest<-zonal(rast_diff,zones=elev_rec_forest,stat="sd",na.rm=TRUE)
716

  
717
################################################
718
#### Figure 11: Spatial lag profiles  
719
#This figure is generated to show the spatial Moran'I for 10 spatial 
720
#for Jan 1 and Sept 1 in 2010 for all models (1 to 7) and methods
721

  
722
index<-1 #index corresponding to Jan 1 #For now create Moran's I for only one date...
723
lf_moran_list_date1 <-lapply(list_raster_obj_files[c("gam_daily","gam_CAI","gam_fss")],
724
                               FUN=function(x){x<-load_obj(x);x$method_mod_obj[[index]][[y_var_name]]})                           
725
index<-244 #index corresponding to Sept 1 #For now create Moran's I for only one date...
726
lf_moran_list_date2 <-lapply(list_raster_obj_files[c("gam_daily","gam_CAI","gam_fss")],
727
                               FUN=function(x){x<-load_obj(x);x$method_mod_obj[[index]][[y_var_name]]})                           
728

  
729
#date_selected <- "20100901"
730
#date_selected <- "20100101"
731
#methods_names <-c("gam","kriging","gwr")
732
methods_names <-c("gam_daily","gam_CAI","gam_FSS")
733

  
734
names_layers<-methods_names
735
#Subset images to eliminate mod_kr
736
lf1 <- list(lf_moran_list_date1[[1]],lf_moran_list_date1[[2]][1:7],lf_moran_list_date1[[3]][1:7])
737
lf2 <- list(lf_moran_list_date2[[1]],lf_moran_list_date2[[2]][1:7],lf_moran_list_date2[[3]][1:7])
738
names(lf1)<-c("gam_daily","gam_CAI","gam_FSS")
739
names(lf2)<-c("gam_daily","gam_CAI","gam_FSS")
740

  
741
### Now extract Moran's I for a range of lag using a list of images
742

  
743
#set maximum lag range
744
nb_lag <-10
745
#Provide list of raster images:
746
list_lf <- list(lf1,lf2)
747

  
748
list_moran_df1 <- calculate_moranI_profile(list_lf[[1]],nb_lag) #for January 1
749
list_moran_df2 <- calculate_moranI_profile(list_lf[[2]],nb_lag) #for September 1
750
names(list_moran_df1)<-c("gam_daily","gam_CAI","gam_FSS")
751
names(list_moran_df2)<-c("gam_daily","gam_CAI","gam_FSS")
752

  
753
#Run accross two dates...
754
#list_moran lapply(list_lf,FUN=calculate_moranI_profile,nb_lag=nb_lag)
755

  
756
### Prepare to plot lag Moran's I profiles
757

  
758
#generate automatic title for exploration if necessary!!
759
list_title_plot<- list(c("Spatial lag profile on January 1, 2010"),
760
                  c("Spatial lag profile on September 1, 2010"))
761
#name used in the panel!!!
762
names_panel_plot <-c("mod1 = lat*long","mod2 = lat*long + LST","mod3 = lat*long + elev","mod4 = lat*long + N_w*E_w",
763
                 "mod5 = lat*long + elev + DISTOC","mod6 = lat*long + elev + LST","mod7 = lat*long + elev + LST*FOREST")
764
layout_m<-c(2,4) # works if set to this?? ok set the resolution...
765
list_moran_df <- list(list_moran_df1,list_moran_df2)
766
list_param_plot_moranI_profile_fun <- list(list_moran_df,list_title_plot,names_panel_plot,layout_m)
767
names(list_param_plot_moranI_profile_fun) <- c("list_moran_df","list_title_plot","names_panel_plot","layout_m")
768

  
769
#debug(plot_moranI_profile_fun)
770
#p<- plot_moranI_profile_fun(1,list_param=list_param_plot_moranI_profile_fun)
667 771
  
668
## CREATE plots
669
#X11()
670
#plot(avg_elev_rec[,1],avg_elev_rec[,2],type="b",ylim=c(-10,1),
671
#       ylab="",xlab="",axes=FALSE)
672
#mtext("tmax difference between FSS and CAI (degree C)",side=2,cex=1.2,line=3,font=2)
673
#mtext("elevation classes (m)",side=1,cex=1.2,line=3,font=2)
674
#lines(avg_elev_rec_forest[,1],avg_elev_rec_forest[,2],col="green",type="b") #Elevation and 100% forest...
675
#breaks_lab<-avg_elev_rec[,1]
676
#elev_lab<-c("0-500","500-1000","1000-1500","1500-2000","2000-4000")
677
#axis(side=1,las=1,
678
#       at=breaks_lab,labels=elev_lab, cex=1.5,font=2) #reduce number of labels to Jan and June
679
#axis(2,cex.axis=1.2,font=2)
680
#legend("bottomleft",legend=c("Elevation", "elev_forest"), 
681
#         cex=1, lwd=1.3,col=c("black","green"),bty="n",
682
#         lty=1)
683
#box()
684
#title(paste("Prediction tmax difference (",mf_selected,"-",mc_selected,") and elevation ",sep=""),cex=1.4,font=2)
685
#savePlot(paste("fig7_diff_prediction_tmax_difference_elevation",mf_selected,mc_selected,date_selected,out_prefix,".png", sep="_"), type="png")
686
#dev.off()
772
list_moranI_plots <- lapply(1:2,FUN=plot_moranI_profile_fun,list_param=list_param_plot_moranI_profile_fun)
773

  
774

  
775
#This function uses list moran_df object from calculate_moranI_profile function!!
776

  
777
#layout_m<-c(2,4) # works if set to this?? ok set the resolution...
778
#layout_m<-c(2*3,4) # works if set to this?? ok set the resolution...
779

  
780
png(paste("Figure11_paper_spatial_correlogram_prediction_models_levelplot_",out_prefix,".png", sep=""),
781
    height=480*layout_m[1],width=480*layout_m[2])
782
    #height=3*480*layout_m[1],width=2*480*layout_m[2])
783
    #height=480*6,width=480*4)
784
#png(paste("Figure11_paper_spatial_correlogram_prediction_models_levelplot_",out_prefix,".png", sep=""),
785
#    height=480,width=480)
786
    #height=480*6,width=480*4)
787

  
788
p1 <- list_moranI_plots[[1]]
789
p2 <- list_moranI_plots[[2]]
790

  
791
grid.arrange(p1,p2,ncol=1)
792
dev.off()
793

  
794
################################################
795
#Figure 12: Monthly climatology, Daily deviation and bias
796

  
797
lf_delta_gam_fss <-extract_list_from_list_obj(load_obj(list_raster_obj_files[["gam_fss"]])$method_mod_obj,"delta") #getting objet
798

  
799
test01 <-stack(lf_delta_gam_fss[[287]]) #Ot.14
800

  
801
test0 <-stack(lf_delta_gam_fss[[288]]) #Ot.15
802
test1 <-stack(lf_delta_gam_fss[[289]]) #Ot.16
803
test2<-stack(lf_delta_gam_fss[[290]]) #Ot.17
804

  
805
plot(test01)
806

  
807
plot(test0)
808
plot(test1)
809
plot(test2)
810

  
687 811

  
688 812
################################################
689
#Figure 9: Tmax and LST averagees
813
#Figure 13: Tmax and LST averagees
690 814

  
691 815
names_tmp<-c("mm_01","mm_02","mm_03","mm_04","mm_05","mm_06","mm_07","mm_08","mm_09","mm_10","mm_11","mm_12")
692 816
LST_s<-subset(s_raster,names_tmp)
693 817
names_tmp<-c("nobs_01","nobs_02","nobs_03","nobs_04","nobs_05","nobs_06","nobs_07","nobs_08",
694 818
             "nobs_09","nobs_10","nobs_11","nobs_12")
695 819
LST_nobs<-subset(s_raster,names_tmp)
696

  
697
#LST_nobs<-mask(LST_nobs,LC_mask,filename="test2.tif")
698
#LST_s<-mask(LST_s,LC_mask,filename="test3.tif")
699
#c("Jan","Feb")
700 820
plot(LST_s)
701 821
plot(LST_nobs)
822
y_range_nobs <- c(0,385)
823
y_range_avg <- c(-15,60)
824
var_name <- "LST"
825

  
826
#It works for any variable with stack of monthly values...!!
827
LST_stat <- plot_mean_nobs_r_stack_by_month(var_s=LST_s,var_nobs=LST_nobs,y_range_nobs,y_range_avg,var_name,out_prefix)
702 828

  
703
#Map 5: LST and TMax
704

  
705
#note differnces in patternin agricultural areas and 
706
min_values<-cellStats(LST_s,"min")
707
max_values<-cellStats(LST_s,"max")
708
mean_values<-cellStats(LST_s,"mean")
709
sd_values<-cellStats(LST_s,"sd")
710
#median_values<-cellStats(molst,"median") Does not extist
711
statistics_LST_s<-cbind(min_values,max_values,mean_values,sd_values) #This shows that some values are extremes...especially in October
712
LST_stat_data<-as.data.frame(statistics_LST_s)
713
rownames(LST_stat_data) <- month.abb
714
names(LST_stat_data)<-c("min","max","mean","sd")
715
# Statistics for number of valid observation stack
716
min_values<-cellStats(LST_nobs,"min")
717
max_values<-cellStats(LST_nobs,"max")
718
mean_values<-cellStats(LST_nobs,"mean")
719
sd_values<-cellStats(LST_nobs,"sd")
720
LST_nobs_stat_data<-cbind(min_values,max_values,mean_values,sd_values) #This shows that some values are extremes...especially in October
721
LST_nobs_stat_data <- as.data.frame(LST_nobs_stat_data)
722
rownames(LST_nobs_stat_data) <- month.abb
723

  
724
#X11(width=12,height=12)
725
#Plot statiscs (mean,min,max) for monthly LST images
726
plot(1:12,LST_stat_data$mean,type="b",ylim=c(-15,60),col="black",xlab="month",ylab="tmax (degree C)")
727
lines(1:12,LST_stat_data$min,type="b",col="blue")
728
lines(1:12,LST_stat_data$max,type="b",col="red")
729
text(1:12,LST_stat_data$mean,rownames(LST_stat_data),cex=1,pos=2)
730

  
731
legend("topleft",legend=c("min","mean","max"), cex=1.5, col=c("blue","black","red"),
732
       lty=1)
733
title(paste("LST statistics for Oregon", "2010",sep=" "))
734
#savePlot("lst_statistics_OR.png",type="png")
735

  
736
#Plot number of valid observations for LST
737
plot(1:12,LST_nobs_stat_data$mean,type="b",ylim=c(0,365),col="black",xlab="month",ylab="tmax (degree C)")
738
lines(1:12,LST_nobs_stat_data$min,type="b",col="blue")
739
lines(1:12,LST_nobs_stat_data$max,type="b",col="red")
740
text(1:12,LST_nobs_stat_data$mean,rownames(LST_stat_data),cex=1,pos=2)
741

  
742
legend("topleft",legend=c("min","mean","max"), cex=1.5, col=c("blue","black","red"),
743
       lty=1)
744
title(paste("LST number of valid observations for Oregon", "2010",sep=" "))
745
#savePlot("lst_nobs_OR.png",type="png")
746

  
747
#met_stations_obj <- load_obj(file.path(in_dir1,met_obj_file_1))
748
#met_stations_obj$monthly_query_ghcn_data
749
#
829
##################################################
830
####### Figure 13 : disussion  
750 831

  
832
##################################################
833
####### Figure 13a: spatial variation -MoranI and  std  
834

  
835
#Use data from accuraccy with no monthly hold out
751 836
raster_obj <- load_obj(list_raster_obj_files$gam_CAI)
752 837
data_month <- (raster_obj$method_mod_obj[[1]]$data_month_s)
753 838
dates<-unique(raster_obj$tb_diagnostic_v$date)
839
#Extract monthly data frame used in fitting 
840
list_data_month_s <-extract_list_from_list_obj(raster_obj$validation_mod_month_obj,"data_s") 
841
year_nbs <- sapply(list_data_month_s,FUN=length) #number of observations per month
842
#Convert sp data.frame and combined them in one unique df, see function define earlier
843
data_month_all <-convert_spdf_to_df_from_list(list_data_month_s) #long rownames
844
#LSTD_bias_avgm<-aggregate(LSTD_bias~month,data=data_month_all,mean)
845
#LSTD_bias_sdm<-aggregate(LSTD_bias~month,data=data_month_all,sd)
846
LST_avgm<-aggregate(LST~month,data=data_month_all,mean)
847
TMax_avgm<-aggregate(TMax~month,data=data_month_all,mean)
848

  
849
#plot(LST_avgm,type="b")
850
#lines(TMax_avgm,col="blue",add=T)
851
plot(data_month_all$month,data_month_all$LST, xlab="month",ylab="LST at station")
852
title(paste("Monthly LST at stations in Oregon", "2010",sep=" "))
853
png(paste("LST_at_stations_per_month_",out_prefix,".png", sep=""), type="png")
854
                  
855

  
856
statistics_LST_s<- LST_stat$avg #extracted earlier!!!
857

  
858
png(paste("Monthly_mean_TMax_LST_at_Stations_",out_prefix,".png", sep=""), type="png")
859

  
860
plot(TMax_avgm,type="b",ylim=c(0,35),col="red",xlab="month",ylab="tmax (degree C)")
861
lines(1:12,LST_avgm$LST,type="b",col="blue")
862
#lines(1:12,statistics_LST_s$mean,type="b",col="darkgreen")
863
text(TMax_avgm[,1],TMax_avgm[,2],rownames(statistics_LST_s),cex=0.8,pos=2)
864
                  
865
legend("topleft",legend=c("TMax","LST"), cex=1, col=c("red","blue"),
866
              lty=1,bty="n")
867
title(paste("Monthly mean tmax and LST at stations in Oregon", "2010",sep=" "))           
868

  
869
##################################################
870
####### Figure 13b: spatial variation -MoranI and  std  
871
#LST for area with FOrest 50> and grass >50
754 872

  
873
#Account for forest
874
data_month_all_forest<-data_month_all[data_month_all$LC1>=50,]
875
data_month_all_grass<-data_month_all[data_month_all$LC6>=50,]
876
#data_month_all_grass<-data_month_all[data_month_all$LC6>=10,]
877
data_month_all_urban<-data_month_all[data_month_all$LC9>=50,]
755 878

  
756
#########################
757
#Figure
758
#LST for area with FOrest 50> and grass >50
879
LST_avgm_forest<-aggregate(LST~month,data=data_month_all_forest,mean)
880
LST_avgm_grass<-aggregate(LST~month,data=data_month_all_grass,mean)
881
LST_avgm_urban<-aggregate(LST~month,data=data_month_all_urban,mean)
882

  
883
plot(TMax_avgm,type="b",ylim=c(-7,42))
884
lines(LST_avgm$LST,col="blue",add=T)
885
lines(LST_avgm_forest,col="green",add=T)
886
lines(LST_avgm_grass,col="red",add=T)
887
lines(LST_avgm_urban,col="pink",add=T)
888
legend("topleft",legend=c("TMax","LST","LST_forest","LST_grass","LST_urban"), cex=0.8, col=c("black","blue","green","red","pink"),
889
       lty=1,bty="n")
890
title(paste("Monthly average tmax for stations in Oregon ", "2010",sep=""))
759 891

  
760
########################
761
#Figure 
892
savePlot(paste("Temporal_profile_res",id[i],out_prefix,".png", sep=""), type="png")  
893

  
894
##################################################
895
####### Figure 13c: spatial variation -MoranI and  std  
762 896

  
763 897
#get the  list of all prediction for a specific method
764 898
#list_lf <-extract_list_from_list_obj(load_obj(list_raster_obj_files$gam_CAI)$method_mod_obj,"dailyTmax") #getting objet
765 899
list_lf_gam_CAI <-extract_list_from_list_obj(load_obj(list_raster_obj_files[["gam_CAI"]])$method_mod_obj,"dailyTmax") #getting objet
766 900
list_lf_gam_fss <-extract_list_from_list_obj(load_obj(list_raster_obj_files[["gam_fss"]])$method_mod_obj,"dailyTmax") #getting objet
767 901

  
768

  
769

  
770 902
#nb_lag <- 10
771 903
list_filters<-lapply(1:nb_lag,FUN=autocor_filter_fun,f_type="queen") #generate lag 10 filters
772 904
list_param_stat_moran <- list(filter=list_filters[[10]],lf_list=list_lf_gam_CAI)
773 905
tt <- stat_moran_std_raster_fun(1,list_param=list_param_stat_moran)
774 906
tt <- mclapply(1:11, list_param=list_param_stat_moran, FUN=stat_moran_std_raster_fun,mc.preschedule=FALSE,mc.cores = 11) #This is the end bracket from mclapply(...) statement
775 907
list_param_stat_moran <- list(filter=list_filters[[10]],lf_list=list_lf_gam_fss)
776
tt_fss <- mclapply(1:365, list_param=list_param_stat_moran, FUN=stat_moran_std_raster_fun,mc.preschedule=FALSE,mc.cores = 11) #This is the end bracket from mclapply(...) statement
777
save(tt_fss,file=file.path(out_dir,"moran_std_tt_fss.RData"))
778 908

  
779
tx<- do.call(rbind,tt_fss)
909
#Takes 1 hour to get the average moran's I for the whole year so load moran_std_tt_fss2.RData
910
#tt_fss <- mclapply(1:365, list_param=list_param_stat_moran, FUN=stat_moran_std_raster_fun,mc.preschedule=FALSE,mc.cores = 11) #This is the end bracket from mclapply(...) statement
911
#save(tt_fss,file=file.path(out_dir,"moran_std_tt_fss.RData"))
912
#tx<- do.call(rbind,tt_fss)
780 913

  
781
mo<-as.integer(strftime(dates, "%m"))          # current month of the date being processed
914
#dates <-strptime(dates, "%Y%m%d")   # interpolation date being processed
915
#mo<-as.integer(strftime(dates, "%m"))          # current month of the date being processed
782 916

  
917
dates<-unique(raster_obj$tb_diagnostic_v$date)
783 918
tt_fss2 <- load_obj("moran_std_tt_fss2.RData")
784 919
tx<- do.call(rbind,tt_fss2)
785 920
dates_proc<-strptime(dates, "%Y%m%d")   # interpolation date being processed
......
807 942
plot(x$month,x$std,type="b",col="red",axes=F)
808 943
  #axis(4,xlab="",ylab="elevation(m)")  
809 944
axis(4,cex=1.2)
945

  
946
##################################################
947
####### Figure 13d: correlation between elevation and LST
948

  
949
LST_s <- subset(s_raster,c("mm_01","mm_02","mm_03","mm_04","mm_05","mm_06","mm_07","mm_08","mm_09","mm_10","mm_11","mm_12"))
950
LST1 <- subset(LST_s,1)
951

  
952
test<-stack(LST_s,elev)
953
names(test) <- c("mm_01","mm_02","mm_03","mm_04","mm_05","mm_06","mm_07","mm_08","mm_09","mm_10","mm_11","mm_12","elev")
954
x_cor <-layerStats(test,stat="pearson",na.rm=T)
955
corr(values(LST1),values(elev))
956
plot(1:12,x_cor[[1]][1:12,13],type="b",ylab="corelation",xlab="month",main="correlation between elevation and LST")
957

  
810 958
#now summarize by month...
811 959

  
812 960
#plot(data_month$TMax,add=TRUE)
......
823 971
#1- (SSRES_m/SST)
824 972
#1- (SSRES_dev/SST)
825 973

  
826
##############
827

  
828
lf_delta_gam_fss <-extract_list_from_list_obj(load_obj(list_raster_obj_files[["gam_fss"]])$method_mod_obj,"delta") #getting objet
829

  
830
test01 <-stack(lf_delta_gam_fss[[287]]) #Ot.14
831

  
832
test0 <-stack(lf_delta_gam_fss[[288]]) #Ot.15
833
test1 <-stack(lf_delta_gam_fss[[289]]) #Ot.16
834
test2<-stack(lf_delta_gam_fss[[290]]) #Ot.17
835

  
836
plot(test01)
837

  
838
plot(test0)
839
plot(test1)
840
plot(test2)
841

  
842
################################################
843
#Figure 10: Spatial lag profiles and stations data  
844

  
845
index<-244 #index corresponding to Sept 1 #For now create Moran's I for only one date...
846
index<-1
847
lf_moran_list<-lapply(list_raster_obj_files[c("gam_daily","gam_CAI","gam_fss")],
848
                               FUN=function(x){x<-load_obj(x);x$method_mod_obj[[index]][[y_var_name]]})                           
849
#lf1 <- raster_prediction_obj_1$method_mod_obj[[index]][[y_var_name]] #select relevant raster images for the given dates
850

  
851
date_selected <- "20100901"
852
date_selected <- "20100101"
853
#methods_names <-c("gam","kriging","gwr")
854
methods_names <-c("gam_daily","gam_CAI","gam_FSS")
855

  
856
names_layers<-methods_names
857
#lf <- (list(lf1,lf4[1:7],lf7[1:7]))
858
lf<-list(lf_moran_list[[1]],lf_moran_list[[2]][1:7],lf_moran_list[[3]][1:7])
859

  
860
names_layers <-c("mod1 = lat*long","mod2 = lat*long + LST","mod3 = lat*long + elev","mod4 = lat*long + N_w*E_w",
861
                 "mod5 = lat*long + elev + DISTOC","mod6 = lat*long + elev + LST","mod7 = lat*long + elev + LST*FOREST")
862 974

  
863 975

  
864
nb_lag <-10
865
#lf
866

  
867
calculate_moranI_profile <- function(nb_lag,lf){
868
  list_filters<-lapply(1:nb_lag,FUN=autocor_filter_fun,f_type="queen") #generate lag 10 filters
869
  #moran_list <- lapply(list_filters,FUN=Moran,x=r)
870
  list_moran_df <- vector("list",length=length(lf))
871
  for (j in 1:length(lf)){
872
    r_stack <- stack(lf[[j]])
873
    list_param_moran <- list(list_filters=list_filters,r_stack=r_stack) #prepare parameters list for function
874
    #moran_r <-moran_multiple_fun(1,list_param=list_param_moran)
875
    nlayers(r_stack) 
876
    moran_I_df <-mclapply(1:nlayers(r_stack), list_param=list_param_moran, FUN=moran_multiple_fun,mc.preschedule=FALSE,mc.cores = 10) #This is the end bracket from mclapply(...) statement
877

  
878
    moran_df <- do.call(cbind,moran_I_df) #bind Moran's I value 10*nlayers data.frame
879
    moran_df$lag <-1:nrow(moran_df)
976
#### Now elev?
977
#LC1<-mask(LC1,mask_ELEV_SRTM)
978
#  cellStats(LC1,"countNA")        #Check that NA have been assigned to water and areas below 0 m
880 979
  
881
    list_moran_df[[j]] <- moran_df
882
  }
883
}
884

  
885
names(list_moran_df)<-c("gam_daily","gam_CAI","gam_fss")
886
list_dd <- vector("list",length=length(list_moran_df))
887

  
888
for(j in 1:length(lf)){
889
  method_name <- names(list_moran_df)[j]
890
  mydata <- list_moran_df[[j]]
891
  dd <- do.call(make.groups, mydata[,-ncol(mydata)]) 
892
  dd$lag <- mydata$lag
893
  dd$method_v <- method_name
894
  list_dd[[j]] <- dd
895
}
896

  
897
dd_combined<- do.call(rbind,list_dd)
898

  
899
layout_m<-c(2,4) #one row two columns
900

  
901
png(paste("Figure_9_spatial_correlogram_tmax_prediction_models_gam_levelplot_",date_selected,out_prefix,".png", sep=""),
902
    height=480*layout_m[1],width=480*layout_m[2])
903
par(mfrow=layout_m)
904
p<-xyplot(data ~ lag | which , data=dd_combined,group=method_v,type="b", as.table=TRUE,
905
          pch=1:3,auto.key=list(columns=3,cex=1.5,font=2),
906
          par.settings = list(
907
          superpose.symbol = list(pch=1:3,col=1:3,pch.cex=1.4),
908
          axis.text = list(font = 2, cex = 1.3),layout=layout_m,
909
                              par.main.text=list(font=2,cex=2),strip.background=list(col="white")),
910
                              par.strip.text=list(font=2,cex=1.5),
911
          strip=strip.custom(factor.levels=names_layers),
912
          xlab=list(label="Spatial lag neighbor", cex=2,font=2),
913
          ylab=list(label="Moran's I", cex=2, font=2))
914

  
915
#Use as.table to reverse order of panel from top to bottom.
916

  
917
#p<-xyplot(data ~ lag | which , dd_combined,group=method_v,type="b",
918
#          par.settings = list(axis.text = list(font = 2, cex = 1.3),layout=layout_m,
919
#                              par.main.text=list(font=2,cex=2),strip.background=list(col="white")),
920
#                              par.strip.text=list(font=2,cex=1.5),
921
#          strip=strip.custom(factor.levels=names_layers),
922
#          xlab=list(label="Spatial lag neighbor", cex=2,font=2),
923
#          ylab=list(label="Moran's I", cex=2, font=2))
924

  
925
print(p)
926

  
927
dev.off()
980
#LC1_50_m<- LC1>50
981
#LC1_100_m<- LC1>=100
982
#LC1_50_m[LC1_50_m==0]<-NA
983
#LC1_100_m[LC1_100_m==0]<-NA
984
#LC1_50<-LC1_50_m*LC1
985
#LC1_100<-LC1_100_m*LC1
986
#avl<-c(0,500,1,500,1000,2,1000,1500,3,1500,2000,4,2000,4000,5)
987
#rclmat<-matrix(avl,ncol=3,byrow=TRUE)
988
#elev_rec<-reclass(ELEV_SRTM,rclmat)  #Loss of layer names when using reclass
989
  
990
#elev_rec_forest<-elev_rec*LC1_100_m
991
#avg_elev_rec<-zonal(rast_diff,zones=elev_rec,stat="mean",na.rm=TRUE)
992
#std_elev_rec<-zonal(rast_diff,zones=elev_rec,stat="sd",na.rm=TRUE)
993
#avg_elev_rec_forest<-zonal(rast_diff,zones=elev_rec_forest,stat="mean",na.rm=TRUE)
994
#std_elev_rec_forest<-zonal(rast_diff,zones=elev_rec_forest,stat="sd",na.rm=TRUE)
995
  
996
## CREATE plots
997
#X11()
998
#plot(avg_elev_rec[,1],avg_elev_rec[,2],type="b",ylim=c(-10,1),
999
#       ylab="",xlab="",axes=FALSE)
1000
#mtext("tmax difference between FSS and CAI (degree C)",side=2,cex=1.2,line=3,font=2)
1001
#mtext("elevation classes (m)",side=1,cex=1.2,line=3,font=2)
1002
#lines(avg_elev_rec_forest[,1],avg_elev_rec_forest[,2],col="green",type="b") #Elevation and 100% forest...
1003
#breaks_lab<-avg_elev_rec[,1]
1004
#elev_lab<-c("0-500","500-1000","1000-1500","1500-2000","2000-4000")
1005
#axis(side=1,las=1,
1006
#       at=breaks_lab,labels=elev_lab, cex=1.5,font=2) #reduce number of labels to Jan and June
1007
#axis(2,cex.axis=1.2,font=2)
1008
#legend("bottomleft",legend=c("Elevation", "elev_forest"), 
1009
#         cex=1, lwd=1.3,col=c("black","green"),bty="n",
1010
#         lty=1)
1011
#box()
1012
#title(paste("Prediction tmax difference (",mf_selected,"-",mc_selected,") and elevation ",sep=""),cex=1.4,font=2)
1013
#savePlot(paste("fig7_diff_prediction_tmax_difference_elevation",mf_selected,mc_selected,date_selected,out_prefix,".png", sep="_"), type="png")
1014
#dev.off()
928 1015

  
929
LST_s <- subset(s_raster,c("mm_01","mm_02","mm_03","mm_04","mm_05","mm_06","mm_07","mm_08","mm_09","mm_10","mm_11","mm_12"))
930
LST1 <- subset(LST_s,1)
931
test<-stack(LST_s,elev)
932
names(test) <- c("mm_01","mm_02","mm_03","mm_04","mm_05","mm_06","mm_07","mm_08","mm_09","mm_10","mm_11","mm_12","elev")
933
x_cor <-layerStats(test,stat="pearson",na.rm=T)
934
corr(values(LST1),values(elev))
935
plot(1:12,x_cor[[1]][1:12,13],type="b",ylab="corelation",xlab="month",main="correlation between elevation and LST")
936 1016

  
937 1017
###################### END OF SCRIPT #######################
938 1018

  
climate/research/oregon/interpolation/multi_timescales_paper_interpolation_functions.R
5 5
#Functions used in the production of figures and data for the multi timescale paper are recorded.
6 6
#AUTHOR: Benoit Parmentier                                                                      #
7 7
#DATE CREATED: 11/25/2013            
8
#DATE MODIFIED: 12/16/2013            
8
#DATE MODIFIED: 12/23/2013            
9 9
#Version: 1
10 10
#PROJECT: Environmental Layers project                                       #
11 11
#################################################################################################
......
32 32

  
33 33
function_analyses_paper <-"multi_timescales_paper_interpolation_functions_12092013.R"
34 34

  
35
calc_stat_by_month_tb <-function(names_mod,tb,month_holdout=F){
36
  #function
37
  add_month_tag<-function(tb){
38
  date<-strptime(tb$date, "%Y%m%d")   # interpolation date being processed
39
  month<-strftime(date, "%m")          # current month of the date being processed
40
  }
41
  
42
  #begin
43
  
44
  tb$month<-add_month_tag(tb)
45
  if(month_holdout==F){
46
    t<-melt(subset(tb,pred_mod==names_mod),
47
          measure=c("mae","rmse","r","me","m50"), 
48
          id=c("pred_mod","prop","month"),
49
          na.rm=T)
50
    avg_tb<-cast(t,pred_mod+prop+month~variable,mean)
51
    sd_tb<-cast(t,pred_mod+prop+month~variable,sd)
52
    n_tb<-cast(t,pred_mod+prop+month~variable,length)
53

  
54
  }
55
  if(month_holdout==T){
56
    t<-melt(subset(tb,pred_mod==names_mod),
57
          measure=c("mae","rmse","r","me","m50"), 
58
          id=c("pred_mod","prop_month","month"),
59
          na.rm=T)
60
    avg_tb<-cast(t,pred_mod+prop_month+month~variable,mean)
61
    sd_tb<-cast(t,pred_mod+prop_month+month~variable,sd)
62
    n_tb<-cast(t,pred_mod+prop_month+month~variable,length)
63

  
64
  }
65
  #n_NA<-cast(t,dst_cat1~variable,is.na)
66
  
67
  #### prepare returning object
68
  prop_obj<-list(tb,avg_tb,sd_tb,n_tb)
69
  names(prop_obj) <-c("tb","avg_tb","sd_tb","n_tb")
70
  
71
  return(prop_obj)
72
}
73

  
35 74
plot_transect_m2<-function(list_trans,r_stack,title_plot,disp=FALSE,m_layers){
36 75
  #This function creates plot of transects for stack of raster images.
37 76
  #Arguments:
......
184 223
  return(moran_v)
185 224
}
186 225

  
226
#Extract moran's I profile from list of images...the list may contain sublist!!! e.g. for diffeferent
227
#methods in interpolation
228
calculate_moranI_profile <- function(lf,nb_lag){
229
  list_filters<-lapply(1:nb_lag,FUN=autocor_filter_fun,f_type="queen") #generate lag 10 filters
230
  #moran_list <- lapply(list_filters,FUN=Moran,x=r)
231
  list_moran_df <- vector("list",length=length(lf))
232
  for (j in 1:length(lf)){
233
    r_stack <- stack(lf[[j]])
234
    list_param_moran <- list(list_filters=list_filters,r_stack=r_stack) #prepare parameters list for function
235
    #moran_r <-moran_multiple_fun(1,list_param=list_param_moran)
236
    nlayers(r_stack) 
237
    moran_I_df <-mclapply(1:nlayers(r_stack), list_param=list_param_moran, FUN=moran_multiple_fun,mc.preschedule=FALSE,mc.cores = 10) #This is the end bracket from mclapply(...) statement
238

  
239
    moran_df <- do.call(cbind,moran_I_df) #bind Moran's I value 10*nlayers data.frame
240
    moran_df$lag <-1:nrow(moran_df)
241
  
242
    list_moran_df[[j]] <- moran_df
243
  }
244
  names(list_moran_df) <- names(lf)
245
  return(list_moran_df)
246
}
247

  
248
plot_moranI_profile_fun <- function(i,list_param){
249
  #extract relevant parameters
250
  list_moran_df <- list_param$list_moran_df[[i]]
251
  title_plot <- list_param$list_title_plot[[i]]
252
  layout_m <- list_param$layout_m
253
  names_panel_plot <- list_param$names_panel_plot
254
  
255
  #date_selected <- list_param$list_date_selected[i]
256
  list_dd <- vector("list",length=length(list_moran_df))
257

  
258
  ## Reorganize data to use xyplot
259
  for(j in 1:length(list_moran_df)){
260
    method_name <- names(list_moran_df)[j]
261
    mydata <- list_moran_df[[j]]
262
    dd <- do.call(make.groups, mydata[,-ncol(mydata)]) 
263
    dd$lag <- mydata$lag
264
    dd$method_v <- method_name
265
    list_dd[[j]] <- dd
266
  }
267
  dd_combined<- do.call(rbind,list_dd)
268
  
269
  ## Set up plot
270
  #layout_m<-c(2,4) #one row two columns
271
  
272
  png(paste("Spatial_correlogram_prediction_models_levelplot_",i,"_",out_prefix,".png", sep=""),
273
      height=480*layout_m[1],width=480*layout_m[2])
274
  par(mfrow=layout_m)
275
  p<-xyplot(data ~ lag | which , data=dd_combined,group=method_v,type="b", as.table=TRUE,
276
            pch=1:3,auto.key=list(columns=3,cex=1.5,font=2),
277
            main=title_plot,
278
            par.settings = list(
279
              superpose.symbol = list(pch=1:3,col=1:3,pch.cex=1.4),
280
              axis.text = list(font = 2, cex = 1.3),layout=layout_m,
281
              par.main.text=list(font=2,cex=2),strip.background=list(col="white")),
282
            par.strip.text=list(font=2,cex=1.5),
283
            strip=strip.custom(factor.levels=names_panel_plot),
284
            xlab=list(label="Spatial lag neighbor", cex=2,font=2),
285
            ylab=list(label="Moran's I", cex=2, font=2))
286
  
287
  print(p)
288
  
289
  dev.off()
290
  
291
  return(p)
292
}
293

  
294
### Calulate Moran's I and std for every image in a year for a specific model...
187 295
#Modfiy...to allow any stat: min, max, mean,sd etc.
188 296
stat_moran_std_raster_fun<-function(i,list_param){
189 297
  f <-list_param$filter
......
210 318
  return(dat_var_stat)
211 319
}
212 320

  
213
plot_accuracy_by_holdout_fun <-function(list_tb,ac_metric){
321
#Add option for plotname
322
plot_accuracy_by_holdout_fun <-function(list_tb,ac_metric,plot_names,names_mod){
214 323
  #
215 324
  list_plots <- vector("list",length=length(list_tb))
216 325
  for (i in 1:length(list_tb)){
217 326
    #i <- i+1
218 327
    tb <-list_tb[[i]]
219
    plot_name <- names(list_tb)[i]
328
    tb_name <- names(list_tb)[i]
329
    plot_name <- plot_names[i]
220 330
    pat_str <- "tb_m"
221
    if(substr(plot_name,start=1,stop=4)== pat_str){
331
    if(substr(tb_name,start=1,stop=4)== pat_str){
222 332
      names_id <- c("pred_mod","prop")
223 333
      plot_formula <- paste(ac_metric,"~prop",sep="",collapse="") 
224 334
    }else{
225 335
      names_id <- c("pred_mod","prop_month")
226 336
      plot_formula <- paste(ac_metric,"~prop_month",collapse="")
227 337
    }
228
    names_mod <-unique(tb$pred_mod)
338
    names_mod <- unique(tb$pred_mod)#should take a closer look later...
339
    #calc_stat_prop_tb_diagnostic <-function(names_mod,names_id,tb){
340

  
229 341
    prop_obj <- calc_stat_prop_tb_diagnostic(names_mod,names_id,tb)
230 342
    avg_tb <- prop_obj$avg_tb
231
  
343
    
344
    avg_tb <- subset(avg_tb,pred_mod!="mod_kr") #removes mod_kr
232 345
    layout_m<-c(1,1) #one row two columns
233 346
    par(mfrow=layout_m)
234 347
    
235 348
    #add option for plot title?
236
    png(paste("Figure__accuracy_",ac_metric,"_prop_month_",plot_name,"_",out_prefix,".png", sep=""),
349
    png(paste("Accuracy_",ac_metric,"_prop_month_",plot_name,"_",out_prefix,".png", sep=""),
237 350
      height=480*layout_m[1],width=480*layout_m[2])
238
    
351
    #correct label:ylabel: monthly holdout proprotion
352
    #ylabel: RMSE (degred C)
239 353
    p <- xyplot(as.formula(plot_formula),group=pred_mod,type="b",
240 354
          data=avg_tb,
241
          main=paste(ac_metric,plot_name,sep=" "),
355
          main=plot_name,
242 356
          pch=1:length(avg_tb$pred_mod),
243 357
          par.settings=list(superpose.symbol = list(
244 358
          pch=1:length(avg_tb$pred_mod))),
245
          auto.key=list(columns=5))
359
          auto.key=list(columns=5),
360
          xlab="Monthly Hold out proportion",
361
          ylab="RMSE (°C)")
362

  
246 363
    print(p)
247 364
  
365
    
248 366
    dev.off()
249 367
    list_plots[[i]] <- p
250 368
  }
......
334 452
  return(unlist(list_raster_name))
335 453
}
336 454

  
455

  
456
#computer averages and nobs per month for a givenvarialbe
457
plot_mean_nobs_r_stack_by_month <-function(var_s,var_nobs,y_range_nobs,y_range_avg,var_name,out_prefix){
458
  
459
  LST_s <- var_s 
460
  LST_nobs <- var_nobs
461

  
462
  #make this more efficient...
463
  min_values<-cellStats(LST_s,"min")
464
  max_values<-cellStats(LST_s,"max")
465
  mean_values<-cellStats(LST_s,"mean")
466
  sd_values<-cellStats(LST_s,"sd")
467
  #median_values<-cellStats(molst,"median") Does not extist
468
  statistics_LST_s<-cbind(min_values,max_values,mean_values,sd_values) #This shows that some values are extremes...especially in October
469
  LST_stat_data<-as.data.frame(statistics_LST_s)
470
  rownames(LST_stat_data) <- month.abb
471
  names(LST_stat_data)<-c("min","max","mean","sd")
472
  # Statistics for number of valid observation stack
473
  min_values<-cellStats(LST_nobs,"min")
474
  max_values<-cellStats(LST_nobs,"max")
475
  mean_values<-cellStats(LST_nobs,"mean")
476
  sd_values<-cellStats(LST_nobs,"sd")
477
  LST_nobs_stat_data<-cbind(min_values,max_values,mean_values,sd_values) #This shows that some values are extremes...especially in October
478
  LST_nobs_stat_data <- as.data.frame(LST_nobs_stat_data)
479
  rownames(LST_nobs_stat_data) <- month.abb
480
  
481
  layout_m <- c(1,1)
482
  png(paste(var_name,"_","mean_by_month","_",out_prefix,".png", sep=""),
483
    height=480*layout_m[1],width=480*layout_m[2])
484
  par(mfrow=layout_m)    
485

  
486
  plot(1:12,LST_stat_data$mean,type="b",ylim=y_range_avg,col="black",xlab="month",ylab="tmax (degree C)")
487
  lines(1:12,LST_stat_data$min,type="b",col="blue")
488
  lines(1:12,LST_stat_data$max,type="b",col="red")
489
  text(1:12,LST_stat_data$mean,rownames(LST_stat_data),cex=1,pos=2)
490
  
491
  legend("topleft",legend=c("min","mean","max"), cex=1.5, col=c("blue","black","red"),
492
         lty=1,bty="n")
493
  
494
  title(paste(var_name,"statistics for Oregon", "2010",sep=" "))
495
  #title(paste("LST statistics for Oregon", "2010",sep=" "))
496
  dev.off()
497
  
498
  layout_m <- c(1,1)
499
  png(paste(var_name,"_","avg_by_month","_",out_prefix,".png", sep=""),
500
    height=480*layout_m[1],width=480*layout_m[2])
501
  par(mfrow=layout_m)    
502

  
503
  #Plot number of valid observations for LST
504
  plot(1:12,LST_nobs_stat_data$mean,type="b",ylim=y_range_nobs,col="black",xlab="month",ylab="tmax (degree C)")
505
  lines(1:12,LST_nobs_stat_data$min,type="b",col="blue")
506
  lines(1:12,LST_nobs_stat_data$max,type="b",col="red")
507
  text(1:12,LST_nobs_stat_data$mean,rownames(LST_stat_data),cex=1,pos=2)
508
  
509
  legend("topleft",legend=c("min","mean","max"), cex=1.5, col=c("blue","black","red"),
510
         lty=1,bty="n")
511
  
512
  title(paste(var_name,"number of valid observations for Oregon", "2010",sep=" "))
513
  #title(paste("LST number of valid observations for Oregon", "2010",sep=" "))
514
  dev.off()
515
  
516
  #Add plot of number of NA per month...
517
  
518
  list_obj <- list(avg=LST_stat_data,nobs=LST_nobs_stat_data)
519
  return(list_obj)
520

  
521
}
522

  
337 523
################### END OF SCRIPT ###################
338 524

  
339 525

  

Also available in: Unified diff