Revision 5445585f
Added by Benoit Parmentier over 8 years ago
| climate/research/oregon/interpolation/global_product_assessment_part0.R | ||
|---|---|---|
| 9 | 9 |
# |
| 10 | 10 |
#AUTHOR: Benoit Parmentier |
| 11 | 11 |
#CREATED ON: 10/27/2016 |
| 12 |
#MODIFIED ON: 11/10/2016
|
|
| 12 |
#MODIFIED ON: 11/14/2016
|
|
| 13 | 13 |
#Version: 1 |
| 14 | 14 |
#PROJECT: Environmental Layers project |
| 15 | 15 |
#COMMENTS: |
| ... | ... | |
| 20 | 20 |
#source /nobackupp6/aguzman4/climateLayers/sharedModules2/etc/environ.sh |
| 21 | 21 |
# |
| 22 | 22 |
#setfacl -Rm u:aguzman4:rwx /nobackupp6/aguzman4/climateLayers/LST_tempSpline/ |
| 23 |
#COMMIT: modifying function generate raster of number of predictions for day with missing tiles
|
|
| 23 |
#COMMIT: adding and debugging function to generate raster of number of predictions for day with missing tiles
|
|
| 24 | 24 |
|
| 25 | 25 |
################################################################################################# |
| 26 | 26 |
|
| ... | ... | |
| 459 | 459 |
out_mosaic_name_overlap_masked <- file.path(out_dir_str,paste("r_overlap_sum_masked_",out_suffix_str_tmp,".tif",sep=""))
|
| 460 | 460 |
|
| 461 | 461 |
r_overlap_m <- mask(r_overlap,r_mask,filename=out_mosaic_name_overlap_masked) |
| 462 |
plot(r_overlap_m) |
|
| 462 |
#plot(r_overlap_m)
|
|
| 463 | 463 |
#plot(spdf_tiles_test,add=T,border="green",usePolypath = FALSE) #added usePolypath following error on brige and NEX |
| 464 | 464 |
|
| 465 | 465 |
r_table <- ratify(r_overlap_m) # build the Raster Attibute table |
| ... | ... | |
| 467 | 467 |
#rat$legend <- paste0("tile_",1:26)
|
| 468 | 468 |
tb_freq <- as.data.frame(freq(r_table)) |
| 469 | 469 |
rat$legend <- tb_freq$value |
| 470 |
|
|
| 471 | 470 |
levels(r_table) <- rat |
| 472 |
#my_col=c('blue','red','green')
|
|
| 471 |
|
|
| 472 |
|
|
| 473 |
res_pix <- 800 |
|
| 474 |
#res_pix <- 480 |
|
| 475 |
col_mfrow <- 1 |
|
| 476 |
row_mfrow <- 1 |
|
| 477 |
|
|
| 478 |
png_filename <- file.path(out_dir,paste("Figure_maximum_overlap_",region_name,"_",out_suffix,".png",sep =""))
|
|
| 479 |
|
|
| 480 |
title_str <- paste("Maximum overlap: Number of predicted pixels for ",variable_name, sep = "")
|
|
| 481 |
|
|
| 482 |
png(filename=png_filename,width = col_mfrow * res_pix,height = row_mfrow * res_pix) |
|
| 483 |
#my_col=c('blue','red','green')
|
|
| 473 | 484 |
my_col <- rainbow(length(tb_freq$value)) |
| 474 | 485 |
|
| 475 |
plot(r_table,col=my_col,legend=F,box=F,axes=F) |
|
| 486 |
plot(r_table,col=my_col,legend=F,box=F,axes=F,main=title_str)
|
|
| 476 | 487 |
legend(x='topright', legend =rat$legend,fill = my_col,cex=0.8) |
| 477 | 488 |
|
| 489 |
dev.off() |
|
| 490 |
|
|
| 478 | 491 |
#r <- raster(matrix(runif(20),5,4)) |
| 479 | 492 |
#r[r>.5] <- NA |
| 480 | 493 |
#g <- as(r, 'SpatialGridDataFrame') |
| ... | ... | |
| 532 | 545 |
|
| 533 | 546 |
list_names_tile_coord |
| 534 | 547 |
df_time_series |
| 535 |
missing_tiles <- df_missing_tiles_day[i] |
|
| 548 |
missing_tiles <- df_missing_tiles_day[i,] |
|
| 549 |
date_str <- missing_tiles$date |
|
| 550 |
|
|
| 536 | 551 |
#r_tiles_s <- list_param$r_tiles_s |
| 537 | 552 |
list_param$list_tiles_predicted_masked |
| 538 | 553 |
|
| 539 |
#df_missing_tiles_day <- subset(df_missing,tot_missing > 0) |
|
| 540 |
#stack() ## all tiles for the day |
|
| 541 |
#df_missing_tiles_day[,-c("tot_missing")]
|
|
| 542 |
df_missing_tiles_day[,!c("tot_missing")]
|
|
| 543 | 554 |
selected_col <- names(list_tiles_predicted_masked) |
| 544 |
df_missing_tiles_day_subset <- subset(df_missing_tiles_day,select=selected_col) |
|
| 545 |
#drops <- c("x","z")
|
|
| 546 |
#DF[ , !(names(DF) %in% drops)] |
|
| 547 |
selected_missing <- df_missing_tiles_day_subset[i,]==1 |
|
| 548 |
names(df_missing_tiles_day)[selected_missing] |
|
| 555 |
missing_tiles_subset <- subset(missing_tiles,select=selected_col) |
|
| 556 |
selected_missing <- missing_tiles_subset==1 |
|
| 557 |
#names(df_missing_tiles_day)[selected_missing] |
|
| 549 | 558 |
|
| 550 | 559 |
#names(list_tiles_predicted_masked)[selected_missing] |
| 551 | 560 |
list_missing_tiles_raster <- list_tiles_predicted_masked[selected_missing] |
| 552 | 561 |
r_tiles_s <- stack(list_missing_tiles_raster) |
| 553 | 562 |
|
| 554 | 563 |
### first sum missing |
| 555 |
datasum <- stackApply(r_tiles_s, 1:nlayers(r_tiles_s), fun = sum)
|
|
| 564 |
datasum <- stackApply(r_tiles_s, 1:nlayers(r_tiles_s), fun = sum) |
|
| 556 | 565 |
|
| 557 | 566 |
### then substract missing tiles... |
| 558 |
r_day_predicted <- r_overlap_m -datasum |
|
| 567 |
r_day_predicted <- r_overlap_m - datasum
|
|
| 559 | 568 |
|
| 569 |
r_table <- ratify(r_day_predicted) # build the Raster Attibute table |
|
| 570 |
rat <- levels(r_table)[[1]]#get the values of the unique cell frot the attribute table |
|
| 571 |
#rat$legend <- paste0("tile_",1:26)
|
|
| 572 |
tb_freq <- as.data.frame(freq(r_table)) |
|
| 573 |
rat$legend <- tb_freq$value |
|
| 574 |
levels(r_table) <- rat |
|
| 575 |
|
|
| 576 |
res_pix <- 800 |
|
| 577 |
#res_pix <- 480 |
|
| 578 |
col_mfrow <- 1 |
|
| 579 |
row_mfrow <- 1 |
|
| 580 |
|
|
| 581 |
png_filename <- file.path(out_dir,paste("Figure_number_of_predictionds_by_pixel_",date_str,"_",region_name,"_",out_suffix,".png",sep =""))
|
|
| 582 |
|
|
| 583 |
title_str <- paste("Number of predicted pixels for ",variable_name," on ",date_str, sep = "")
|
|
| 584 |
|
|
| 585 |
png(filename=png_filename,width = col_mfrow * res_pix,height = row_mfrow * res_pix) |
|
| 586 |
#my_col=c('blue','red','green')
|
|
| 587 |
my_col <- rainbow(length(tb_freq$value)) |
|
| 588 |
plot(r_table,col=my_col,legend=F,box=F,axes=F,main=title_str) |
|
| 589 |
legend(x='topright', legend =rat$legend,fill = my_col,cex=0.8) |
|
| 590 |
|
|
| 591 |
dev.off() |
|
| 592 |
|
|
| 593 |
### Day missing reclass above |
|
| 594 |
|
|
| 595 |
## do this in gdalcalc? |
|
| 596 |
r_missing_day <- r_day_predicted == 0 |
|
| 597 |
|
|
| 598 |
res_pix <- 800 |
|
| 599 |
#res_pix <- 480 |
|
| 600 |
col_mfrow <- 1 |
|
| 601 |
row_mfrow <- 1 |
|
| 602 |
|
|
| 603 |
png_filename <- file.path(out_dir,paste("Figure_missing_predictionds_by_pixel_",date_str,"_",region_name,"_",out_suffix,".png",sep =""))
|
|
| 604 |
|
|
| 605 |
title_str <- paste("Number of predicted pixels for ",variable_name," on ",date_str, sep = "")
|
|
| 606 |
|
|
| 607 |
png(filename=png_filename,width = col_mfrow * res_pix,height = row_mfrow * res_pix) |
|
| 608 |
#my_col=c('blue','red','green')
|
|
| 609 |
my_col <- c("black","red")
|
|
| 610 |
plot(r_missing_day,col=my_col,legend=F,box=F,axes=F,main=title_str) |
|
| 611 |
legend(x='topright', legend =c("prediced","missing"),fill = my_col,cex=0.8)
|
|
| 612 |
|
|
| 613 |
dev.off() |
|
| 614 |
|
|
| 560 | 615 |
### generate retunr object |
| 561 | 616 |
returm(r_day_predicted) |
| 562 | 617 |
} |
Also available in: Unified diff
adding and debugging function to generate raster of number of predictions for day with missing tiles