Revision 18837c79
Added by Benoit Parmentier over 11 years ago
| climate/research/oregon/interpolation/results_interpolation_date_output_analyses.R | ||
|---|---|---|
| 4 | 4 |
#Part 1: Script produces plots for every selected date |
| 5 | 5 |
#Part 2: Examine |
| 6 | 6 |
#AUTHOR: Benoit Parmentier |
| 7 |
#DATE: 05/10/2013
|
|
| 7 |
#DATE: 06/10/2013
|
|
| 8 | 8 |
|
| 9 | 9 |
#PROJECT: NCEAS INPLANT: Environment and Organisms --TASK#???-- |
| 10 | 10 |
|
| 11 | 11 |
################################################################################################## |
| 12 | 12 |
|
| 13 |
## Function(s) used in script |
|
| 13 | 14 |
|
| 14 |
### Parameters and arguments |
|
| 15 |
##Paths to inputs and output |
|
| 16 |
#Select relevant dates and load R objects created during the interpolation step |
|
| 17 |
|
|
| 18 |
##Paths to inputs and output |
|
| 19 |
#script_path<-"/home/parmentier/Data/IPLANT_project/Venezuela_interpolation/Venezuela_01142013/" |
|
| 20 |
#in_path <- "/home/parmentier/Data/IPLANT_project/Venezuela_interpolation/Venezuela_01142013/input_data/" |
|
| 21 |
#out_path<- "/home/parmentier/Data/IPLANT_project/Venezuela_interpolation/Venezuela_01142013/output_data/" |
|
| 22 |
#infile_covar<-"covariates__venezuela_region__VE_01292013.tif" #this is an output from covariate script |
|
| 23 |
#date_selected<-c("20000101") ##This is for year 2000!!!
|
|
| 24 |
#raster_prediction_obj<-load_obj("raster_prediction_obj_dailyTmin_365d_GAM_fus5_all_lstd_03292013.RData")
|
|
| 25 |
#out_prefix<-"_365d_GAM_fus5_all_lstd_03132013" |
|
| 26 |
#out_prefix<-"_365d_GAM_fus5_all_lstd_03142013" #User defined output prefix |
|
| 27 |
#out_prefix<-"_365d_GAM_fus5_all_lstd_03292013" #User defined output prefix |
|
| 28 |
#var<-"TMIN" |
|
| 29 |
#gam_fus_mod<-load_obj("gam_fus_mod_365d_GAM_fus5_all_lstd_02202013.RData")
|
|
| 30 |
#validation_mod_obj<-load_obj("gam_fus_validation_mod_365d_GAM_fus5_all_lstd_02202013.RData")
|
|
| 31 |
#clim_method_mod_obj<-load_obj("gamclim_fus_mod_365d_GAM_fus5_all_lstd_02202013.RData")
|
|
| 15 |
load_obj <- function(f) |
|
| 16 |
{
|
|
| 17 |
env <- new.env() |
|
| 18 |
nm <- load(f, env)[1] |
|
| 19 |
env[[nm]] |
|
| 20 |
} |
|
| 32 | 21 |
|
| 33 |
#rnames<-c("x","y","lon","lat","N","E","N_w","E_w","elev","slope","aspect","CANHEIGHT","DISTOC")
|
|
| 34 |
#lc_names<-c("LC1","LC2","LC3","LC4","LC5","LC6","LC7","LC8","LC9","LC10","LC11","LC12")
|
|
| 35 |
#lst_names<-c("mm_01","mm_02","mm_03","mm_04","mm_05","mm_06","mm_07","mm_08","mm_09","mm_10","mm_11","mm_12",
|
|
| 36 |
# "nobs_01","nobs_02","nobs_03","nobs_04","nobs_05","nobs_06","nobs_07","nobs_08", |
|
| 37 |
# "nobs_09","nobs_10","nobs_11","nobs_12") |
|
| 38 |
#covar_names<-c(rnames,lc_names,lst_names) |
|
| 22 |
#Also used in validation script...place new files for common functions used in interpolation |
|
| 39 | 23 |
|
| 40 |
#list_param_results_analyses<-list(in_path,out_path,script_path,raster_prediction_obj,interpolation_method,infile_covar,covar_names,date_selected,var,out_prefix) |
|
| 41 |
#names(list_param_results_analyses)<-c("in_path","out_path","script_path","raster_prediction_obj", "interpolation_method",
|
|
| 42 |
# "infile_covar","covar_names","date_selected","var","out_prefix") |
|
| 24 |
extract_from_list_obj<-function(obj_list,list_name){
|
|
| 25 |
list_tmp<-vector("list",length(obj_list))
|
|
| 26 |
for (i in 1:length(obj_list)){
|
|
| 27 |
tmp<-obj_list[[i]][[list_name]] #double bracket to return data.frame |
|
| 28 |
list_tmp[[i]]<-tmp |
|
| 29 |
} |
|
| 30 |
tb_list_tmp<-do.call(rbind,list_tmp) #long rownames |
|
| 31 |
return(tb_list_tmp) #this is a data.frame |
|
| 32 |
} |
|
| 43 | 33 |
|
| 44 | 34 |
plots_assessment_by_date<-function(j,list_param){
|
| 45 | 35 |
###Function to assess results from interpolation predictions |
| ... | ... | |
| 70 | 60 |
library(grid) |
| 71 | 61 |
library(lattice) |
| 72 | 62 |
|
| 73 |
## Function(s) used in script |
|
| 74 |
|
|
| 75 |
load_obj <- function(f) |
|
| 76 |
{
|
|
| 77 |
env <- new.env() |
|
| 78 |
nm <- load(f, env)[1] |
|
| 79 |
env[[nm]] |
|
| 80 |
} |
|
| 81 |
|
|
| 82 | 63 |
### BEGIN SCRIPT |
| 83 | 64 |
#Parse input parameters |
| 84 | 65 |
|
| ... | ... | |
| 92 | 73 |
raster_prediction_obj<-list_param$raster_prediction_obj |
| 93 | 74 |
method_mod_obj<-raster_prediction_obj$method_mod_obj |
| 94 | 75 |
validation_mod_obj<-raster_prediction_obj$validation_mod_obj |
| 95 |
clim_method_mod_obj <- raster_prediction_obj$clim_method_mod_obj |
|
| 96 | 76 |
|
| 97 | 77 |
if (var=="TMAX"){
|
| 98 | 78 |
y_var_name<-"dailyTmax" |
| ... | ... | |
| 141 | 121 |
metrics_s<-validation_mod_obj[[index]]$metrics_s |
| 142 | 122 |
data_v<-validation_mod_obj[[index]]$data_v |
| 143 | 123 |
data_s<-validation_mod_obj[[index]]$data_s |
| 144 |
data_month<-clim_method_mod_obj[[index]]$data_month |
|
| 145 |
formulas<-clim_method_mod_obj[[index]]$formulas |
|
| 124 |
formulas<-method_mod_obj[[index]]$formulas |
|
| 146 | 125 |
|
| 126 |
|
|
| 147 | 127 |
#Adding layer LST to the raster stack of covariates |
| 148 | 128 |
#The names of covariates can be changed... |
| 149 | 129 |
|
| ... | ... | |
| 166 | 146 |
rmse<-metrics_v$rmse[nrow(metrics_v)] |
| 167 | 147 |
rmse_f<-metrics_s$rmse[nrow(metrics_s)] |
| 168 | 148 |
|
| 169 |
png(file.path(out_path,paste("LST_",y_var_month,"_scatterplot_",sampling_dat$date,"_",sampling_dat$prop,"_",sampling_dat$run_samp,
|
|
| 170 |
out_prefix,".png", sep=""))) |
|
| 171 |
plot(data_month[[y_var_month]],data_month$LST,xlab=paste("Station mo ",y_var_month,sep=""),ylab=paste("LST mo ",y_var_month,sep=""))
|
|
| 172 |
title(paste("LST vs ", y_var_month,"for",datelabel,sep=" "))
|
|
| 173 |
abline(0,1) |
|
| 174 |
nb_point<-paste("n=",length(data_month[[y_var_month]]),sep="")
|
|
| 175 |
mean_bias<-paste("Mean LST bias= ",format(mean(data_month$LSTD_bias,na.rm=TRUE),digits=3),sep="")
|
|
| 176 |
#Add the number of data points on the plot |
|
| 177 |
legend("topleft",legend=c(mean_bias,nb_point),bty="n")
|
|
| 178 |
dev.off() |
|
| 179 |
|
|
| 180 |
## Figure 2: Daily_tmax_monthly_TMax_scatterplot, modify for TMin!! |
|
| 181 |
|
|
| 182 |
png(file.path(out_path,paste("Month_day_scatterplot_",y_var_name,"_",y_var_month,"_",sampling_dat$date,"_",sampling_dat$prop,"_",sampling_dat$run_samp,
|
|
| 183 |
out_prefix,".png", sep=""))) |
|
| 184 |
plot(data_s[[y_var_name]]~data_s[[y_var_month]],xlab=paste("Month") ,ylab=paste("Daily for",datelabel),main="across stations in VE")
|
|
| 185 |
nb_point<-paste("ns=",length(data_s[[y_var_month]]),sep="")
|
|
| 186 |
nb_point2<-paste("ns_obs=",length(data_s[[y_var_month]])-sum(is.na(data_s[[y_var_name]])),sep="")
|
|
| 187 |
nb_point3<-paste("n_month=",length(data_month[[y_var_month]]),sep="")
|
|
| 188 |
#Add the number of data points on the plot |
|
| 189 |
legend("topleft",legend=c(nb_point,nb_point2,nb_point3),bty="n",cex=0.8)
|
|
| 190 |
dev.off() |
|
| 191 |
|
|
| 192 |
## Figure 3: Predicted_tmax_versus_observed_scatterplot |
|
| 193 |
|
|
| 194 |
#This is for mod_kr!! add other models later... |
|
| 195 |
model_name<-"mod_kr" #can be looped through models later on... |
|
| 196 |
|
|
| 197 |
png(file.path(out_path,paste("Predicted_versus_observed_scatterplot_",y_var_name,"_",model_name,"_",sampling_dat$date,"_",sampling_dat$prop,"_",
|
|
| 198 |
sampling_dat$run_samp,out_prefix,".png", sep=""))) |
|
| 199 |
y_range<-range(c(data_s[[model_name]],data_v[[model_name]]),na.rm=T) |
|
| 200 |
x_range<-range(c(data_s[[y_var_name]],data_v[[y_var_name]]),na.rm=T) |
|
| 201 |
col_t<- c("black","red")
|
|
| 202 |
pch_t<- c(1,2) |
|
| 203 |
plot(data_s[[model_name]],data_s[[y_var_name]], |
|
| 204 |
xlab=paste("Actual daily for",datelabel),ylab="Pred daily",
|
|
| 205 |
ylim=y_range,xlim=x_range,col=col_t[1],pch=pch_t[1]) |
|
| 206 |
points(data_v[[model_name]],data_v[[y_var_name]],col=col_t[2],pch=pch_t[2]) |
|
| 207 |
grid(lwd=0.5, col="black") |
|
| 208 |
abline(0,1) |
|
| 209 |
legend("topleft",legend=c("training","testing"),pch=pch_t,col=col_t,bty="n",cex=0.8)
|
|
| 210 |
title(paste("Predicted_versus_observed_",y_var_name,"_",model_name,"_",datelabel,sep=" "))
|
|
| 211 |
nb_point1<-paste("ns_obs=",length(data_s[[y_var_name]])-sum(is.na(data_s[[model_name]])),sep="")
|
|
| 212 |
nb_point2<-paste("nv_obs=",length(data_v[[y_var_name]])-sum(is.na(data_v[[model_name]])),sep="")
|
|
| 213 |
|
|
| 214 |
rmse_str1<-paste("RMSE= ",format(rmse,digits=3),sep="")
|
|
| 215 |
rmse_str2<-paste("RMSE_f= ",format(rmse_f,digits=3),sep="")
|
|
| 149 |
if (interpolation_method=="gam_CAI" | interpolation_method=="gam_fusion"){
|
|
| 150 |
clim_method_mod_obj <- raster_prediction_obj$clim_method_mod_obj |
|
| 151 |
data_month<-clim_method_mod_obj[[index]]$data_month |
|
| 152 |
|
|
| 153 |
png(file.path(out_path,paste("LST_",y_var_month,"_scatterplot_",sampling_dat$date,"_",sampling_dat$prop,"_",sampling_dat$run_samp,
|
|
| 154 |
out_prefix,".png", sep=""))) |
|
| 155 |
plot(data_month[[y_var_month]],data_month$LST,xlab=paste("Station mo ",y_var_month,sep=""),ylab=paste("LST mo ",y_var_month,sep=""))
|
|
| 156 |
title(paste("LST vs ", y_var_month,"for",datelabel,sep=" "))
|
|
| 157 |
abline(0,1) |
|
| 158 |
nb_point<-paste("n=",length(data_month[[y_var_month]]),sep="")
|
|
| 159 |
mean_bias<-paste("Mean LST bias= ",format(mean(data_month$LSTD_bias,na.rm=TRUE),digits=3),sep="")
|
|
| 160 |
#Add the number of data points on the plot |
|
| 161 |
legend("topleft",legend=c(mean_bias,nb_point),bty="n")
|
|
| 162 |
dev.off() |
|
| 163 |
|
|
| 164 |
## Figure 2: Daily_tmax_monthly_TMax_scatterplot, modify for TMin!! |
|
| 165 |
|
|
| 166 |
png(file.path(out_path,paste("Month_day_scatterplot_",y_var_name,"_",y_var_month,"_",sampling_dat$date,"_",sampling_dat$prop,"_",sampling_dat$run_samp,
|
|
| 167 |
out_prefix,".png", sep=""))) |
|
| 168 |
plot(data_s[[y_var_name]]~data_s[[y_var_month]],xlab=paste("Month") ,ylab=paste("Daily for",datelabel),main="across stations in VE")
|
|
| 169 |
nb_point<-paste("ns=",length(data_s[[y_var_month]]),sep="")
|
|
| 170 |
nb_point2<-paste("ns_obs=",length(data_s[[y_var_month]])-sum(is.na(data_s[[y_var_name]])),sep="")
|
|
| 171 |
nb_point3<-paste("n_month=",length(data_month[[y_var_month]]),sep="")
|
|
| 172 |
#Add the number of data points on the plot |
|
| 173 |
legend("topleft",legend=c(nb_point,nb_point2,nb_point3),bty="n",cex=0.8)
|
|
| 174 |
dev.off() |
|
| 175 |
|
|
| 176 |
## Figure 3: monthly stations used |
|
| 177 |
|
|
| 178 |
png(file.path(out_path,paste("Monthly_data_study_area_", y_var_name,
|
|
| 179 |
out_prefix,".png", sep=""))) |
|
| 180 |
plot(raster(rast_pred_temp,layer=5)) |
|
| 181 |
plot(data_month,col="black",cex=1.2,pch=4,add=TRUE) |
|
| 182 |
title("Monthly ghcn station in Venezuela for January")
|
|
| 183 |
dev.off() |
|
| 216 | 184 |
|
| 217 |
#Add the number of data points on the plot |
|
| 218 |
legend("bottomright",legend=c(nb_point1,nb_point2,rmse_str1,rmse_str2),bty="n",cex=0.8)
|
|
| 219 |
dev.off() |
|
| 185 |
} |
|
| 186 |
## Figure 4: Predicted_tmax_versus_observed_scatterplot |
|
| 187 |
|
|
| 188 |
names_mod <- names(method_mod_obj[[index]][[y_var_name]]) #names of models to plot |
|
| 189 |
#model_name<-"mod_kr" #can be looped through models later on... |
|
| 190 |
|
|
| 191 |
for (k in 1:length(names_mod)){
|
|
| 192 |
model_name <- names_mod[k] |
|
| 193 |
png(file.path(out_path,paste("Predicted_versus_observed_scatterplot_",y_var_name,"_",model_name,"_",sampling_dat$date,"_",sampling_dat$prop,"_",
|
|
| 194 |
sampling_dat$run_samp,out_prefix,".png", sep=""))) |
|
| 195 |
y_range<-range(c(data_s[[model_name]],data_v[[model_name]]),na.rm=T) |
|
| 196 |
x_range<-range(c(data_s[[y_var_name]],data_v[[y_var_name]]),na.rm=T) |
|
| 197 |
col_t<- c("black","red")
|
|
| 198 |
pch_t<- c(1,2) |
|
| 199 |
plot(data_s[[model_name]],data_s[[y_var_name]], |
|
| 200 |
xlab=paste("Actual daily for",datelabel),ylab="Pred daily",
|
|
| 201 |
ylim=y_range,xlim=x_range,col=col_t[1],pch=pch_t[1]) |
|
| 202 |
points(data_v[[model_name]],data_v[[y_var_name]],col=col_t[2],pch=pch_t[2]) |
|
| 203 |
grid(lwd=0.5, col="black") |
|
| 204 |
abline(0,1) |
|
| 205 |
legend("topleft",legend=c("training","testing"),pch=pch_t,col=col_t,bty="n",cex=0.8)
|
|
| 206 |
title(paste("Predicted_versus_observed_",y_var_name,"_",model_name,"_",datelabel,sep=" "))
|
|
| 207 |
nb_point1<-paste("ns_obs=",length(data_s[[y_var_name]])-sum(is.na(data_s[[model_name]])),sep="")
|
|
| 208 |
nb_point2<-paste("nv_obs=",length(data_v[[y_var_name]])-sum(is.na(data_v[[model_name]])),sep="")
|
|
| 209 |
|
|
| 210 |
rmse_str1<-paste("RMSE= ",format(rmse,digits=3),sep="")
|
|
| 211 |
rmse_str2<-paste("RMSE_f= ",format(rmse_f,digits=3),sep="")
|
|
| 212 |
|
|
| 213 |
#Add the number of data points on the plot |
|
| 214 |
legend("bottomright",legend=c(nb_point1,nb_point2,rmse_str1,rmse_str2),bty="n",cex=0.8)
|
|
| 215 |
dev.off() |
|
| 216 |
} |
|
| 220 | 217 |
|
| 221 |
## Figure 4a: prediction raster images
|
|
| 218 |
## Figure 5a: prediction raster images
|
|
| 222 | 219 |
png(file.path(out_path,paste("Raster_prediction_",y_var_name,"_",sampling_dat$date,"_",sampling_dat$prop,"_",sampling_dat$run_samp,
|
| 223 | 220 |
out_prefix,".png", sep=""))) |
| 224 | 221 |
#paste(metrics_v$pred_mod,format(metrics_v$rmse,digits=3),sep=":") |
| ... | ... | |
| 227 | 224 |
levelplot(rast_pred_temp) |
| 228 | 225 |
dev.off() |
| 229 | 226 |
|
| 230 |
## Figure 4b: prediction raster images
|
|
| 227 |
## Figure 5b: prediction raster images
|
|
| 231 | 228 |
png(file.path(out_path,paste("Raster_prediction_plot",sampling_dat$date,"_",sampling_dat$prop,"_",sampling_dat$run_samp,
|
| 232 | 229 |
out_prefix,".png", sep=""))) |
| 233 | 230 |
#paste(metrics_v$pred_mod,format(metrics_v$rmse,digits=3),sep=":") |
| ... | ... | |
| 235 | 232 |
plot(rast_pred_temp) |
| 236 | 233 |
dev.off() |
| 237 | 234 |
|
| 238 |
## Figure 5: training and testing stations used
|
|
| 235 |
## Figure 6: training and testing stations used
|
|
| 239 | 236 |
png(file.path(out_path,paste("Training_testing_stations_map_",y_var_name,"_",sampling_dat$date,"_",sampling_dat$prop,"_",sampling_dat$run_samp,
|
| 240 | 237 |
out_prefix,".png", sep=""))) |
| 241 | 238 |
plot(raster(rast_pred_temp,layer=5)) |
| ... | ... | |
| 246 | 243 |
pch=c(2,1),bty="n") |
| 247 | 244 |
dev.off() |
| 248 | 245 |
|
| 249 |
## Figure 6: monthly stations used |
|
| 250 |
|
|
| 251 |
png(file.path(out_path,paste("Monthly_data_study_area_", y_var_name,
|
|
| 252 |
out_prefix,".png", sep=""))) |
|
| 253 |
plot(raster(rast_pred_temp,layer=5)) |
|
| 254 |
plot(data_month,col="black",cex=1.2,pch=4,add=TRUE) |
|
| 255 |
title("Monthly ghcn station in Venezuela for January")
|
|
| 256 |
dev.off() |
|
| 257 |
|
|
| 258 | 246 |
## Figure 7: delta surface and bias |
| 259 | 247 |
|
| 260 |
if (interpolation_method=="gam_fus"){
|
|
| 248 |
if (interpolation_method=="gam_fusion"){
|
|
| 261 | 249 |
png(file.path(out_path,paste("Bias_delta_surface_",y_var_name,"_",sampling_dat$date[i],"_",sampling_dat$prop[i],
|
| 262 | 250 |
"_",sampling_dat$run_samp[i],out_prefix,".png", sep=""))) |
| 263 | 251 |
|
| ... | ... | |
| 288 | 276 |
|
| 289 | 277 |
#Figure 9: histogram for all images... |
| 290 | 278 |
|
| 279 |
#Write out accuracy information: |
|
| 280 |
|
|
| 281 |
#add sd later... |
|
| 282 |
write.table(tb_diagnostic_v,) |
|
| 283 |
tb_diagnostic_v <- raster_prediction_obj$tb_diagnostic_v |
|
| 284 |
raster_prediction_obj$summary_metrics_v |
|
| 285 |
raster_prediction_obj$summary_month_metrics_v$metric_month_avg |
|
| 286 |
raster_prediction_obj$summary_month_metrics_v$metric_month_sd |
|
| 287 |
|
|
| 288 |
#write.table(tb_diagnostic_v, file= file.path(out_path,interpolation_method,"_tb_diagnostic_v",out_prefix,".txt",sep=""), sep=",",overwrite=FALSE) |
|
| 289 |
write.table(tb, file= paste(path,"/","results2_gwr_Assessment_measure_all",out_prefix,".txt",sep=""), sep=",") |
|
| 290 |
|
|
| 291 | 291 |
#histogram(rast_pred_temp) |
| 292 | 292 |
list_output_analyses<-list(metrics_s,metrics_v) |
| 293 | 293 |
names(list_output_analyses) <- c("metrics_s", "metrics_v")
|
Also available in: Unified diff
results output script, modifications to output graphs for added methods: gam_daily, gwr_daily, kriging_daily