1
|
######################################## IBS 2013 POSTER #######################################
|
2
|
############################ Scripts for figures and analyses for the the IBS poster #####################################
|
3
|
#This script creates the figures used in the IBS 2013 poster.
|
4
|
#It uses inputs from interpolation objects created at earlier stages... #
|
5
|
#AUTHOR: Benoit Parmentier #
|
6
|
#DATE: 12/27/2012 #
|
7
|
#PROJECT: NCEAS INPLANT: Environment and Organisms --TASK#491-- #
|
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(maptools)
|
20
|
library(graphics)
|
21
|
library(parallel) # Urbanek S. and Ripley B., package for multi cores & parralel processing
|
22
|
library(raster)
|
23
|
library(rasterVis)
|
24
|
library(plotrix) # Draw circle on graph and additional plotting options
|
25
|
library(reshape) # Data format and type transformation
|
26
|
## Functions
|
27
|
#loading R objects that might have similar names
|
28
|
load_obj <- function(f)
|
29
|
{
|
30
|
env <- new.env()
|
31
|
nm <- load(f, env)[1]
|
32
|
env[[nm]]
|
33
|
}
|
34
|
|
35
|
###Parameters and arguments
|
36
|
|
37
|
infile1<- "ghcn_or_tmax_covariates_06262012_OR83M.shp" #GHCN shapefile containing variables for modeling 2010
|
38
|
#infile2<-"list_10_dates_04212012.txt" #List of 10 dates for the regression
|
39
|
infile2<-"list_365_dates_04212012.txt" #list of dates
|
40
|
infile3<-"LST_dates_var_names.txt" #LST dates name
|
41
|
#infile4<-"models_interpolation_05142012.txt" #Interpolation model names
|
42
|
infile5<-"mean_day244_rescaled.rst" #mean LST for day 244
|
43
|
inlistf<-"list_files_05032012.txt" #list of raster images containing the Covariates
|
44
|
infile6<-"OR83M_state_outline.shp"
|
45
|
#stat_loc<-read.table(paste(path,"/","location_study_area_OR_0602012.txt",sep=""),sep=",", header=TRUE)
|
46
|
|
47
|
obj_list<-"list_obj_12272012.txt" #Results of fusion from the run on ATLAS
|
48
|
#obj_list<-"list_obj_08262012.txt" #Results of fusion from the run on ATLAS
|
49
|
path<-"/home/parmentier/Data/IPLANT_project/methods_interpolation_comparison_10242012" #Jupiter LOCATION on Atlas for kriging #Jupiter Location on XANDERS
|
50
|
setwd(path)
|
51
|
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";
|
52
|
#Number of kriging model
|
53
|
out_prefix<-"methods_comp_12272012_" #User defined output prefix
|
54
|
|
55
|
filename<-sub(".shp","",infile1) #Removing the extension from file.
|
56
|
ghcn<-readOGR(".", filename) #reading shapefile
|
57
|
|
58
|
### PREPARING RASTER COVARIATES STACK #######
|
59
|
|
60
|
#CRS<-proj4string(ghcn) #Storing projection information (ellipsoid, datum,etc.)
|
61
|
lines<-read.table(paste(path,"/",inlistf,sep=""), sep="") #Column 1 contains the names of raster files
|
62
|
inlistvar<-lines[,1] #column 3 the list of models to use...?
|
63
|
|
64
|
inlistvar<-paste(path,"/",as.character(inlistvar),sep="")
|
65
|
covar_names<-as.character(lines[,2]) #Column two contains short names for covaraites
|
66
|
|
67
|
s_raster<- stack(inlistvar) #Creating a stack of raster images from the list of variables.
|
68
|
layerNames(s_raster)<-covar_names #Assigning names to the raster layers
|
69
|
projection(s_raster)<-proj_str
|
70
|
|
71
|
#Create mask
|
72
|
pos<-match("LC10",layerNames(s_raster))
|
73
|
LC10<-subset(s_raster,pos)
|
74
|
LC10[is.na(LC10)]<-0 #Since NA values are 0, we assign all zero to NA
|
75
|
mask_land<-LC10<100
|
76
|
mask_land_NA<-mask_land
|
77
|
mask_land_NA[mask_land_NA==0]<-NA
|
78
|
|
79
|
data_name<-"mask_land_OR"
|
80
|
raster_name<-paste(data_name,".rst", sep="")
|
81
|
writeRaster(mask_land, filename=raster_name,overwrite=TRUE) #Writing the data in a raster file format...(IDRISI)
|
82
|
#writeRaster(r2, filename=raster_name,overwrite=TRUE) #Writing the data in a raster file format...(IDRISI)
|
83
|
|
84
|
pos<-match("ELEV_SRTM",layerNames(s_raster))
|
85
|
ELEV_SRTM<-raster(s_raster,pos)
|
86
|
elev<-ELEV_SRTM
|
87
|
elev[-0.050<elev]<-NA #Remove all negative elevation lower than 50 meters...
|
88
|
|
89
|
mask_elev_NA<-elev
|
90
|
|
91
|
pos<-match("mm_01",layerNames(s_raster))
|
92
|
mm_01<-subset(s_raster,pos)
|
93
|
mm_01<-mm_01-273.15
|
94
|
mm_01<-mask(mm_01,mask_land_NA)
|
95
|
#mention this is the last... files
|
96
|
|
97
|
##################### METHODS COMPARISON ###########################
|
98
|
|
99
|
######################################################################
|
100
|
# PART 1 : USING ACCURACY METRICS FOR FIVE METHODS COMPARISON
|
101
|
# Boxplots and histograms
|
102
|
#start function here...
|
103
|
|
104
|
lines<-read.table(paste(path,"/",obj_list,sep=""), sep=",") #Column 1 contains the names RData objects
|
105
|
inlistobj<-lines[,1]
|
106
|
tinlistobj<-paste(path,"/",as.character(inlistobj),sep="")
|
107
|
obj_names<-as.character(lines[,2]) #Column two contains short names for obj. model
|
108
|
|
109
|
tb_metrics_fun<-function(list_obj,path_data,names_obj){
|
110
|
nel<-length(inlistobj)
|
111
|
#method_mod <-vector("list",nel) #list of one row data.frame
|
112
|
method_tb <-vector("list",nel) #list of one row data.frame
|
113
|
for (k in 1:length(inlistobj)){
|
114
|
#obj_tmp<-load_obj(as.character(inlistobj[i]))
|
115
|
#method_mod[[i]]<-obj_tmp
|
116
|
#names(method_mod[[i]])<-obj_names[i]
|
117
|
mod_tmp<-load_obj(as.character(inlistobj[k]))
|
118
|
tb<-mod_tmp[[1]][[3]][0,] #copy of the data.frame structure that holds the acuary metrics
|
119
|
#mod_tmp<-method_mod[[k]]
|
120
|
for (i in 1:365){ # Assuming 365 days of prediction
|
121
|
tmp<-mod_tmp[[i]][[3]]
|
122
|
tb<-rbind(tb,tmp)
|
123
|
}
|
124
|
rm(mod_tmp)
|
125
|
for(i in 4:(ncol(tb))){ # start of the for loop #1
|
126
|
tb[,i]<-as.numeric(as.character(tb[,i]))
|
127
|
}
|
128
|
method_tb[[k]]<-tb
|
129
|
}
|
130
|
names(method_tb)<-names_obj
|
131
|
return(method_tb)
|
132
|
}
|
133
|
|
134
|
tmp44<-tb_metrics_fun(as.character(inlistobj),path,obj_names)
|
135
|
#Condensed, and added other comparison, monthly comparison...:ok
|
136
|
|
137
|
plot_model_boxplot_combined_fun<-function(tb_list,path_data,obj_names,mod_selected,out_prefix,layout_m){
|
138
|
|
139
|
method_stat<-vector("list",length(obj_names)) #This contains summary information based on accuracy metrics (MAE,RMSE)
|
140
|
names_method<-obj_names
|
141
|
metrics<-c("MAE","RMSE")
|
142
|
tb_metric_list<-vector("list",length(metrics))
|
143
|
tb_metric_list_na<-vector("list",length(metrics))
|
144
|
mean_list<-vector("list",length(metrics))
|
145
|
sd_list<-vector("list",length(metrics))
|
146
|
na_mod_list<-vector("list",length(metrics))
|
147
|
|
148
|
for(i in 1:length(metrics)){ # Reorganizing information in terms of metrics
|
149
|
#for(k in 1:length(tb_list)){ # start of the for main loop to all methods
|
150
|
#tb<-tb_list[[k]]
|
151
|
#metrics<-as.character(unique(tb$metric)) #Name of accuracy metrics (RMSE,MAE etc.)
|
152
|
metric_name<-paste("tb_t_",metrics[i],sep="")
|
153
|
png(paste("boxplot",metric_name,out_prefix,"_combined.png", sep="_"),height=480*layout_m[1],width=480*layout_m[2])
|
154
|
par(mfrow=layout_m)
|
155
|
for(k in 1:length(tb_list)){ # start of the for main loop to all methods
|
156
|
#}#for(i in 1:length(metrics)){ # Reorganizing information in terms of metrics
|
157
|
tb<-tb_list[[k]]
|
158
|
#metric_name<-paste("tb_t_",metrics[i],sep="")
|
159
|
tb_metric<-subset(tb, metric==metrics[i])
|
160
|
assign(metric_name,tb_metric)
|
161
|
tb_metric_list[[i]]<-tb_metric
|
162
|
tb_processed<-tb_metric
|
163
|
mod_pat<-glob2rx("mod*")
|
164
|
var_pat<-grep(mod_pat,names(tb_processed),value=FALSE) # using grep with "value" extracts the matching names
|
165
|
#mod_pat<-mod_selected
|
166
|
#var_pat<-grep(mod_pat,names(tb_processed),value=FALSE) # using grep with "value" extracts the matching names
|
167
|
na_mod<-colSums(!is.na(tb_processed[,var_pat]))
|
168
|
for (j in 1:length(na_mod)){
|
169
|
if (na_mod[j]<183){
|
170
|
tmp_name<-names(na_mod)[j]
|
171
|
pos<-match(tmp_name,names(tb_processed))
|
172
|
tb_processed<-tb_processed[,-pos] #Remove columns that have too many missing values!!!
|
173
|
}
|
174
|
}
|
175
|
tb_metric_list_na[[i]]<-tb_processed
|
176
|
mod_pat<-glob2rx("mod*")
|
177
|
var_pat<-grep(mod_pat,names(tb_processed),value=FALSE)
|
178
|
#Plotting box plots
|
179
|
|
180
|
#png(paste("boxplot",metric_name,names_methods[k],out_prefix,".png", sep="_"))
|
181
|
boxplot(tb_processed[,var_pat],main=names_methods[k], ylim=c(1,5),
|
182
|
ylab= metrics[i], outline=FALSE) #ADD TITLE RELATED TO METHODS...
|
183
|
|
184
|
#Add assessment of missing prediction over the year.
|
185
|
mean_metric<-sapply(tb_processed[,var_pat],mean,na.rm=T)
|
186
|
sd_metric<-sapply(tb_processed[,var_pat],sd,na.rm=T)
|
187
|
mean_list[[i]]<-mean_metric
|
188
|
sd_list[[i]]<-sd_metric
|
189
|
na_mod_list[[i]]<-na_mod_list
|
190
|
#Now calculate monthly averages and overall averages over full year
|
191
|
method_stat<-list(mean_list,sd_list,na_mod_list)
|
192
|
method_stat[[k]]<-list(mean_list,sd_list,na_mod_list)
|
193
|
names(method_stat[[k]])<-c("mean_metrics","sd_metrics","na_metrics")
|
194
|
names(mean_list)<-metrics
|
195
|
method_mean[[k]]<-mean_list
|
196
|
names_methods<-obj_names
|
197
|
#names(method_stat)<-obj_names
|
198
|
}
|
199
|
dev.off() #Close file where figures are drawn
|
200
|
}
|
201
|
return(method_stat)
|
202
|
}
|
203
|
|
204
|
tb_list<-tmp44
|
205
|
mod_selected<-""
|
206
|
layout_plot<-c(1,5)
|
207
|
mean_methods<-plot_model_boxplot_fun(tb_list,path,obj_names,mod_selected,out_prefix)
|
208
|
mean_methods_2<-plot_model_boxplot_combined_fun(tb_list,path,obj_names,mod_selected,out_prefix,layout_m=layout_plot)
|
209
|
|
210
|
##################### PART II #######################
|
211
|
|
212
|
##PLOTTING OF ONE DATE TO COMPARE METHODS!!!
|
213
|
|
214
|
lf_raster_fus<-"_365d_GAM_fusion_all_lstd_12272012.rst"
|
215
|
lf_raster_cai<-"_365d_GAM_CAI4_all_12272012.rst"
|
216
|
date_selected<-"20100103"
|
217
|
titles<-list(cai=c("cai mod1","cai mod4","cai mod7"),
|
218
|
fusion=c("fusion mod1"," fusion mod4"," fusion mod7"))
|
219
|
|
220
|
mask_rast<-mask_elev_NA
|
221
|
mod_selected1<-c(1,4,7)
|
222
|
mod_selected2<-c(1,4,7)
|
223
|
#lf_raster_fus<-file_pat1
|
224
|
#lf_raster_cai<-file_pat2
|
225
|
file_pat1<-glob2rx(paste("*tmax_predicted*",date_selected,"*",lf_raster_cai,sep="")) #Search for files in relation to fusion
|
226
|
#lf_cai<-list.files(pattern=file_pat) #Search for files in relation to fusion
|
227
|
file_pat2<-glob2rx(paste("*tmax_predicted*",date_selected,"*",lf_raster_fus,sep="")) #Search for files in relation to fusion
|
228
|
#lf_fus<-list.files(pattern=file_pat) #Search for files in relation to fusion
|
229
|
layout_plot<-c(2,3)
|
230
|
raster_plots_interpolation_fun<-function(file_pat1,file_pat2,mod_selected1,mod_selected2,titles,mask_rast,
|
231
|
layout_m,out_suffix){
|
232
|
layout_m<-layout_plot
|
233
|
lf_cai<-list.files(pattern=file_pat1) #Search for files in relation to fusion
|
234
|
lf_fus<-list.files(pattern=file_pat2) #Search for files in relation to fusion
|
235
|
|
236
|
r1<-stack(lf_cai[mod_selected1]) #CAI
|
237
|
r2<-stack(lf_fus[mod_selected2])#FUS
|
238
|
predictions<-stack(r1,r2)
|
239
|
predictions<-mask(predictions,mask_rast)
|
240
|
layerNames(predictions)<-unlist(titles)
|
241
|
|
242
|
s.range <- c(min(minValue(predictions)), max(maxValue(predictions)))
|
243
|
col.breaks <- pretty(s.range, n=50)
|
244
|
lab.breaks <- pretty(s.range, n=5)
|
245
|
temp.colors <- colorRampPalette(c('blue', 'white', 'red'))
|
246
|
X11(height=6,width=12)
|
247
|
#plot(predictions, breaks=col.breaks, col=rev(heat.colors(length(col.breaks)-1)),
|
248
|
# axis=list(at=lab.breaks, labels=lab.breaks))
|
249
|
plot(predictions, breaks=col.breaks, col=temp.colors(length(col.breaks)-1),
|
250
|
axis=list(at=lab.breaks, labels=lab.breaks))
|
251
|
#plot(reg_outline, add=TRUE)
|
252
|
savePlot(paste("comparison_one_date_CAI_fusion_tmax_prediction_",date_selected,out_prefix,".png", sep=""),type="png")
|
253
|
#png(paste("boxplot",metric_name,out_prefix,"_combined.png", sep="_"),height=480*layout_m[1],width=480*layout_m[2])
|
254
|
#par(mfrow=layout_m)
|
255
|
png(paste("comparison_one_date_CAI_fusion_tmax_prediction_levelplot_",date_selected,out_prefix,".png", sep=""),
|
256
|
height=480*layout_m[1],width=480*layout_m[2])
|
257
|
levelplot(predictions,main="comparison", ylab=NULL,xlab=NULL,par.settings = list(axis.text = list(font = 2, cex = 1.5),
|
258
|
par.main.text=list(font=2,cex=2),strip.background=list(col="white")),par.strip.text=list(font=2,cex=1.5),
|
259
|
#col.regions=temp.colors,at=seq(-1,1,by=0.02))
|
260
|
col.regions=temp.colors(25))
|
261
|
dev.off()
|
262
|
#savePlot(paste("comparison_one_date_CAI_fusion_tmax_prediction_levelplot_",date_selected,out_prefix,".png", sep=""),type="png")
|
263
|
}
|
264
|
|
265
|
raster_plots_interpolation_fun(file_pat1,file_pat2,
|
266
|
mod_selected1,mod_selected2,titles,mask_rast,layout_plot,out_prefix)
|
267
|
|
268
|
|
269
|
#### FIGURE 3: Transect map
|
270
|
|
271
|
### FIGURE 4: transect plot
|
272
|
|
273
|
|
274
|
|
275
|
#### END OF THE SCRIPT #########
|
276
|
|
277
|
|
278
|
#This can be entered as textfile or option later...ok for running now on 12/07/2012
|
279
|
|
280
|
|
281
|
#Figure 1: Boxplots for all methods and models...
|