Revision 0d15aebd
Added by Benoit Parmentier about 8 years ago
| climate/research/oregon/interpolation/global_product_assessment_part1_functions.R | ||
|---|---|---|
| 4 | 4 |
#Combining tables and figures for individual runs for years and tiles. |
| 5 | 5 |
#AUTHOR: Benoit Parmentier |
| 6 | 6 |
#CREATED ON: 05/24/2016 |
| 7 |
#MODIFIED ON: 09/19/2016
|
|
| 7 |
#MODIFIED ON: 09/20/2016
|
|
| 8 | 8 |
#Version: 1 |
| 9 | 9 |
#PROJECT: Environmental Layers project |
| 10 | 10 |
#COMMENTS: Initial commit, script based on part NASA biodiversity conference |
| ... | ... | |
| 19 | 19 |
# |
| 20 | 20 |
#setfacl -Rmd user:aguzman4:rwx /nobackupp8/bparmen1/output_run10_1500x4500_global_analyses_pred_1992_10052015 |
| 21 | 21 |
|
| 22 |
##COMMIT: plotting function time series, fixing bugs
|
|
| 22 |
##COMMIT: adding windowing and more options to plotting function time series
|
|
| 23 | 23 |
|
| 24 | 24 |
################################################################################################# |
| 25 | 25 |
|
| ... | ... | |
| 789 | 789 |
} |
| 790 | 790 |
|
| 791 | 791 |
|
| 792 |
plot_observation_predictions_time_series <- function(df_pix_time_series,var_name1,var_name2,station_id,scaling=NULL,out_dir=".",out_suffix=""){
|
|
| 792 |
plot_observation_predictions_time_series <- function(df_pix_time_series,var_name1,var_name2,list_windows=NULL,time_series_id=NULL,scaling_factors=NULL,out_dir=".",out_suffix=""){
|
|
| 793 | 793 |
#This function plots time series for observed and predicted points. |
| 794 | 794 |
#This assumes that time series with residuals were extracted from raster stacks. |
| 795 | 795 |
# |
| ... | ... | |
| 797 | 797 |
#1) df_pix_time_series |
| 798 | 798 |
#2) var_name1: name of variables use |
| 799 | 799 |
#3) var_name2: |
| 800 |
#4) station_id |
|
| 801 |
#5) scaling: if null no scaling |
|
| 802 |
#6) out_dir |
|
| 803 |
#7) out_suffix |
|
| 800 |
#4) list_windows: |
|
| 801 |
#5) time_series_id |
|
| 802 |
#6) scaling_factors: if null no scaling |
|
| 803 |
#7) out_dir |
|
| 804 |
#8) out_suffix |
|
| 804 | 805 |
# |
| 805 | 806 |
|
| 806 | 807 |
########## Start script ######### |
| ... | ... | |
| 808 | 809 |
## Screen of NaN and make sure we have NA |
| 809 | 810 |
df_pix_ts[is.na(df_pix_ts)] <- NA |
| 810 | 811 |
|
| 811 |
#Scale if necessary |
|
| 812 |
if(!is.null(scaling)){
|
|
| 812 |
#Scale if necessary, this should be a list
|
|
| 813 |
if(!is.null(scaling_factors)){
|
|
| 813 | 814 |
|
| 814 |
df_pix_ts[[var_name2]] <- df_pix_ts[[var_name2]]*scaling |
|
| 815 |
df_pix_ts[[var_name2]] <- df_pix_ts[[var_name2]]*scaling_factors[2]
|
|
| 815 | 816 |
} |
| 816 | 817 |
|
| 818 |
if(is.null(time_series_id)){
|
|
| 819 |
time_series_id <- unique(na.omit(df_pix_ts$id)) |
|
| 820 |
} |
|
| 821 |
|
|
| 817 | 822 |
#### var_name 1 |
| 818 |
d_z_obs <- zoo(as.numeric(df_pix_ts[[y_var_name]]),as.Date(df_pix_ts$date))
|
|
| 823 |
d_z_obs <- zoo(as.numeric(df_pix_ts[[var_name1]]),as.Date(df_pix_ts$date))
|
|
| 819 | 824 |
#plot(d_z_obs) |
| 820 | 825 |
|
| 821 | 826 |
#### var_name 2 |
| 822 |
d_z_var <- zoo(as.numeric(df_pix_ts[[var_pred_mosaic]]),as.Date(df_pix_ts$date)) #make sure date is a date object !!!
|
|
| 827 |
d_z_var <- zoo(as.numeric(df_pix_ts[[var_name2]]),as.Date(df_pix_ts$date)) #make sure date is a date object !!!
|
|
| 823 | 828 |
#names(d_z_var) <- var_pred_mosaic |
| 824 | 829 |
#plot(d_z_var) |
| 825 | 830 |
|
| ... | ... | |
| 829 | 834 |
#d_z_res <- zoo(as.numeric(df_pix_ts[[paste0("res_",var_pred_mosaic)]]),as.Date(df_pix_ts$date))
|
| 830 | 835 |
#plot(d_z_diff) |
| 831 | 836 |
|
| 832 |
d_z_all <- merge(d_z_var,d_z_obs,d_z_diff)
|
|
| 833 |
names(d_z_all) <- c("pred","obs","diff")
|
|
| 834 |
d_z <- merge(d_z_var,d_z_obs)
|
|
| 835 |
names(d_z) <- c("pred","obs")
|
|
| 837 |
d_z_all <- merge(d_z_obs,d_z_var,d_z_diff)
|
|
| 838 |
names(d_z_all) <- c("obs","pred","diff")
|
|
| 839 |
d_z <- merge(d_z_obs,d_z_var)
|
|
| 840 |
names(d_z) <- c("obs","pred")
|
|
| 836 | 841 |
|
| 837 | 842 |
range(d_z_diff,na.rm=T) |
| 838 | 843 |
mean(d_z_diff,na.rm=T) |
| 839 | 844 |
quantile(d_z_diff,na.rm=T) |
| 845 |
|
|
| 846 |
##Make data.frame with dates for later use!! |
|
| 847 |
#from libary mosaic |
|
| 848 |
|
|
| 849 |
rmse_val <- rmse_fun(d_z_diff) |
|
| 850 |
df_basic_stat <- c(fav_stats(d_z_diff),rmse_val=rmse_val) |
|
| 851 |
|
|
| 852 |
#df_basic_stat <- fav_stats(d_z_all) #does not work on multiple ts |
|
| 853 |
#df_basic_stat$date <- date_processed |
|
| 854 |
#df_basic_stat$reg <- reg |
|
| 855 |
#quantile(data_stations_var_pred$res_mod1,c(1,5,10,90,95,99)) |
|
| 856 |
df_quantile_val <- quantile(na.omit(d_z_diff),c(0,0.01,0.05,0.10,0.45,0.50,0.55,0.90,0.95,0.99,1)) |
|
| 857 |
#names(df_quantile_val) |
|
| 858 |
#as.list(df_quantile_val) |
|
| 859 |
#df_test <- data.frame(names(df_quantile_val))[numeric(0), ] |
|
| 860 |
|
|
| 840 | 861 |
|
| 841 | 862 |
range_dates <- range(as.Date(df_pix_time_series$date)) |
| 842 | 863 |
day_start <- range_dates[1] |
| ... | ... | |
| 848 | 869 |
#range_year <- range(df_pix_ts$year) |
| 849 | 870 |
#start_year <- range_year[1] |
| 850 | 871 |
#end_year <- range_year[2] |
| 851 |
var="tmax" |
|
| 872 |
#var="tmax"
|
|
| 852 | 873 |
|
| 853 | 874 |
res_pix <- 1000 |
| 854 | 875 |
#res_pix <- 480 |
| 855 | 876 |
col_mfrow <- 2 |
| 856 | 877 |
row_mfrow <- 1 |
| 857 | 878 |
|
| 858 |
png_filename <- file.path(out_dir,paste("Figure_5a_time_series_profile_",region_name,"_",id_name,"_",
|
|
| 859 |
y_var_name,"_",start_date,"_",end_date,"_",out_suffix,".png",sep ="")) |
|
| 860 |
title_str <- paste("Daily", var_name1,"and", var_name2,"for", station_id," for the ", start_date,"-",end_date," time period",sep=" ")
|
|
| 861 |
|
|
| 862 |
png(filename=png_filename,width = col_mfrow * res_pix,height = row_mfrow * res_pix) |
|
| 863 |
#this is the whole time series |
|
| 864 |
#plot(d_z_var,ylab="tmax in deg C",xlab="Daily time steps", |
|
| 865 |
# main=title_str,cex=3,font=2, |
|
| 866 |
# cex.main=1.5,cex.lab=1.5,font.lab=2, |
|
| 867 |
# lty=3) |
|
| 868 |
#range(df_pix_ts[[var_pred_mosaic]],na.rm=T) |
|
| 869 |
|
|
| 870 |
#plot(d_z,plot.type="single",col=c("blue","red"))
|
|
| 879 |
png_filename1 <- file.path(out_dir,paste("Figure_5a_time_series_profile_combined_",region_name,"_",time_series_id,"_",
|
|
| 880 |
var_name1,"_",var_name2, |
|
| 881 |
"_",start_date,"_",end_date,"_",out_suffix,".png",sep ="")) |
|
| 882 |
png_filename2 <- file.path(out_dir,paste("Figure_5a_time_series_profile_separate_",region_name,"_",time_series_id,"_",
|
|
| 883 |
var_name1,"_",var_name2, |
|
| 884 |
"_",start_date,"_",end_date,"_",out_suffix,".png",sep ="")) |
|
| 885 |
png_filename3 <- file.path(out_dir,paste("Figure_5a_time_series_profile_separate_with_difference_",region_name,"_",
|
|
| 886 |
time_series_id,"_", var_name1,"_",var_name2, |
|
| 887 |
"_",start_date,"_",end_date,"_",out_suffix,".png",sep ="")) |
|
| 888 |
|
|
| 889 |
title_str <- paste("Daily", var_name1,"and", var_name2,"for", time_series_id,"over the", start_date,"-",end_date,"time period",sep=" ")
|
|
| 871 | 890 |
col_str <- c("blue","red")
|
| 891 |
title_str2 <- paste("Daily", var_name1,"-", var_name2,"and difference","for", time_series_id,"over the", start_date,"-",end_date,"time period",sep=" ")
|
|
| 892 |
col_str2 <- c("blue","red","darkgreen")
|
|
| 893 |
|
|
| 894 |
png(filename=png_filename1,width = col_mfrow * res_pix,height = row_mfrow * res_pix) |
|
| 895 |
#this is the whole time series |
|
| 872 | 896 |
plot(d_z,plot.type="single",col=col_str, |
| 873 | 897 |
ylab="tmax in deg C",xlab="Daily time steps", |
| 874 | 898 |
main=title_str,cex=1,font=2,type="l", |
| 875 | 899 |
cex.main=1.5,cex.lab=1.5,font.lab=2) |
| 876 |
|
|
| 877 |
#plot(df_pix_ts[[var_pred_mosaic]] ~ as.Date(df_pix_ts$date),ylab="tmax in deg C",xlab="Daily time steps", |
|
| 878 |
# main=title_str,cex=1,font=2,type="l", |
|
| 879 |
# cex.main=1.5,cex.lab=1.5,font.lab=2, |
|
| 880 |
# lty=3) |
|
| 881 |
#lines(as.numeric(df_pix_ts[[y_var_name]]) ~ as.Date(df_pix_ts$date),ylab="tmax in deg C",xlab="Daily time steps", |
|
| 882 |
# main=title_str,cex=1,font=2,col="red",type="l", |
|
| 883 |
# cex.main=1.5,cex.lab=1.5,font.lab=2, |
|
| 884 |
# lty=3) |
|
| 900 |
legend("topleft",legend=c(var_name1,var_name2),
|
|
| 901 |
col=col_str,lty=c(1,1), |
|
| 902 |
cex=1.2,bty="n") |
|
| 903 |
dev.off() |
|
| 904 |
|
|
| 905 |
png(filename=png_filename2,width = col_mfrow * res_pix,height = row_mfrow * res_pix) |
|
| 906 |
#this is the whole time series |
|
| 907 |
plot(d_z,plot.type="multiple",col=col_str, |
|
| 908 |
ylab="tmax in deg C",xlab="Daily time steps", |
|
| 909 |
main=title_str,cex=1,font=2,type="l", |
|
| 910 |
cex.main=1.5,cex.lab=1.5,font.lab=2) |
|
| 885 | 911 |
|
| 912 |
legend("topleft",legend=c(var_name1,var_name2),
|
|
| 913 |
col=col_str,lty=c(1,1), |
|
| 914 |
cex=1.2,bty="n") |
|
| 886 | 915 |
dev.off() |
| 916 |
|
|
| 917 |
png(filename=png_filename3,width = col_mfrow * res_pix,height = row_mfrow * res_pix) |
|
| 918 |
#this is the whole time series |
|
| 919 |
plot(d_z_all,plot.type="multiple",col=col_str2, |
|
| 920 |
ylab="tmax in deg C",xlab="Daily time steps", |
|
| 921 |
main=title_str2,cex=1,font=2,type="l", |
|
| 922 |
cex.main=1.5,cex.lab=1.5,font.lab=2) |
|
| 887 | 923 |
|
| 924 |
legend("topleft",legend=c(var_name1,var_name2,"diff"),
|
|
| 925 |
col=col_str2,lty=c(1,1), |
|
| 926 |
cex=1.2,bty="n") |
|
| 927 |
dev.off() |
|
| 928 |
|
|
| 888 | 929 |
#for list_windows go in a loop |
| 889 | 930 |
if(!is.null(list_windows)){
|
| 890 | 931 |
|
| ... | ... | |
| 892 | 933 |
|
| 893 | 934 |
### get smaller window |
| 894 | 935 |
|
| 895 |
window_range <- list_windows[i]
|
|
| 936 |
window_range <- list_windows[[i]]
|
|
| 896 | 937 |
#day_start <- "1999-01-01" #PARAM 12 arg 12 |
| 897 | 938 |
#day_end <- "1999-12-31" #PARAM 13 arg 13 |
| 898 | 939 |
day_start <- window_range[1] |
| 899 |
day_end <- windwow_range[2]
|
|
| 940 |
day_end <- window_range[2] |
|
| 900 | 941 |
|
| 901 | 942 |
start_date <- as.Date(day_start) |
| 902 | 943 |
end_date <- as.Date(day_end) |
| 903 |
start_year <- year(start_date) |
|
| 904 |
end_year <- year(end_date) |
|
| 944 |
#start_year <- year(start_date)
|
|
| 945 |
#end_year <- year(end_date)
|
|
| 905 | 946 |
|
| 906 | 947 |
### now select from the original time series |
| 907 | 948 |
# |
| 908 | 949 |
d_z_all_w <- window(d_z_all,start=start_date,end=end_date) |
| 909 | 950 |
d_z_w <- window(d_z,start=start_date,end=end_date) |
| 910 | 951 |
|
| 952 |
|
|
| 953 |
#### prepare inputs for plots |
|
| 954 |
png_filename4 <- file.path(out_dir,paste("Figure_5b_time_series_profile_window_combined_",region_name,"_",time_series_id,"_",
|
|
| 955 |
var_name1,"_",var_name2, |
|
| 956 |
"_",start_date,"_",end_date,"_",out_suffix,".png",sep ="")) |
|
| 957 |
png_filename5 <- file.path(out_dir,paste("Figure_5b_time_series_profile_window_separate_",region_name,"_",time_series_id,"_",
|
|
| 958 |
var_name1,"_",var_name2, |
|
| 959 |
"_",start_date,"_",end_date,"_",out_suffix,".png",sep ="")) |
|
| 960 |
png_filename6 <- file.path(out_dir,paste("Figure_5b_time_series_profile_window_separate_with_difference_",region_name,"_",
|
|
| 961 |
time_series_id,"_", var_name1,"_",var_name2, |
|
| 962 |
"_",start_date,"_",end_date,"_",out_suffix,".png",sep ="")) |
|
| 963 |
|
|
| 964 |
title_str <- paste("Daily", var_name1,"and", var_name2,"for", time_series_id,"over the", start_date,"-",end_date,"time period",sep=" ")
|
|
| 965 |
col_str <- c("blue","red")
|
|
| 966 |
title_str2 <- paste("Daily", var_name1,"-", var_name2,"and difference","for", time_series_id,"over the", start_date,"-",end_date,"time period",sep=" ")
|
|
| 967 |
col_str2 <- c("blue","red","darkgreen")
|
|
| 968 |
|
|
| 911 | 969 |
res_pix <- 1000 |
| 912 | 970 |
#res_pix <- 480 |
| 913 | 971 |
col_mfrow <- 2 |
| 914 | 972 |
row_mfrow <- 1 |
| 915 |
|
|
| 916 |
#png_filename <- file.path(out_dir,paste("Figure5_b_time_series_profile_",region_name,"_",out_suffix,".png",sep =""))
|
|
| 917 |
#title_str <- paste("Predicted daily ", station_id," ",var," for the ", start_year,"-",end_year," time period",sep="")
|
|
| 918 |
png_filename <- file.path(out_dir,paste("Figure_5_b_time_series_profile_",region_name,"_",id_name,"_",y_var_name,
|
|
| 919 |
"_",start_date,"_",end_date,"_",out_suffix,".png",sep ="")) |
|
| 920 |
title_str <- paste("Daily", var_name1,"and", var_name2,"for", station_id," for the ", start_date,"-",end_date," time period",sep=" ")
|
|
| 921 |
|
|
| 922 |
png(filename=png_filename,width = col_mfrow * res_pix,height = row_mfrow * res_pix) |
|
| 923 | 973 |
|
| 924 |
#add legend later!!
|
|
| 925 |
col_str <- c("blue","red")
|
|
| 974 |
png(filename=png_filename4,width = col_mfrow * res_pix,height = row_mfrow * res_pix)
|
|
| 975 |
#this is the windowed time series
|
|
| 926 | 976 |
plot(d_z_w,plot.type="single",col=col_str, |
| 927 | 977 |
ylab="tmax in deg C",xlab="Daily time steps", |
| 928 | 978 |
main=title_str,cex=1,font=2,type="l", |
| 929 | 979 |
cex.main=1.5,cex.lab=1.5,font.lab=2) |
| 930 |
|
|
| 980 |
legend("topleft",legend=c(var_name1,var_name2),
|
|
| 981 |
col=col_str,lty=c(1,1), |
|
| 982 |
cex=1.2,bty="n") |
|
| 983 |
dev.off() |
|
| 984 |
|
|
| 985 |
png(filename=png_filename5,width = col_mfrow * res_pix,height = row_mfrow * res_pix) |
|
| 986 |
#this is the windowed time series |
|
| 987 |
plot(d_z_w,plot.type="multiple",col=col_str, |
|
| 988 |
ylab="tmax in deg C",xlab="Daily time steps", |
|
| 989 |
main=title_str,cex=1,font=2,type="l", |
|
| 990 |
cex.main=1.5,cex.lab=1.5,font.lab=2) |
|
| 991 |
|
|
| 992 |
legend("topleft",legend=c(var_name1,var_name2),
|
|
| 993 |
col=col_str,lty=c(1,1), |
|
| 994 |
cex=1.2,bty="n") |
|
| 995 |
dev.off() |
|
| 996 |
|
|
| 997 |
png(filename=png_filename6,width = col_mfrow * res_pix,height = row_mfrow * res_pix) |
|
| 998 |
#this is the windowed time series |
|
| 999 |
plot(d_z_all_w,plot.type="multiple",col=col_str2, |
|
| 1000 |
ylab="tmax in deg C",xlab="Daily time steps", |
|
| 1001 |
main=title_str2,cex=1,font=2,type="l", |
|
| 1002 |
cex.main=1.5,cex.lab=1.5,font.lab=2) |
|
| 1003 |
|
|
| 1004 |
legend("topleft",legend=c(var_name1,var_name2,"diff"),
|
|
| 1005 |
col=col_str2,lty=c(1,1), |
|
| 1006 |
cex=1.2,bty="n") |
|
| 931 | 1007 |
dev.off() |
| 932 | 1008 |
|
| 1009 |
list_png_files <- list(png_filename4,png_filename5,png_filename5) |
|
| 933 | 1010 |
} |
| 934 | 1011 |
|
| 935 | 1012 |
} |
| 936 | 1013 |
|
| 937 |
#return() |
|
| 1014 |
plot_time_seris_obj <- list(df_basic_stat,df_quantile_val) |
|
| 1015 |
names(plot_time_seris_obj) <- c("df_basic_stat","df_quantile_val")
|
|
| 1016 |
return(plot_time_seris_obj) |
|
| 938 | 1017 |
} |
| 939 | 1018 |
|
| 1019 |
rmse_fun <- function(x){
|
|
| 1020 |
#x: residuals vector |
|
| 1021 |
rmse_val <-sqrt(mean(x^2,na.rm=T)) |
|
| 1022 |
return(rmse_val) |
|
| 1023 |
} |
|
| 940 | 1024 |
|
| 941 | 1025 |
############################ END OF SCRIPT ################################## |
Also available in: Unified diff
adding windowing and more options to plotting function time series