Revision 20935590
Added by Benoit Parmentier about 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/01/2016
|
|
| 12 |
#MODIFIED ON: 11/03/2016
|
|
| 13 | 13 |
#Version: 1 |
| 14 | 14 |
#PROJECT: Environmental Layers project |
| 15 | 15 |
#COMMENTS: testing of files by tiles and combining listing |
| ... | ... | |
| 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: combining tif tiles and shapefiles to examine potential gaps
|
|
| 23 |
#COMMIT: experimenting with detection of gaps, testing rasterize
|
|
| 24 | 24 |
|
| 25 | 25 |
################################################################################################# |
| 26 | 26 |
|
| ... | ... | |
| 378 | 378 |
centroids_pts <- remove_errors_list(centroids_pts) #[[!inherits(shps_tiles,"try-error")]] |
| 379 | 379 |
#centroids_pts <- tmp_pts |
| 380 | 380 |
|
| 381 |
r <- raster(infile_mask) |
|
| 381 |
r_mask <- raster(infile_mask)
|
|
| 382 | 382 |
plot(r) |
| 383 |
plot(shps_tiles[[26]],add=T,border="blue",usePolypath = FALSE) #added usePolypath following error on brige and NEX
|
|
| 383 |
plot(shps_tiles[[1]],add=T,border="blue",usePolypath = FALSE) #added usePolypath following error on brige and NEX
|
|
| 384 | 384 |
|
| 385 | 385 |
## find overlap |
| 386 | 386 |
#http://gis.stackexchange.com/questions/156660/loop-to-check-multiple-polygons-overlap-in-r |
| ... | ... | |
| 394 | 394 |
# |
| 395 | 395 |
} |
| 396 | 396 |
|
| 397 |
names(shps_tiles) <- list_names_tile_coord
|
|
| 397 |
names(shps_tiles) <- basename(unlist(in_dir_reg))
|
|
| 398 | 398 |
r_ref <- raster(list_lf_raster_tif_tiles[[1]][1]) |
| 399 |
list_r_ref <- lapply(1:length(in_dir_reg), function(x){raster(list_lf_raster_tif_tiles[[i]][1])})
|
|
| 399 | 400 |
tile_spdf <- shps_tiles[[1]] |
| 400 | 401 |
tile_coord <- basename(in_dir_reg[1]) |
| 401 | 402 |
date_val <- df_missing$date[1] |
| 402 | 403 |
|
| 403 |
mclapply(1:length(shps_tiles), |
|
| 404 |
debug(rasterize_tile_day) |
|
| 405 |
list_predicted <- rasterize_tile_day(1, |
|
| 406 |
list_spdf=shps_tiles, |
|
| 407 |
df_missing=df_missing, |
|
| 408 |
list_r_ref=list_r_ref, |
|
| 409 |
col_name="overlap", |
|
| 410 |
date_val=df_missing$date[1]) |
|
| 411 |
list_predicted <- mclapply(1:length(shps_tiles), |
|
| 404 | 412 |
FUN=rasterize_tile_day, |
| 405 | 413 |
list_spdf=shps_tiles, |
| 406 | 414 |
df_missing=df_missing, |
| 407 |
r_ref=list_r |
|
| 408 |
) |
|
| 409 |
rasterize_tile_day <- function(i,list_spdf,df_missing,list_r_ref,date_val){
|
|
| 415 |
col_name = "overlap", |
|
| 416 |
list_r_ref=list_r_ref, |
|
| 417 |
date_val=df_missing$date[1], |
|
| 418 |
mc.preschedule=FALSE, |
|
| 419 |
mc.cores = num_cores) |
|
| 420 |
|
|
| 421 |
#http://stackoverflow.com/questions/19586945/how-to-legend-a-raster-using-directly-the-raster-attribute-table-and-displaying |
|
| 422 |
# |
|
| 423 |
#http://gis.stackexchange.com/questions/148398/how-does-spatial-polygon-over-polygon-work-to-when-aggregating-values-in-r |
|
| 424 |
#ok other way of doing this: |
|
| 425 |
#1. find overlap assuming all predictions! |
|
| 426 |
#2. Us raster image with number of overlaps in the mosaic tile |
|
| 427 |
#3. for every pixel generate and ID (tile ID) as integer, there should be 26 layers at the mosaic extent |
|
| 428 |
#4. generate a table? for each pixel it can say if is part of a specific tile |
|
| 429 |
#5. workout a formula to generate the number of predictions for each pixel based on tile predicted for each date!! |
|
| 430 |
rasterize_tile_day <- function(i,list_spdf,df_missing,list_r_ref,col_name,date_val){
|
|
| 410 | 431 |
# |
| 411 | 432 |
tile_spdf <- list_spdf[[i]] |
| 412 |
tile_coord <- names(tile_spdf)[i]
|
|
| 433 |
tile_coord <- names(list_spdf)[i]
|
|
| 413 | 434 |
r_ref <- list_r_ref[[i]] |
| 414 | 435 |
|
| 415 | 436 |
df_tmp <- subset(df_missing,date==date_val,select=tile_coord) |
| ... | ... | |
| 425 | 446 |
tile_spdf$tile_coord <- tile_coord |
| 426 | 447 |
tile_spdf$overlap <- 1 |
| 427 | 448 |
|
| 428 |
r <- rasterize(tile_spdf,r_ref,"predicted") |
|
| 429 |
|
|
| 449 |
#r <- rasterize(tile_spdf,r_ref,"predicted") |
|
| 450 |
#r <- rasterize(tile_spdf,r_ref,col_name) |
|
| 451 |
r <- raster(r_ref,crs=projection(r_ref)) #new layer without values |
|
| 452 |
if(col_name=="overlap"){
|
|
| 453 |
set1f <- function(x){rep(1, x)}
|
|
| 454 |
r <- init(r, fun=set1f, overwrite=TRUE) |
|
| 455 |
} |
|
| 456 |
if(col_name=="predicted"){
|
|
| 457 |
set1f <- function(x){rep(val, x)}
|
|
| 458 |
r <- init(r, fun=set1f, overwrite=TRUE) |
|
| 459 |
} |
|
| 460 |
|
|
| 430 | 461 |
return(r) |
| 431 | 462 |
} |
| 432 | 463 |
|
Also available in: Unified diff
experimenting with detection of gaps, testing rasterize