| 1 |
|
######################################## Paper Methods_comparison: Analyses part 5 #######################################
|
| 2 |
|
############################ Scripts for figures and analyses for paper 2 #####################################
|
| 3 |
|
#This script performs analyses and create figures for the FSS paper.
|
| 4 |
|
#It uses inputs from interpolation objects created at earlier stages...
|
| 5 |
|
#Note that this is exploratory code i.e. not part of the worklfow.
|
| 6 |
|
#AUTHOR: Benoit Parmentier #
|
| 7 |
|
#DATE: 09/13/2013 #
|
| 8 |
|
#PROJECT: NCEAS INPLANT: Environment and Organisms --TASK#491-- #
|
| 9 |
|
###################################################################################################
|
| 10 |
|
|
| 11 |
|
###Loading R library and packages
|
| 12 |
|
#library(gtools) # loading some useful tools
|
| 13 |
|
library(mgcv) # GAM package by Wood 2006 (version 2012)
|
| 14 |
|
library(sp) # Spatial pacakge with class definition by Bivand et al. 2008
|
| 15 |
|
library(spdep) # Spatial package with methods and spatial stat. by Bivand et al. 2012
|
| 16 |
|
library(rgdal) # GDAL wrapper for R, spatial utilities (Keitt et al. 2012)
|
| 17 |
|
library(gstat) # Kriging and co-kriging by Pebesma et al. 2004
|
| 18 |
|
library(automap) # Automated Kriging based on gstat module by Hiemstra et al. 2008
|
| 19 |
|
library(spgwr)
|
| 20 |
|
library(maptools)
|
| 21 |
|
library(graphics)
|
| 22 |
|
library(parallel) # Urbanek S. and Ripley B., package for multi cores & parralel processing
|
| 23 |
|
library(raster)
|
| 24 |
|
library(rasterVis)
|
| 25 |
|
library(plotrix) # Draw circle on graph and additional plotting options
|
| 26 |
|
library(reshape) # Data format and type transformation
|
| 27 |
|
|
| 28 |
|
##################### Function used in the script ##############
|
| 29 |
|
|
| 30 |
|
## Extract a list of object from an object: Useful to extract information from
|
| 31 |
|
## RData objects saved in the interpolation phase.
|
|
1 |
### Analyses and exploration of results for single time scale methods
|
|
2 |
|
|
3 |
### Loading R library and packages
|
|
4 |
#library used in the workflow production:
|
|
5 |
library(gtools) # loading some useful tools
|
|
6 |
library(mgcv) # GAM package by Simon Wood
|
|
7 |
library(sp) # Spatial pacakge with class definition by Bivand et al.
|
|
8 |
library(spdep) # Spatial pacakge with methods and spatial stat. by Bivand et al.
|
|
9 |
library(rgdal) # GDAL wrapper for R, spatial utilities
|
|
10 |
library(gstat) # Kriging and co-kriging by Pebesma et al.
|
|
11 |
library(fields) # NCAR Spatial Interpolation methods such as kriging, splines
|
|
12 |
library(raster) # Hijmans et al. package for raster processing
|
|
13 |
library(gdata) # various tools with xls reading, cbindX
|
|
14 |
library(rasterVis) # Raster plotting functions
|
|
15 |
library(parallel) # Parallelization of processes with multiple cores
|
|
16 |
library(maptools) # Tools and functions for sp and other spatial objects e.g. spCbind
|
|
17 |
library(maps) # Tools and data for spatial/geographic objects
|
|
18 |
library(reshape) # Change shape of object, summarize results
|
|
19 |
library(plotrix) # Additional plotting functions
|
|
20 |
library(plyr) # Various tools including rbind.fill
|
|
21 |
library(spgwr) # GWR method
|
|
22 |
library(automap) # Kriging automatic fitting of variogram using gstat
|
|
23 |
library(rgeos) # Geometric, topologic library of functions
|
|
24 |
#RPostgreSQL # Interface R and Postgres, not used in this script
|
|
25 |
|
|
26 |
#Additional libraries not used in workflow
|
|
27 |
library(pgirmess) # Krusall Wallis test with mulitple options, Kruskalmc {pgirmess}
|
|
28 |
library(ncf)
|
|
29 |
|
|
30 |
#### FUNCTION USED IN SCRIPT
|
|
31 |
|
|
32 |
function_analyses_paper <-"contribution_of_covariates_paper_interpolation_functions_10152013.R"
|
| 32 |
33 |
|
| 33 |
34 |
load_obj <- function(f)
|
| 34 |
35 |
{
|
| ... | ... | |
| 37 |
38 |
env[[nm]]
|
| 38 |
39 |
}
|
| 39 |
40 |
|
| 40 |
|
### Need to improve this function!!!
|
| 41 |
|
calc_stat_prop_tb_diagnostic <-function(names_mod,names_id,tb){
|
| 42 |
|
|
| 43 |
|
t<-melt(subset(tb,pred_mod==names_mod),
|
| 44 |
|
measure=c("mae","rmse","r","me","m50"),
|
| 45 |
|
id=names_id,
|
| 46 |
|
na.rm=T)
|
| 47 |
|
char_tmp <-rep("+",length=length(names_id)-1)
|
| 48 |
|
var_summary <-paste(names_id,sep="",collapse=char_tmp)
|
| 49 |
|
var_summary_formula <-paste(var_summary,collpase="~variable")
|
| 50 |
|
avg_tb<-cast(t,var_summary_formula,mean)
|
| 51 |
|
sd_tb<-cast(t,var_summary_formula,sd)
|
| 52 |
|
n_tb<-cast(t,var_summary_formula,length)
|
| 53 |
|
#n_NA<-cast(t,dst_cat1~variable,is.na)
|
| 54 |
|
|
| 55 |
|
#### prepare returning object
|
| 56 |
|
prop_obj<-list(tb,avg_tb,sd_tb,n_tb)
|
| 57 |
|
names(prop_obj) <-c("tb","avg_tb","sd_tb","n_tb")
|
| 58 |
|
|
| 59 |
|
return(prop_obj)
|
| 60 |
|
}
|
| 61 |
41 |
|
| 62 |
|
#Calculate the difference between training and testing in two different data.frames. Columns to substract are provided.
|
| 63 |
|
diff_df<-function(tb_s,tb_v,list_metric_names){
|
| 64 |
|
tb_diff<-vector("list", length(list_metric_names))
|
| 65 |
|
for (i in 1:length(list_metric_names)){
|
| 66 |
|
metric_name<-list_metric_names[i]
|
| 67 |
|
tb_diff[[i]] <-tb_s[,c(metric_name)] - tb_v[,c(metric_name)]
|
| 68 |
|
}
|
| 69 |
|
names(tb_diff)<-list_metric_names
|
| 70 |
|
tb_diff<-as.data.frame(do.call(cbind,tb_diff))
|
| 71 |
|
return(tb_diff)
|
| 72 |
|
}
|
|
42 |
##############################
|
|
43 |
#### Parameters and constants
|
| 73 |
44 |
|
| 74 |
|
################## PARAMETERS ##########
|
| 75 |
|
|
| 76 |
|
|
| 77 |
|
#path to gam CAI and kriging analyes with hold-out
|
| 78 |
|
in_dir1 <- "/home/parmentier/Data/IPLANT_project/Oregon_interpolation/Oregon_03142013/output_data_365d_gam_CAI_lst_comb3_08312013/"
|
| 79 |
|
in_dir2 <- "/home/parmentier/Data/IPLANT_project/Oregon_interpolation/Oregon_03142013/output_data_365d_gam_CAI_lst_comb3_09012013"
|
| 80 |
|
in_dir3 <-"/home/parmentier/Data/IPLANT_project/Oregon_interpolation/Oregon_03142013/output_data_365d_gam_CAI_lst_comb3_09032013"
|
| 81 |
|
in_dir4 <- "/home/parmentier/Data/IPLANT_project/Oregon_interpolation/Oregon_03142013/output_data_365d_kriging_CAI_lst_comb3_09042013"
|
| 82 |
|
|
| 83 |
|
#path to gam fusion and kriging fusion analyes with hold-out
|
| 84 |
|
in_dir5 <- "/data/project/layers/commons/Oregon_interpolation/output_data_365d_gam_fus_lst_comb3_09092013"
|
| 85 |
|
in_dir6 <- "/data/project/layers/commons/Oregon_interpolation/output_data_365d_gam_fus_lst_comb3_09102013"
|
| 86 |
|
in_dir7 <-"/data/project/layers/commons/Oregon_interpolation/output_data_365d_gam_fus_lst_comb3_09112013"
|
| 87 |
|
in_dir8 <-"/data/project/layers/commons/Oregon_interpolation/output_data_365d_kriging_fus_lst_comb3_09122013"
|
| 88 |
|
in_dir9 <-"/data/project/layers/commons/Oregon_interpolation/output_data_365d_kriging_fus_lst_comb3_09132013"
|
| 89 |
|
in_dir10 <-"/data/project/layers/commons/Oregon_interpolation/output_data_365d_kriging_fus_lst_comb3_09142013"
|
| 90 |
|
in_dir11 <-"/data/project/layers/commons/Oregon_interpolation/output_data_365d_gwr_CAI_lst_comb3_09162013"
|
| 91 |
|
in_dir12 <-"/data/project/layers/commons/Oregon_interpolation/output_data_365d_gwr_CAI_lst_comb3_09172013"
|
| 92 |
|
in_dir13 <-"/data/project/layers/commons/Oregon_interpolation/output_data_365d_gwr_cai_lst_comb3_09282013"
|
| 93 |
|
in_dir14 <-"/data/project/layers/commons/Oregon_interpolation/output_data_365d_gwr_fus_lst_comb3_09232013"
|
| 94 |
|
in_dir15 <-"/data/project/layers/commons/Oregon_interpolation/output_data_365d_gwr_fus_lst_comb3_09262013"
|
| 95 |
|
in_dir16 <-"/data/project/layers/commons/Oregon_interpolation/output_data_365d_gam_cai_lst_comb3_10042013"
|
| 96 |
|
|
| 97 |
|
|
| 98 |
|
#better as list and load one by one specific element from the object
|
| 99 |
|
raster_prediction_obj1 <-load_obj(file.path(in_dir1,"raster_prediction_obj_gam_CAI_dailyTmax_365d_gam_CAI_lst_comb3_08312013.RData"))
|
| 100 |
|
raster_prediction_obj2 <-load_obj(file.path(in_dir2,"raster_prediction_obj_gam_CAI_dailyTmax_365d_gam_CAI_lst_comb3_09012013.RData"))
|
| 101 |
|
raster_prediction_obj3 <-load_obj(file.path(in_dir3,"raster_prediction_obj_gam_CAI_dailyTmax_365d_gam_CAI_lst_comb3_09032013.RData"))
|
| 102 |
|
raster_prediction_obj4 <-load_obj(file.path(in_dir4,"raster_prediction_obj_kriging_CAI_dailyTmax_365d_kriging_CAI_lst_comb3_09042013.RData"))
|
| 103 |
|
|
| 104 |
|
raster_prediction_obj5 <-load_obj(file.path(in_dir5,"raster_prediction_obj_gam_fusion_dailyTmax_365d_gam_fus_lst_comb3_09092013.RData"))
|
| 105 |
|
raster_prediction_obj6 <-load_obj(file.path(in_dir6,"raster_prediction_obj_gam_fusion_dailyTmax_365d_gam_fus_lst_comb3_09102013.RData"))
|
| 106 |
|
raster_prediction_obj7 <-load_obj(file.path(in_dir7,"raster_prediction_obj_gam_fusion_dailyTmax_365d_gam_fus_lst_comb3_09112013.RData"))
|
| 107 |
|
raster_prediction_obj8 <-load_obj(file.path(in_dir8,"raster_prediction_obj_kriging_fusion_dailyTmax_365d_kriging_fus_lst_comb3_09122013.RData"))
|
| 108 |
|
raster_prediction_obj9 <-load_obj(file.path(in_dir9,"raster_prediction_obj_kriging_fusion_dailyTmax_365d_kriging_fus_lst_comb3_09132013.RData"))
|
| 109 |
|
raster_prediction_obj10 <-load_obj(file.path(in_dir10,"raster_prediction_obj_kriging_fusion_dailyTmax_365d_kriging_fus_lst_comb3_09142013.RData"))
|
| 110 |
|
raster_prediction_obj11 <-load_obj(file.path(in_dir11,"raster_prediction_obj_gwr_CAI_dailyTmax_365d_gwr_CAI_lst_comb3_09162013.RData"))
|
| 111 |
|
raster_prediction_obj12 <-load_obj(file.path(in_dir12,"raster_prediction_obj_gwr_CAI_dailyTmax_365d_gwr_CAI_lst_comb3_09172013.RData"))
|
| 112 |
|
raster_prediction_obj13 <-load_obj(file.path(in_dir13,"raster_prediction_obj_gwr_CAI_dailyTmax_365d_gwr_cai_lst_comb3_09282013.RData"))
|
| 113 |
|
raster_prediction_obj14 <-load_obj(file.path(in_dir14,"raster_prediction_obj_gwr_fusion_dailyTmax_365d_gwr_fus_lst_comb3_09232013.RData"))
|
| 114 |
|
raster_prediction_obj15 <-load_obj(file.path(in_dir15,"raster_prediction_obj_gwr_fusion_dailyTmax_365d_gwr_fus_lst_comb3_09262013.RData"))
|
| 115 |
|
|
| 116 |
|
raster_prediction_obj16 <-load_obj(file.path(in_dir16,"raster_prediction_obj_gam_CAI_dailyTmax_365d_gam_cai_lst_comb3_10042013.RData"))
|
| 117 |
|
|
| 118 |
|
out_dir<-"/home/parmentier/Data/IPLANT_project/paper_multitime_scale__analyses_tables_fig_09032013"
|
|
45 |
script_path<-"/home/parmentier/Data/IPLANT_project/env_layers_scripts/" #path to script
|
|
46 |
source(file.path(script_path,function_analyses_paper)) #source all functions used in this script.
|
|
47 |
|
|
48 |
in_dir1 <-"/home/parmentier/Data/IPLANT_project/Oregon_interpolation/Oregon_03142013/output_data_365d_gam_day_lst_comb3_08132013"
|
|
49 |
in_dir2 <-"/home/parmentier/Data/IPLANT_project/Oregon_interpolation/Oregon_03142013/output_data_365d_gam_day_lst_comb4_08152013"
|
|
50 |
#kriging results:
|
|
51 |
in_dir3 <-"/home/parmentier/Data/IPLANT_project/Oregon_interpolation/Oregon_03142013/output_data_365d_kriging_day_lst_comb3_07112013"
|
|
52 |
#gwr results:
|
|
53 |
in_dir4 <-"/home/parmentier/Data/IPLANT_project/Oregon_interpolation/Oregon_03142013/output_data_365d_gwr_day_lst_comb3_part1_07122013"
|
|
54 |
#multisampling results (gam)
|
|
55 |
#in_dir5<- "/home/parmentier/Data/IPLANT_project/Oregon_interpolation/Oregon_03142013/output_data_365d_gam_daily_mults10_lst_comb3_08082013"
|
|
56 |
#in_dir6<- "/home/parmentier/Data/IPLANT_project/Oregon_interpolation/Oregon_03142013/output_data_365d_kriging_daily_mults10_lst_comb3_08062013"
|
|
57 |
#in_dir7<- "/home/parmentier/Data/IPLANT_project/Oregon_interpolation/Oregon_03142013/output_data_365d_gwr_daily_mults10_lst_comb3_08072013"
|
|
58 |
#Hold-out every two days over 365 days
|
|
59 |
in_dir5 <-"/data/project/layers/commons/Oregon_interpolation/output_data_365d_gam_daily_mults1_lst_comb3_10122013"
|
|
60 |
in_dir6 <- "/data/project/layers/commons/Oregon_interpolation/output_data_365d_kriging_daily_mults1_lst_comb3_10112013"
|
|
61 |
in_dir7 <-"/data/project/layers/commons/Oregon_interpolation/output_data_365d_gwr_daily_mults1_lst_comb3_10132013"
|
|
62 |
|
|
63 |
out_dir<-"/home/parmentier/Data/IPLANT_project/paper_analyses_tables_fig_08032013"
|
| 119 |
64 |
setwd(out_dir)
|
|
65 |
|
|
66 |
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 |
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 |
CRS_locs_WGS84<-CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +towgs84=0,0,0") #Station coords WGS84
|
|
69 |
y_var_name <- "dailyTmax"
|
|
70 |
out_prefix<-"analyses_10152013"
|
|
71 |
|
|
72 |
#method_interpolation <- "gam_daily"
|
|
73 |
covar_obj_file_1 <- "covar_obj__365d_gam_day_lst_comb3_08132013.RData"
|
|
74 |
met_obj_file_1 <- "met_stations_outfiles_obj_gam_daily__365d_gam_day_lst_comb3_08132013.RData"
|
|
75 |
#met_stations_outfiles_obj_gam_daily__365d_gam_day_lst_comb3_08132013.RData
|
|
76 |
|
|
77 |
#raster_prediciton object for baseline 1 () s(lat,lon) + s(elev)) and baseline 2 (slat,lon))
|
|
78 |
raster_obj_file_1 <- "raster_prediction_obj_gam_daily_dailyTmax_365d_gam_day_lst_comb3_08132013.RData"
|
|
79 |
raster_obj_file_2 <- "raster_prediction_obj_gam_daily_dailyTmax_365d_gam_day_lst_comb4_08152013.RData"
|
|
80 |
|
|
81 |
raster_obj_file_3 <- "raster_prediction_obj_kriging_daily_dailyTmax_365d_kriging_day_lst_comb3_07112013.RData"
|
|
82 |
raster_obj_file_4 <- "raster_prediction_obj_gwr_daily_dailyTmax_365d_gwr_day_lst_comb3_part1_07122013.RData"
|
|
83 |
#multisampling using baseline lat,lon + elev
|
|
84 |
#raster_obj_file_5 <- "raster_prediction_obj_gam_daily_dailyTmax_365d_gam_daily_mults10_lst_comb3_08082013.RData"
|
|
85 |
#raster_obj_file_6 <- "raster_prediction_obj_kriging_daily_dailyTmax_365d_kriging_daily_mults10_lst_comb3_08062013.RData"
|
|
86 |
#raster_obj_file_7 <- "raster_prediction_obj_gwr_daily_dailyTmax_365d_gwr_daily_mults10_lst_comb3_08072013.RData"
|
|
87 |
|
|
88 |
raster_obj_file_5 <- "raster_prediction_obj_gam_daily_dailyTmax_365d_gam_daily_mults1_lst_comb3_10122013.RData"
|
|
89 |
raster_obj_file_6 <- "raster_prediction_obj_kriging_daily_dailyTmax_365d_kriging_daily_mults1_lst_comb3_10112013.RData"
|
|
90 |
raster_obj_file_7 <- "raster_prediction_obj_gwr_daily_dailyTmax_365d_gwr_daily_mults1_lst_comb3_10132013.RData"
|
|
91 |
|
|
92 |
#Load objects containing training, testing, models objects
|
|
93 |
|
|
94 |
met_stations_obj <- load_obj(met_stations_outfiles_obj_file)
|
|
95 |
covar_obj <-load_obj(file.path(in_dir1,covar_obj_file_1)) #Reading covariates object for GAM daily method
|
|
96 |
infile_covariates <- covar_obj$infile_covariates
|
|
97 |
infile_reg_outline <- covar_obj$infile_reg_outline
|
|
98 |
covar_names<- covar_obj$covar_names
|
|
99 |
#####
|
|
100 |
s_raster <- brick(infile_covariates)
|
|
101 |
names(s_raster)<-covar_names
|
|
102 |
|
|
103 |
raster_prediction_obj_1 <-load_obj(file.path(in_dir1,raster_obj_file_1)) #comb3 (baseline 2)
|
|
104 |
raster_prediction_obj_2 <-load_obj(file.path(in_dir2,raster_obj_file_2)) #comb4 (baseline 1)
|
|
105 |
raster_prediction_obj_3 <-load_obj(file.path(in_dir3,raster_obj_file_3)) #comb3/mod1 baseline 2, kriging
|
|
106 |
raster_prediction_obj_4 <-load_obj(file.path(in_dir4,raster_obj_file_4)) #comb3/mod1 baseline 2, gwr
|
|
107 |
raster_prediction_obj_5 <-load_obj(file.path(in_dir5,raster_obj_file_5)) #gam daily multisampling 10 to 70%
|
|
108 |
raster_prediction_obj_6 <-load_obj(file.path(in_dir6,raster_obj_file_6)) #kriging daily multisampling 10 to 70%
|
|
109 |
raster_prediction_obj_7 <-load_obj(file.path(in_dir7,raster_obj_file_7)) #gwr daily multisampling 10 to 70%
|
|
110 |
|
|
111 |
|
|
112 |
############### BEGIN SCRIPT #################
|
|
113 |
|
|
114 |
############ PART 1: Exploration of surfaces bias, delta and climatology surfaces ###########
|
|
115 |
|
|
116 |
#"/home/parmentier/Data/IPLANT_project/Oregon_interpolation/Oregon_03142013/output_data_365d_gam_day_lst_comb4_07152013/"
|
|
117 |
in_dir<-"/home/parmentier/Data/IPLANT_project/Oregon_interpolation/Oregon_03142013/output_data_365d_kriging_cai_lst_comb3_07312013"
|
|
118 |
|
|
119 |
out_dir<-""
|
|
120 |
setwd(in_dir)
|
| 120 |
121 |
y_var_name <- "dailyTmax"
|
| 121 |
122 |
y_var_month <- "TMax"
|
| 122 |
123 |
#y_var_month <- "LSTD_bias"
|
| 123 |
|
out_suffix <- "_OR_10102013"
|
| 124 |
|
#script_path<-"/data/project/layers/commons/data_workflow/env_layers_scripts/"
|
| 125 |
|
#### FUNCTION USED IN SCRIPT
|
| 126 |
124 |
|
| 127 |
|
function_analyses_paper <-"contribution_of_covariates_paper_interpolation_functions_09232013.R"
|
|
125 |
out_prefix<-"/home/parmentier/Data/IPLANT_project/Oregon_interpolation/Oregon_03142013/output_data_365d_GAM_fus_all_lst_05312013"
|
|
126 |
method_interpolation <- "kriging_CAI"
|
|
127 |
covar_obj_file <- "covar_obj__365d_kriging_cai_lst_comb3_07312013.RData"
|
|
128 |
raster_obj_file <- "raster_prediction_obj_kriging_CAI_dailyTmax_365d_kriging_cai_lst_comb3_07312013.RData"
|
|
129 |
script_path<-"/data/project/layers/commons/data_workflow/env_layers_scripts/"
|
| 128 |
130 |
|
| 129 |
|
script_path<-"/home/parmentier/Data/IPLANT_project/env_layers_scripts/" #path to script
|
| 130 |
|
source(file.path(script_path,function_analyses_paper)) #source all functions used in this script.
|
|
131 |
source(file.path(script_path,"interpolation_method_day_function_multisampling_07052013.R")) #Include GAM_day
|
| 131 |
132 |
|
| 132 |
|
#################################################################################
|
| 133 |
|
############ ANALYSES 1: Average accuracy per proportion for monthly hold out in muli-timescale mehtods... #######
|
|
133 |
#Load objects containing training, testing, models objects
|
|
134 |
covar_obj <-load_obj(covar_obj_file)
|
|
135 |
infile_covariates <- covar_obj$infile_covariates
|
|
136 |
infile_reg_outline <- covar_obj$infile_reg_outline
|
|
137 |
covar_names<- covar_obj$covar_names
|
| 134 |
138 |
|
| 135 |
|
tb_mv_gam_CAI <-rbind(raster_prediction_obj1$tb_month_diagnostic_v,raster_prediction_obj2$tb_month_diagnostic_v,raster_prediction_obj3$tb_month_diagnostic_v)
|
| 136 |
|
tb_ms_gam_CAI <-rbind(raster_prediction_obj1$tb_month_diagnostic_s,raster_prediction_obj2$tb_month_diagnostic_s,raster_prediction_obj3$tb_month_diagnostic_s)
|
|
139 |
raster_prediction_obj <-load_obj(raster_obj_file)
|
| 137 |
140 |
|
| 138 |
|
tb_v_gam_CAI <-rbind(raster_prediction_obj1$tb_diagnostic_v,raster_prediction_obj2$tb_diagnostic_v,raster_prediction_obj3$tb_diagnostic_v)
|
| 139 |
|
tb_s_gam_CAI <-rbind(raster_prediction_obj1$tb_diagnostic_s,raster_prediction_obj2$tb_diagnostic_s,raster_prediction_obj3$tb_diagnostic_s)
|
| 140 |
|
#prop_obj_gam_CAI_v <- calc_stat_prop_tb_diagnostic(names_mod,names_id,tb_v)
|
|
141 |
names(raster_prediction_obj) #list of two objects
|
| 141 |
142 |
|
| 142 |
|
tb_mv_gwr_CAI <-rbind(raster_prediction_obj11$tb_month_diagnostic_v,raster_prediction_obj12$tb_month_diagnostic_v,raster_prediction_obj13$tb_month_diagnostic_v)
|
| 143 |
|
tb_ms_gwr_CAI <-rbind(raster_prediction_obj11$tb_month_diagnostic_s,raster_prediction_obj12$tb_month_diagnostic_s,raster_prediction_obj13$tb_month_diagnostic_s)
|
|
143 |
raster_prediction_obj$summary_metrics_v
|
| 144 |
144 |
|
| 145 |
|
tb_v_gwr_CAI <-rbind(raster_prediction_obj11$tb_diagnostic_v,raster_prediction_obj12$tb_diagnostic_v,raster_prediction_obj13$tb_diagnostic_v)
|
| 146 |
|
tb_s_gwr_CAI <-rbind(raster_prediction_obj11$tb_diagnostic_s,raster_prediction_obj12$tb_diagnostic_s,raster_prediction_obj13$tb_diagnostic_s)
|
|
145 |
j<-1 #selected month for climatology
|
|
146 |
i<-1
|
|
147 |
data_s <- raster_prediction_obj$method_mod_obj[[i]]$data_s #training data
|
|
148 |
data_v <- raster_prediction_obj$method_mod_obj[[i]]$data_v #testing data
|
| 147 |
149 |
|
| 148 |
|
tb_mv_kriging_CAI <- raster_prediction_obj4$tb_month_diagnostic_v
|
| 149 |
|
tb_ms_kriging_CAI <- raster_prediction_obj4$tb_month_diagnostic_s
|
|
150 |
method_mod_obj <- raster_prediction_obj$method_mod_obj #this object contains daily information, training, testing and images
|
| 150 |
151 |
|
| 151 |
|
tb_v_kriging_CAI <- raster_prediction_obj4$tb_diagnostic_v
|
| 152 |
|
tb_s_kriging_CAI <- raster_prediction_obj4$tb_diagnostic_s
|
|
152 |
clim_method_mod_obj <- raster_prediction_obj$clim_method_mod_obj
|
|
153 |
data_month1 <- clim_method_mod_obj[[1]]$data_month #monthly data
|
|
154 |
data_month <- clim_method_mod_obj[[j]]$data_month #monthly data
|
|
155 |
clim_mod_obj_month <- clim_method_mod_obj[[j]]
|
|
156 |
names(clim_mod_obj_month)
|
| 153 |
157 |
|
| 154 |
|
### SAME for gam fusion
|
|
158 |
#####
|
|
159 |
s_raster <- brick(infile_covariates)
|
|
160 |
names(s_raster)<-covar_names
|
|
161 |
#data_s$y_var <- data_s[[y_var_name]]
|
|
162 |
#formula<-"y_var ~ s(lat,lon,elev_s)"
|
|
163 |
#date<-strptime(sampling_dat$date[i], "%Y%m%d") # interpolation date being processed
|
|
164 |
month<-strftime(date, "%m") # current month of the date being processed
|
|
165 |
month<-"01"
|
|
166 |
LST_month<-paste("mm_",month,sep="") # name of LST month to be matched
|
|
167 |
pos<-match("LST",names(s_raster)) #Find the position of the layer with name "LST", if not present pos=NA
|
|
168 |
s_raster<-dropLayer(s_raster,pos) # If it exists drop layer
|
|
169 |
LST<-subset(s_raster,LST_month)
|
|
170 |
names(LST)<-"LST"
|
|
171 |
s_raster<-addLayer(s_raster,LST) #Adding current month
|
| 155 |
172 |
|
| 156 |
|
tb_mv_gam_fus <-rbind(raster_prediction_obj5$tb_month_diagnostic_v,raster_prediction_obj6$tb_month_diagnostic_v,raster_prediction_obj7$tb_month_diagnostic_v)
|
| 157 |
|
tb_ms_gam_fus <-rbind(raster_prediction_obj5$tb_month_diagnostic_s,raster_prediction_obj6$tb_month_diagnostic_s,raster_prediction_obj7$tb_month_diagnostic_s)
|
| 158 |
173 |
|
| 159 |
|
tb_v_gam_fus <-rbind(raster_prediction_obj5$tb_diagnostic_v,raster_prediction_obj6$tb_diagnostic_v,raster_prediction_obj7$tb_diagnostic_v)
|
| 160 |
|
tb_s_gam_fus <-rbind(raster_prediction_obj5$tb_diagnostic_s,raster_prediction_obj6$tb_diagnostic_s,raster_prediction_obj7$tb_diagnostic_s)
|
|
174 |
### MONTH MODELS
|
|
175 |
index<-1
|
|
176 |
data_month$y_var<-data_month[[y_var_month]]
|
|
177 |
mod1 <- clim_method_mod_obj[[1]]$mod[[1]]
|
| 161 |
178 |
|
| 162 |
|
tb_mv_gwr_fus <-rbind(raster_prediction_obj14$tb_month_diagnostic_v,raster_prediction_obj15$tb_month_diagnostic_v)
|
| 163 |
|
tb_ms_gwr_fus <-rbind(raster_prediction_obj14$tb_month_diagnostic_s,raster_prediction_obj15$tb_month_diagnostic_s)
|
| 164 |
179 |
|
| 165 |
|
tb_v_gwr_fus <-rbind(raster_prediction_obj14$tb_diagnostic_v,raster_prediction_obj15$tb_diagnostic_v)
|
| 166 |
|
tb_s_gwr_fus <-rbind(raster_prediction_obj14$tb_diagnostic_s,raster_prediction_obj15$tb_diagnostic_s)
|
|
180 |
clim_method_mod_obj[[1]]$clim #list of files containing model predictions...
|
|
181 |
clim_rast<-stack(clim_method_mod_obj[[index]]$clim)
|
|
182 |
delta_rast<-raster(method_mod_obj[[index]]$delta) #only one delta image!!!
|
|
183 |
pred_temp<-as.character(method_mod_obj[[index]][[y_var_name]]) #list of files with path included
|
|
184 |
rast_pred_temp_s <-stack(pred_temp) #stack of temperature predictions from models (daily)
|
| 167 |
185 |
|
| 168 |
|
tb_mv_kriging_fus <-rbind(raster_prediction_obj8$tb_month_diagnostic_v,raster_prediction_obj9$tb_month_diagnostic_v,raster_prediction_obj10$tb_month_diagnostic_v)
|
| 169 |
|
tb_ms_kriging_fus <-rbind(raster_prediction_obj8$tb_month_diagnostic_s,raster_prediction_obj9$tb_month_diagnostic_s,raster_prediction_obj10$tb_month_diagnostic_s)
|
|
186 |
names(delta_rast)<-"delta"
|
|
187 |
rast_temp_date<-stack(clim_rast,delta_rast)
|
|
188 |
#rast_temp_date<-mask(rast_temp_date,LC_mask,file=file.path(out_path,"test.tif"),overwrite=TRUE)
|
|
189 |
#bias_d_rast<-raster("fusion_bias_LST_20100103_30_1_10d_GAM_fus5_all_lstd_02082013.rst")
|
| 170 |
190 |
|
| 171 |
|
tb_v_kriging_fus <-rbind(raster_prediction_obj8$tb_diagnostic_v,raster_prediction_obj9$tb_diagnostic_v,raster_prediction_obj10$tb_diagnostic_v)
|
| 172 |
|
tb_s_kriging_fus <-rbind(raster_prediction_obj8$tb_diagnostic_s,raster_prediction_obj9$tb_diagnostic_s,raster_prediction_obj10$tb_diagnostic_s)
|
|
191 |
plot(rast_temp_date)
|
| 173 |
192 |
|
| 174 |
|
list_tb <- list(tb_v_gam_CAI,tb_v_kriging_CAI,tb_v_gwr_CAI,tb_s_gam_CAI,tb_s_kriging_CAI,tb_s_gwr_CAI,
|
| 175 |
|
tb_mv_gam_CAI,tb_mv_kriging_CAI,tb_mv_gwr_CAI,tb_ms_gam_CAI,tb_ms_kriging_CAI,tb_ms_gwr_CAI,
|
| 176 |
|
tb_v_gam_fus,tb_v_kriging_fus,tb_v_gwr_fus,tb_s_gam_fus,tb_s_kriging_fus,tb_s_gwr_fus,
|
| 177 |
|
tb_mv_gam_fus,tb_mv_kriging_fus,tb_mv_gwr_fus,tb_ms_gam_fus,tb_ms_kriging_fus,tb_ms_gwr_fus)
|
|
193 |
month_m_rast<-subset(clim_rast,"mod1")
|
|
194 |
day_m_rast<-subset(rast_pred_temp_s,1)
|
| 178 |
195 |
|
| 179 |
|
names(list_tb) <- c("tb_v_gam_CAI","tb_v_kriging_CAI","tb_v_gwr_CAI","tb_s_gam_CAI","tb_s_kriging_CAI","tb_s_gwr_CAI",
|
| 180 |
|
"tb_mv_gam_CAI","tb_mv_kriging_CAI","tb_mv_gwr_CAI","tb_ms_gam_CAI","tb_ms_kriging_CAI","tb_ms_gwr_CAI",
|
| 181 |
|
"tb_v_gam_fus","tb_v_kriging_fus","tb_v_gwr_fus","tb_s_gam_fus","tb_s_kriging_fus","tb_s_gwr_fus",
|
| 182 |
|
"tb_mv_gam_fus","tb_mv_kriging_fus","tb_mv_gwr_fus","tb_ms_gam_fus","tb_ms_kriging_fus","tb_ms_gwr_fus")
|
|
196 |
test<- month_m_rast + delta_rast
|
|
197 |
diff <- test - day_m_rast #this is equal to zeor roughly...
|
| 183 |
198 |
|
| 184 |
|
#list_tb <-list(tb_v_gam_CAI,tb_v_kriging_CAI,tb_v_gwr_CAI,tb_s_gam_CAI,tb_s_kriging_CAI,tb_s_gwr_CAI,tb_mv_gam_CAI,tb_mv_kriging_CAI,tb_ms_gam_CAI,tb_ms_kriging_CAI,tb_ms_gwr_CAI,tb_ms_gwr_CAI #Add fusion here
|
| 185 |
|
# tb_v_gam_fus,tb_v_kriging_fus,tb_s_gam_fus,tb_s_kriging_fus,tb_mv_gam_fus,tb_mv_kriging_fus,tb_ms_gam_fus,tb_ms_kriging_fus) #Add fusion here
|
| 186 |
|
#names(list_tb) <- c("tb_v_gam_CAI","tb_v_kriging_CAI","tb_v_gwr_CAI","tb_s_gam_CAI","tb_s_kriging_CAI","tb_s_gwr_CAI","tb_mv_gam_CAI","tb_mv_kriging_CAI","tb_ms_gam_CAI","tb_ms_kriging_CAI","tb_ms_gwr_CAI","tb_ms_gwr_CAI" #Add fusion here
|
| 187 |
|
# "tb_v_gam_fus","tb_v_kriging_fus","tb_s_gam_fus","tb_s_kriging_fus","tb_mv_gam_fus","tb_mv_kriging_fus","tb_ms_gam_fus","tb_ms_kriging_fus") #Add fusion here
|
|
199 |
s_sgdf<-as(day_m_rast,"SpatialGridDataFrame") #Conversion to spatial grid data frame, only convert the necessary layers!!
|
|
200 |
#s_spdf<-as(day_m_rast,"SpatialPointsDataFrame") #Conversion to spatial grid data frame, only convert the necessary layers!!
|
| 188 |
201 |
|
| 189 |
|
##### DAILY AVERAGE ACCURACY : PLOT AND DIFFERENCES...Cd
|
|
202 |
names(s_sgdf)<-"var1.pred"
|
| 190 |
203 |
|
| 191 |
|
for(i in 1:length(list_tb)){
|
| 192 |
|
#i <- i+1
|
| 193 |
|
tb <-list_tb[[i]]
|
| 194 |
|
plot_name <- names(list_tb)[i]
|
| 195 |
|
pat_str <- "tb_m"
|
| 196 |
|
if(substr(plot_name,start=1,stop=4)== pat_str){
|
| 197 |
|
names_id <- c("pred_mod","prop")
|
| 198 |
|
plot_formula <- paste("rmse","~prop",sep="",collapse="")
|
| 199 |
|
|
| 200 |
|
}else{
|
| 201 |
|
names_id <- c("pred_mod","prop_month")
|
| 202 |
|
plot_formula <- paste("rmse","~prop_month",collapse="")
|
| 203 |
|
}
|
| 204 |
|
names_mod <-unique(tb$pred_mod)
|
| 205 |
|
|
| 206 |
|
prop_obj <- calc_stat_prop_tb_diagnostic(names_mod,names_id,tb)
|
| 207 |
|
|
| 208 |
|
avg_tb <- prop_obj$avg_tb
|
| 209 |
|
|
| 210 |
|
layout_m<-c(1,1) #one row two columns
|
| 211 |
|
par(mfrow=layout_m)
|
| 212 |
|
|
| 213 |
|
png(paste("Figure__accuracy_rmse_prop_month_",plot_name,out_suffix,".png", sep=""),
|
| 214 |
|
height=480*layout_m[1],width=480*layout_m[2])
|
| 215 |
|
|
| 216 |
|
p<- xyplot(as.formula(plot_formula),group=pred_mod,type="b",
|
| 217 |
|
data=avg_tb,
|
| 218 |
|
main=paste("rmse ",plot_name,sep=" "),
|
| 219 |
|
pch=1:length(avg_tb$pred_mod),
|
| 220 |
|
par.settings=list(superpose.symbol = list(
|
| 221 |
|
pch=1:length(avg_tb$pred_mod))),
|
| 222 |
|
auto.key=list(columns=5))
|
| 223 |
|
print(p)
|
|
204 |
mod1$krige_output<-s_sgdf
|
|
205 |
|
|
206 |
plot(mod1) #does not work because we set krige_output to null!!!
|
|
207 |
formula_mod<-formula("y_var~lat*lon + elev_s")
|
|
208 |
col_names<-all.vars(formula_mod) #extract terms names from formula object
|
|
209 |
if (length(col_names)==1){
|
|
210 |
data_fit <-data_month
|
|
211 |
}else{
|
|
212 |
data_fit <- remove_na_spdf(col_names,data_month)
|
|
213 |
}
|
|
214 |
ref_rast<-as(subset(s_raster,1),"SpatialGridDataFrame")
|
|
215 |
s_spdf<-select_var_stack(s_raster,formula_mod,spdf=TRUE) #This only works if s_raster is in memory!!! need to be modified
|
|
216 |
|
|
217 |
proj4string(data_fit)<-proj4string(s_spdf)
|
|
218 |
test_mod <- autoKrige(formula_mod, input_data=data_fit,new_data=s_spdf,data_variogram=data_fit)
|
|
219 |
plot(test_mod)
|
|
220 |
prediction_spdf = test_mod$krige_output
|
|
221 |
sample_variogram = test_mod$exp_var
|
|
222 |
variogram_model = test_mod$var_model
|
|
223 |
|
|
224 |
#### CHECKING THE INPUTS FROM COVARIATES
|
|
225 |
LC1 <- subset(s_raster,"LC1")
|
|
226 |
plot(LC1,colNA=c("red"))
|
|
227 |
|
|
228 |
LC2 <- subset(s_raster,"LC2")
|
|
229 |
plot(LC2,colNA=c("red"))
|
|
230 |
|
|
231 |
LC_names<- paste("LC",1:10,sep="")
|
|
232 |
lc_reg_s <-subset(s_raster,LC_names)
|
|
233 |
plot(lc_reg_s,colNA="red")
|
|
234 |
|
|
235 |
plot(subset(s_raster,"CANHGHT"),colNA="red")
|
|
236 |
|
|
237 |
#Now create mask based on water areas
|
|
238 |
|
|
239 |
LC12<-raster(lc_reg_s,layer=nlayers(lc_reg_s)) #this is open water
|
|
240 |
LC_mask<-LC12
|
|
241 |
LC_mask[LC_mask==100]<-NA
|
|
242 |
LC_mask <- LC_mask > 100
|
|
243 |
|
|
244 |
CANHGHT <- subset(s_raster,"CANHGHT")
|
|
245 |
|
|
246 |
lc_path<-"/data/project/layers/commons/data_workflow/inputs/lc-consensus-global"
|
|
247 |
infile_modis_grid<-"/data/project/layers/commons/data_workflow/inputs/modis_grid/modis_sinusoidal_grid_world.shp" #modis grid tiling system, global
|
|
248 |
infile_elev<-"/data/project/layers/commons/data_workflow/inputs/dem-cgiar-srtm-1km-tif/srtm_1km.tif" #elevation at 1km, global extent to be replaced by the new fused product
|
|
249 |
infile_canheight<-"/data/project/layers/commons/data_workflow/inputs/treeheight-simard2011/Simard_Pinto_3DGlobalVeg_JGR.tif" #Canopy height, global extent
|
|
250 |
infile_distoc <- "/data/project/layers/commons/data_workflow/inputs/distance_to_coast/GMT_intermediate_coast_distance_01d_rev.tif" #distance to coast, global extent at 0.01 deg
|
|
251 |
|
|
252 |
CANH<-raster(infile_canheight)
|
|
253 |
|
|
254 |
LC1_W<-raster(list.files(path=lc_path,full.names=T)[4])
|
|
255 |
|
|
256 |
#Correlation matrix for a subset
|
|
257 |
r<-subset(s_raster,5:8)
|
|
258 |
t44<-layerStats(r,"pearson",na.rm=T)
|
|
259 |
image(t44[[1]])
|
|
260 |
|
|
261 |
###########################################################################################
|
|
262 |
############ PART 2: Granularity-Autocorrelation analyses of predicted surfaces ###########
|
|
263 |
#####
|
|
264 |
|
|
265 |
lf2 <- raster_prediction_obj_2$method_mod_obj[[index]][[y_var_name]]
|
|
266 |
lf2 #contains the models for gam
|
|
267 |
|
|
268 |
pred_temp_s <-stack(lf2)
|
|
269 |
date_selected <- "20109101"
|
|
270 |
#names_layers <-c("mod1=s(lat,long)+s(elev)","mod4=s(lat,long)+s(LST)","diff=mod1-mod4")
|
|
271 |
names_layers <-c("mod1 = s(lat,long)","mod2 = s(lat,long)+s(elev)","mod3 = s(lat,long)+s(N_w)","mod4 = s(lat,long)+s(E_w)",
|
|
272 |
"mod5 = s(lat,long)+s(LST)","mod6 = s(lat,long)+s(DISTOC)","mod7 = s(lat,long)+s(LC1)",
|
|
273 |
"mod8 = s(lat,long)+s(LC1,LST)","mod9 = s(lat,long)+s(CANHGHT)","mod10 = s(lat,long)+s(LST,CANHGHT)")
|
|
274 |
|
|
275 |
#names_layers<-names(pred_temp_s)
|
|
276 |
#names(pred_temp_s)<-names_layers
|
|
277 |
|
|
278 |
s.range <- c(min(minValue(pred_temp_s)), max(maxValue(pred_temp_s)))
|
|
279 |
#s.range <- s.range+c(5,-5)
|
|
280 |
col.breaks <- pretty(s.range, n=200)
|
|
281 |
lab.breaks <- pretty(s.range, n=100)
|
|
282 |
temp.colors <- colorRampPalette(c('blue', 'white', 'red'))
|
|
283 |
max_val<-s.range[2]
|
|
284 |
min_val <-s.range[1]
|
|
285 |
#max_val<- -10
|
|
286 |
#min_val <- 0
|
|
287 |
layout_m<-c(4,3) #one row two columns
|
|
288 |
|
|
289 |
p<-levelplot(pred_temp_s,main="Interpolated Surfaces Model Comparison baseline 1", ylab=NULL,xlab=NULL,
|
|
290 |
par.settings = list(axis.text = list(font = 2, cex = 1.3),layout=layout_m,
|
|
291 |
par.main.text=list(font=2,cex=2),strip.background=list(col="white")),par.strip.text=list(font=2,cex=1.5),
|
|
292 |
names.attr=names_layers,col.regions=temp.colors,at=seq(max_val,min_val,by=0.01))
|
|
293 |
#col.regions=temp.colors(25))
|
|
294 |
print(p)
|
|
295 |
|
|
296 |
#####################################
|
|
297 |
### Create spatial correlogram ####
|
|
298 |
|
|
299 |
r_mod5 <- subset(pred_temp_s,"mod5") #wiht LST
|
|
300 |
r_mod1 <- subset(pred_temp_s,"mod1") #wiht lat,long
|
|
301 |
r_mod2 <- subset(pred_temp_s,"mod2") #wiht elev
|
|
302 |
|
|
303 |
df_mod5 <- as(r_mod5,"SpatialPointsDataFrame")
|
|
304 |
df_mod1 <- as(r_mod1,"SpatialPointsDataFrame")
|
|
305 |
df_mod2 <- as(r_mod2,"SpatialPointsDataFrame")
|
|
306 |
|
|
307 |
r_stack <-stack(subset(s_raster,"mm_01"),pred_temp_s)
|
|
308 |
df_rs <- as(r_stack,"SpatialPointsDataFrame")
|
|
309 |
|
|
310 |
|
|
311 |
correg_t<-correlog(coordinates(df_mod5),df_mod5$mod5)
|
|
312 |
correg_t<-correlog(coordinates(df_mod5)[,1],coordinates(df_mod5)[,2],df_mod5$mod5)
|
|
313 |
|
|
314 |
data_s<-(as.data.frame((list_data_s[[1]])))
|
|
315 |
data_v<-(list_data_v[[1]])
|
|
316 |
|
|
317 |
data_s <-na.omit(data_s[,c("x","y","LST",y_var_name,"elev")])
|
|
318 |
correg_t1 <- correlog(as.matrix(cbind(data_s$x,data_s$y)),z=data_s$LST)
|
|
319 |
correg_t2 <- correlog(as.matrix(cbind(data_s$x,data_s$y)),z=data_s$elev)
|
|
320 |
correg_t3 <- correlog(as.matrix(cbind(data_s$x,data_s$y)),z=data_s[[y_var_name]])
|
|
321 |
correg_t4 <- correlog(as.matrix(cbind(data_s$x,data_s$y)),z=data_s[,c("LST",y_var_name,"elev")])
|
|
322 |
|
|
323 |
plot(correg_t1)
|
|
324 |
plot(correg_t2,add=T)
|
|
325 |
|
|
326 |
correg_t1[,1]
|
|
327 |
correg_t2[,1]
|
|
328 |
|
|
329 |
df_x <- as.data.frame(cbind(correg_t1[,1],correg_t1[,2],correg_t2[,2],correg_t3[,2]))
|
|
330 |
names(df_x)<-c("dist","LST",y_var_name,"elev")
|
|
331 |
|
|
332 |
xyplot(LST+dailyTmax+elev~dist,df_x,type="b",
|
|
333 |
auto.key=list(title="Var", space = "right", cex=1.0),
|
|
334 |
par.settings = list(superpose.symbol=list(pch = 0:3, cex=1)),
|
|
335 |
)
|
|
336 |
|
|
337 |
### For the whole image:
|
|
338 |
|
|
339 |
sp_correlogram_fun <- function(i,list_param){
|
| 224 |
340 |
|
| 225 |
|
dev.off()
|
|
341 |
df <- list_param$list_df[[i]]
|
|
342 |
var_zname <- list_param$var_zname[i]
|
|
343 |
order_lag <- list_param$order_lag[i]
|
|
344 |
method_cor <- list_param$method_cor
|
|
345 |
nb_obj <- list_param$nb_obj
|
|
346 |
randomisation_par <- list_param$randomisation_par
|
|
347 |
sp.cor <- sp.correlogram(nb_obj, df[[var_zname]], order=order_lag,
|
|
348 |
method=method_cor, randomisation=randomisation_par)
|
|
349 |
return(sp.cor)
|
| 226 |
350 |
|
| 227 |
351 |
}
|
| 228 |
352 |
|
| 229 |
|
#xyplot( rmse ~ prop_month | pred_mod,type="b",data=as.data.frame(avg_tb))
|
| 230 |
|
|
| 231 |
|
##### Calculate differences
|
| 232 |
|
|
| 233 |
|
metric_names <- c("mae","rmse","me","r")
|
| 234 |
|
diff_kriging_CAI <- diff_df(tb_s_kriging_CAI,tb_v_kriging_CAI,metric_names)
|
| 235 |
|
diff_gam_CAI <- diff_df(tb_s_gam_CAI[tb_s_gam_CAI$pred_mod!="mod_kr"],tb_v_gam_CAI,metric_names)
|
| 236 |
|
diff_gwr_CAI <- diff_df(tb_s_gwr_CAI,tb_v_gwr_CAI,metric_names)
|
| 237 |
|
|
| 238 |
|
layout_m<-c(1,1) #one row two columns
|
| 239 |
|
par(mfrow=layout_m)
|
| 240 |
|
|
| 241 |
|
png(paste("Figure__accuracy_rmse_prop_month_",plot_name,out_suffix,".png", sep=""),
|
| 242 |
|
height=480*layout_m[1],width=480*layout_m[2])
|
| 243 |
|
boxplot(diff_kriging_CAI$rmse,diff_gam_CAI$rmse,diff_gwr_CAI$rmse,names=c("kriging_CAI","gam_CAI","gwr_CAI"),
|
| 244 |
|
main="Difference between training and testing daily rmse")
|
| 245 |
|
dev.off()
|
| 246 |
|
|
| 247 |
|
#remove prop 0,
|
| 248 |
|
diff_kriging_CAI <- diff_df(tb_s_kriging_CAI[tb_s_kriging_CAI$prop_month!=0,],tb_v_kriging_CAI[tb_v_kriging_CAI$prop_month!=0,],metric_names)
|
| 249 |
|
diff_gam_CAI <- diff_df(tb_s_gam_CAI[tb_s_gam_CAI$prop_month!=0,],tb_v_gam_CAI[tb_v_gam_CAI$prop_month!=0,],metric_names)
|
| 250 |
|
diff_gwr_CAI <- diff_df(tb_s_gwr_CAI[tb_s_gwr_CAI$prop_month!=0,],tb_v_gwr_CAI[tb_v_gwr_CAI$prop_month!=0,],metric_names)
|
| 251 |
|
boxplot(diff_kriging_CAI$rmse,diff_gam_CAI$rmse,diff_gwr_CAI$rmse,names=c("kriging_CAI","gam_CAI","gwr_CAI"),
|
| 252 |
|
main="Difference between training and testing daily rmse")
|
| 253 |
|
|
| 254 |
|
#now monthly accuracy
|
| 255 |
|
metric_names <- c("mae","rmse","me","r")
|
| 256 |
|
diff_kriging_m_CAI <- diff_df(tb_ms_kriging_CAI[tb_ms_kriging_CAI$prop!=0,],tb_mv_kriging_CAI,metric_names)
|
| 257 |
|
diff_gam_m_CAI <- diff_df(tb_ms_gam_CAI[tb_ms_gam_CAI$prop!=0,],tb_mv_gam_CAI,metric_names)
|
| 258 |
|
diff_gwr_m_CAI <- diff_df(tb_ms_gwr_CAI[tb_ms_gwr_CAI$prop!=0,],tb_mv_gwr_CAI,metric_names)
|
| 259 |
|
|
| 260 |
|
layout_m<-c(1,1) #one row two columns
|
| 261 |
|
par(mfrow=layout_m)
|
| 262 |
|
|
| 263 |
|
png(paste("Figure__accuracy_rmse_prop_month_",plot_name,out_suffix,".png", sep=""),
|
| 264 |
|
height=480*layout_m[1],width=480*layout_m[2])
|
| 265 |
|
boxplot(diff_kriging_m_CAI$rmse,diff_gam_m_CAI$rmse,diff_gwr_m_CAI$rmse,names=c("kriging_CAI","gam_CAI","gwr_CAI"),
|
| 266 |
|
main="Difference between training and monhtly testing rmse")
|
| 267 |
|
dev.off()
|
| 268 |
|
|
| 269 |
|
#boxplot(diff_kriging_m_CAI$rmse,diff_gam_m_CAI$rmse,diff_gwr_CAI,names=c("kriging_CAI","gam_CAI","gwr_CAI"),
|
| 270 |
|
# main="Difference between training and monhtly testing rmse")
|
| 271 |
|
|
| 272 |
|
### For fusion
|
| 273 |
|
|
| 274 |
|
metric_names <- c("mae","rmse","me","r")
|
| 275 |
|
diff_kriging_fus <- diff_df(tb_s_kriging_fus,tb_v_kriging_fus,metric_names)
|
| 276 |
|
diff_gam_fus <- diff_df(tb_s_gam_fus,tb_v_gam_fus,metric_names)
|
| 277 |
|
diff_gwr_fus <- diff_df(tb_s_gwr_fus,tb_v_gwr_fus,metric_names)
|
| 278 |
|
|
| 279 |
|
layout_m<-c(1,1) #one row two columns
|
| 280 |
|
par(mfrow=layout_m)
|
| 281 |
|
|
| 282 |
|
png(paste("Figure__accuracy_rmse_prop_month_",plot_name,out_suffix,".png", sep=""),
|
| 283 |
|
height=480*layout_m[1],width=480*layout_m[2])
|
| 284 |
|
boxplot(diff_kriging_fus$rmse,diff_gam_fus$rmse,diff_gwr_fus$rmse,names=c("kriging_fus","gam_fus","gwr_fus"),
|
| 285 |
|
main="Difference between training and testing daily rmse")
|
| 286 |
|
dev.off()
|
| 287 |
|
|
| 288 |
|
metric_names <- c("mae","rmse","me","r")
|
| 289 |
|
diff_kriging_m_fus <- diff_df(tb_ms_kriging_fus[tb_ms_kriging_fus$prop!=0,],tb_mv_kriging_fus[tb_mv_kriging_fus$prop!=0,],metric_names)
|
| 290 |
|
diff_gam_m_fus <- diff_df(tb_ms_gam_fus[tb_ms_gam_fus$prop!=0,],tb_mv_gam_fus[tb_mv_gam_fus$prop!=0,],metric_names)
|
| 291 |
|
diff_gwr_m_fus <- diff_df(tb_ms_gwr_fus[tb_ms_gwr_fus$prop!=0,],tb_mv_gwr_fus[tb_mv_gwr_fus$prop!=0,],metric_names)
|
| 292 |
|
|
| 293 |
|
layout_m<-c(1,1) #one row two columns
|
| 294 |
|
par(mfrow=layout_m)
|
| 295 |
|
|
| 296 |
|
png(paste("Figure__accuracy_rmse_prop_month_",plot_name,out_suffix,".png", sep=""),
|
| 297 |
|
height=480*layout_m[1],width=480*layout_m[2])
|
| 298 |
|
boxplot(diff_kriging_m_fus$rmse,diff_gam_m_fus$rmse,diff_gwr_m_fus$rmse, names=c("kriging_fus","gam_fus","gwr_fus"),
|
| 299 |
|
main="Difference between training and testing FUS rmse")
|
| 300 |
|
dev.off()
|
| 301 |
|
|
| 302 |
|
### NOW PLOT OF COMPARISON BETWEEN Kriging and GAM
|
| 303 |
|
|
| 304 |
|
#Now get variance and range for holdout an dmethods.
|
| 305 |
|
|
| 306 |
|
tb_v_gam_CAI
|
| 307 |
|
tb_v_gam_fus
|
| 308 |
|
tb_v_kriging_CAI
|
| 309 |
|
tb_v_kriging_fus
|
| 310 |
|
|
| 311 |
|
methods_names <- c("tb_v_gam_CAI","tb_v_gam_fus","tb_v_kriging_CAI","tb_v_kriging_fus","tb_v_gwr_CAI","tb_v_gwr_fus")
|
| 312 |
|
list_prop_obj <- vector("list",length=length(methods_names))
|
| 313 |
|
for(i in 1:length(methods_names)){
|
| 314 |
|
tb <- list_tb[[methods_names[i]]]
|
| 315 |
|
names_id <- c("pred_mod","prop_month")
|
| 316 |
|
names_mod <-unique(tb$pred_mod)
|
| 317 |
|
list_prop_obj[[i]] <- calc_stat_prop_tb_diagnostic(names_mod,names_id,tb)
|
| 318 |
|
names(list_prop_obj)[i] <- methods_names[i]
|
| 319 |
|
#avg_tb <- prop_obj$avg_tb
|
|
353 |
# 'nb' - neighbourhood of each cell
|
|
354 |
#r.nb <- dnearneigh(as.matrix(xy), d1=0.5, d2=1.5)
|
|
355 |
# 'nb' - an alternative way to specify the neighbourhood
|
|
356 |
# r.nb <- cell2nb(nrow=side, ncol=side, type="queen")
|
|
357 |
#sp.cor <- sp.correlogram(r.nb, df_mod5$mod5, order=15,
|
|
358 |
# method="I", randomisation=FALSE)
|
|
359 |
|
|
360 |
r_stack <-stack(subset(s_raster,c("mm_01","mm_07")),pred_temp_s)
|
|
361 |
names(r_stack)[1:2]<-c("mm01","mm_07")
|
|
362 |
df_rs <- as(r_stack,"SpatialPointsDataFrame")
|
|
363 |
r.nb <- dnearneigh(coordinates(df_rs), d1=res(s_raster)[1]/2, d2=1.5*res(s_raster)[1]) #lag1
|
|
364 |
|
|
365 |
#Do not run... slow
|
|
366 |
rk.nb14 <- knearneigh(coordinates(df_rs), k=14) #lag1
|
|
367 |
#rk_nb14 <- knearneigh(coordinates(df_rs), k=14) #lag1
|
|
368 |
#save(rk_nb14, file = "rk_nb14.RData")
|
|
369 |
|
|
370 |
|
|
371 |
#rk_nb7 <- knearneigh(coordinates(df_rs), k=7) #lag1
|
|
372 |
#save(rk_nb7, file = "rk_nb7.RData")
|
|
373 |
|
|
374 |
#lrk_nb7 <- knn2nb(rk_nb7)
|
|
375 |
|
|
376 |
#m_LST1 <- moran.test(df_rs$mm01,nb2listw(lrk_nb7),na.action=na.omit,zero.policy=TRUE)
|
|
377 |
#sp.cor <- sp.correlogram(lrk_nb7, df_rs$mm01, order=7,
|
|
378 |
# method="I", randomisation=FALSE)
|
|
379 |
|
|
380 |
|
|
381 |
list_df <- list(df_rs,df_rs,df_rs,df_rs,df_rs)
|
|
382 |
var_zname <- c("mm01","mm_07","mod1","mod2","mod5")
|
|
383 |
order_lag <- c(14,14,14,14,14)
|
|
384 |
method_cor <- "I"
|
|
385 |
nb_obj <- r.nb
|
|
386 |
randomisation_par <- "FALSE"
|
|
387 |
|
|
388 |
list_param_spat_correlog <- list(list_df,var_zname,order_lag,method_cor,nb_obj,randomisation_par)
|
|
389 |
names(list_param_spat_correlog) <- c("list_df","var_zname","order_lag","method_cor","nb_obj","randomisation_par")
|
|
390 |
|
|
391 |
#debug(sp_correlogram_fun)
|
|
392 |
#list_sp_correlog <-sp_correlogram_fun(2,list_param_spat_correlog)
|
|
393 |
#r_qc_s <- lapply(1:length(infile_var),FUN=import_list_modis_layers_fun,list_param=list_param_import_modis)
|
|
394 |
#r_qc_s <-mclapply(1:11,FUN=import_list_modis_layers_fun,list_param=list_param_import_modis,mc.preschedule=FALSE,mc.cores = 11) #This is the end bracket from mclapply(...) statement
|
|
395 |
list_sp_correlog <-mclapply(1:length(list_df),FUN=sp_correlogram_fun,list_param=list_param_spat_correlog ,mc.preschedule=FALSE,mc.cores = 5) #This is the end bracket from mclapply(...) statement
|
|
396 |
#does not work...
|
|
397 |
print(list_sp_correlog[[1]])
|
|
398 |
plot(list_sp_correlog[[1]])
|
|
399 |
print(list_sp_correlog[[2]])
|
|
400 |
plot(list_sp_correlog[[2]])
|
|
401 |
|
|
402 |
##### Use filter option to compute lag Moran's I
|
|
403 |
|
|
404 |
#Queen's case for 5 lags...: should do this in a function to generate filters...
|
|
405 |
#lag 1: 2*1+1 rows
|
|
406 |
f1 <- matrix(c(1,1,1,
|
|
407 |
1,0,1,
|
|
408 |
1,1,1), nrow=3)
|
|
409 |
#lag 2: 2*2+1 rows
|
|
410 |
f2 <- matrix(c(1,1,1,1,1, #filter for lag 2
|
|
411 |
1,0,0,0,1,
|
|
412 |
1,0,0,0,1,
|
|
413 |
1,0,0,0,1,
|
|
414 |
1,1,1,1,1),nrow=5)
|
|
415 |
f3 <- matrix(c(1,1,1,1,1,1,1,
|
|
416 |
1,0,0,0,0,0,1,
|
|
417 |
1,0,0,0,0,0,1,
|
|
418 |
1,0,0,0,0,0,1,
|
|
419 |
1,0,0,0,0,0,1,
|
|
420 |
1,0,0,0,0,0,1,
|
|
421 |
1,1,1,1,1,1,1),nrow=7)
|
|
422 |
f4 <- matrix(c(1,1,1,1,1,1,1,1,1,
|
|
423 |
1,0,0,0,0,0,0,0,1,
|
|
424 |
1,0,0,0,0,0,0,0,1,
|
|
425 |
1,0,0,0,0,0,0,0,1,
|
|
426 |
1,0,0,0,0,0,0,0,1,
|
|
427 |
1,0,0,0,0,0,0,0,1,
|
|
428 |
1,0,0,0,0,0,0,0,1,
|
|
429 |
1,0,0,0,0,0,0,0,1,
|
|
430 |
1,1,1,1,1,1,1,1,1),nrow=9)
|
|
431 |
r<- subset(s_raster,"mm_07")
|
|
432 |
Moran(r,f1)
|
|
433 |
Moran(r,f2)
|
|
434 |
Moran(r,f3)
|
|
435 |
|
|
436 |
#generate automatically filters for MORAN's I in the image...
|
|
437 |
|
|
438 |
autocor_filter_fun <-function(no_lag=1,f_type="queen"){
|
|
439 |
if(f_type=="queen"){
|
|
440 |
no_rows <- 2*no_lag +1
|
|
441 |
border_row <-rep(1,no_rows)
|
|
442 |
other_row <- c(1,rep(0,no_rows-2),1)
|
|
443 |
other_rows <- rep(other_row,no_rows-2)
|
|
444 |
mat_data<- c(border_row,other_rows,border_row)
|
|
445 |
autocor_filter<-matrix(mat_data,nrow=no_rows)
|
|
446 |
}
|
|
447 |
#if(f_type=="rook){} #add later
|
|
448 |
return(autocor_filter)
|
| 320 |
449 |
}
|
| 321 |
450 |
|
| 322 |
|
ac_prop_tb_list <- extract_list_from_list_obj(list_prop_obj,"avg_tb")
|
| 323 |
|
nb_rows <- sapply(ac_prop_tb_list,FUN=nrow)
|
| 324 |
|
method_interp_names<-c("gam_CAI","gam_fus","kriging_CAI","kriging_fus","gwr_CAI","gwr_fus")
|
| 325 |
|
for(i in 1:length(methods_names)){
|
| 326 |
|
avg_tb<-ac_prop_tb_list[[i]]
|
| 327 |
|
avg_tb$method_interp <-rep(x=method_interp_names[i],times=nb_rows[i])
|
| 328 |
|
ac_prop_tb_list[[i]] <- avg_tb
|
|
451 |
#moran_multipe_fun<-function(i,list_param)
|
|
452 |
# lapply(list_filters,FUN=Moran,x=r)
|
|
453 |
|
|
454 |
r<- subset(r_stack,"mod1")
|
|
455 |
Moran(r,f1)
|
|
456 |
Moran(r,f2)
|
|
457 |
list_filters<-lapply(1:5,FUN=autocor_filter_fun,f_type="queen")
|
|
458 |
|
|
459 |
Moran(r,list_filters[[1]])
|
|
460 |
Moran(r,list_filters[[2]])
|
|
461 |
|
|
462 |
plot(subset(s_raster,"mm_09"))
|
|
463 |
|
|
464 |
r_stack <-stack(subset(s_raster,c("mm_09")),pred_temp_s)
|
|
465 |
names(r_stack)[1]<-c("mm_09")
|
|
466 |
|
|
467 |
r<- subset(r_stack,"mod1")
|
|
468 |
Moran(r) #with lag 1 and default rooks lag correlation
|
|
469 |
|
|
470 |
list_filters<-lapply(1:5,FUN=autocor_filter_fun,f_type="queen")
|
|
471 |
|
|
472 |
|
|
473 |
#cacluate Moran's I for 5 lags for one layer
|
|
474 |
moran_list <- lapply(list_filters,FUN=Moran,x=r)
|
|
475 |
|
|
476 |
moran_multiple_fun<-function(i,list_param){
|
|
477 |
#un
|
|
478 |
list_filters <-list_param$list_filters
|
|
479 |
r <- subset(list_param$r_stack,i)
|
|
480 |
moran_list <- lapply(list_filters,FUN=Moran,x=r)
|
|
481 |
moran_v <-as.data.frame(unlist(moran_list))
|
|
482 |
names(moran_v)<-names(r)
|
|
483 |
return(moran_v)
|
| 329 |
484 |
}
|
| 330 |
|
#lapply(methods_names,function(i) {rep(x[i],nrows[i],times[i])},times=nb_rows)
|
| 331 |
|
|
| 332 |
|
#names(ac_prop_tb_list) <- names(list_prop_obj)
|
| 333 |
|
|
| 334 |
|
t44 <- do.call(rbind,ac_prop_tb_list) #contains all accuracy by method, proportion, model, sample etc.
|
| 335 |
|
View(t44)
|
| 336 |
|
|
| 337 |
|
t44[which.min(t44$rmse),] #Find the mimum rmse across all models and methods...
|
| 338 |
|
|
| 339 |
|
test <- t44[order(t44$rmse),]
|
| 340 |
|
test[1:24,]
|
| 341 |
|
|
| 342 |
|
test2<-test[test$method_interp%in% c("gam_fus","gam_CAI"),]
|
| 343 |
|
test2[1:24,]
|
| 344 |
|
#head(ac_prop_tb)
|
| 345 |
|
test3<-subset(test,prop_month==0 & method_interp%in%c("gam_CAI"))
|
| 346 |
|
#test3<-test[test$method_interp%in% c("gam_CAI"),]
|
| 347 |
|
test3[1:24,]
|
| 348 |
|
|
| 349 |
|
#list_prop_obj$avg_tb
|
| 350 |
|
#xyplot(as.formula(plot_formula),group=pred_mod,type="b",
|
| 351 |
|
# data=avg_tb,
|
| 352 |
|
# main=paste("rmse ",plot_name,sep=" "),
|
| 353 |
|
# pch=1:length(avg_tb$pred_mod),
|
| 354 |
|
# par.settings=list(superpose.symbol = list(
|
| 355 |
|
# pch=1:length(avg_tb$pred_mod))),
|
| 356 |
|
# auto.key=list(columns=5))
|
| 357 |
|
|
| 358 |
|
## DAILY DEVIATIONS WITH MODELS
|
| 359 |
|
### Examining results with models for daily deviation called 2: need to subset the models
|
| 360 |
|
|
| 361 |
|
#c("mod1.dev_mod1","mod1.dev_mod2","mod2.dev_mod1","mod2.dev_mod2","mod3.dev_mod1","mod3.dev_mod2","mod4.dev_mod1",
|
| 362 |
|
# "mod4.dev_mod2","mod5.dev_mod1","mod5.dev_mod2","mod6.dev_mod1","mod6.dev_mod2","mod7.dev_mod1","mod7.dev_mod2",
|
| 363 |
|
# "mod_kr.dev_mod1" "mod_kr.dev_mod2")
|
| 364 |
|
|
| 365 |
|
tb_s_gam_CAI_selected <-subset(tb_s_gam_CAI,prop_month==0 & pred_mod%in%c("mod1","mod4","mod5","mod6","mod7","mod8"))
|
| 366 |
|
tb_v_gam_CAI_selected <-subset(tb_v_gam_CAI,prop_month==0 & pred_mod%in%c("mod1","mod4","mod5","mod6","mod7","mod8"))
|
| 367 |
|
tb_s_gwr_CAI_selected <-subset(tb_s_gwr_CAI,prop_month==0 & pred_mod%in%c("mod1","mod4","mod5","mod6","mod7","mod8"))
|
| 368 |
|
tb_v_gwr_CAI_selected <-subset(tb_v_gwr_CAI,prop_month==0 & pred_mod%in%c("mod1","mod4","mod5","mod6","mod7","mod8"))
|
| 369 |
|
tb_s_kriging_CAI_selected <-subset(tb_s_kriging_CAI,prop_month==0 & pred_mod%in%c("mod1","mod4","mod5","mod6","mod7","mod8"))
|
| 370 |
|
tb_v_kriging_CAI_selected <-subset(tb_v_kriging_CAI,prop_month==0 & pred_mod%in%c("mod1","mod4","mod5","mod6","mod7","mod8"))
|
| 371 |
|
|
| 372 |
|
tb_s_gam_CAI2_selected <- subset(raster_prediction_obj16$tb_diagnostic_s,pred_mod%in%c("mod1.dev_mod1","mod1.dev_mod2","mod3.dev_mod1",
|
| 373 |
|
"mod3.dev_mod2","mod4.dev_mod1", "mod4.dev_mod2",
|
| 374 |
|
"mod5.dev_mod1","mod5.dev_mod2","mod6.dev_mod1",
|
| 375 |
|
"mod6.dev_mod2","mod7.dev_mod1","mod7.dev_mod2"))
|
| 376 |
|
tb_v_gam_CAI2_selected <- subset(raster_prediction_obj16$tb_diagnostic_v,pred_mod%in%c("mod1.dev_mod1","mod1.dev_mod2","mod3.dev_mod1",
|
| 377 |
|
"mod3.dev_mod2","mod4.dev_mod1", "mod4.dev_mod2",
|
| 378 |
|
"mod5.dev_mod1","mod5.dev_mod2","mod6.dev_mod1",
|
| 379 |
|
"mod6.dev_mod2","mod7.dev_mod1","mod7.dev_mod2"))
|
| 380 |
|
|
| 381 |
|
metric_names <- c("mae","rmse","me","r")
|
| 382 |
|
|
| 383 |
|
diff_gam_cai2 <- diff_df(tb_s_gam_CAI2_selected,tb_v_gam_CAI2_selected,metric_names)
|
| 384 |
|
diff_kriging_CAI_selected <- diff_df(tb_s_kriging_CAI_selected,tb_v_kriging_CAI_selected,metric_names)
|
| 385 |
|
diff_gam_CAI_selected <- diff_df(tb_s_gam_CAI_selected,tb_v_gam_CAI_selected,metric_names)
|
| 386 |
|
diff_gwr_CAI_selected <- diff_df(tb_s_gwr_CAI_selected,tb_v_gwr_CAI_selected,metric_names)
|
| 387 |
|
|
| 388 |
|
layout_m<-c(1,1) #one row two columns
|
| 389 |
|
par(mfrow=layout_m)
|
| 390 |
|
|
| 391 |
|
png(paste("Figure__accuracy_rmse_prop_month_",plot_name,out_suffix,".png", sep=""),
|
| 392 |
|
height=480*layout_m[1],width=480*layout_m[2])
|
| 393 |
|
boxplot(diff_gam_cai2$rmse,diff_gam_CAI_selected$rmse,
|
| 394 |
|
diff_kriging_CAI_selected$rmse,diff_gwr_CAI_selected$rmse,names=c("gam_cai2","gam_CAI","kriging","gwr"))
|
| 395 |
|
|
| 396 |
|
dev.off()
|
|
485 |
|
|
486 |
list_filters<-lapply(1:10,FUN=autocor_filter_fun,f_type="queen")
|
|
487 |
|
|
488 |
list_param_moran <- list(list_filters=list_filters,r_stack=r_stack)
|
|
489 |
#moran_r <-moran_multiple_fun(1,list_param=list_param_moran)
|
|
490 |
nlayers(r_stack)
|
|
491 |
moran_I_df <-mclapply(1:nlayers(r_stack), list_param=list_param_moran, FUN=moran_multiple_fun,mc.preschedule=FALSE,mc.cores = 11) #This is the end bracket from mclapply(...) statement
|
|
492 |
|
|
493 |
moran_df <- do.call(cbind,moran_I_df)
|
|
494 |
moran_df$lag <-1:nrow(moran_df)
|
|
495 |
|
|
496 |
#melt(moran_df,id=names(moran_df))
|
|
497 |
#moran_df <- do.call(rbind,moran_I_df)
|
|
498 |
mydata<-moran_df
|
|
499 |
dd <- do.call(make.groups, mydata[,-ncol(mydata)])
|
|
500 |
dd$lag <- mydata$lag
|
|
501 |
#names(dd)[2]<-"models"
|
|
502 |
names_layers <-c("LST",names_layers)
|
|
503 |
|
|
504 |
xyplot(data ~ lag | which, dd,type="b",strip=strip.custom(factor.levels=names_layers))
|
|
505 |
|
|
506 |
#solve problem wiht name
|
analyses paperm, experimentation with correlograms in predictions