Revision e6ad0547
Added by Benoit Parmentier about 8 years ago
| climate/research/oregon/interpolation/global_product_assessment_part0_functions.R | ||
|---|---|---|
| 197 | 197 |
return(shp_obj) |
| 198 | 198 |
} |
| 199 | 199 |
|
| 200 |
rasterize_tile_day <- function(i,list_spdf,df_missing,list_r_ref,col_name,date_val){
|
|
| 200 |
rasterize_tile_day <- function(i,list_spdf,df_missing,list_r_ref,col_name,date_val,out_dir=".",out_suffix=""){
|
|
| 201 | 201 |
# |
| 202 |
#This function creates a raster image from tiles and missing information data.frame. |
|
| 202 | 203 |
# |
| 204 |
#INPUTS: |
|
| 205 |
#1) i : counter to consider tile being processed |
|
| 206 |
#2) list_spdf: list of spatial polygon data.frame to convert in raster |
|
| 207 |
#3) df_missing: data.frame with information used in rasterization |
|
| 208 |
#4) list_r_ref: reference raster to use |
|
| 209 |
#5) col_name: name of column containing value used in the rasteriation, the column is stored in the df_missing data.frame |
|
| 210 |
#6) date_val: |
|
| 211 |
#7) out_dir: |
|
| 212 |
#8) out_suffix |
|
| 213 |
#OUTPUTS |
|
| 214 |
#1) raster_name: output raster_name |
|
| 215 |
# |
|
| 216 |
|
|
| 217 |
#### Start script #### |
|
| 218 |
|
|
| 203 | 219 |
tile_spdf <- list_spdf[[i]] |
| 204 | 220 |
tile_coord <- names(list_spdf)[i] |
| 205 | 221 |
r_ref <- list_r_ref[[i]] |
| ... | ... | |
| 233 | 249 |
r <- init(r, fun=set1f, overwrite=TRUE) |
| 234 | 250 |
} |
| 235 | 251 |
|
| 236 |
out_dir <- "." #can set this up later |
|
| 237 |
out_suffix_str <- paste0(col_name,"") # can set this parameter later
|
|
| 252 |
#out_dir <- "." #can set this up later
|
|
| 253 |
out_suffix_str <- paste0(col_name,out_suffix) # can set this parameter later
|
|
| 238 | 254 |
raster_name <- file.path(out_dir,paste("r_",tile_id,"_",tile_coord,"_",out_suffix_str,".tif",sep=""))
|
| 239 | 255 |
#raster_name <- |
| 240 | 256 |
writeRaster(r, NAflag=NA_flag_val,filename=raster_name,overwrite=TRUE) #unweighted mean |
| ... | ... | |
| 243 | 259 |
} |
| 244 | 260 |
|
| 245 | 261 |
generate_raster_number_of_prediction_by_day <- function(i,list_param){
|
| 246 |
|
|
| 262 |
# |
|
| 247 | 263 |
##This function generates raster of missing pixels and number of predictions for days with missing tiles for a given region. |
| 264 |
# |
|
| 248 | 265 |
#INPUTS |
| 249 | 266 |
#1) list_tiles_predicted_masked |
| 250 | 267 |
#2) df_missing_tiles_day |
| ... | ... | |
| 255 | 272 |
#7) scaling <- list_param$scaling |
| 256 | 273 |
#8) python_bin <- list_param$python_bin |
| 257 | 274 |
#9) data_type <- list_param$data_type |
| 258 |
#10) out_suffix |
|
| 259 |
#11) out_dir |
|
| 275 |
#10) plotting figures: if True plot png for missing day predictions |
|
| 276 |
#11) out_suffix |
|
| 277 |
#12) out_dir |
|
| 278 |
#OUTPUTS |
|
| 260 | 279 |
|
| 261 | 280 |
|
| 262 | 281 |
###### Start script ##### |
| ... | ... | |
| 272 | 291 |
scaling <- list_param$scaling |
| 273 | 292 |
python_bin <- list_param$python_bin |
| 274 | 293 |
data_type <- list_param$data_type |
| 294 |
plotting_figures <- list_param$plotting_figures |
|
| 275 | 295 |
out_suffix <- list_param$out_suffix |
| 276 | 296 |
out_dir <- list_param$out_dir |
| 277 | 297 |
|
| ... | ... | |
| 300 | 320 |
##### Sum missing tiles in the stack and generate number of predictions by pixels |
| 301 | 321 |
## This stores files in the temp dir |
| 302 | 322 |
raster_name_data_sum <- file.path(out_dir,paste("r_data_sum","_",region_name,"_masked_",date_str,file_format,sep=""))
|
| 303 |
r_data_sum <- stackApply(r_tiles_s, 1:nlayers(r_tiles_s), fun = sum,filename=raster_name_data_sum) |
|
| 323 |
r_data_sum <- stackApply(r_tiles_s, 1:nlayers(r_tiles_s), fun = sum,filename=raster_name_data_sum,overwrite=TRUE)
|
|
| 304 | 324 |
|
| 305 | 325 |
### then substract missing tiles... |
| 306 | 326 |
raster_name_number_prediction <- file.path(out_dir,paste("r_day_number_of_prediction_sum_day_mosaiced","_",region_name,"_masked_",date_str,file_format,sep=""))
|
| ... | ... | |
| 335 | 355 |
#rat$legend <- tb_freq$value |
| 336 | 356 |
#levels(r_table) <- rat |
| 337 | 357 |
|
| 338 |
res_pix <- 800 |
|
| 339 |
#res_pix <- 480 |
|
| 340 |
col_mfrow <- 1 |
|
| 341 |
row_mfrow <- 1 |
|
| 342 |
|
|
| 343 |
png_filename_number_of_predictions <- file.path(out_dir,paste("Figure_number_of_predictions_by_pixel_",date_str,"_",region_name,"_",out_suffix,".png",sep =""))
|
|
| 358 |
if(plotting_figures==TRUE){
|
|
| 344 | 359 |
|
| 345 |
title_str <- paste("Number of predicted pixels for ",variable_name," on ",date_str, sep = "")
|
|
| 360 |
res_pix <- 800 |
|
| 361 |
#res_pix <- 480 |
|
| 362 |
col_mfrow <- 1 |
|
| 363 |
row_mfrow <- 1 |
|
| 346 | 364 |
|
| 347 |
png(filename=png_filename_number_of_predictions,width = col_mfrow * res_pix,height = row_mfrow * res_pix) |
|
| 348 |
#my_col=c('blue','red','green')
|
|
| 349 |
my_col <- rainbow(length(tb_freq$value)) |
|
| 350 |
plot(r_day_predicted,col=my_col,legend=F,box=F,axes=F,main=title_str) |
|
| 351 |
legend(x='topright', legend =tb_freq$value,fill = my_col,cex=0.8) |
|
| 352 |
dev.off() |
|
| 365 |
png_filename_number_of_predictions <- file.path(out_dir,paste("Figure_number_of_predictions_by_pixel_",date_str,"_",region_name,"_",out_suffix,".png",sep =""))
|
|
| 366 |
title_str <- paste("Number of predicted pixels for ",variable_name," on ",date_str, sep = "")
|
|
| 367 |
|
|
| 368 |
png(filename=png_filename_number_of_predictions,width = col_mfrow * res_pix,height = row_mfrow * res_pix) |
|
| 369 |
#my_col=c('blue','red','green')
|
|
| 370 |
my_col <- rainbow(length(tb_freq$value)) |
|
| 371 |
plot(r_day_predicted,col=my_col,legend=F,box=F,axes=F,main=title_str) |
|
| 372 |
legend(x='topright', legend =tb_freq$value,fill = my_col,cex=0.8) |
|
| 373 |
dev.off() |
|
| 374 |
}else{
|
|
| 375 |
png_filename_number_of_predictions <- NULL |
|
| 376 |
} |
|
| 353 | 377 |
|
| 354 | 378 |
### Day missing reclass above |
| 355 | 379 |
## change here |
| ... | ... | |
| 370 | 394 |
|
| 371 | 395 |
r_missing_day <- raster(raster_name_missing) |
| 372 | 396 |
|
| 373 |
res_pix <- 800 |
|
| 374 |
#res_pix <- 480 |
|
| 375 |
col_mfrow <- 1 |
|
| 376 |
row_mfrow <- 1 |
|
| 377 |
|
|
| 378 |
png_filename_missing_predictions <- file.path(out_dir,paste("Figure_missing_predictions_by_pixel_",date_str,"_",region_name,"_",out_suffix,".png",sep =""))
|
|
| 397 |
if(plotting_figures==TRUE){
|
|
| 379 | 398 |
|
| 380 |
title_str <- paste("Number of predicted pixels for ",variable_name," on ",date_str, sep = "")
|
|
| 399 |
res_pix <- 800 |
|
| 400 |
#res_pix <- 480 |
|
| 401 |
col_mfrow <- 1 |
|
| 402 |
row_mfrow <- 1 |
|
| 381 | 403 |
|
| 382 |
png(filename=png_filename_missing_predictions,width = col_mfrow * res_pix,height = row_mfrow * res_pix) |
|
| 383 |
#my_col=c('blue','red','green')
|
|
| 384 |
my_col <- c("black","red")
|
|
| 385 |
plot(r_missing_day,col=my_col,legend=F,box=F,axes=F,main=title_str) |
|
| 386 |
legend(x='topright', legend =c("prediced","missing"),fill = my_col,cex=0.8)
|
|
| 404 |
png_filename_missing_predictions <- file.path(out_dir,paste("Figure_missing_predictions_by_pixel_",date_str,"_",region_name,"_",out_suffix,".png",sep =""))
|
|
| 405 |
title_str <- paste("Number of predicted pixels for ",variable_name," on ",date_str, sep = "")
|
|
| 387 | 406 |
|
| 388 |
dev.off() |
|
| 407 |
png(filename=png_filename_missing_predictions,width = col_mfrow * res_pix,height = row_mfrow * res_pix) |
|
| 408 |
#my_col=c('blue','red','green')
|
|
| 409 |
my_col <- c("black","red")
|
|
| 410 |
plot(r_missing_day,col=my_col,legend=F,box=F,axes=F,main=title_str) |
|
| 411 |
legend(x='topright', legend =c("prediced","missing"),fill = my_col,cex=0.8)
|
|
| 412 |
|
|
| 413 |
dev.off() |
|
| 414 |
|
|
| 415 |
}else{
|
|
| 416 |
png_filename_missing_predictions <- NULL |
|
| 417 |
} |
|
| 389 | 418 |
|
| 390 | 419 |
### writeout data |
| 391 | 420 |
#extension_str <- extension(lf_files) |
| ... | ... | |
| 441 | 470 |
|
| 442 | 471 |
########################## START SCRIPT ######################################### |
| 443 | 472 |
|
| 473 |
#browser() |
|
| 444 | 474 |
#system("ls /nobackup/bparmen1"
|
| 445 |
out_dir <- in_dir
|
|
| 475 |
#out_dir <- in_dir #use directory of out_dir
|
|
| 446 | 476 |
if(create_out_dir_param==TRUE){
|
| 447 | 477 |
out_dir <- create_dir_fun(out_dir,out_suffix) |
| 448 | 478 |
setwd(out_dir) |
| ... | ... | |
| 539 | 569 |
#undebug(check_missing) |
| 540 | 570 |
|
| 541 | 571 |
test_missing <- try(lapply(1:length(list_lf_raster_tif_tiles),function(i){check_missing(lf=list_lf_raster_tif_tiles[[i]],
|
| 542 |
pattern_str=NULL, |
|
| 543 |
in_dir=out_dir,
|
|
| 544 |
date_start=start_date, |
|
| 545 |
date_end=end_date, |
|
| 546 |
item_no=item_no, #9 for predicted tiles |
|
| 547 |
out_suffix=out_suffix, |
|
| 548 |
num_cores=num_cores, |
|
| 549 |
out_dir=".")}))
|
|
| 572 |
pattern_str=NULL,
|
|
| 573 |
in_dir=in_dir,
|
|
| 574 |
date_start=start_date,
|
|
| 575 |
date_end=end_date,
|
|
| 576 |
item_no=item_no, #9 for predicted tiles
|
|
| 577 |
out_suffix=out_suffix,
|
|
| 578 |
num_cores=num_cores,
|
|
| 579 |
out_dir=out_dir)}))
|
|
| 550 | 580 |
|
| 551 | 581 |
|
| 552 | 582 |
#test_missing <- try(lapply(1:1,function(i){check_missing(lf=list_lf_raster_tif_tiles[[i]],
|
| ... | ... | |
| 576 | 606 |
#http://stackoverflow.com/questions/26220913/replace-na-with-na |
| 577 | 607 |
#Use dfr[dfr=="<NA>"]=NA where dfr is your dataframe. |
| 578 | 608 |
names(df_lf_tiles_time_series) <- unlist(basename(in_dir_reg)) |
| 579 |
filename_df_lf_tiles <- paste0("df_files_by_tiles_predicted_tif_",region_name,"_",pred_mod_name,"_",out_suffix)
|
|
| 609 |
filename_df_lf_tiles <- file.path(out_dir,paste0("df_files_by_tiles_predicted_tif_",region_name,"_",pred_mod_name,"_",out_suffix))
|
|
| 580 | 610 |
write.table(df_lf_tiles_time_series,file=filename_df_lf_tiles) |
| 581 | 611 |
|
| 582 | 612 |
###Now combined missing in one table? |
| ... | ... | |
| 590 | 620 |
df_missing$reg <- region_name |
| 591 | 621 |
df_missing$date <- day_to_mosaic |
| 592 | 622 |
|
| 593 |
filename_df_missing <- paste0("df_missing_by_tiles_predicted_tif_",region_name,"_",pred_mod_name,"_",out_suffix)
|
|
| 623 |
filename_df_missing <- file.path(out_dir,paste0("df_missing_by_tiles_predicted_tif_",region_name,"_",pred_mod_name,"_",out_suffix))
|
|
| 594 | 624 |
write.table(df_missing,file=filename_df_missing) |
| 595 | 625 |
|
| 596 | 626 |
######################## |
| ... | ... | |
| 678 | 708 |
list_r_ref=list_r_ref, |
| 679 | 709 |
col_name = "overlap", |
| 680 | 710 |
date_val=df_missing$date[1], |
| 711 |
out_dir = out_dir, |
|
| 712 |
out_suffix = "", |
|
| 681 | 713 |
mc.preschedule=FALSE, |
| 682 | 714 |
mc.cores = num_cores) |
| 683 | 715 |
|
| ... | ... | |
| 808 | 840 |
#r_tiles_s <- r_tiles_stack |
| 809 | 841 |
names_tiles <- basename(in_dir_reg) |
| 810 | 842 |
|
| 811 |
|
|
| 843 |
browser() |
|
| 812 | 844 |
list_param_generate_raster_number_pred <- list(list_tiles_predicted_masked,df_missing_tiles_day,r_overlap_m, |
| 813 | 845 |
num_cores,region_name,data_type,scaling,python_bin, |
| 846 |
plotting_figures, |
|
| 814 | 847 |
NA_flag_val,out_suffix,out_dir) |
| 815 | 848 |
|
| 816 | 849 |
names(list_param_generate_raster_number_pred) <- c("list_tiles_predicted_masked","df_missing_tiles_day","r_overlap_m",
|
| 817 | 850 |
"num_cores","region_name","data_type","scaling","python_bin", |
| 851 |
"plotting_figures", |
|
| 818 | 852 |
"NA_flag_val","out_suffix","out_dir") |
| 819 | 853 |
|
| 820 | 854 |
#function_product_assessment_part0_functions <- "global_product_assessment_part0_functions_11152016b.R" |
| ... | ... | |
| 825 | 859 |
#browser() |
| 826 | 860 |
#5.10pm |
| 827 | 861 |
#test_number_pix_predictions <- generate_raster_number_of_prediction_by_day(1,list_param=list_param_generate_raster_number_pred) |
| 862 |
|
|
| 828 | 863 |
if(nrow(df_missing_tiles_day)>0){
|
| 829 | 864 |
|
| 830 | 865 |
obj_number_pix_predictions <- mclapply(1:nrow(df_missing_tiles_day), |
| ... | ... | |
| 850 | 885 |
tb_freq_overlap,png_filename_maximum_overlap,obj_number_pix_predictions) |
| 851 | 886 |
names(predictions_tiles_missing_obj) <- c("df_lf_tiles_time_series","df_missing_tiles_day","raster_name_overlap",
|
| 852 | 887 |
"tb_freq_overlap","png_maximum_overlap","obj_number_pix_predictions") |
| 853 |
browser() |
|
| 888 |
|
|
| 854 | 889 |
|
| 855 | 890 |
return(predictions_tiles_missing_obj) |
| 856 | 891 |
} |
Also available in: Unified diff
adding plotting_figures parameter to deal with error in png creation in forks on pfe NEX