Project

General

Profile

« Previous | Next » 

Revision 215cad8e

Added by Benoit Parmentier over 10 years ago

contributions of covariates OR paper, slight modifications, adding LST climatologies figures

View differences:

climate/research/oregon/interpolation/contribution_of_covariates_paper_interpolation.R
4 4 
#different covariates using two baselines. Accuracy methods are added in the the function script to evaluate results.
5 5 
#Figures, tables and data for the contribution of covariate paper are also produced in the script.
6 6 
#AUTHOR: Benoit Parmentier
7 
#MMODIFIED ON: 10/15/2013
8 
#Version: 4
7 
#MODIFIED ON: 05/21/2014
8 
#Version: 5
9 9 
#PROJECT: Environmental Layers project
10 10 
#################################################################################################
11 11 

  		




  ......




  37
  37
  
    
library(ncf)


  		




  38
  38
  
    
#### FUNCTION USED IN SCRIPT


  		




  39
  39
  
    

  		




  40
  
  
    
function_analyses_paper <-"contribution_of_covariates_paper_interpolation_functions_10152013.R"


  		




  
  40
  
    
function_analyses_paper <-"contribution_of_covariates_paper_interpolation_functions_05212014.R"


  		




  41
  41
  
    

  		




  42
  42
  
    
##############################


  		




  43
  43
  
    
#### Parameters and constants  


  		




  ......




  60
  60
  
    
in_dir6 <- "/data/project/layers/commons/Oregon_interpolation/output_data_365d_kriging_daily_mults1_lst_comb3_10112013"


  		




  61
  61
  
    
in_dir7 <-"/data/project/layers/commons/Oregon_interpolation/output_data_365d_gwr_daily_mults1_lst_comb3_10132013"


  		




  62
  62
  
    

  		




  63
  
  
    
out_dir<-"/home/parmentier/Data/IPLANT_project/paper_analyses_tables_fig_08032013"


  		




  64
  
  
    
setwd(out_dir)


  		




  65
  63
  
    

  		




  66
  64
  
    
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


  		




  67
  65
  
    
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"


  		




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


  		




  69
  67
  
    
y_var_name <- "dailyTmax"


  		




  70
  
  
    
out_prefix<-"analyses_10152013"


  		




  
  68
  
    
out_prefix<-"_05252014"


  		




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


  		




  
  70
  
    
setwd(out_dir)


  		




  
  71
  
    

  		




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


  		




  
  73
  
    
create_out_dir_param = TRUE


  		




  
  74
  
    

  		




  
  75
  
    
#Create output directory


  		




  
  76
  
    

  		




  
  77
  
    
if(create_out_dir_param==TRUE){


  		




  
  78
  
    
  out_dir <- create_dir_fun(out_dir,out_prefix)


  		




  
  79
  
    
  setwd(out_dir)


  		




  
  80
  
    
}else{


  		




  
  81
  
    
  setwd(out_dir) #use previoulsy defined directory


  		




  
  82
  
    
}


  		




  71
  83
  
    

  		




  72
  84
  
    
#method_interpolation <- "gam_daily"


  		




  73
  85
  
    
covar_obj_file_1 <- "covar_obj__365d_gam_day_lst_comb3_08132013.RData"


  		




  ......




  99
  111
  
    
#####


  		




  100
  112
  
    
s_raster <- brick(infile_covariates)


  		




  101
  113
  
    
names(s_raster)<-covar_names


  		




  
  114
  
    
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


  		




  
  115
  
    
ref_rast_d001 <-"/data/project/layers/commons/data_workflow/inputs/region_outlines_ref_files/mean_day001_rescaled.rst"


  		




  102
  116
  
    

  		




  103
  117
  
    
raster_prediction_obj_1 <-load_obj(file.path(in_dir1,raster_obj_file_1)) #comb3 (baseline 2)


  		




  104
  118
  
    
raster_prediction_obj_2 <-load_obj(file.path(in_dir2,raster_obj_file_2)) #comb4 (baseline 1)


  		




  ......




  223
  237
  
    
####### Now create figures #############


  		




  224
  238
  
    

  		




  225
  239
  
    
#figure 1: study area


  		




  226
  
  
    
#figure 2: methodological worklfow


  		




  227
  
  
    
#figure 3:Figure 3. MAE/RMSE and distance to closest fitting station.


  		




  228
  
  
    
#Figure 4. RMSE and MAE, mulitisampling and hold out for FSS and GAM.


  		




  229
  
  
    
#Figure 5. Overtraining tendency


  		




  230
  
  
    
#Figure 6: Spatial pattern of prediction for one day


  		




  
  240
  
    
#figure 2: methodological worklfow (generated outside R)


  		




  
  241
  
    
#figure 3: LST daily and monthly climatology


  		




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


  		




  
  243
  
    
#Figure 5. RMSE and MAE, mulitisampling and hold out for FSS and GAM.


  		




  
  244
  
    
#Figure 6. Overtraining tendency


  		




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


  		




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


  		




  
  247
  
    
#Figure 8: Figure 8. (a) Monthly MAE averages for the three interpolation methods: GAM, GWR and Kriging.(b) Monthly MAE boxplot for GAM.


  		




  
  248
  
    
#Figure 9: difference image


  		




  231
  249
  
    

  		




  232
  250
  
    
### Figure 1: Oregon study area


  		




  233
  251
  
    
#3 parameters from output


  		




  ......




  286
  304
  
    
box()


  		




  287
  305
  
    
dev.off()


  		




  288
  306
  
    

  		




  
  307
  
    
#########################


  		




  289
  308
  
    
### Figure 2:  Method comparison workflow 


  		




  290
  309
  
    

  		




  291
  310
  
    
# Workflow figure is not generated in R


  		




  292
  311
  
    

  		




  
  312
  
    
##########################


  		




  
  313
  
    
### Figure 3: LST averaging: daily mean compared to monthly mean


  		




  
  314
  
    

  		




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


  		




  
  316
  
    

  		




  
  317
  
    
lst_md <- raster(ref_rast_name)


  		




  
  318
  
    
projection(lst_md)<-projection(s_raster)


  		




  
  319
  
    
lst_mm_09<-subset(s_raster,"mm_09")


  		




  
  320
  
    

  		




  
  321
  
    
lst_md<-raster(ref_rast_d001)


  		




  
  322
  
    
lst_md<- lst_md - 273.16


  		




  
  323
  
    
lst_mm_01<-subset(s_raster,"mm_01")


  		




  
  324
  
    

  		




  
  325
  
    
no_brks <- length(seq(min_val,max_val,by=0.1))-1


  		




  
  326
  
    
#temp.colors <- colorRampPalette(c('blue', 'white', 'red'))


  		




  
  327
  
    
#temp.colors <- colorRampPalette(c('blue', 'lightgoldenrodyellow', 'red'))


  		




  
  328
  
    
#temp.colors <- matlab.like(no_brks)


  		




  
  329
  
    
temp.colors <- colorRampPalette(c('blue', 'khaki', 'red'))


  		




  
  330
  
    

  		




  
  331
  
    
png(filename=paste("Figure_3_paper_Comparison_daily_monthly_mean_lst",out_prefix,".png",sep=""),width=960,height=480)


  		




  
  332
  
    
par(mfrow=c(1,2))


  		




  
  333
  
    
plot(lst_md,col=temp.colors(25))


  		




  
  334
  
    
plot(interp_area,add=TRUE)


  		




  
  335
  
    
title("Mean for January 1")


  		




  
  336
  
    
plot(lst_mm_01,col=temp.colors(25))


  		




  
  337
  
    
plot(interp_area,add=TRUE)


  		




  
  338
  
    
title("Mean for month of January")


  		




  
  339
  
    
dev.off()


  		




  
  340
  
    

  		




  293
  341
  
    
################################################


  		




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


  		




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


  		




  295
  343
  
    

  		




  296
  344
  
    
#Analysis accuracy in term of distance to closest station


  		




  297
  345
  
    
#Assign model's names


  		












  

    
      climate/research/oregon/interpolation/contribution_of_covariates_paper_interpolation_functions.R
    
  





  4
  4
  
    
# interpolation code.


  		




  5
  5
  
    
#Figures and data for the contribution of covariate paper are also produced.


  		




  6
  6
  
    
#AUTHOR: Benoit Parmentier                                                                      #


  		




  7
  
  
    
#DATE: 10/15/2013            


  		




  
  7
  
    
#DATE: 05/21/2014            


  		




  8
  8
  
    
#Version: 2


  		




  9
  9
  
    
#PROJECT: Environmental Layers project                                       #


  		




  10
  10
  
    
#################################################################################################


  		




  ......




  29
  29
  
    

  		




  30
  30
  
    
#### FUNCTION USED IN SCRIPT


  		




  31
  31
  
    

  		




  32
  
  
    
function_analyses_paper <-"contribution_of_covariates_paper_interpolation_functions_08152013.R"


  		




  
  32
  
    
function_analyses_paper <-"contribution_of_covariates_paper_interpolation_functions_05252014.R"


  		




  33
  33
  
    

  		




  34
  34
  
    
load_obj <- function(f)


  		




  35
  35
  
    
{


  		




  ......




  38
  38
  
    
  env[[nm]]


  		




  39
  39
  
    
}


  		




  40
  40
  
    

  		




  
  41
  
    
create_dir_fun <- function(out_dir,out_suffix){


  		




  
  42
  
    
  if(!is.null(out_suffix)){


  		




  
  43
  
    
    out_dir <- paste(out_dir,"_",out_suffix,sep="")


  		




  
  44
  
    
    #out_dir <- file.path(out_dir,out_name)


  		




  
  45
  
    
  }


  		




  
  46
  
    
  #create if does not exists


  		




  
  47
  
    
  if(!file.exists(out_dir)){


  		




  
  48
  
    
    dir.create(out_dir)


  		




  
  49
  
    
  }


  		




  
  50
  
    
  return(out_dir)


  		




  
  51
  
    
}


  		




  
  52
  
    

  		




  41
  53
  
    
extract_list_from_list_obj<-function(obj_list,list_name){


  		




  42
  54
  
    
  #Create a list of an object from a given list of object using a name prodived as input


  		




  43
  55
  
    
  


  		




  ......




  57
  69
  
    
  list_tmp<-vector("list",length(obj_list))


  		




  58
  70
  
    
  for (i in 1:length(obj_list)){


  		




  59
  71
  
    
    tmp<-obj_list[[i]][[list_name]] #double bracket to return data.frame


  		




  60
  
  
    
    list_tmp[[i]]<-tmp


  		




  
  72
  
    
    list_tmp[[i]]<- as.data.frame(tmp) #deal with spdf cases


  		




  61
  73
  
    
  }


  		




  62
  74
  
    
  tb_list_tmp<-do.call(rbind.fill,list_tmp) #long rownames


  		




  63
  75
  
    
  #tb_list_tmp<-do.call(rbind,list_tmp) #long rownames


  		












Also available in:  Unified diff