Revision 0060900c
Added by Benoit Parmentier over 8 years ago
| climate/research/oregon/interpolation/global_product_assessment_part1.R | ||
|---|---|---|
| 79 | 79 |
source(file.path(script_path,function_assessment_part2_functions)) #source all functions used in this script |
| 80 | 80 |
source(file.path(script_path,function_assessment_part3)) #source all functions used in this script |
| 81 | 81 |
|
| 82 |
#Product assessment |
|
| 83 |
function_product_assessment_part1_functions <- "global_product_assessment_part1_05242016.R" |
|
| 84 |
source(file.path(script_path,function_product_assessment_part1_functions)) #source all functions used in this script |
|
| 85 |
|
|
| 82 | 86 |
############################### |
| 83 | 87 |
####### Parameters, constants and arguments ### |
| 84 | 88 |
|
| ... | ... | |
| 301 | 305 |
dim(data_stations) #one station only but repitition of records because of tiles and dates!!! |
| 302 | 306 |
#> dim(data_stations) |
| 303 | 307 |
#[1] 100458 90 |
| 304 |
#This is a lot of replication!!! okay cut it down |
|
| 305 | 308 |
|
| 306 | 309 |
#data_stations_temp <- aggregate(id ~ date, data = data_stations, min) |
| 307 | 310 |
#data_stations_temp <- aggregate(id ~ x + y + date + dailyTmax + mod1 + res_mod1 , data = data_stations, min) |
| ... | ... | |
| 310 | 313 |
#> dim(data_stations_temp) |
| 311 | 314 |
#[1] 11171 5 |
| 312 | 315 |
|
| 316 |
data_stations_temp$date_str <- data_stations_temp$date |
|
| 317 |
data_stations_temp$date <- as.Date(strptime(data_stations_temp$date_str,"%Y%m%d")) |
|
| 318 |
|
|
| 319 |
##Find stations used for training and testing |
|
| 320 |
#data_stations_temp2 <- aggregate(id ~ training,data = data_stations, sum ) #+ mod1 + res_mod1 , data = data_stations, min) |
|
| 321 |
#data_stations_temp2 <- aggregate(date ~ training,data = data_stations, sum ) #+ mod1 + res_mod1 , data = data_stations, min) |
|
| 322 |
|
|
| 313 | 323 |
############### PART3: Make raster stack and display maps ############# |
| 314 | 324 |
#### Extract corresponding raster for given dates and plot stations used |
| 315 | 325 |
|
| ... | ... | |
| 493 | 503 |
#scaling |
| 494 | 504 |
#start_date: subset dates |
| 495 | 505 |
#end_date: subset dates |
| 496 |
|
|
| 506 |
df2 <- data_stations_temp |
|
| 507 |
|
|
| 508 |
|
|
| 497 | 509 |
df_selected <- subset(df,select=station_id) |
| 498 | 510 |
#df_selected <- subset(pix_ts,select=station_id) |
| 499 | 511 |
names(df_selected) <- station_id |
| ... | ... | |
| 502 | 514 |
if(!is.null(scaling)){
|
| 503 | 515 |
df_selected[[station_id]]<- df_selected[[station_id]]*scaling |
| 504 | 516 |
} |
| 517 |
if(!is.null(df2)){
|
|
| 518 |
df_selected2 <- df2 |
|
| 519 |
rm(df2) |
|
| 520 |
d_z_tmp2 <- zoo(df_selected2$dailyTmax,df_selected2$date) |
|
| 521 |
names(d_z_tmp2)<-station_id |
|
| 522 |
}else{
|
|
| 523 |
df_selected2 <- NULL |
|
| 524 |
} |
|
| 505 | 525 |
|
| 506 | 526 |
d_z_tmp <- zoo(df_selected[[station_id]],df_selected$date) |
| 507 | 527 |
#d_z_tmp <- zoo(df[[station_id]],df$date) |
| ... | ... | |
| 529 | 549 |
main=title_str,cex=3,font=2, |
| 530 | 550 |
cex.main=1.5,cex.lab=1.5,font.lab=2, |
| 531 | 551 |
lty=3) |
| 532 |
|
|
| 552 |
if(!is.null(df_selected2)){
|
|
| 553 |
lines(d_z_tmp2,ylab="tmax in deg C",xlab="Daily time steps", |
|
| 554 |
main=title_str,cex=3,font=2, |
|
| 555 |
col="red", |
|
| 556 |
cex.main=1.5,cex.lab=1.5,font.lab=2, |
|
| 557 |
lty=3) |
|
| 558 |
} |
|
| 533 | 559 |
dev.off() |
| 534 | 560 |
|
| 535 | 561 |
day_start <- "1984-01-01" #PARAM 12 arg 12 |
| ... | ... | |
| 541 | 567 |
end_year <- year(end_date) |
| 542 | 568 |
|
| 543 | 569 |
d_z <- window(d_z_tmp,start=start_date,end=end_date) |
| 544 |
#d_z2 <- window(d_z_tmp2,start=start_date,end=end_date) |
|
| 570 |
if(!is.null(df_selected2)){
|
|
| 571 |
d_z2 <- window(d_z_tmp2,start=start_date,end=end_date) |
|
| 572 |
} |
|
| 545 | 573 |
|
| 546 | 574 |
res_pix <- 1000 |
| 547 | 575 |
#res_pix <- 480 |
| ... | ... | |
| 557 | 585 |
main=title_str,cex=3,font=2, |
| 558 | 586 |
cex.main=1.5,cex.lab=1.5,font.lab=2, |
| 559 | 587 |
lty=3) |
| 588 |
if(!is.null(df_selected2)){
|
|
| 589 |
lines(d_z2,ylab="tmax in deg C",xlab="Daily time steps", |
|
| 590 |
main=title_str,cex=3,font=2, |
|
| 591 |
col="red", |
|
| 592 |
cex.main=1.5,cex.lab=1.5,font.lab=2, |
|
| 593 |
lty=3) |
|
| 594 |
} |
|
| 560 | 595 |
|
| 561 | 596 |
dev.off() |
| 562 | 597 |
|
| ... | ... | |
| 570 | 605 |
start_year <- year(start_date) |
| 571 | 606 |
end_year <- year(end_date) |
| 572 | 607 |
d_z <- window(d_z_tmp,start=start_date,end=end_date) |
| 573 |
#d_z2 <- window(d_z_tmp2,start=start_date,end=end_date) |
|
| 608 |
if(!is.null(df_selected2)){
|
|
| 609 |
d_z2 <- window(d_z_tmp2,start=start_date,end=end_date) |
|
| 610 |
} |
|
| 574 | 611 |
|
| 575 | 612 |
res_pix <- 1000 |
| 576 | 613 |
#res_pix <- 480 |
| ... | ... | |
| 586 | 623 |
main=title_str,cex=3,font=2, |
| 587 | 624 |
cex.main=1.5,cex.lab=1.5,font.lab=2, |
| 588 | 625 |
lty=3) |
| 626 |
if(!is.null(df_selected2)){
|
|
| 627 |
lines(d_z2,ylab="tmax in deg C",xlab="Daily time steps", |
|
| 628 |
main=title_str,cex=3,font=2, |
|
| 629 |
col="red", |
|
| 630 |
cex.main=1.5,cex.lab=1.5,font.lab=2, |
|
| 631 |
lty=3) |
|
| 632 |
} |
|
| 589 | 633 |
|
| 590 | 634 |
dev.off() |
| 591 | 635 |
|
| 636 |
#################### |
|
| 592 | 637 |
###### Now add figures with additional met stations? |
| 593 | 638 |
|
| 594 |
############################ END OF SCRIPT ################################## |
|
| 639 |
#df_selected2 <- data_stations_temp |
|
| 640 |
|
|
| 641 |
#d_z_tmp <- zoo(df_selected[[station_id]],df_selected$date) |
|
| 642 |
#d_z_tmp2 <- zoo(df_selected2$dailyTmax,df_selected2$date) |
|
| 643 |
|
|
| 644 |
#d_z_tmp <- zoo(df[[station_id]],df$date) |
|
| 645 |
#names(d_z_tmp)<-station_id |
|
| 646 |
#names(d_z_tmp2)<-station_id |
|
| 647 |
|
|
| 648 |
#min(d_z_tmp$ID_8) |
|
| 649 |
#max(d_z_tmp$ID_8) |
|
| 650 |
#day_start <- "1984-01-01" #PARAM 12 arg 12 |
|
| 651 |
#day_end <- "2014-12-31" #PARAM 13 arg 13 |
|
| 652 |
#day_start <- "1991-01-01" #PARAM 12 arg 12 |
|
| 653 |
#day_end <- "1992-12-31" #PARAM 13 arg 13 |
|
| 654 |
|
|
| 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 |
|
|
| 660 |
#res_pix <- 1000 |
|
| 661 |
#res_pix <- 480 |
|
| 662 |
#col_mfrow <- 2 |
|
| 663 |
#row_mfrow <- 1 |
|
| 664 |
|
|
| 665 |
#png_filename <- file.path(out_dir,paste("Figure5a_time_series_profile_",region_name,"_",out_suffix,".png",sep =""))
|
|
| 666 |
#title_str <- paste("Predicted daily ", station_id," ",var," for the ", start_year,"-",end_year," time period",sep="")
|
|
| 667 |
|
|
| 668 |
#png(filename=png_filename,width = col_mfrow * res_pix,height = row_mfrow * res_pix) |
|
| 669 |
#this is the whole time series |
|
| 670 |
|
|
| 671 |
#plot(d_z_tmp,ylab="tmax in deg C",xlab="Daily time steps", |
|
| 672 |
# main=title_str,cex=3,font=2, |
|
| 673 |
# cex.main=1.5,cex.lab=1.5,font.lab=2, |
|
| 674 |
# lty=3) |
|
| 675 |
#lines(d_z_tmp2,ylab="tmax in deg C",xlab="Daily time steps", |
|
| 676 |
# main=title_str,cex=3,font=2, |
|
| 677 |
# col="red", |
|
| 678 |
# cex.main=1.5,cex.lab=1.5,font.lab=2, |
|
| 679 |
# lty=3) |
|
| 680 |
# |
|
| 681 |
#dev.off() |
|
| 682 |
|
|
| 683 |
#This is a lot of replication!!! okay cut it down |
|
| 684 |
data_stations$date <- as.Date(strptime(data_stations_temp$date_str,"%Y%m%d")) |
|
| 685 |
data_stations_tmp <- data_stations |
|
| 686 |
data_stations_tmp <- data_stations |
|
| 687 |
|
|
| 688 |
data_stations_tmp$date <- as.Date(strptime(data_stations_tmp$date,"%Y%m%d")) |
|
| 689 |
#data_stations_tmp$date <- as.Date(strptime(data_stations_tmp$date_str,"%Y%m%d")) |
|
| 690 |
d_z4 <- |
|
| 691 |
d_z_tmp4 <- zoo(data_stations_tmp$dailyTmax,data_stations_tmp$date) |
|
| 692 |
plot(d_z_tmp4,cex.main=1.5,cex.lab=1.5,font.lab=2, |
|
| 693 |
lty=3) |
|
| 694 |
lines(d_z_tmp,ylab="tmax in deg C",xlab="Daily time steps", |
|
| 695 |
main=title_str,cex=3,font=2, |
|
| 696 |
col="red", |
|
| 697 |
cex.main=1.5,cex.lab=1.5,font.lab=2, |
|
| 698 |
lty=3) |
|
| 699 |
|
|
| 700 |
day_start <- "1991-01-01" #PARAM 12 arg 12 |
|
| 701 |
day_end <- "1992-12-31" #PARAM 13 arg 13 |
|
| 702 |
|
|
| 703 |
start_date <- as.Date(day_start) |
|
| 704 |
end_date <- as.Date(day_end) |
|
| 705 |
start_year <- year(start_date) |
|
| 706 |
end_year <- year(end_date) |
|
| 707 |
d_z4 <- window(d_z_tmp4,start=start_date,end=end_date) |
|
| 708 |
|
|
| 709 |
### |
|
| 710 |
############################ END OF SCRIPT ################################## |
|
Also available in: Unified diff
aggregating training and testing ifnormation by day and quick first plot as test