Revision f23e6e04
Added by Benoit Parmentier over 11 years ago
| climate/research/oregon/interpolation/GAM_fusion_function_multisampling_validation_metrics.R | ||
|---|---|---|
| 2 | 2 |
|
| 3 | 3 |
#The interpolation is done first at the monthly add delta. |
| 4 | 4 |
#AUTHOR: Benoit Parmentier |
| 5 |
#DATE: 05/01/2013
|
|
| 5 |
#DATE: 05/06/2013
|
|
| 6 | 6 |
|
| 7 | 7 |
#Change this to allow explicitly arguments... |
| 8 | 8 |
#Arguments: |
| ... | ... | |
| 71 | 71 |
############### BEGIN SCRIPT ########### |
| 72 | 72 |
|
| 73 | 73 |
#PARSING INPUT PARAMETERS |
| 74 |
out_path <- list_param$out_path |
|
| 74 | 75 |
day_list<- list_param$rast_day_year_list[[i]] |
| 75 | 76 |
#day_list <-rast_day_yearlist[[i]] #list of prediction for the current date... |
| 76 | 77 |
names_mod<-names(day_list) |
| ... | ... | |
| 103 | 104 |
run_info<-cbind(sampling_dat_day,n=nv) |
| 104 | 105 |
run_info[rep(seq_len(nrow(run_info)), each=N),] #repeating same row n times |
| 105 | 106 |
metrics_v_df<-cbind(metrics_v_obj$metrics,run_info) |
| 107 |
metrics_v_df["var_interp"]<-rep(y_var_name,times=nrow(metrics_v_df)) |
|
| 108 |
#Name of the variable interpolated, useful for cross-comparison between methods at later stages |
|
| 106 | 109 |
|
| 107 | 110 |
run_info<-cbind(sampling_dat_day,n=ns) |
| 108 | 111 |
run_info[rep(seq_len(nrow(run_info)), each=N),] |
| 109 | 112 |
metrics_s_df<-cbind(metrics_s_obj$metrics,run_info) |
| 113 |
metrics_s_df["var_interp"]<-rep(y_var_name,times=nrow(metrics_s_df)) |
|
| 114 |
#Name of the variable interpolated, useful for cross-comparison between methods at later stages |
|
| 110 | 115 |
|
| 111 | 116 |
data_v<-spCbind(data_v,metrics_v_obj$residuals) |
| 112 | 117 |
data_s<-spCbind(data_s,metrics_s_obj$residuals) |
| ... | ... | |
| 131 | 136 |
|
| 132 | 137 |
#### Function to plot boxplot from data.frame table of accuracy metrics |
| 133 | 138 |
|
| 134 |
boxplot_from_tb <-function(tb_diagnostic,metric_names,out_prefix){
|
|
| 139 |
boxplot_from_tb <-function(tb_diagnostic,metric_names,out_prefix,out_path){
|
|
| 135 | 140 |
#now boxplots and mean per models |
| 136 | 141 |
library(gdata) #Nesssary to use cbindX |
| 137 | 142 |
|
| 138 | 143 |
### Start script |
| 144 |
y_var_name<-unique(tb_diagnostic$var_interp) #extract the name of interpolated variable: dailyTmax, dailyTmin |
|
| 139 | 145 |
|
| 140 | 146 |
mod_names<-sort(unique(tb_diagnostic$pred_mod)) #models that have accuracy metrics |
| 141 | 147 |
t<-melt(tb_diagnostic, |
| 142 | 148 |
#measure=mod_var, |
| 143 | 149 |
id=c("date","pred_mod","prop"),
|
| 144 | 150 |
na.rm=F) |
| 151 |
t$value<-as.numeric(t$value) #problem with char!!! |
|
| 145 | 152 |
avg_tb<-cast(t,pred_mod~variable,mean) |
| 146 |
|
|
| 153 |
avg_tb$var_interp<-rep(y_var_name,times=nrow(avg_tb)) |
|
| 147 | 154 |
median_tb<-cast(t,pred_mod~variable,median) |
| 155 |
|
|
| 156 |
#avg_tb<-cast(t,pred_mod~variable,mean) |
|
| 148 | 157 |
tb<-tb_diagnostic |
| 149 | 158 |
|
| 150 | 159 |
#mod_names<-sort(unique(tb$pred_mod)) #kept for clarity |
| ... | ... | |
| 162 | 171 |
mod_var<-grep(mod_pat,names(mod_metrics),value=TRUE) # using grep with "value" extracts the matching names |
| 163 | 172 |
#browser() |
| 164 | 173 |
test<-mod_metrics[mod_var] |
| 165 |
png(paste("boxplot_metric_",metric_ac, out_prefix,".png", sep=""))
|
|
| 174 |
png(file.path(out_path,paste("boxplot_metric_",metric_ac, out_prefix,".png", sep="")))
|
|
| 175 |
#boxplot(test,outline=FALSE,horizontal=FALSE,cex=0.5, |
|
| 176 |
# ylab=paste(metric_ac,"in degree C",sep=" ")) |
|
| 177 |
|
|
| 166 | 178 |
boxplot(test,outline=FALSE,horizontal=FALSE,cex=0.5, |
| 167 |
ylab=paste(metric_ac,"in degree C",sep=" ")) |
|
| 179 |
ylab=paste(metric_ac,"in degree C",sep=" "),axisnames=FALSE,axes=FALSE) |
|
| 180 |
axis(1, labels = FALSE) |
|
| 181 |
## Create some text labels |
|
| 182 |
labels <- labels<- names(test) |
|
| 183 |
## Plot x axis labels at default tick marks |
|
| 184 |
text(1:ncol(test), par("usr")[3] - 0.25, srt = 45, adj = 1,
|
|
| 185 |
labels = labels, xpd = TRUE) |
|
| 186 |
axis(2) |
|
| 187 |
box() |
|
| 168 | 188 |
#legend("bottomleft",legend=paste(names(rows_total),":",rows_total,sep=""),cex=0.7,bty="n")
|
| 169 | 189 |
#title(as.character(t(paste(t(names(rows_total)),":",rows_total,sep=""))),cex=0.8) |
| 170 | 190 |
title(paste(metric_ac,"for",y_var_name,sep=" "),cex=0.8) |
| ... | ... | |
| 174 | 194 |
avg_tb$n<-rows_total #total number of predictions on which the mean is based |
| 175 | 195 |
median_tb$n<-rows_total |
| 176 | 196 |
summary_obj<-list(avg_tb,median_tb) |
| 197 |
names(summary_obj)<-c("avg","median")
|
|
| 177 | 198 |
return(summary_obj) |
| 178 | 199 |
} |
| 179 |
#boxplot_month_from_tb(tb_diagnostic,metric_names,out_prefix) |
|
| 200 |
#boxplot_month_from_tb(tb_diagnostic,metric_names,out_prefix,out_path)
|
|
| 180 | 201 |
## Function to display metrics by months/seasons |
| 181 |
boxplot_month_from_tb <-function(tb_diagnostic,metric_names,out_prefix){
|
|
| 202 |
boxplot_month_from_tb <-function(tb_diagnostic,metric_names,out_prefix,out_path){
|
|
| 182 | 203 |
|
| 183 | 204 |
#Generate boxplot per month for models and accuracy metrics |
| 184 | 205 |
#Input parameters: |
| ... | ... | |
| 189 | 210 |
|
| 190 | 211 |
################# |
| 191 | 212 |
## BEGIN |
| 192 |
|
|
| 213 |
y_var_name<-unique(tb_diagnostic$var_interp) #extract the name of interpolated variable: dailyTmax, dailyTmin |
|
| 193 | 214 |
date_f<-strptime(tb_diagnostic$date, "%Y%m%d") # interpolation date being processed |
| 194 | 215 |
tb_diagnostic$month<-strftime(date_f, "%m") # current month of the date being processed |
| 195 | 216 |
mod_names<-sort(unique(tb_diagnostic$pred_mod)) #models that have accuracy metrics |
| ... | ... | |
| 199 | 220 |
#measure=mod_var, |
| 200 | 221 |
id=c("date","pred_mod","prop","month"),
|
| 201 | 222 |
na.rm=F) |
| 223 |
t$value<-as.numeric(t$value) #problem with char!!! |
|
| 202 | 224 |
tb_mod_m_avg <-cast(t,pred_mod+month~variable,mean) #monthly mean for every model |
| 225 |
tb_mod_m_avg$var_interp<-rep(y_var_name,times=nrow(tb_mod_m_avg)) |
|
| 226 |
|
|
| 203 | 227 |
tb_mod_m_sd <-cast(t,pred_mod+month~variable,sd) #monthly sd for every model |
| 204 | 228 |
|
| 205 |
tb_mod_m_list <-lapply(mod_names, function(k) subset(tb_mod_m, pred_mod==k)) #this creates a list of 5 based on models names |
|
| 206 |
|
|
| 229 |
tb_mod_m_list <-lapply(mod_names, function(k) subset(tb_mod_m_avg, pred_mod==k)) #this creates a list of 5 based on models names
|
|
| 230 |
|
|
| 207 | 231 |
for (k in 1:length(mod_names)){
|
| 208 | 232 |
mod_metrics <-tb_mod_list[[k]] |
| 209 | 233 |
current_mod_name<- mod_names[k] |
| ... | ... | |
| 211 | 235 |
metric_ac<-metric_names[j] |
| 212 | 236 |
col_selected<-c(metric_ac,"month") |
| 213 | 237 |
test<-mod_metrics[col_selected] |
| 214 |
png(paste("boxplot_metric_",metric_ac,"_",current_mod_name,"_by_month_",out_prefix,".png", sep=""))
|
|
| 238 |
png(file.path(out_path,paste("boxplot_metric_",metric_ac,"_",current_mod_name,"_by_month_",out_prefix,".png", sep="")))
|
|
| 215 | 239 |
boxplot(test[[metric_ac]]~test[[c("month")]],outline=FALSE,horizontal=FALSE,cex=0.5,
|
| 216 |
ylab=paste(metric_ac,"in degree C",sep=" ")) |
|
| 240 |
ylab=paste(metric_ac,"in degree C",sep=" "),,axisnames=FALSE,axes=FALSE) |
|
| 241 |
#boxplot(test[[metric_ac]]~test[[c("month")]],outline=FALSE,horizontal=FALSE,cex=0.5,
|
|
| 242 |
# ylab=paste(metric_ac,"in degree C",sep=" ")) |
|
| 243 |
axis(1, labels = FALSE) |
|
| 244 |
## Create some text labels |
|
| 245 |
labels <- month.abb # abbreviated names for each month |
|
| 246 |
## Plot x axis labels at default tick marks |
|
| 247 |
text(1:length(labels), par("usr")[3] - 0.25, srt = 45, adj = 1,
|
|
| 248 |
labels = labels, xpd = TRUE) |
|
| 249 |
axis(2) |
|
| 250 |
box() |
|
| 217 | 251 |
#legend("bottomleft",legend=paste(names(rows_total),":",rows_total,sep=""),cex=0.7,bty="n")
|
| 218 | 252 |
title(paste(metric_ac,"for",current_mod_name,"by month",sep=" ")) |
| 219 | 253 |
dev.off() |
Also available in: Unified diff
validation script, modifications to store results in output folder and additional minor changes