Revision a5f071ae
Added by Benoit Parmentier over 8 years ago
| climate/research/oregon/interpolation/global_product_assessment_part1.R | ||
|---|---|---|
| 144 | 144 |
#run_figure_by_year <- TRUE # param 10, arg 7 |
| 145 | 145 |
list_year_predicted <- "1984,2014" |
| 146 | 146 |
scaling <- 0.01 #was scaled on 100 |
| 147 |
#if scaling is null then perform no scaling!! |
|
| 147 | 148 |
|
| 148 | 149 |
df_centroids_fname <- "/data/project/layers/commons/NEX_data/climateLayers/out/reg4/mosaic/output_reg4_1999/df_centroids_19990701_reg4_1999.txt" |
| 149 | 150 |
|
| ... | ... | |
| 420 | 421 |
df_points$year <- year_str |
| 421 | 422 |
df_points$day <- day_str |
| 422 | 423 |
|
| 423 |
|
|
| 424 | 424 |
############### PART5: Extraction of temporal profile ############# |
| 425 | 425 |
#### Extract time series from raster stack |
| 426 | 426 |
### Add dates to mkae it a date object? |
| ... | ... | |
| 448 | 448 |
df_points_day_extracted <- read.table(df_points_extracted_fname,sep=",") |
| 449 | 449 |
} |
| 450 | 450 |
|
| 451 |
pix_ts <- as.data.frame(t(df_points_day_extracted))
|
|
| 452 |
rownames(df_points_day_extracted)
|
|
| 451 |
head(df_points_day) #contains ID, dates and coordinates
|
|
| 452 |
df_points_day$id
|
|
| 453 | 453 |
|
| 454 |
lf_var <- names(r_mosaic) |
|
| 454 |
pix_ts <- as.data.frame(t(df_points_day_extracted)) |
|
| 455 |
pix_ts <- pix_ts[-1,] |
|
| 456 |
#var_names <- rownames(df_points_day_extracted) #same as lf_mosaic_list but without (*.tif) |
|
| 457 |
var_names <- rownames(pix_ts) #same as lf_mosaic_list but without (*.tif) |
|
| 458 |
var_id <- df_points_day$id |
|
| 459 |
df_points_day_extracted$id <- var_id |
|
| 460 |
|
|
| 461 |
#lf_var <- names(r_mosaic) |
|
| 455 | 462 |
### Not get the data from the time series |
| 456 |
data_pixel <- df_ts_pix[id_selected,] |
|
| 457 |
data_pixel <- as.data.frame(data_pixel) |
|
| 458 |
|
|
| 459 |
pix_ts <- t(as.data.frame(subset(data_pixel,select=r_ts_name))) #can subset to range later |
|
| 463 |
#data_pixel <- df_ts_pix[id_selected,] |
|
| 464 |
#data_pixel <- as.data.frame(data_pixel) |
|
| 465 |
#pix_ts <- t(as.data.frame(subset(data_pixel,select=r_ts_name))) #can subset to range later |
|
| 460 | 466 |
#pix_ts <- subset(as.data.frame(pix_ts),select=r_ts_name) |
| 461 |
pix_ts <- (as.data.frame(pix_ts)) |
|
| 462 |
## Process the coliform data |
|
| 463 |
|
|
| 464 |
#there are several measurements per day for some stations !!! |
|
| 465 |
#id_name <- data_pixel[[var_ID]] |
|
| 466 |
|
|
| 467 |
#df_tmp <-data_var[data_var$LOCATION_ID==id_name,] |
|
| 468 |
df_tmp <- subset(data_var,data_var$ID_stat==id_name) |
|
| 469 |
|
|
| 470 |
#df_points_extracted <- extract(r_mosaic,data_stations,df=T) |
|
| 471 |
#df_points_extracted_fname <- paste0("df_points_extracted_fname",".txt")
|
|
| 472 |
#write.table(df_points_extracted,file=df_points_extracted_fname) |
|
| 473 | 467 |
|
| 474 |
#df_points <- read.table(df_points_extracted_fname,sep=",") |
|
| 475 |
#df_points_tmp <- df_points |
|
| 476 |
#df_points <- as.data.frame(t(df_points)) |
|
| 477 |
#names(df_points) <- paste0("ID_",1:ncol(df_points))
|
|
| 468 |
list_dates_produced <- unlist(mclapply(1:length(var_names), |
|
| 469 |
FUN=extract_date, |
|
| 470 |
x=var_names, |
|
| 471 |
item_no=12, |
|
| 472 |
mc.preschedule=FALSE, |
|
| 473 |
mc.cores = num_cores)) |
|
| 478 | 474 |
|
| 479 |
#df_centroids <- read.table(df_centroids_fname,sep=",") |
|
| 475 |
list_dates_produced_date_val <- as.Date(strptime(list_dates_produced,"%Y%m%d")) |
|
| 476 |
pix_ts$date <- list_dates_produced_date_val |
|
| 477 |
pix_ts[,1]*scaling #scale? |
|
| 480 | 478 |
|
| 481 |
unique_date_tb <-table(df_points$date) |
|
| 482 |
unique_date <- unique(df_points$date) |
|
| 479 |
names(pix_ts) <- var_id |
|
| 483 | 480 |
|
| 484 |
station_id <- 8 |
|
| 485 |
var_name <-paste0("ID_",station_id)
|
|
| 481 |
pix_ts <- read.table("pix_ts2.txt",sep=",")
|
|
| 482 |
names(pix_ts) <- c(var_id,"files","date_str") |
|
| 483 |
pix_ts$date <- as.Date(strptime(pix_ts$date_str,"%Y%m%d")) |
|
| 484 |
|
|
| 485 |
#head(pix_ts) |
|
| 486 | 486 |
|
| 487 |
##Screen for unique date values |
|
| 488 |
if(max(unique_date_tb)>1){
|
|
| 489 |
# formula_str <- paste(var_name," ~ ","TRIP_START_DATE_f",sep="") |
|
| 490 |
var_pix <- aggregate(ID_8 ~ date, data = df_points, mean) #aggregate by date |
|
| 487 |
###start of plotting |
|
| 488 |
### makes this a function: |
|
| 489 |
id_selected <- "82111099999" |
|
| 490 |
#station_id <- 8 |
|
| 491 |
station_id <- id_selected |
|
| 492 |
df <- pix_ts |
|
| 493 |
#scaling |
|
| 494 |
#start_date: subset dates |
|
| 495 |
#end_date: subset dates |
|
| 496 |
|
|
| 497 |
df_selected <- subset(df,select=station_id) |
|
| 498 |
#df_selected <- subset(pix_ts,select=station_id) |
|
| 499 |
names(df_selected) <- station_id |
|
| 500 |
df_selected$date <- df$date |
|
| 501 |
|
|
| 502 |
if(!is.null(scaling)){
|
|
| 503 |
df_selected[[station_id]]<- df_selected[[station_id]]*scaling |
|
| 491 | 504 |
} |
| 492 | 505 |
|
| 493 |
var_pix$ID_8 <- var_pix$ID_8*scaling |
|
| 506 |
d_z_tmp <- zoo(df_selected[[station_id]],df_selected$date) |
|
| 507 |
#d_z_tmp <- zoo(df[[station_id]],df$date) |
|
| 508 |
names(d_z_tmp)<-station_id |
|
| 509 |
#min(d_z_tmp$ID_8) |
|
| 510 |
#max(d_z_tmp$ID_8) |
|
| 511 |
day_start <- "1984-01-01" #PARAM 12 arg 12 |
|
| 512 |
day_end <- "2014-12-31" #PARAM 13 arg 13 |
|
| 513 |
start_date <- as.Date(day_start) |
|
| 514 |
end_date <- as.Date(day_end) |
|
| 515 |
start_year <- year(start_date) |
|
| 516 |
end_year <- year(end_date) |
|
| 494 | 517 |
|
| 495 |
d_z_tmp <- zoo(var_pix$ID_8,var_pix$date) |
|
| 496 |
names(d_z_tmp)<-"ID_8" |
|
| 497 |
min(d_z_tmp$ID_8) |
|
| 498 |
max(d_z_tmp$ID_8) |
|
| 518 |
res_pix <- 1000 |
|
| 519 |
#res_pix <- 480 |
|
| 520 |
col_mfrow <- 2 |
|
| 521 |
row_mfrow <- 1 |
|
| 522 |
|
|
| 523 |
png_filename <- file.path(out_dir,paste("Figure5a_time_series_profile_",region_name,"_",out_suffix,".png",sep =""))
|
|
| 524 |
title_str <- paste("Predicted daily ", station_id," ",var," for the ", start_year,"-",end_year," time period",sep="")
|
|
| 525 |
|
|
| 526 |
png(filename=png_filename,width = col_mfrow * res_pix,height = row_mfrow * res_pix) |
|
| 527 |
#this is the whole time series |
|
| 528 |
plot(d_z_tmp,ylab="tmax in deg C",xlab="Daily time steps", |
|
| 529 |
main=title_str,cex=3,font=2, |
|
| 530 |
cex.main=1.5,cex.lab=1.5,font.lab=2, |
|
| 531 |
lty=3) |
|
| 499 | 532 |
|
| 500 |
plot(d_z_tmp) #this is the whole time series
|
|
| 533 |
dev.off()
|
|
| 501 | 534 |
|
| 502 |
day_start <- "1986-01-01" #PARAM 12 arg 12
|
|
| 535 |
day_start <- "1984-01-01" #PARAM 12 arg 12
|
|
| 503 | 536 |
day_end <- "1998-12-31" #PARAM 13 arg 13 |
| 504 | 537 |
|
| 505 | 538 |
start_date <- as.Date(day_start) |
| ... | ... | |
| 515 | 548 |
col_mfrow <- 2 |
| 516 | 549 |
row_mfrow <- 1 |
| 517 | 550 |
|
| 518 |
png_filename <- file.path(out_dir,paste("Figure5a_time_series_profile_",region_name,"_",out_suffix,".png",sep =""))
|
|
| 519 |
title_str <- paste("Predicted daily ",variable_name," for the ", start_year,"-",end_year," time period",sep="")
|
|
| 551 |
png_filename <- file.path(out_dir,paste("Figure5b_time_series_profile_",region_name,"_",out_suffix,".png",sep =""))
|
|
| 552 |
title_str <- paste("Predicted daily ", station_id," ",var," for the ", start_year,"-",end_year," time period",sep="")
|
|
| 520 | 553 |
|
| 521 | 554 |
png(filename=png_filename,width = col_mfrow * res_pix,height = row_mfrow * res_pix) |
| 522 | 555 |
|
| ... | ... | |
| 527 | 560 |
|
| 528 | 561 |
dev.off() |
| 529 | 562 |
|
| 530 |
#### Subset for 5b
|
|
| 563 |
#### Subset for 5c: zoom in
|
|
| 531 | 564 |
|
| 532 | 565 |
day_start <- "1991-01-01" #PARAM 12 arg 12 |
| 533 | 566 |
day_end <- "1992-12-31" #PARAM 13 arg 13 |
| ... | ... | |
| 544 | 577 |
col_mfrow <- 2 |
| 545 | 578 |
row_mfrow <- 1 |
| 546 | 579 |
|
| 547 |
png_filename <- file.path(out_dir,paste("Figure5b_subset_time_series_profile_",region_name,"_",out_suffix,".png",sep =""))
|
|
| 548 |
title_str <- paste("Predicted daily ",variable_name," for the ", start_year,"-",end_year," time period",sep="")
|
|
| 580 |
png_filename <- file.path(out_dir,paste("Figure5c_subset_time_series_profile_",region_name,"_",out_suffix,".png",sep =""))
|
|
| 581 |
title_str <- paste("Predicted daily ", station_id," ",var," for the ", start_year,"-",end_year," time period",sep="")
|
|
| 549 | 582 |
|
| 550 | 583 |
png(filename=png_filename,width = col_mfrow * res_pix,height = row_mfrow * res_pix) |
| 551 | 584 |
|
| ... | ... | |
| 556 | 589 |
|
| 557 | 590 |
dev.off() |
| 558 | 591 |
|
| 559 |
#data_pixel <- data_df[id_selected,] |
|
| 560 |
#data_pixel$rainfall <- as.numeric(data_pixel$rainfall) |
|
| 561 |
#d_z_tmp <-zoo(data_pixel$rainfall,as.Date(data_pixel$date)) |
|
| 562 |
#names(d_z_tmp)<- "rainfall" |
|
| 563 |
#data_pixel <- as.data.frame(data_pixel) |
|
| 564 |
#d_z_tmp2 <- zoo(data_pixel[[var_name]],as.Date(data_pixel$date)) |
|
| 565 |
|
|
| 566 |
#df_tmp <- subset(data_var,data_var$ID_stat==id_name) |
|
| 567 |
#if(da) |
|
| 568 |
|
|
| 592 |
###### Now add figures with additional met stations? |
|
| 569 | 593 |
|
| 570 | 594 |
############################ END OF SCRIPT ################################## |
Also available in: Unified diff
adding time series object for plotting extracted time series profiles from new stack of 30 years for region 4