Revision 57893002
Added by Benoit Parmentier about 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: 06/10/2013
|
|
| 7 |
#DATE: 08/05/2013
|
|
| 8 | 8 |
|
| 9 | 9 |
#PROJECT: NCEAS INPLANT: Environment and Organisms --TASK#???-- |
| 10 | 10 |
|
| ... | ... | |
| 34 | 34 |
plots_assessment_by_date<-function(j,list_param){
|
| 35 | 35 |
###Function to assess results from interpolation predictions |
| 36 | 36 |
#AUTHOR: Benoit Parmentier |
| 37 |
#DATE: 05/10/2013
|
|
| 37 |
#DATE: 08/05/2013
|
|
| 38 | 38 |
#PROJECT: NCEAS INPLANT: Environment and Organisms --TASK#363-- |
| 39 | 39 |
|
| 40 | 40 |
#1) in_path |
| ... | ... | |
| 42 | 42 |
#3) script_path |
| 43 | 43 |
#4) raster_prediction_obj |
| 44 | 44 |
#5) interpolation_method |
| 45 |
#6) infile_covariates |
|
| 46 |
#7) covar_names |
|
| 45 |
#7) covar_obj: covariates object contains file name and names of covariates |
|
| 47 | 46 |
#8) date_selected_results |
| 48 | 47 |
#9) var |
| 49 | 48 |
#10) out_prefix |
| ... | ... | |
| 67 | 66 |
var<-list_param$var #variable being interpolated |
| 68 | 67 |
out_path <- list_param$out_path |
| 69 | 68 |
interpolation_method <- list_param$interpolation_method |
| 70 |
infile_covariates <- list_param$infile_covariates |
|
| 71 |
covar_names<-list_param$covar_names |
|
| 69 |
infile_covariates <- list_param$covar_obj$infile_covariates
|
|
| 70 |
covar_names<-list_param$covar_obj$covar_names
|
|
| 72 | 71 |
|
| 73 | 72 |
raster_prediction_obj<-list_param$raster_prediction_obj |
| 74 | 73 |
method_mod_obj<-raster_prediction_obj$method_mod_obj |
| 75 | 74 |
validation_mod_obj<-raster_prediction_obj$validation_mod_obj |
| 76 | 75 |
|
| 76 |
#This should not be set here...? master script |
|
| 77 | 77 |
if (var=="TMAX"){
|
| 78 | 78 |
y_var_name<-"dailyTmax" |
| 79 | 79 |
y_var_month<-"TMax" |
| ... | ... | |
| 111 | 111 |
|
| 112 | 112 |
#Get raster stack of interpolated surfaces |
| 113 | 113 |
index<-i_dates[[j]] |
| 114 |
pred_temp<-as.character(method_mod_obj[[index]][[y_var_name]]) #list of files with path included |
|
| 114 |
pred_temp<-as.character(method_mod_obj[[index]][[y_var_name]]) #list of daily prediction files with path included
|
|
| 115 | 115 |
rast_pred_temp_s <-stack(pred_temp) #stack of temperature predictions from models (daily) |
| 116 | 116 |
rast_pred_temp <-mask(rast_pred_temp_s,LC_mask,file=file.path(out_path,"test.tif"),overwrite=TRUE) |
| 117 | 117 |
|
| ... | ... | |
| 146 | 146 |
rmse<-metrics_v$rmse[nrow(metrics_v)] |
| 147 | 147 |
rmse_f<-metrics_s$rmse[nrow(metrics_s)] |
| 148 | 148 |
|
| 149 |
if (interpolation_method=="gam_CAI" | interpolation_method=="gam_fusion"){
|
|
| 149 |
#Set as constant in master script ?? : c("gam_CAI","kriging_CAI","gwr_CAI","gam_fusion","kriging_fusion","gwr_fusion")
|
|
| 150 |
if (interpolation_method %in% c("gam_CAI","kriging_CAI","gwr_CAI","gam_fusion","kriging_fusion","gwr_fusion")){
|
|
| 151 |
#if multi-time scale method then the raster prediction object contains a "clim_method_mod_obj" |
|
| 150 | 152 |
clim_method_mod_obj <- raster_prediction_obj$clim_method_mod_obj |
| 151 | 153 |
data_month<-clim_method_mod_obj[[index]]$data_month |
| 152 | 154 |
|
| ... | ... | |
| 156 | 158 |
title(paste("LST vs ", y_var_month,"for",datelabel,sep=" "))
|
| 157 | 159 |
abline(0,1) |
| 158 | 160 |
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="")
|
|
| 161 |
LSTD_bias <- data_month$TMax - data_month$LST #in case it is a CAI method, calculate bias |
|
| 162 |
mean_bias<-paste("Mean LST bias= ",format(mean(LSTD_bias,na.rm=TRUE),digits=3),sep="")
|
|
| 160 | 163 |
#Add the number of data points on the plot |
| 161 | 164 |
legend("topleft",legend=c(mean_bias,nb_point),bty="n")
|
| 162 | 165 |
dev.off() |
| ... | ... | |
| 221 | 224 |
#paste(metrics_v$pred_mod,format(metrics_v$rmse,digits=3),sep=":") |
| 222 | 225 |
names(rast_pred_temp)<-paste(metrics_v$pred_mod,format(metrics_v$rmse,digits=3),sep=":") |
| 223 | 226 |
#plot(rast_pred_temp) |
| 224 |
levelplot(rast_pred_temp) |
|
| 227 |
levelplot(rast_pred_temp) #not working...takes too long to plot?
|
|
| 225 | 228 |
dev.off() |
| 226 | 229 |
|
| 227 | 230 |
## Figure 5b: prediction raster images |
| ... | ... | |
| 245 | 248 |
|
| 246 | 249 |
## Figure 7: delta surface and bias |
| 247 | 250 |
|
| 248 |
if (interpolation_method=="gam_fusion"){
|
|
| 251 |
if (interpolation_method%in% c("gam_fusion","kriging_fusion","gwr_fusion")){
|
|
| 249 | 252 |
png(file.path(out_path,paste("Bias_delta_surface_",y_var_name,"_",sampling_dat$date[i],"_",sampling_dat$prop[i],
|
| 250 | 253 |
"_",sampling_dat$run_samp[i],out_prefix,".png", sep=""))) |
| 251 | 254 |
|
| ... | ... | |
| 253 | 256 |
delta_rast<-raster(method_mod_obj[[index]]$delta) #only one delta image!!! |
| 254 | 257 |
names(delta_rast)<-"delta" |
| 255 | 258 |
rast_temp_date<-stack(bias_rast,delta_rast) |
| 259 |
layers_names <- names(rast_temp_date) |
|
| 256 | 260 |
rast_temp_date<-mask(rast_temp_date,LC_mask,file=file.path(out_path,"test.tif"),overwrite=TRUE) |
| 257 | 261 |
#bias_d_rast<-raster("fusion_bias_LST_20100103_30_1_10d_GAM_fus5_all_lstd_02082013.rst")
|
| 262 |
names(rast_temp_date) <-layers_names |
|
| 258 | 263 |
plot(rast_temp_date) |
| 259 | 264 |
dev.off() |
| 260 | 265 |
} |
| 261 | 266 |
|
| 262 |
if (interpolation_method=="gam_CAI"){
|
|
| 267 |
if (interpolation_method %in% c("gam_CAI","kriging_CAI","gwr_CAI","gam_fusion","kriging_fusion","gwr_fusion")){
|
|
| 263 | 268 |
png(file.path(out_path,paste("clim_surface_",y_var_name,"_",sampling_dat$date[i],"_",sampling_dat$prop[i],
|
| 264 | 269 |
"_",sampling_dat$run_samp[i],out_prefix,".png", sep=""))) |
| 265 | 270 |
|
| 266 | 271 |
clim_rast<-stack(clim_method_mod_obj[[index]]$clim) |
| 267 | 272 |
delta_rast<-raster(method_mod_obj[[index]]$delta) #only one delta image!!! |
| 268 | 273 |
names(delta_rast)<-"delta" |
| 274 |
layers_names <- c(names(clim_rast),"delta") |
|
| 269 | 275 |
rast_temp_date<-stack(clim_rast,delta_rast) |
| 270 |
rast_temp_date<-mask(rast_temp_date,LC_mask,file=file.path(out_path,"test.tif"),overwrite=TRUE) |
|
| 276 |
rast_temp_date<-mask(rast_temp_date,LC_mask,file=file.path(out_path,"test.tif"),overwrite=TRUE) #loosing names here
|
|
| 271 | 277 |
#bias_d_rast<-raster("fusion_bias_LST_20100103_30_1_10d_GAM_fus5_all_lstd_02082013.rst")
|
| 278 |
names(rast_temp_date) <-layers_names |
|
| 272 | 279 |
plot(rast_temp_date) |
| 273 | 280 |
|
| 274 | 281 |
dev.off() |
| 275 | 282 |
} |
| 276 | 283 |
|
| 277 | 284 |
#Figure 9: histogram for all images... |
| 285 |
## Add later...? distance to closest fitting station? |
|
| 286 |
|
|
| 287 |
#tb_diagnostic_v <- raster_prediction_obj$tb_diagnostic_v |
|
| 288 |
#raster_prediction_obj$summary_metrics_v |
|
| 289 |
#raster_prediction_obj$summary_month_metrics_v$metric_month_avg |
|
| 290 |
#raster_prediction_obj$summary_month_metrics_v$metric_month_sd |
|
| 278 | 291 |
|
| 279 | 292 |
#Write out accuracy information: |
| 280 | 293 |
|
| 281 | 294 |
#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 |
|
|
| 295 |
#write.table(tb_diagnostic_v,) |
|
| 288 | 296 |
#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=",") |
|
| 297 |
#write.table(tb, file= paste(path,"/","results2_gwr_Assessment_measure_all",out_prefix,".txt",sep=""), sep=",")
|
|
| 290 | 298 |
|
| 291 | 299 |
#histogram(rast_pred_temp) |
| 292 | 300 |
list_output_analyses<-list(metrics_s,metrics_v) |
Also available in: Unified diff
results interpolation script, fixing bugs and slight changes