1 |
170edade
|
Benoit Parmentier
|
##################################### 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
|