Revision fabd94a4
Added by Benoit Parmentier over 8 years ago
| climate/research/oregon/interpolation/global_product_assessment_part1.R | ||
|---|---|---|
| 593 | 593 |
i<-1 |
| 594 | 594 |
|
| 595 | 595 |
#Product assessment |
| 596 |
function_product_assessment_part1_functions <- "global_product_assessment_part1_functions_09192016b.R" |
|
| 597 |
source(file.path(script_path,function_product_assessment_part1_functions)) #source all functions used in this script |
|
| 598 |
#undebug(combine_measurements_and_predictions_df) |
|
| 599 | 596 |
out_suffix_str <- paste0(region_name,"_",out_suffix) |
| 600 | 597 |
#this can be run with mclapply, very fast right now: |
| 601 | 598 |
station_summary_obj <- combine_measurements_and_predictions_df(i=i, |
| ... | ... | |
| 617 | 614 |
#df_pix_ts_filename <- file.path(out_dir,paste0("df_pix_ts_",id_name,out_suffix,y_var_name,".txt"))
|
| 618 | 615 |
#write.table(df_pix_ts,df_pix_ts_filename,sep=",") |
| 619 | 616 |
|
| 620 |
|
|
| 621 | 617 |
##################### plot time series: make this a function!!! #### |
| 622 | 618 |
###start of plotting |
| 623 | 619 |
### makes this a function: |
| ... | ... | |
| 627 | 623 |
var_pred_mosaic <- paste0(var_pred,"_mosaic") |
| 628 | 624 |
id_name <- list_selected_ID[i] |
| 629 | 625 |
station_id <- id_name |
| 630 |
#> myDF[ is.na(myDF) ] <- NA |
|
| 631 | 626 |
|
| 632 |
df_pix_ts[is.na(df_pix_ts)] <- NA |
|
| 633 |
if(!is.null(scaling)){
|
|
| 634 |
df_pix_ts[[var_pred_mosaic]] <- df_pix_ts[[var_pred_mosaic]]*scaling |
|
| 635 |
} |
|
| 627 |
var_name1 <- y_var_name |
|
| 628 |
var_name2 <- "mod1_mosaic" |
|
| 629 |
out_suffix_str <- paste(region_name,out_suffix,sep="_") |
|
| 630 |
scaling_factors <- c("NA",scaling) #no scaling for y_var_name but scaling for var_name2
|
|
| 631 |
list_windows <- list(c("1999-01-01","1999-12-31"))
|
|
| 636 | 632 |
|
| 637 |
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 !!! |
|
| 638 |
names(d_z_var) <- var_pred_mosaic |
|
| 639 |
plot(d_z_var) |
|
| 640 |
d_z_obs <- zoo(as.numeric(df_pix_ts[[y_var_name]]),as.Date(df_pix_ts$date)) |
|
| 641 |
plot(d_z_obs) |
|
| 642 |
|
|
| 643 |
d_z_res <- zoo(as.numeric(df_pix_ts[[paste0("res_",var_pred_mosaic)]]),as.Date(df_pix_ts$date))
|
|
| 644 |
plot(d_z_res) |
|
| 645 |
|
|
| 646 |
d_z <- merge(d_z_var,d_z_obs) |
|
| 647 |
range(d_z_res,na.rm=T) |
|
| 648 |
mean(d_z_res,na.rm=T) |
|
| 649 |
quantile(d_z_res,na.rm=T) |
|
| 650 |
plot(d_z,plot.type="single",col=c("blue","red"))
|
|
| 651 |
names(d_z) <- c("pred","obs")
|
|
| 652 |
|
|
| 653 |
#day_start <- "1984-01-01" #PARAM 12 arg 12 |
|
| 654 |
#day_end <- "2014-12-31" #PARAM 13 arg 13 |
|
| 655 |
#start_date <- as.Date(day_start) |
|
| 656 |
#end_date <- as.Date(day_end) |
|
| 657 |
#start_year <- year(start_date) |
|
| 658 |
#end_year <- year(end_date) |
|
| 659 |
range_year <- range(df_pix_ts$year) |
|
| 660 |
start_year <- range_year[1] |
|
| 661 |
end_year <- range_year[2] |
|
| 662 |
var="tmax" |
|
| 663 |
|
|
| 664 |
res_pix <- 1000 |
|
| 665 |
#res_pix <- 480 |
|
| 666 |
col_mfrow <- 2 |
|
| 667 |
row_mfrow <- 1 |
|
| 668 |
|
|
| 669 |
png_filename <- file.path(out_dir,paste("Figure5a_time_series_profile_",region_name,"_",id_name,y_var_name,out_suffix,".png",sep =""))
|
|
| 670 |
title_str <- paste("Predicted daily ", station_id," ",var," for the ", start_year,"-",end_year," time period",sep="")
|
|
| 671 |
|
|
| 672 |
png(filename=png_filename,width = col_mfrow * res_pix,height = row_mfrow * res_pix) |
|
| 673 |
#this is the whole time series |
|
| 674 |
#plot(d_z_var,ylab="tmax in deg C",xlab="Daily time steps", |
|
| 675 |
# main=title_str,cex=3,font=2, |
|
| 676 |
# cex.main=1.5,cex.lab=1.5,font.lab=2, |
|
| 677 |
# lty=3) |
|
| 678 |
#range(df_pix_ts[[var_pred_mosaic]],na.rm=T) |
|
| 679 |
|
|
| 680 |
plot(df_pix_ts[[var_pred_mosaic]] ~ as.Date(df_pix_ts$date),ylab="tmax in deg C",xlab="Daily time steps", |
|
| 681 |
main=title_str,cex=1,font=2,type="l", |
|
| 682 |
cex.main=1.5,cex.lab=1.5,font.lab=2, |
|
| 683 |
lty=3) |
|
| 684 |
lines(as.numeric(df_pix_ts[[y_var_name]]) ~ as.Date(df_pix_ts$date),ylab="tmax in deg C",xlab="Daily time steps", |
|
| 685 |
main=title_str,cex=1,font=2,col="red",type="l", |
|
| 686 |
cex.main=1.5,cex.lab=1.5,font.lab=2, |
|
| 687 |
lty=3) |
|
| 688 |
|
|
| 689 |
|
|
| 690 |
dev.off() |
|
| 691 |
|
|
| 692 |
### get smaller window |
|
| 693 |
day_start <- "1994-01-01" #PARAM 12 arg 12 |
|
| 694 |
day_end <- "1999-12-31" #PARAM 13 arg 13 |
|
| 695 |
|
|
| 696 |
start_date <- as.Date(day_start) |
|
| 697 |
end_date <- as.Date(day_end) |
|
| 698 |
start_year <- year(start_date) |
|
| 699 |
end_year <- year(end_date) |
|
| 700 |
d_z_tmp <- d_z_var |
|
| 701 |
d_z <- window(d_z_tmp,start=start_date,end=end_date) |
|
| 702 |
names(d_z) <- var_pred_mosaic |
|
| 703 |
df_var <- as.data.frame(d_z) |
|
| 704 |
|
|
| 705 |
d_z_tmp <- d_z_obs |
|
| 706 |
d_z_tmp <- d_z |
|
| 707 |
d_z_tmp <- window(d_z_tmp,start=start_date,end=end_date) |
|
| 708 |
plot(d_z_tmp,col=c("blue","red"),plot.type="single")
|
|
| 709 |
df_z <- as.data.frame(d_z) |
|
| 710 |
|
|
| 711 |
names(d_z) <- y_var_name |
|
| 712 |
df_obs <- as.data.frame(d_z) |
|
| 713 |
names(d_z) <- y_var_name |
|
| 714 |
|
|
| 715 |
if(!is.null(df_selected2)){
|
|
| 716 |
d_z2 <- window(d_z_tmp2,start=start_date,end=end_date) |
|
| 717 |
} |
|
| 718 |
|
|
| 719 |
res_pix <- 1000 |
|
| 720 |
#res_pix <- 480 |
|
| 721 |
col_mfrow <- 2 |
|
| 722 |
row_mfrow <- 1 |
|
| 723 |
|
|
| 724 |
png_filename <- file.path(out_dir,paste("Figure5b_time_series_profile_",region_name,"_",out_suffix,".png",sep =""))
|
|
| 725 |
title_str <- paste("Predicted daily ", station_id," ",var," for the ", start_year,"-",end_year," time period",sep="")
|
|
| 726 |
|
|
| 727 |
png(filename=png_filename,width = col_mfrow * res_pix,height = row_mfrow * res_pix) |
|
| 728 |
|
|
| 729 |
plot(df_pix_ts[[var_pred_mosaic]] ~ as.Date(df_pix_ts$date),ylab="tmax in deg C",xlab="Daily time steps", |
|
| 730 |
main=title_str,cex=1,font=2,type="l", |
|
| 731 |
cex.main=1.5,cex.lab=1.5,font.lab=2, |
|
| 732 |
lty=3) |
|
| 733 |
lines(as.numeric(df_pix_ts[[y_var_name]]) ~ as.Date(df_pix_ts$date),ylab="tmax in deg C",xlab="Daily time steps", |
|
| 734 |
main=title_str,cex=1,font=2,col="red",type="l", |
|
| 735 |
cex.main=1.5,cex.lab=1.5,font.lab=2, |
|
| 736 |
lty=3) |
|
| 737 |
|
|
| 738 |
dev.off() |
|
| 739 |
|
|
| 740 |
#### Subset for 5c: zoom in |
|
| 741 |
|
|
| 742 |
day_start <- "1991-01-01" #PARAM 12 arg 12 |
|
| 743 |
day_end <- "1992-12-31" #PARAM 13 arg 13 |
|
| 744 |
|
|
| 745 |
start_date <- as.Date(day_start) |
|
| 746 |
end_date <- as.Date(day_end) |
|
| 747 |
start_year <- year(start_date) |
|
| 748 |
end_year <- year(end_date) |
|
| 749 |
d_z <- window(d_z_tmp,start=start_date,end=end_date) |
|
| 750 |
if(!is.null(df_selected2)){
|
|
| 751 |
d_z2 <- window(d_z_tmp2,start=start_date,end=end_date) |
|
| 752 |
} |
|
| 753 |
|
|
| 754 |
res_pix <- 1000 |
|
| 755 |
#res_pix <- 480 |
|
| 756 |
col_mfrow <- 2 |
|
| 757 |
row_mfrow <- 1 |
|
| 758 |
|
|
| 759 |
png_filename <- file.path(out_dir,paste("Figure5c_subset_time_series_profile_",region_name,"_",out_suffix,".png",sep =""))
|
|
| 760 |
title_str <- paste("Predicted daily ", station_id," ",var," for the ", start_year,"-",end_year," time period",sep="")
|
|
| 761 |
|
|
| 762 |
png(filename=png_filename,width = col_mfrow * res_pix,height = row_mfrow * res_pix) |
|
| 763 |
|
|
| 764 |
plot(d_z,ylab="tmax in deg C",xlab="Daily time steps", |
|
| 765 |
main=title_str,cex=3,font=2, |
|
| 766 |
cex.main=1.5,cex.lab=1.5,font.lab=2, |
|
| 767 |
lty=3) |
|
| 768 |
if(!is.null(df_selected2)){
|
|
| 769 |
lines(d_z2,ylab="tmax in deg C",xlab="Daily time steps", |
|
| 770 |
main=title_str,cex=3,font=2, |
|
| 771 |
col="red", |
|
| 772 |
cex.main=1.5,cex.lab=1.5,font.lab=2, |
|
| 773 |
lty=3) |
|
| 774 |
} |
|
| 633 |
function_product_assessment_part1_functions <- "global_product_assessment_part1_functions_09202016.R" |
|
| 634 |
source(file.path(script_path,function_product_assessment_part1_functions)) #source all functions used in this script |
|
| 775 | 635 |
|
| 776 |
dev.off() |
|
| 636 |
###TO DO: |
|
| 637 |
#1) compute correlation r, RMSE, MAE etc by station and profile |
|
| 638 |
#2) Compute temporal autocorrelation by station and profile |
|
| 639 |
#3) compute range, min, max etc by stations |
|
| 640 |
|
|
| 641 |
undebug(plot_observation_predictions_time_series) |
|
| 642 |
#test_plot <- plot_observation_predictions_time_series(df_pix_time_series=df_pix_ts, |
|
| 643 |
# var_name1=var_name1, |
|
| 644 |
# var_name2=var_name2, |
|
| 645 |
# list_windows=list_windows, |
|
| 646 |
# time_series_id=id_name, |
|
| 647 |
# scaling_factors=scaling_factors, |
|
| 648 |
# out_dir=out_dir, |
|
| 649 |
# out_suffix=out_suffix_str) |
|
| 650 |
out_suffix_str2 <- paste(region_name,out_suffix,"test",sep="_") |
|
| 651 |
|
|
| 652 |
test_plot <- plot_observation_predictions_time_series(df_pix_time_series=df_pix_ts, |
|
| 653 |
var_name1=var_name1, |
|
| 654 |
var_name2=var_name2, |
|
| 655 |
list_windows=list_windows, |
|
| 656 |
time_series_id=id_name, |
|
| 657 |
scaling_factors=NULL, |
|
| 658 |
out_dir=out_dir, |
|
| 659 |
out_suffix=out_suffix_str2) |
|
| 777 | 660 |
|
| 778 | 661 |
############### PART5: Make raster stack and display maps ############# |
| 779 | 662 |
#### Extract corresponding raster for given dates and plot stations used |
| 780 | 663 |
|
| 664 |
## TODO: make movies from time series in png |
|
| 781 | 665 |
|
| 782 | 666 |
#start_date <- day_to_mosaic_range[1] |
| 783 | 667 |
#end_date <- day_to_mosaic_range[2] |
| ... | ... | |
| 827 | 711 |
|
| 828 | 712 |
lf_mosaic_plot_fig <- mclapply(1:length(lf_mosaic_scaled),FUN=plot_raster_mosaic,list_param=list_param_plot_raster_mosaic,mc.preschedule=FALSE,mc.cores = num_cores) |
| 829 | 713 |
|
| 830 |
#################### |
|
| 831 |
###### Now add figures with additional met stations? |
|
| 832 |
|
|
| 833 |
#df_selected2 <- data_stations_var_pred |
|
| 834 |
|
|
| 835 |
#d_z_tmp <- zoo(df_selected[[station_id]],df_selected$date) |
|
| 836 |
#d_z_tmp2 <- zoo(df_selected2$dailyTmax,df_selected2$date) |
|
| 837 |
|
|
| 838 |
#d_z_tmp <- zoo(df[[station_id]],df$date) |
|
| 839 |
#names(d_z_tmp)<-station_id |
|
| 840 |
#names(d_z_tmp2)<-station_id |
|
| 841 |
|
|
| 842 |
#min(d_z_tmp$ID_8) |
|
| 843 |
#max(d_z_tmp$ID_8) |
|
| 844 |
#day_start <- "1984-01-01" #PARAM 12 arg 12 |
|
| 845 |
#day_end <- "2014-12-31" #PARAM 13 arg 13 |
|
| 846 |
#day_start <- "1991-01-01" #PARAM 12 arg 12 |
|
| 847 |
#day_end <- "1992-12-31" #PARAM 13 arg 13 |
|
| 848 |
|
|
| 849 |
#start_date <- as.Date(day_start) |
|
| 850 |
#end_date <- as.Date(day_end) |
|
| 851 |
#start_year <- year(start_date) |
|
| 852 |
#end_year <- year(end_date) |
|
| 853 |
|
|
| 854 |
#res_pix <- 1000 |
|
| 855 |
#res_pix <- 480 |
|
| 856 |
#col_mfrow <- 2 |
|
| 857 |
#row_mfrow <- 1 |
|
| 858 |
|
|
| 859 |
#png_filename <- file.path(out_dir,paste("Figure5a_time_series_profile_",region_name,"_",out_suffix,".png",sep =""))
|
|
| 860 |
#title_str <- paste("Predicted daily ", station_id," ",var," for the ", start_year,"-",end_year," 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 |
|
|
| 865 |
#plot(d_z_tmp,ylab="tmax in deg C",xlab="Daily time steps", |
|
| 866 |
# main=title_str,cex=3,font=2, |
|
| 867 |
# cex.main=1.5,cex.lab=1.5,font.lab=2, |
|
| 868 |
# lty=3) |
|
| 869 |
#lines(d_z_tmp2,ylab="tmax in deg C",xlab="Daily time steps", |
|
| 870 |
# main=title_str,cex=3,font=2, |
|
| 871 |
# col="red", |
|
| 872 |
# cex.main=1.5,cex.lab=1.5,font.lab=2, |
|
| 873 |
# lty=3) |
|
| 874 |
# |
|
| 875 |
#dev.off() |
|
| 876 |
|
|
| 877 |
#This is a lot of replication!!! okay cut it down |
|
| 878 |
data_stations$date <- as.Date(strptime(data_stations_var_pred$date_str,"%Y%m%d")) |
|
| 879 |
data_stations_tmp <- data_stations |
|
| 880 |
data_stations_tmp <- data_stations |
|
| 881 |
|
|
| 882 |
data_stations_tmp$date <- as.Date(strptime(data_stations_tmp$date,"%Y%m%d")) |
|
| 883 |
#data_stations_tmp$date <- as.Date(strptime(data_stations_tmp$date_str,"%Y%m%d")) |
|
| 884 |
#d_z4 <- |
|
| 885 |
d_z_tmp4 <- zoo(data_stations_tmp$dailyTmax,data_stations_tmp$date) |
|
| 886 |
plot(d_z_tmp,cex.main=1.5,cex.lab=1.5,font.lab=2, |
|
| 887 |
lty=3) |
|
| 888 |
lines(d_z_tmp4,ylab="tmax in deg C",xlab="Daily time steps", |
|
| 889 |
main=title_str,cex=3,font=2, |
|
| 890 |
col="red", |
|
| 891 |
cex.main=1.5,cex.lab=1.5,font.lab=2, |
|
| 892 |
lty=3) |
|
| 893 |
|
|
| 894 |
day_start <- "1991-01-01" #PARAM 12 arg 12 |
|
| 895 |
day_end <- "1992-12-31" #PARAM 13 arg 13 |
|
| 896 |
|
|
| 897 |
start_date <- as.Date(day_start) |
|
| 898 |
end_date <- as.Date(day_end) |
|
| 899 |
start_year <- year(start_date) |
|
| 900 |
end_year <- year(end_date) |
|
| 901 |
d_z4 <- window(d_z_tmp4,start=start_date,end=end_date) |
|
| 902 |
|
|
| 903 |
data_stations_tmp$date[7190] |
|
| 904 |
|
|
| 905 |
###TO DO: |
|
| 906 |
#1) compute correlation r, RMSE, MAE etc by station and profile |
|
| 907 |
#2) Compute temporal autocorrelation by station and profile |
|
| 908 |
#3) |
|
| 909 |
|
|
| 910 | 714 |
#### PLOT ACCURACY METRICS: First test #### |
| 911 | 715 |
##this will be cleaned up later: |
| 912 | 716 |
|
Also available in: Unified diff
testing ploting time series profile function and clean up of main code for product assessment part1