Revision f101c997
Added by Benoit Parmentier over 8 years ago
| climate/research/oregon/interpolation/global_product_assessment_part1.R | ||
|---|---|---|
| 376 | 376 |
c( -115.330122,36.677139), |
| 377 | 377 |
c( -78.396011,35.696143)) |
| 378 | 378 |
|
| 379 |
|
|
| 380 | 379 |
test_day_query2 <- lapply(list_lat_long,FUN=query_for_station_lat_long,df_points_spdf=df_point_day,step_x=1,step_y=1) |
| 381 | 380 |
#test_day_query <-query_for_station_lat_long(c(-72,42),df_points_spdf=df_point_day,step_x=1,step_y=0.25) |
| 382 | 381 |
df_stations_selected <- do.call(rbind,test_day_query2) |
| 383 | 382 |
proj4string(df_stations_selected) <- proj_str |
| 384 | 383 |
#debug(query_for_station_lat_long) |
| 385 | 384 |
|
| 385 |
df_stations_locations <- do.call(rbind,list_lat_long) #these are the stations to match |
|
| 386 |
##layer to be clipped |
|
| 387 |
if(class(countries_shp)!="SpatialPolygonsDataFrame"){
|
|
| 388 |
reg_layer <- readOGR(dsn=dirname(countries_shp),sub(".shp","",basename(countries_shp)))
|
|
| 389 |
}else{
|
|
| 390 |
reg_layer <- countries_shp |
|
| 391 |
} |
|
| 392 |
df_stations_selected$lab_id <- 1:nrow(df_stations_selected) |
|
| 393 |
#now clip |
|
| 394 |
#use points that were used to extract |
|
| 395 |
reg_layer_clipped <- gClip(shp=reg_layer, bb=df_points , keep.attribs=TRUE,outDir=NULL,outSuffix=NULL) |
|
| 396 |
|
|
| 397 |
plot(reg_layer_clipped) |
|
| 398 |
plot(df_stations_selected,add=T,col="blue",cex=2) |
|
| 399 |
text(df_stations_selected$x,df_stations_selected$y, df_stations_selected$id, col="red", cex=0.8) |
|
| 400 |
|
|
| 386 | 401 |
##Next use the day results and select a mean station, quartile and min and max? |
| 387 | 402 |
|
| 388 | 403 |
#list_id_data_v <- unique(data_stations_var_pred_data_v$id) |
| ... | ... | |
| 595 | 610 |
#Product assessment |
| 596 | 611 |
out_suffix_str <- paste0(region_name,"_",out_suffix) |
| 597 | 612 |
#this can be run with mclapply, very fast right now: |
| 598 |
station_summary_obj <- combine_measurements_and_predictions_df(i=i, |
|
| 613 |
#station_summary_obj <- combine_measurements_and_predictions_df(i=i, |
|
| 614 |
# df_raster=df_raster, |
|
| 615 |
# df_time_series=df_time_series, |
|
| 616 |
# df_ts_pix=df_ts_pix, |
|
| 617 |
# data_var=data_var, |
|
| 618 |
# list_selected_ID=list_selected_ID, |
|
| 619 |
# r_ts_name=r_ts_name, |
|
| 620 |
# var_name=var_name, |
|
| 621 |
# var_pred = var_pred, |
|
| 622 |
# out_dir =out_dir, |
|
| 623 |
# out_suffix=out_suffix_str, |
|
| 624 |
# plot_fig=F) |
|
| 625 |
df_pix_ts <- station_summary_obj$df_pix_ts |
|
| 626 |
|
|
| 627 |
##### |
|
| 628 |
##combine information for the 11 selected stations |
|
| 629 |
|
|
| 630 |
list_station_summary_obj <- mclapply(1:length(list_selected_ID), |
|
| 631 |
FUN=combine_measurements_and_predictions_df, |
|
| 599 | 632 |
df_raster=df_raster, |
| 600 | 633 |
df_time_series=df_time_series, |
| 601 | 634 |
df_ts_pix=df_ts_pix, |
| ... | ... | |
| 606 | 639 |
var_pred = var_pred, |
| 607 | 640 |
out_dir =out_dir, |
| 608 | 641 |
out_suffix=out_suffix_str, |
| 609 |
plot_fig=F) |
|
| 610 |
df_pix_ts <- station_summary_obj$df_pix_ts |
|
| 642 |
plot_fig=F, |
|
| 643 |
mc.preschedule=FALSE, |
|
| 644 |
mc.cores = num_cores) |
|
| 645 |
|
|
| 611 | 646 |
#station_summary_obj <- list(nb_zero,nb_NA, df_pix_ts) |
| 612 | 647 |
|
| 613 | 648 |
#id_name <- list_selected_ID[i] |
| ... | ... | |
| 627 | 662 |
var_name1 <- y_var_name |
| 628 | 663 |
var_name2 <- "mod1_mosaic" |
| 629 | 664 |
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
|
|
| 665 |
scaling_factors <- c(NA,scaling) #no scaling for y_var_name but scaling for var_name2
|
|
| 631 | 666 |
list_windows <- list(c("1999-01-01","1999-12-31"))
|
| 632 | 667 |
|
| 633 |
function_product_assessment_part1_functions <- "global_product_assessment_part1_functions_09202016.R" |
|
| 668 |
function_product_assessment_part1_functions <- "global_product_assessment_part1_functions_09202016b.R"
|
|
| 634 | 669 |
source(file.path(script_path,function_product_assessment_part1_functions)) #source all functions used in this script |
| 635 | 670 |
|
| 636 | 671 |
###TO DO: |
| ... | ... | |
| 638 | 673 |
#2) Compute temporal autocorrelation by station and profile |
| 639 | 674 |
#3) compute range, min, max etc by stations |
| 640 | 675 |
|
| 641 |
undebug(plot_observation_predictions_time_series) |
|
| 676 |
#undebug(plot_observation_predictions_time_series)
|
|
| 642 | 677 |
#test_plot <- plot_observation_predictions_time_series(df_pix_time_series=df_pix_ts, |
| 643 | 678 |
# var_name1=var_name1, |
| 644 | 679 |
# var_name2=var_name2, |
| ... | ... | |
| 649 | 684 |
# out_suffix=out_suffix_str) |
| 650 | 685 |
out_suffix_str2 <- paste(region_name,out_suffix,"test",sep="_") |
| 651 | 686 |
|
| 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) |
|
| 687 |
#list_df_pix_ts <-extract_from_list_obj(list_station_summary_obj,"df_pix_ts") |
|
| 688 |
list_df_pix_ts <- lapply(1:length(list_station_summary_obj),FUN=function(i){list_station_summary_obj[[i]]$df_pix_ts})
|
|
| 689 |
|
|
| 690 |
debug(plot_observation_predictions_time_series) |
|
| 691 |
test_plot <- mclapply(list_df_pix_ts, |
|
| 692 |
FUN=plot_observation_predictions_time_series, |
|
| 693 |
var_name1=var_name1, |
|
| 694 |
var_name2=var_name2, |
|
| 695 |
list_windows=list_windows, |
|
| 696 |
time_series_id=NULL,#read in from data |
|
| 697 |
scaling_factors=scaling_factors, |
|
| 698 |
out_dir=out_dir, |
|
| 699 |
out_suffix=out_suffix_str2, |
|
| 700 |
mc.preschedule=FALSE, |
|
| 701 |
mc.cores = num_cores) |
|
| 660 | 702 |
|
| 661 | 703 |
############### PART5: Make raster stack and display maps ############# |
| 662 | 704 |
#### Extract corresponding raster for given dates and plot stations used |
Also available in: Unified diff
running combining function and time series plotting for 11 selected stations for env layers meeting