Revision e450523b
Added by Benoit Parmentier about 10 years ago
| climate/research/oregon/interpolation/multi_timescales_paper_interpolation_functions.R | ||
|---|---|---|
| 613 | 613 |
|
| 614 | 614 |
### Plotting and computing average MAE per station for different methods |
| 615 | 615 |
plot_MAE_per_station_fun <- function(list_data_v,names_var,interp_method,var_background,stat_loc,out_suffix){
|
| 616 |
#Function to create a series of residuals MAE plots... |
|
| 616 |
#Function to compute residuals MAE per stations |
|
| 617 |
#Plots of maps of MAE per stations with raster background. |
|
| 618 |
#Plots of histograms of MAE per model and method |
|
| 617 | 619 |
|
| 618 | 620 |
mae_fun<-function(x){mean(abs(x))} #Mean Absolute Error give a residuals vector
|
| 619 | 621 |
sd_abs_fun<-function(x){sd(abs(x))} #sd Absolute Error give a residuals vector
|
| ... | ... | |
| 625 | 627 |
data_v_combined <-convert_spdf_to_df_from_list(list_data_v) #long rownames |
| 626 | 628 |
|
| 627 | 629 |
#names_var_all<-c("res_mod1","res_mod2","res_mod3","res_mod4","res_mod5","res_mod6","res_mod7")#,"res_mod8","res_mod9","res_mod10")
|
| 628 |
names_var_all <- res_model_name <- paste("res",names_var,sep="_")
|
|
| 629 |
|
|
| 630 |
res_model_name <- paste("res",names_var,sep="_")
|
|
| 631 |
covar_names_mod <- c("elev_s", "lat", "lon", "E_w", "N_w", "DISTOC", "LC1")
|
|
| 632 |
names_var_all <- c(res_model_name,covar_names_mod) |
|
| 630 | 633 |
t<-melt(data_v_combined, |
| 631 | 634 |
measure=names_var_all, |
| 632 | 635 |
id=c("id"),
|
| ... | ... | |
| 644 | 647 |
data_v_mae<-spTransform(data_v_mae,CRS(CRS_interp)) #Project from WGS84 to new coord. system |
| 645 | 648 |
|
| 646 | 649 |
list_p_mae <- vector("list", length(names_var_all))
|
| 647 |
#names_var <- c("mod1","mod2","mod3","mod7")
|
|
| 650 |
list_p_hist <- vector("list", length(names_var_all))
|
|
| 648 | 651 |
|
| 649 | 652 |
for (k in 1:length(names_var)){
|
| 650 | 653 |
model_name <- names_var[k] |
| ... | ... | |
| 657 | 660 |
na.rm=TRUE) |
| 658 | 661 |
p3 <- p2 + p1 + p2 #to force legend... |
| 659 | 662 |
list_p_mae[[k]] <- p3 |
| 663 |
|
|
| 664 |
#Now add histogram!! |
|
| 665 |
p_hist <-histogram(df_tmp[[res_model_name]], |
|
| 666 |
col=c("grey"),
|
|
| 667 |
ylab=list(label="Percent of total",cex=1.5), |
|
| 668 |
xlab=list(label="Residuals",cex=1.5), |
|
| 669 |
main=list(label=paste("MAE per station for ",model_name," ",interp_method, sep=""),
|
|
| 670 |
cex=1.8), |
|
| 671 |
par.settings = list(axis.text = list(font = 2, cex = 1.3), |
|
| 672 |
par.main.text=list(font=2,cex=2),strip.background=list(col="white")), |
|
| 673 |
par.strip.text=list(font=2,cex=1.5) |
|
| 674 |
) |
|
| 675 |
list_p_hist[[k]] <- p_hist |
|
| 676 |
|
|
| 660 | 677 |
} |
| 661 | 678 |
|
| 662 |
data_mae_obj <- list(list_p_mae,data_v_mae)
|
|
| 663 |
names(data_mae_obj) <- c("list_p_mae","data_v_mae")
|
|
| 679 |
data_mae_obj <- list(list_p_mae,list_p_hist,data_v_mae,data_v_combined)
|
|
| 680 |
names(data_mae_obj) <- c("list_p_mae","list_p_hist","data_v_mae","data_v_combined")
|
|
| 664 | 681 |
return(data_mae_obj) |
| 665 | 682 |
} |
| 666 | 683 |
|
| 684 |
### Plotting and computing average MAE per station for different methods |
|
| 685 |
plot_residuals_map_fun <- function(list_data_v,date_selected,index,names_var,interp_method,var_background){
|
|
| 686 |
|
|
| 687 |
#Function create residuals map from list of station data |
|
| 688 |
#index <- 244 |
|
| 689 |
data_v <- list_data_v[[index]] |
|
| 690 |
names_mod <- names_var |
|
| 691 |
date_proc<-strptime(date_selected, "%Y%m%d") # interpolation date being processed |
|
| 692 |
mo<-as.integer(strftime(date_proc, "%m")) # current month of the date being processed |
|
| 693 |
day<-as.integer(strftime(date_proc, "%d")) |
|
| 694 |
year<-as.integer(strftime(date_proc, "%Y")) |
|
| 695 |
datelabel=format(ISOdate(year,mo,day),"%b %d, %Y") |
|
| 696 |
list_p <- vector("list", length(names_mod))
|
|
| 697 |
for (k in 1:length(names_mod)){
|
|
| 698 |
model_name <- names_mod[k] |
|
| 699 |
res_model_name <- paste("res",model_name,sep="_")
|
|
| 700 |
#cx <- ((data_v[[res_model_name]])*2) |
|
| 701 |
p1 <- levelplot(var_background,#margin=F, |
|
| 702 |
scales = list(draw = FALSE), colorkey = FALSE,par.settings = GrTheme) |
|
| 703 |
p2 <- bubble(data_v,res_model_name, |
|
| 704 |
main=paste("Residuals ",model_name," ",interp_method," ",datelabel,sep=""))
|
|
| 705 |
p3 <- p2 + p1 + p2 #to force legend... |
|
| 706 |
list_p[[k]] <- p3 |
|
| 707 |
} |
|
| 708 |
return(list_p) |
|
| 709 |
} |
|
| 667 | 710 |
|
| 668 | 711 |
################### END OF SCRIPT ################### |
| 669 | 712 |
|
Also available in: Unified diff
revisions2 multi-timescale paper modifications of functions for plotting histogram and bowplots of residuals