Revision 48cbb42e
Added by Benoit Parmentier over 8 years ago
| climate/research/oregon/interpolation/global_product_assessment_part1.R | ||
|---|---|---|
| 19 | 19 |
# |
| 20 | 20 |
#setfacl -Rmd user:aguzman4:rwx /nobackupp8/bparmen1/output_run10_1500x4500_global_analyses_pred_1992_10052015 |
| 21 | 21 |
|
| 22 |
#COMMIT: testing function combine extraction and assessment data, also moved to function script
|
|
| 22 |
#COMMIT: plotting extracted predicted values and measured tmax
|
|
| 23 | 23 |
|
| 24 | 24 |
################################################################################################# |
| 25 | 25 |
|
| ... | ... | |
| 85 | 85 |
source(file.path(script_path,function_assessment_part3)) #source all functions used in this script |
| 86 | 86 |
|
| 87 | 87 |
#Product assessment |
| 88 |
function_product_assessment_part1_functions <- "global_product_assessment_part1_functions_09192016.R" |
|
| 88 |
function_product_assessment_part1_functions <- "global_product_assessment_part1_functions_09192016b.R"
|
|
| 89 | 89 |
source(file.path(script_path,function_product_assessment_part1_functions)) #source all functions used in this script |
| 90 | 90 |
|
| 91 | 91 |
############################### |
| ... | ... | |
| 593 | 593 |
i<-1 |
| 594 | 594 |
|
| 595 | 595 |
#Product assessment |
| 596 |
#function_product_assessment_part1_functions <- "global_product_assessment_part1_functions_09192016.R"
|
|
| 597 |
#source(file.path(script_path,function_product_assessment_part1_functions)) #source all functions used in this script
|
|
| 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 | 598 |
#undebug(combine_measurements_and_predictions_df) |
| 599 |
|
|
| 599 |
out_suffix_str <- paste0(region_name,"_",out_suffix) |
|
| 600 | 600 |
#this can be run with mclapply, very fast right now: |
| 601 | 601 |
station_summary_obj <- combine_measurements_and_predictions_df(i=i, |
| 602 | 602 |
df_raster=df_raster, |
| ... | ... | |
| 607 | 607 |
r_ts_name=r_ts_name, |
| 608 | 608 |
var_name=var_name, |
| 609 | 609 |
var_pred = var_pred, |
| 610 |
plot_fig=T) |
|
| 610 |
out_dir =out_dir, |
|
| 611 |
out_suffix=out_suffix_str, |
|
| 612 |
plot_fig=F) |
|
| 611 | 613 |
df_pix_ts <- station_summary_obj$df_pix_ts |
| 612 | 614 |
#station_summary_obj <- list(nb_zero,nb_NA, df_pix_ts) |
| 613 | 615 |
|
| 614 |
id_name <- list_selected_ID[i] |
|
| 615 |
df_pix_ts_filename <- file.path(out_dir,paste0("df_pix_ts_",id_name,out_suffix,y_var_name,".txt"))
|
|
| 616 |
write.table(df_pix_ts,df_pix_ts_filename,sep=",") |
|
| 616 |
#id_name <- list_selected_ID[i]
|
|
| 617 |
#df_pix_ts_filename <- file.path(out_dir,paste0("df_pix_ts_",id_name,out_suffix,y_var_name,".txt"))
|
|
| 618 |
#write.table(df_pix_ts,df_pix_ts_filename,sep=",")
|
|
| 617 | 619 |
|
| 618 | 620 |
|
| 619 | 621 |
##################### plot time series: make this a function!!! #### |
| 620 | 622 |
###start of plotting |
| 621 | 623 |
### makes this a function: |
| 622 | 624 |
|
| 623 |
pix_ts$date <- list_dates_produced_date_val |
|
| 624 |
pix_ts[,1]*scaling #scale? |
|
| 625 |
|
|
| 626 |
names(pix_ts) <- var_id |
|
| 627 |
|
|
| 628 |
pix_ts <- read.table("pix_ts2.txt",sep=",")
|
|
| 629 |
names(pix_ts) <- c(var_id,"files","date_str") |
|
| 630 |
pix_ts$date <- as.Date(strptime(pix_ts$date_str,"%Y%m%d")) |
|
| 631 |
|
|
| 625 |
#pix_ts$date <- as.Date(strptime(pix_ts$date_str,"%Y%m%d")) |
|
| 632 | 626 |
#head(pix_ts) |
| 627 |
var_pred_mosaic <- paste0(var_pred,"_mosaic") |
|
| 628 |
id_name <- list_selected_ID[i] |
|
| 629 |
station_id <- id_name |
|
| 630 |
#> myDF[ is.na(myDF) ] <- NA |
|
| 633 | 631 |
|
| 634 |
id_selected <- "82111099999" |
|
| 635 |
#station_id <- 8 |
|
| 636 |
station_id <- id_selected |
|
| 637 |
df <- pix_ts |
|
| 638 |
#scaling |
|
| 639 |
#start_date: subset dates |
|
| 640 |
#end_date: subset dates |
|
| 641 |
df2 <- data_stations_var_pred |
|
| 642 |
|
|
| 643 |
|
|
| 644 |
df_selected <- subset(df,select=station_id) |
|
| 645 |
#df_selected <- subset(pix_ts,select=station_id) |
|
| 646 |
names(df_selected) <- station_id |
|
| 647 |
df_selected$date <- df$date |
|
| 648 |
|
|
| 632 |
df_pix_ts[is.na(df_pix_ts)] <- NA |
|
| 649 | 633 |
if(!is.null(scaling)){
|
| 650 |
df_selected[[station_id]]<- df_selected[[station_id]]*scaling |
|
| 651 |
} |
|
| 652 |
if(!is.null(df2)){
|
|
| 653 |
df_selected2 <- df2 |
|
| 654 |
rm(df2) |
|
| 655 |
d_z_tmp2 <- zoo(df_selected2$dailyTmax,df_selected2$date) |
|
| 656 |
names(d_z_tmp2)<-station_id |
|
| 657 |
}else{
|
|
| 658 |
df_selected2 <- NULL |
|
| 634 |
df_pix_ts[[var_pred_mosaic]] <- df_pix_ts[[var_pred_mosaic]]*scaling |
|
| 659 | 635 |
} |
| 660 | 636 |
|
| 661 |
d_z_tmp <- zoo(df_selected[[station_id]],df_selected$date) |
|
| 662 |
#d_z_tmp <- zoo(df[[station_id]],df$date) |
|
| 663 |
names(d_z_tmp)<-station_id |
|
| 664 |
#min(d_z_tmp$ID_8) |
|
| 665 |
#max(d_z_tmp$ID_8) |
|
| 666 |
day_start <- "1984-01-01" #PARAM 12 arg 12 |
|
| 667 |
day_end <- "2014-12-31" #PARAM 13 arg 13 |
|
| 668 |
start_date <- as.Date(day_start) |
|
| 669 |
end_date <- as.Date(day_end) |
|
| 670 |
start_year <- year(start_date) |
|
| 671 |
end_year <- year(end_date) |
|
| 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" |
|
| 672 | 663 |
|
| 673 | 664 |
res_pix <- 1000 |
| 674 | 665 |
#res_pix <- 480 |
| 675 | 666 |
col_mfrow <- 2 |
| 676 | 667 |
row_mfrow <- 1 |
| 677 | 668 |
|
| 678 |
png_filename <- file.path(out_dir,paste("Figure5a_time_series_profile_",region_name,"_",out_suffix,".png",sep =""))
|
|
| 669 |
png_filename <- file.path(out_dir,paste("Figure5a_time_series_profile_",region_name,"_",id_name,y_var_name,out_suffix,".png",sep =""))
|
|
| 679 | 670 |
title_str <- paste("Predicted daily ", station_id," ",var," for the ", start_year,"-",end_year," time period",sep="")
|
| 680 | 671 |
|
| 681 | 672 |
png(filename=png_filename,width = col_mfrow * res_pix,height = row_mfrow * res_pix) |
| 682 | 673 |
#this is the whole time series |
| 683 |
plot(d_z_tmp,ylab="tmax in deg C",xlab="Daily time steps", |
|
| 684 |
main=title_str,cex=3,font=2, |
|
| 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", |
|
| 685 | 682 |
cex.main=1.5,cex.lab=1.5,font.lab=2, |
| 686 | 683 |
lty=3) |
| 687 |
if(!is.null(df_selected2)){
|
|
| 688 |
lines(d_z_tmp2,ylab="tmax in deg C",xlab="Daily time steps", |
|
| 689 |
main=title_str,cex=3,font=2, |
|
| 690 |
col="red", |
|
| 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", |
|
| 691 | 686 |
cex.main=1.5,cex.lab=1.5,font.lab=2, |
| 692 | 687 |
lty=3) |
| 693 |
} |
|
| 688 |
|
|
| 689 |
|
|
| 694 | 690 |
dev.off() |
| 695 | 691 |
|
| 696 |
day_start <- "1984-01-01" #PARAM 12 arg 12 |
|
| 697 |
day_end <- "1998-12-31" #PARAM 13 arg 13 |
|
| 692 |
### get smaller window |
|
| 693 |
day_start <- "1994-01-01" #PARAM 12 arg 12 |
|
| 694 |
day_end <- "1999-12-31" #PARAM 13 arg 13 |
|
| 698 | 695 |
|
| 699 | 696 |
start_date <- as.Date(day_start) |
| 700 | 697 |
end_date <- as.Date(day_end) |
| 701 | 698 |
start_year <- year(start_date) |
| 702 | 699 |
end_year <- year(end_date) |
| 703 |
|
|
| 700 |
d_z_tmp <- d_z_var |
|
| 704 | 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 |
|
|
| 705 | 715 |
if(!is.null(df_selected2)){
|
| 706 | 716 |
d_z2 <- window(d_z_tmp2,start=start_date,end=end_date) |
| 707 | 717 |
} |
| ... | ... | |
| 716 | 726 |
|
| 717 | 727 |
png(filename=png_filename,width = col_mfrow * res_pix,height = row_mfrow * res_pix) |
| 718 | 728 |
|
| 719 |
plot(d_z,ylab="tmax in deg C",xlab="Daily time steps", |
|
| 720 |
main=title_str,cex=3,font=2, |
|
| 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", |
|
| 721 | 735 |
cex.main=1.5,cex.lab=1.5,font.lab=2, |
| 722 | 736 |
lty=3) |
| 723 |
if(!is.null(df_selected2)){
|
|
| 724 |
lines(d_z2,ylab="tmax in deg C",xlab="Daily time steps", |
|
| 725 |
main=title_str,cex=3,font=2, |
|
| 726 |
col="red", |
|
| 727 |
cex.main=1.5,cex.lab=1.5,font.lab=2, |
|
| 728 |
lty=3) |
|
| 729 |
} |
|
| 730 | 737 |
|
| 731 | 738 |
dev.off() |
| 732 | 739 |
|
Also available in: Unified diff
plotting extracted predicted values and measured tmax