Revision 03b5035b
Added by Benoit Parmentier over 8 years ago
| climate/research/oregon/interpolation/global_product_assessment_part1_functions.R | ||
|---|---|---|
| 19 | 19 |
# |
| 20 | 20 |
#setfacl -Rmd user:aguzman4:rwx /nobackupp8/bparmen1/output_run10_1500x4500_global_analyses_pred_1992_10052015 |
| 21 | 21 |
|
| 22 |
##COMMIT: function for combining extracted values and |
|
| 22 |
##COMMIT: plotting function time series, fixing bugs |
|
| 23 |
|
|
| 23 | 24 |
################################################################################################# |
| 24 | 25 |
|
| 25 | 26 |
### Loading R library and packages |
| ... | ... | |
| 810 | 811 |
#Scale if necessary |
| 811 | 812 |
if(!is.null(scaling)){
|
| 812 | 813 |
|
| 813 |
df_pix_ts[[var_pred_mosaic]] <- df_pix_ts[[var_pred_mosaic]]*scaling
|
|
| 814 |
df_pix_ts[[var_name2]] <- df_pix_ts[[var_name2]]*scaling
|
|
| 814 | 815 |
} |
| 815 | 816 |
|
| 816 | 817 |
#### var_name 1 |
| 817 | 818 |
d_z_obs <- zoo(as.numeric(df_pix_ts[[y_var_name]]),as.Date(df_pix_ts$date)) |
| 818 |
plot(d_z_obs) |
|
| 819 |
#plot(d_z_obs)
|
|
| 819 | 820 |
|
| 820 | 821 |
#### var_name 2 |
| 821 | 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 !!! |
| 822 | 823 |
#names(d_z_var) <- var_pred_mosaic |
| 823 |
plot(d_z_var) |
|
| 824 |
#plot(d_z_var)
|
|
| 824 | 825 |
|
| 825 |
d_z_diff <- d_z_var - d_z_obs |
|
| 826 |
d_z_diff <- d_z_var - d_z_obs #this is residuals if var is predicted
|
|
| 826 | 827 |
#d_z_res <- d_z_var - d_z_obs |
| 827 | 828 |
|
| 828 | 829 |
#d_z_res <- zoo(as.numeric(df_pix_ts[[paste0("res_",var_pred_mosaic)]]),as.Date(df_pix_ts$date))
|
| 829 |
plot(d_z_res)
|
|
| 830 |
#plot(d_z_diff)
|
|
| 830 | 831 |
|
| 831 |
d_z_all <- merge(d_z_var,d_z_obs,d_z_res)
|
|
| 832 |
names(d_z_all) <- c("pred","obs","res")
|
|
| 832 |
d_z_all <- merge(d_z_var,d_z_obs,d_z_diff)
|
|
| 833 |
names(d_z_all) <- c("pred","obs","diff")
|
|
| 833 | 834 |
d_z <- merge(d_z_var,d_z_obs) |
| 834 | 835 |
names(d_z) <- c("pred","obs")
|
| 835 | 836 |
|
| 836 |
range(d_z_res,na.rm=T)
|
|
| 837 |
mean(d_z_res,na.rm=T)
|
|
| 838 |
quantile(d_z_res,na.rm=T)
|
|
| 837 |
range(d_z_diff,na.rm=T)
|
|
| 838 |
mean(d_z_diff,na.rm=T)
|
|
| 839 |
quantile(d_z_diff,na.rm=T)
|
|
| 839 | 840 |
|
| 840 | 841 |
range_dates <- range(as.Date(df_pix_time_series$date)) |
| 842 |
day_start <- range_dates[1] |
|
| 843 |
day_end <- range_dates[2] |
|
| 841 | 844 |
#day_start <- "1984-01-01" #PARAM 12 arg 12 |
| 842 | 845 |
#day_end <- "2014-12-31" #PARAM 13 arg 13 |
| 843 |
#start_date <- as.Date(day_start) |
|
| 844 |
#end_date <- as.Date(day_end) |
|
| 845 |
start_year <- year(range_dates[1]) |
|
| 846 |
end_year <- year(range_dates[2]) |
|
| 846 |
start_date <- as.Date(day_start) |
|
| 847 |
end_date <- as.Date(day_end) |
|
| 847 | 848 |
#range_year <- range(df_pix_ts$year) |
| 848 | 849 |
#start_year <- range_year[1] |
| 849 | 850 |
#end_year <- range_year[2] |
| ... | ... | |
| 854 | 855 |
col_mfrow <- 2 |
| 855 | 856 |
row_mfrow <- 1 |
| 856 | 857 |
|
| 857 |
png_filename <- file.path(out_dir,paste("Figure_5a_time_series_profile_",region_name,"_",id_name,y_var_name,out_suffix,".png",sep =""))
|
|
| 858 |
title_str <- paste("Daily", var_name1,"and", var_name2,"for", station_id," for the ", start_year,"-",end_year," time period",sep=" ")
|
|
| 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=" ")
|
|
| 859 | 861 |
|
| 860 | 862 |
png(filename=png_filename,width = col_mfrow * res_pix,height = row_mfrow * res_pix) |
| 861 | 863 |
#this is the whole time series |
| ... | ... | |
| 888 | 890 |
|
| 889 | 891 |
for(i in 1:length(list_windows)){
|
| 890 | 892 |
|
| 891 |
|
|
| 892 | 893 |
### get smaller window |
| 893 | 894 |
|
| 894 | 895 |
window_range <- list_windows[i] |
| 895 |
#day_start <- "1994-01-01" #PARAM 12 arg 12
|
|
| 896 |
#day_start <- "1999-01-01" #PARAM 12 arg 12
|
|
| 896 | 897 |
#day_end <- "1999-12-31" #PARAM 13 arg 13 |
| 897 | 898 |
day_start <- window_range[1] |
| 898 | 899 |
day_end <- windwow_range[2] |
| ... | ... | |
| 904 | 905 |
|
| 905 | 906 |
### now select from the original time series |
| 906 | 907 |
# |
| 907 |
d_z_w <- window(d_z_tmp,start=start_date,end=end_date) |
|
| 908 |
names(d_z) <- var_pred_mosaic |
|
| 909 |
df_var <- as.data.frame(d_z) |
|
| 910 |
|
|
| 911 |
d_z_tmp <- d_z_obs |
|
| 912 |
d_z_tmp <- d_z |
|
| 913 |
d_z_tmp <- window(d_z_tmp,start=start_date,end=end_date) |
|
| 914 |
plot(d_z_tmp,col=c("blue","red"),plot.type="single")
|
|
| 915 |
df_z <- as.data.frame(d_z) |
|
| 916 |
|
|
| 917 |
names(d_z) <- y_var_name |
|
| 918 |
df_obs <- as.data.frame(d_z) |
|
| 919 |
names(d_z) <- y_var_name |
|
| 908 |
d_z_all_w <- window(d_z_all,start=start_date,end=end_date) |
|
| 909 |
d_z_w <- window(d_z,start=start_date,end=end_date) |
|
| 920 | 910 |
|
| 921 | 911 |
res_pix <- 1000 |
| 922 | 912 |
#res_pix <- 480 |
| 923 | 913 |
col_mfrow <- 2 |
| 924 | 914 |
row_mfrow <- 1 |
| 925 | 915 |
|
| 926 |
png_filename <- file.path(out_dir,paste("Figure5b_time_series_profile_",region_name,"_",out_suffix,".png",sep =""))
|
|
| 927 |
title_str <- paste("Predicted daily ", station_id," ",var," for the ", start_year,"-",end_year," time period",sep="")
|
|
| 928 |
|
|
| 929 |
png(filename=png_filename,width = col_mfrow * res_pix,height = row_mfrow * res_pix) |
|
| 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=" ")
|
|
| 930 | 921 |
|
| 931 |
plot(df_pix_ts[[var_pred_mosaic]] ~ as.Date(df_pix_ts$date),ylab="tmax in deg C",xlab="Daily time steps", |
|
| 932 |
main=title_str,cex=1,font=2,type="l", |
|
| 933 |
cex.main=1.5,cex.lab=1.5,font.lab=2, |
|
| 934 |
lty=3) |
|
| 935 |
lines(as.numeric(df_pix_ts[[y_var_name]]) ~ as.Date(df_pix_ts$date),ylab="tmax in deg C",xlab="Daily time steps", |
|
| 936 |
main=title_str,cex=1,font=2,col="red",type="l", |
|
| 937 |
cex.main=1.5,cex.lab=1.5,font.lab=2, |
|
| 938 |
lty=3) |
|
| 922 |
png(filename=png_filename,width = col_mfrow * res_pix,height = row_mfrow * res_pix) |
|
| 939 | 923 |
|
| 924 |
#add legend later!! |
|
| 925 |
col_str <- c("blue","red")
|
|
| 926 |
plot(d_z_w,plot.type="single",col=col_str, |
|
| 927 |
ylab="tmax in deg C",xlab="Daily time steps", |
|
| 928 |
main=title_str,cex=1,font=2,type="l", |
|
| 929 |
cex.main=1.5,cex.lab=1.5,font.lab=2) |
|
| 930 |
|
|
| 940 | 931 |
dev.off() |
| 941 | 932 |
|
| 942 | 933 |
} |
Also available in: Unified diff
plotting function time series, fixing bugs