Revision f3ed9122
Added by Benoit Parmentier about 8 years ago
| climate/research/oregon/interpolation/global_product_assessment_part0.R | ||
|---|---|---|
| 12 | 12 |
#MODIFIED ON: 11/04/2016 |
| 13 | 13 |
#Version: 1 |
| 14 | 14 |
#PROJECT: Environmental Layers project |
| 15 |
#COMMENTS: moving functions to function script global product assessment part0
|
|
| 15 |
#COMMENTS: |
|
| 16 | 16 |
#TODO: |
| 17 | 17 |
#1) |
| 18 | 18 |
#First source these files: |
| ... | ... | |
| 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: experimenting with detection of gaps, testing rasterize
|
|
| 23 |
#COMMIT: computing maximum overlap for each pixel in a region
|
|
| 24 | 24 |
|
| 25 | 25 |
################################################################################################# |
| 26 | 26 |
|
| ... | ... | |
| 54 | 54 |
library(zoo) |
| 55 | 55 |
library(xts) |
| 56 | 56 |
library(lubridate) |
| 57 |
library(mosaic)
|
|
| 57 |
#library(mosaic) #not installed on NEX
|
|
| 58 | 58 |
|
| 59 | 59 |
###### Function used in the script ####### |
| 60 | 60 |
|
| ... | ... | |
| 86 | 86 |
source(file.path(script_path,function_assessment_part3)) #source all functions used in this script |
| 87 | 87 |
|
| 88 | 88 |
#Product assessment |
| 89 |
function_product_assessment_part0_functions <- "global_product_assessment_part0_functions_10312016.R"
|
|
| 89 |
function_product_assessment_part0_functions <- "global_product_assessment_part0_functions_11042016.R"
|
|
| 90 | 90 |
source(file.path(script_path,function_product_assessment_part0_functions)) #source all functions used in this script |
| 91 | 91 |
##Don't load part 1 and part2, mosaic package does not work on NEX |
| 92 | 92 |
#function_product_assessment_part1_functions <- "global_product_assessment_part1_functions_09192016b.R" |
| ... | ... | |
| 331 | 331 |
|
| 332 | 332 |
#collect info: read in all shapefiles |
| 333 | 333 |
|
| 334 |
|
|
| 335 | 334 |
obj_centroids_shp <- centroids_shp_fun(1,list_shp_reg_files=in_dir_shp_list) |
| 336 | 335 |
|
| 337 |
|
|
| 338 | 336 |
obj_centroids_shp <- mclapply(1:length(in_dir_shp_list), |
| 339 | 337 |
FUN=centroids_shp_fun, |
| 340 | 338 |
list_shp_reg_files=in_dir_shp_list, |
| ... | ... | |
| 372 | 370 |
|
| 373 | 371 |
names(shps_tiles) <- basename(unlist(in_dir_reg)) |
| 374 | 372 |
r_ref <- raster(list_lf_raster_tif_tiles[[1]][1]) |
| 375 |
list_r_ref <- lapply(1:length(in_dir_reg), function(x){raster(list_lf_raster_tif_tiles[[i]][1])})
|
|
| 373 |
list_r_ref <- lapply(1:length(in_dir_reg), function(i){raster(list_lf_raster_tif_tiles[[i]][1])})
|
|
| 376 | 374 |
tile_spdf <- shps_tiles[[1]] |
| 377 | 375 |
tile_coord <- basename(in_dir_reg[1]) |
| 378 | 376 |
date_val <- df_missing$date[1] |
| 379 | 377 |
|
| 380 |
debug(rasterize_tile_day) |
|
| 378 |
function_product_assessment_part0_functions <- "global_product_assessment_part0_functions_11042016.R" |
|
| 379 |
source(file.path(script_path,function_product_assessment_part0_functions)) #source all functions used in this script |
|
| 380 |
|
|
| 381 |
undebug(rasterize_tile_day) |
|
| 381 | 382 |
list_predicted <- rasterize_tile_day(1, |
| 382 | 383 |
list_spdf=shps_tiles, |
| 383 | 384 |
df_missing=df_missing, |
| 384 | 385 |
list_r_ref=list_r_ref, |
| 385 | 386 |
col_name="overlap", |
| 386 | 387 |
date_val=df_missing$date[1]) |
| 387 |
list_predicted <- mclapply(1:length(shps_tiles), |
|
| 388 |
##check that everything is correct: |
|
| 389 |
plot(r_mask) |
|
| 390 |
plot(list_predicted,add=T) |
|
| 391 |
plot(spdf_tiles_test,add=T,border="green",usePolypath = FALSE) #added usePolypath following error on brige and NEX |
|
| 392 |
|
|
| 393 |
list_predicted <- mclapply(1:6, |
|
| 388 | 394 |
FUN=rasterize_tile_day, |
| 389 | 395 |
list_spdf=shps_tiles, |
| 390 | 396 |
df_missing=df_missing, |
| 397 |
list_r_ref=list_r_ref, |
|
| 391 | 398 |
col_name = "overlap", |
| 399 |
date_val=df_missing$date[1], |
|
| 400 |
mc.preschedule=FALSE, |
|
| 401 |
mc.cores = num_cores) |
|
| 402 |
|
|
| 403 |
list_predicted <- mclapply(1:length(shps_tiles), |
|
| 404 |
FUN=rasterize_tile_day, |
|
| 405 |
list_spdf=shps_tiles, |
|
| 406 |
df_missing=df_missing, |
|
| 392 | 407 |
list_r_ref=list_r_ref, |
| 408 |
col_name = "overlap", |
|
| 393 | 409 |
date_val=df_missing$date[1], |
| 394 | 410 |
mc.preschedule=FALSE, |
| 395 | 411 |
mc.cores = num_cores) |
| 396 | 412 |
|
| 413 |
### Make a list of file |
|
| 414 |
filename_list_mosaics_weights_m <- file.path(out_dir_str,paste("list_to_mosaics_","weights_",mosaic_method,"_",out_suffix_str_tmp,".txt",sep=""))
|
|
| 415 |
|
|
| 416 |
#writeLines(unlist(list_weights_m),con=filename_list_mosaics_weights_m) #weights files to mosaic |
|
| 417 |
#writeLines(unlist(list_weights_prod_m),con=filename_list_mosaics_prod_weights_m) #prod weights files to mosaic |
|
| 418 |
|
|
| 419 |
writeLines(unlist(list_weights_m),con=filename_list_mosaics_weights_m) #weights files to mosaic |
|
| 420 |
|
|
| 421 |
#out_mosaic_name_weights_m <- r_weights_sum_raster_name <- file.path(out_dir,paste("r_weights_sum_m_",mosaic_method,"_weighted_mean_",out_suffix,".tif",sep=""))
|
|
| 422 |
#out_mosaic_name_prod_weights_m <- r_weights_sum_raster_name <- file.path(out_dir,paste("r_prod_weights_sum_m_",mosaic_method,"_weighted_mean_",out_suffix,".tif",sep=""))
|
|
| 423 |
out_mosaic_name_weights_m <- file.path(out_dir_str,paste("r_weights_sum_m_",mosaic_method,"_weighted_mean_",out_suffix_str_tmp,".tif",sep=""))
|
|
| 424 |
|
|
| 425 |
rast_ref_name <- infile_mask |
|
| 426 |
mostaic_python <- "/nobackupp6/aguzman4/climateLayers/sharedCode/" |
|
| 427 |
|
|
| 428 |
#python /nobackupp6/aguzman4/climateLayers/sharedCode//gdal_merge_sum.py --config GDAL_CACHEMAX=1500 --overwrite=TRUE -o /nobackupp8/bparmen1/climateLayers/out |
|
| 429 |
mosaic_overlap_tiles_obj <- mosaic_python_merge(NA_flag_val=NA_flag_val, |
|
| 430 |
module_path=mosaic_python, |
|
| 431 |
module_name="gdal_merge_sum.py", |
|
| 432 |
input_file=filename_list_mosaics_weights_m, |
|
| 433 |
out_mosaic_name=out_mosaic_name_weights_m, |
|
| 434 |
raster_ref_name = rast_ref_name) ##if NA, not take into account |
|
| 435 |
r_weights_sum_raster_name <- mosaic_weights_obj$out_mosaic_name |
|
| 436 |
cmd_str1 <- mosaic_weights_obj$cmd_str |
|
| 437 |
|
|
| 397 | 438 |
#http://stackoverflow.com/questions/19586945/how-to-legend-a-raster-using-directly-the-raster-attribute-table-and-displaying |
| 398 | 439 |
# |
| 399 | 440 |
#http://gis.stackexchange.com/questions/148398/how-does-spatial-polygon-over-polygon-work-to-when-aggregating-values-in-r |
| ... | ... | |
| 403 | 444 |
#3. for every pixel generate and ID (tile ID) as integer, there should be 26 layers at the mosaic extent |
| 404 | 445 |
#4. generate a table? for each pixel it can say if is part of a specific tile |
| 405 | 446 |
#5. workout a formula to generate the number of predictions for each pixel based on tile predicted for each date!! |
| 406 |
|
|
| 407 |
|
|
| 447 |
|
|
| 408 | 448 |
### use rasterize |
| 409 | 449 |
spdf_tiles <- do.call(bind, shps_tiles) #bind all tiles together in one shapefile |
| 410 | 450 |
spdf_tiles <- do.call(intersect, shps_tiles) |
Also available in: Unified diff
computing maximum overlap for each pixel in a region