Revision d0d1abc2
Added by Benoit Parmentier over 8 years ago
| climate/research/oregon/interpolation/global_product_assessment_part0_functions.R | ||
|---|---|---|
| 9 | 9 |
# |
| 10 | 10 |
#AUTHOR: Benoit Parmentier |
| 11 | 11 |
#CREATED ON: 10/31/2016 |
| 12 |
#MODIFIED ON: 11/27/2016
|
|
| 12 |
#MODIFIED ON: 11/28/2016
|
|
| 13 | 13 |
#Version: 1 |
| 14 | 14 |
#PROJECT: Environmental Layers project |
| 15 | 15 |
#COMMENTS: removing unused functions and clean up for part0 global prodduct assessment part0 |
| ... | ... | |
| 26 | 26 |
# |
| 27 | 27 |
#setfacl -Rmd user:aguzman4:rwx /nobackupp8/bparmen1/output_run10_1500x4500_global_analyses_pred_1992_10052015 |
| 28 | 28 |
|
| 29 |
##COMMIT: modification of inputs for checking missing
|
|
| 29 |
##COMMIT: cleaning outputs generated by R from temporary dir
|
|
| 30 | 30 |
|
| 31 | 31 |
################################################################################################# |
| 32 | 32 |
|
| ... | ... | |
| 286 | 286 |
r_tiles_s <- stack(list_missing_tiles_raster) |
| 287 | 287 |
|
| 288 | 288 |
##### Sum missing tiles in the stack and generate number of predictions by pixels |
| 289 |
|
|
| 290 |
datasum <- stackApply(r_tiles_s, 1:nlayers(r_tiles_s), fun = sum) |
|
| 289 |
## This stores files in the temp dir |
|
| 290 |
raster_name_data_sum <- file.path(out_dir,paste("r_data_sum","_",region_name,"_masked_",date_str,file_format,sep=""))
|
|
| 291 |
r_data_sum <- stackApply(r_tiles_s, 1:nlayers(r_tiles_s), fun = sum,filename=raster_name_data_sum) |
|
| 291 | 292 |
|
| 292 | 293 |
### then substract missing tiles... |
| 293 |
r_day_predicted <- r_overlap_m - datasum
|
|
| 294 |
#r_day_predicted <- r_overlap_m - r_data_sum
|
|
| 294 | 295 |
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=""))
|
| 295 |
writeRaster(r_day_predicted, NAflag=NA_flag_val,filename=raster_name_number_prediction,overwrite=TRUE) |
|
| 296 | 296 |
|
| 297 |
r_table <- ratify(r_day_predicted) # build the Raster Attibute table |
|
| 298 |
rat <- levels(r_table)[[1]]#get the values of the unique cell frot the attribute table |
|
| 297 |
#raster_name_day_predicted <- file.path(out_dir,paste("r_day_predicted","_",region_name,"_masked_",date_str,file_format,sep=""))
|
|
| 298 |
## do this in gdalcalc or overlay function to go faster? |
|
| 299 |
python_cmd <- file.path(python_bin,"gdal_calc.py") |
|
| 300 |
cmd_str1 <- paste(python_cmd, |
|
| 301 |
paste("-A ", filename(r_overlap_m),sep=""),
|
|
| 302 |
paste("-B ", raster_name_data_sum,sep=""),
|
|
| 303 |
paste("--outfile=", raster_name_number_prediction,sep=""),
|
|
| 304 |
paste("--type=",data_type,sep=""),
|
|
| 305 |
"--co='COMPRESS=LZW'", |
|
| 306 |
paste("--NoDataValue=",NA_flag_val,sep=""),
|
|
| 307 |
paste("--calc='(A-B)*",scaling,"'",sep=""),
|
|
| 308 |
"--overwrite",sep=" ") #division by zero is problematic... |
|
| 309 |
system(cmd_str1) |
|
| 310 |
|
|
| 311 |
#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=""))
|
|
| 312 |
#writeRaster(r_day_predicted, NAflag=NA_flag_val,filename=raster_name_number_prediction,overwrite=TRUE) |
|
| 313 |
r_day_predicted <- raster(raster_name_number_prediction) |
|
| 314 |
### Uses too much resources, change here... |
|
| 315 |
#r_table <- ratify(r_day_predicted) # build the Raster Attibute table |
|
| 316 |
#rat <- levels(r_table)[[1]]#get the values of the unique cell frot the attribute table |
|
| 317 |
tb_freq <- as.data.frame(freq(r_day_predicted)) |
|
| 318 |
|
|
| 299 | 319 |
#rat$legend <- paste0("tile_",1:26)
|
| 300 |
tb_freq <- as.data.frame(freq(r_table)) |
|
| 301 |
rat$legend <- tb_freq$value |
|
| 302 |
levels(r_table) <- rat |
|
| 320 |
#tb_freq <- as.data.frame(freq(r_table))
|
|
| 321 |
#rat$legend <- tb_freq$value
|
|
| 322 |
#levels(r_table) <- rat
|
|
| 303 | 323 |
|
| 304 | 324 |
res_pix <- 800 |
| 305 | 325 |
#res_pix <- 480 |
| ... | ... | |
| 313 | 333 |
png(filename=png_filename,width = col_mfrow * res_pix,height = row_mfrow * res_pix) |
| 314 | 334 |
#my_col=c('blue','red','green')
|
| 315 | 335 |
my_col <- rainbow(length(tb_freq$value)) |
| 316 |
plot(r_table,col=my_col,legend=F,box=F,axes=F,main=title_str)
|
|
| 336 |
plot(r_day_predicted,col=my_col,legend=F,box=F,axes=F,main=title_str)
|
|
| 317 | 337 |
legend(x='topright', legend =rat$legend,fill = my_col,cex=0.8) |
| 318 | 338 |
|
| 319 | 339 |
dev.off() |
| 320 | 340 |
|
| 321 | 341 |
### Day missing reclass above |
| 342 |
## change here |
|
| 343 |
r_missing_day <- r_day_predicted == 0 |
|
| 322 | 344 |
|
| 323 | 345 |
## do this in gdalcalc or overlay function to go faster? |
| 324 |
r_missing_day <- r_day_predicted == 0 |
|
| 346 |
python_cmd <- file.path(python_bin,"gdal_calc.py") |
|
| 347 |
cmd_str2 <- paste(python_cmd, |
|
| 348 |
paste("-A ", raster_name_number_prediction,sep=""),
|
|
| 349 |
paste("--outfile=",r_m_weighted_mean_raster_name,sep=""),
|
|
| 350 |
paste("--type=",data_type,sep=""),
|
|
| 351 |
"--co='COMPRESS=LZW'", |
|
| 352 |
paste("--NoDataValue=",NA_flag_val,sep=""),
|
|
| 353 |
paste("--calc='(A=0)*",scaling,"'",sep=""),
|
|
| 354 |
"--overwrite",sep=" ") #division by zero is problematic... |
|
| 355 |
system(cmd_str2) |
|
| 325 | 356 |
|
| 326 | 357 |
res_pix <- 800 |
| 327 | 358 |
#res_pix <- 480 |
| ... | ... | |
| 347 | 378 |
raster_name_missing <- file.path(out_dir,paste("r_missing_day_mosaiced","_",region_name,"_masked_",date_str,file_format,sep=""))
|
| 348 | 379 |
writeRaster(r_missing_day, NAflag=NA_flag_val,filename=raster_name_missing,overwrite=TRUE) |
| 349 | 380 |
|
| 381 |
##Took 13-14 minutes for 28 tiles and one date...!!! |
|
| 382 |
|
|
| 383 |
#if(tmp_files==F){ #if false...delete all files with "_tmp"
|
|
| 384 |
# lf_tmp <- unlist(lf_accuracy_training_raster) |
|
| 385 |
# ##now delete temporary files... |
|
| 386 |
# file.remove(lf_tmp) |
|
| 387 |
#} |
|
| 388 |
|
|
| 350 | 389 |
### generate return object |
| 351 | 390 |
obj_number_day_predicted <- list(raster_name_number_prediction,raster_name_missing) |
| 352 | 391 |
names(obj_number_day_predicted) <- c("raster_name_number_prediction","raster_name_missing")
|
| ... | ... | |
| 522 | 561 |
|
| 523 | 562 |
list_missing <- lapply(1:length(test_missing),FUN=function(i){df_time_series <- test_missing[[i]]$df_time_series$missing})
|
| 524 | 563 |
|
| 564 |
browser() |
|
| 525 | 565 |
df_missing <- as.data.frame(do.call(cbind,list_missing)) |
| 526 | 566 |
names(df_missing) <- unlist(basename(in_dir_reg)) |
| 527 | 567 |
df_missing$tot_missing <- rowSums (df_missing, na.rm = FALSE, dims = 1) |
| ... | ... | |
| 746 | 786 |
|
| 747 | 787 |
#undebug(generate_raster_number_of_prediction_by_day) |
| 748 | 788 |
#4.51pm |
| 749 |
#browser()
|
|
| 789 |
browser() |
|
| 750 | 790 |
#5.10pm |
| 751 | 791 |
#test_number_pix_predictions <- generate_raster_number_of_prediction_by_day(1,list_param=list_param_generate_raster_number_pred) |
| 752 | 792 |
if(nrow(df_missing_tiles_day)>0){
|
Also available in: Unified diff
cleaning outputs generated by R from temporary dir