1
|
##################################### METHOD COMPARISON ##########################################
|
2
|
#################################### CLIMATE INTERPOLATION ########################################
|
3
|
#This script utilizes the R ojbects created during the interpolation phase. #
|
4
|
#At this stage the stcrip produce figures of various accuracy metrics. #
|
5
|
#AUTHOR: Benoit Parmentier #
|
6
|
#DATE: 08/26/2012 #
|
7
|
#PROJECT: NCEAS INPLANT: Environment and Organisms --TASK#??-- #
|
8
|
###################################################################################################
|
9
|
|
10
|
###Loading R library and packages
|
11
|
#library(gtools) # loading some useful tools
|
12
|
library(mgcv) # GAM package by Wood 2006 (version 2012)
|
13
|
library(sp) # Spatial pacakge with class definition by Bivand et al. 2008
|
14
|
library(spdep) # Spatial package with methods and spatial stat. by Bivand et al. 2012
|
15
|
library(rgdal) # GDAL wrapper for R, spatial utilities (Keitt et al. 2012)
|
16
|
library(gstat) # Kriging and co-kriging by Pebesma et al. 2004
|
17
|
library(automap) # Automated Kriging based on gstat module by Hiemstra et al. 2008
|
18
|
library(spgwr)
|
19
|
library(gpclib)
|
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
|
26
|
|
27
|
## Functions
|
28
|
#loading R objects that might have similar names
|
29
|
load_obj <- function(f)
|
30
|
{
|
31
|
env <- new.env()
|
32
|
nm <- load(f, env)[1]
|
33
|
env[[nm]]
|
34
|
}
|
35
|
|
36
|
###Parameters and arguments
|
37
|
|
38
|
infile1<- "ghcn_or_tmax_covariates_06262012_OR83M.shp" #GHCN shapefile containing variables for modeling 2010
|
39
|
#infile2<-"list_10_dates_04212012.txt" #List of 10 dates for the regression
|
40
|
infile2<-"list_365_dates_04212012.txt" #list of dates
|
41
|
infile3<-"LST_dates_var_names.txt" #LST dates name
|
42
|
infile4<-"models_interpolation_05142012.txt" #Interpolation model names
|
43
|
infile5<-"mean_day244_rescaled.rst" #mean LST for day 244
|
44
|
inlistf<-"list_files_05032012.txt" #list of raster images containing the Covariates
|
45
|
|
46
|
obj_list<-"list_obj_08262012.txt" #Results of fusion from the run on ATLAS
|
47
|
path<-"/home/parmentier/Data/IPLANT_project/methods_interpolation_comparison" #Jupiter LOCATION on Atlas for kriging #Jupiter Location on XANDERS
|
48
|
#path<-"/Users/benoitparmentier/Dropbox/Data/NCEAS/Oregon_covariates/" #Local dropbox folder on Benoit's laptop
|
49
|
setwd(path)
|
50
|
proj_str="+proj=lcc +lat_1=43 +lat_2=45.5 +lat_0=41.75 +lon_0=-120.5 +x_0=400000 +y_0=0 +ellps=GRS80 +units=m +no_defs";
|
51
|
#Number of kriging model
|
52
|
out_prefix<-"methods_08262012_" #User defined output prefix
|
53
|
|
54
|
lines<-read.table(paste(path,"/",inlistf,sep=""), sep="") #Column 1 contains the names of raster files
|
55
|
inlistvar<-lines[,1]
|
56
|
inlistvar<-paste(path,"/",as.character(inlistvar),sep="")
|
57
|
|
58
|
#mention this is the last... files
|
59
|
|
60
|
### RESULTS COMPARISON
|
61
|
|
62
|
# PART 1 : using R object created during the interpolation phase
|
63
|
|
64
|
lines<-read.table(paste(path,"/",obj_list,sep=""), sep=",") #Column 1 contains the names RData objects
|
65
|
inlistobj<-lines[,1]
|
66
|
inlistobj<-paste(path,"/",as.character(inlistobj),sep="")
|
67
|
obj_names<-as.character(lines[,2]) #Column two contains short names for obj. model
|
68
|
|
69
|
nel<-length(inlistobj)
|
70
|
method_mod <-vector("list",nel) #list of one row data.frame
|
71
|
method_tb <-vector("list",nel) #list of one row data.frame
|
72
|
method_mean<-vector("list",nel)
|
73
|
|
74
|
for (i in 1:length(inlistobj)){
|
75
|
obj_tmp<-load_obj(inlistobj[i])
|
76
|
method_mod[[i]]<-obj_tmp
|
77
|
#names(method_mod[[i]])<-obj_names[i]
|
78
|
}
|
79
|
obj_tmp<-load_obj(inlistobj[i])
|
80
|
|
81
|
names(method_mod)<-obj_names
|
82
|
|
83
|
#Condense and add other comparison
|
84
|
|
85
|
for(k in 1:length(method_mod)){ # start of the for main loop to all methods
|
86
|
|
87
|
tb<-method_mod[[k]][[1]][[3]][0,] #copy
|
88
|
mod_tmp<-method_mod[[k]]
|
89
|
|
90
|
for (i in 1:365){ # Assuming 365 days of prediction
|
91
|
tmp<-mod_tmp[[i]][[3]]
|
92
|
tb<-rbind(tb,tmp)
|
93
|
}
|
94
|
rm(mod_tmp)
|
95
|
|
96
|
for(i in 4:(ncol(tb))){ # start of the for loop #1
|
97
|
tb[,i]<-as.numeric(as.character(tb[,i]))
|
98
|
}
|
99
|
|
100
|
method_tb[[k]]<-tb
|
101
|
tb_RMSE<-subset(tb, metric=="RMSE")
|
102
|
tb_MAE<-subset(tb,metric=="MAE")
|
103
|
tb_ME<-subset(tb,metric=="ME")
|
104
|
tb_R2<-subset(tb,metric=="R2")
|
105
|
tb_RMSE_f<-subset(tb, metric=="RMSE_f")
|
106
|
tb_MAE_f<-subset(tb,metric=="MAE_f")
|
107
|
|
108
|
tb_diagnostic1<-rbind(tb_RMSE,tb_MAE,tb_ME,tb_R2)
|
109
|
|
110
|
na_mod<-colSums(!is.na(tb_RMSE[,4:ncol(tb)]))
|
111
|
for (j in 4:ncol(tb)){
|
112
|
|
113
|
if (na_mod[j-3]<183){
|
114
|
tb_RMSE<-tb_RMSE[,-j] #Remove columns that has too many missing values!!!
|
115
|
}
|
116
|
}
|
117
|
|
118
|
na_mod<-colSums(!is.na(tb_MAE[,4:ncol(tb)]))
|
119
|
for (j in 4:ncol(tb)){
|
120
|
|
121
|
if (na_mod[j-3]<183){
|
122
|
tb_MAE<-tb_MAE[,-j] #Remove columns that has too many missing values!!!
|
123
|
}
|
124
|
}
|
125
|
|
126
|
na_mod<-colSums(!is.na(tb_MAE_f[,4:ncol(tb)]))
|
127
|
for (j in 4:ncol(tb)){
|
128
|
|
129
|
if (na_mod[j-3]<183){
|
130
|
tb_MAE_f<-tb_MAE_f[,-j] #Remove columns that has too many missing values!!!
|
131
|
}
|
132
|
}
|
133
|
|
134
|
na_mod<-colSums(!is.na(tb_ME[,4:ncol(tb)]))
|
135
|
for (j in 4:ncol(tb)){
|
136
|
|
137
|
if (na_mod[j-3]<183){
|
138
|
tb_ME<-tb_ME[,-j] #Remove columns that has too many missing values!!!
|
139
|
}
|
140
|
}
|
141
|
|
142
|
#Add assessment of missing prediction over the year.
|
143
|
|
144
|
mean_RMSE<-sapply(tb_RMSE[,4:ncol(tb_RMSE)],mean,na.rm=T)
|
145
|
mean_MAE<-sapply(tb_MAE[,4:ncol(tb_MAE)],mean,na.rm=T)
|
146
|
mean_R2<-sapply(tb_R2[,4:ncol(tb_R2)],mean, n.rm=T)
|
147
|
mean_ME<-sapply(tb_ME[,4:ncol(tb_ME)],mean,na.rm=T)
|
148
|
mean_MAE_f<-sapply(tb_MAE[,4:ncol(tb_MAE_f)],mean,na.rm=T)
|
149
|
mean_RMSE_f<-sapply(tb_RMSE_f[,4:ncol(tb_RMSE_f)],mean,na.rm=T)
|
150
|
mean_list<-list(mean_RMSE,mean_MAE,mean_R2,mean_ME,mean_MAE_f,mean_RMSE_f)
|
151
|
names(mean_list)<-c("RMSE","MAE","R2","ME","MAE_f","RMSE_f")
|
152
|
method_mean[[k]]<-mean_list
|
153
|
names_methods<-obj_names
|
154
|
|
155
|
# Now create plots
|
156
|
|
157
|
png(paste("RMSE_for_",names_methods[k],out_prefix,".png", sep=""))
|
158
|
boxplot(tb_RMSE[,4:ncol(tb_RMSE)],main=names_methods[k],
|
159
|
ylab= "RMSE", outline=FALSE) #ADD TITLE RELATED TO METHODS...
|
160
|
dev.off()
|
161
|
|
162
|
#boxplot(tb_RMSE[,4:ncol(tb_RMSE)],main=names_methods[k],outline=FALSE) #ADD TITLE RELATED TO METHODS...
|
163
|
png(paste("MAE_for_",names_methods[k],out_prefix,".png", sep=""))
|
164
|
boxplot(tb_MAE[,4:ncol(tb_MAE)],main=names_methods[k],
|
165
|
ylab= "MAE", outline=FALSE) #ADD TITLE RELATED TO METHODS...
|
166
|
dev.off()
|
167
|
|
168
|
#boxplot(tb_RMSE[,4:ncol(tb_RMSE)],main=names_methods[k],outline=FALSE) #ADD TITLE RELATED TO METHODS...
|
169
|
png(paste("ME_for_",names_methods[k],out_prefix,".png", sep=""))
|
170
|
boxplot(tb_ME[,4:ncol(tb_MAE)],main=names_methods[k],
|
171
|
ylab= "ME", outline=FALSE) #ADD TITLE RELATED TO METHODS...
|
172
|
dev.off()
|
173
|
|
174
|
}
|
175
|
|
176
|
names(method_mean)<-obj_names
|
177
|
#Add summary mean graphs!! HERE
|
178
|
|
179
|
write.table(as.data.frame(method_mean$gam_fus_mod1$MAE), "methods_mean_gam_MAE_test1.txt", sep=",")
|
180
|
write.table(as.data.frame(method_mean$fus_CAI$MAE), "methods_mean_fus_CAI_MAE_test1.txt", sep=",")
|
181
|
|
182
|
#### END OF THE SCRIPT
|