Revision b8203902
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/03/2016
|
|
| 12 |
#MODIFIED ON: 11/04/2016
|
|
| 13 | 13 |
#Version: 1 |
| 14 | 14 |
#PROJECT: Environmental Layers project |
| 15 |
#COMMENTS: testing of files by tiles and combining listing
|
|
| 15 |
#COMMENTS: moving functions to function script global product assessment part0
|
|
| 16 | 16 |
#TODO: |
| 17 | 17 |
#1) |
| 18 | 18 |
#First source these files: |
| ... | ... | |
| 331 | 331 |
|
| 332 | 332 |
#collect info: read in all shapefiles |
| 333 | 333 |
|
| 334 |
centroids_shp_fun <- function(i,list_shp_reg_files){
|
|
| 335 |
# |
|
| 336 |
shp_filename <- list_shp_reg_files[[i]] |
|
| 337 |
layer_name <- sub(".shp","",basename(shp_filename))
|
|
| 338 |
path_to_shp <- dirname(shp_filename) |
|
| 339 |
shp1 <- try(readOGR(path_to_shp, layer_name)) #use try to resolve error below |
|
| 340 |
#shp_61.0_-160.0 |
|
| 341 |
#Geographical CRS given to non-conformant data: -186.331747678 |
|
| 342 |
|
|
| 343 |
#shp1<-readOGR(dirname(list_shp_reg_files[[i]]),sub(".shp","",basename(list_shp_reg_files[[i]])))
|
|
| 344 |
if (!inherits(shp1,"try-error")) {
|
|
| 345 |
pt <- gCentroid(shp1) |
|
| 346 |
#centroids_pts[[i]] <- pt |
|
| 347 |
}else{
|
|
| 348 |
pt <- shp1 |
|
| 349 |
#centroids_pts[[i]] <- pt |
|
| 350 |
} |
|
| 351 |
#shps_tiles[[i]] <- shp1 |
|
| 352 |
#centroids_pts[[i]] <- centroids |
|
| 353 |
|
|
| 354 |
shp_obj <- list(shp1,pt) |
|
| 355 |
names(shp_obj) <- c("spdf","centroid")
|
|
| 356 |
return(shp_obj) |
|
| 357 |
} |
|
| 358 | 334 |
|
| 359 | 335 |
obj_centroids_shp <- centroids_shp_fun(1,list_shp_reg_files=in_dir_shp_list) |
| 360 | 336 |
|
| ... | ... | |
| 427 | 403 |
#3. for every pixel generate and ID (tile ID) as integer, there should be 26 layers at the mosaic extent |
| 428 | 404 |
#4. generate a table? for each pixel it can say if is part of a specific tile |
| 429 | 405 |
#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){
|
|
| 431 |
# |
|
| 432 |
tile_spdf <- list_spdf[[i]] |
|
| 433 |
tile_coord <- names(list_spdf)[i] |
|
| 434 |
r_ref <- list_r_ref[[i]] |
|
| 435 |
|
|
| 436 |
df_tmp <- subset(df_missing,date==date_val,select=tile_coord) |
|
| 437 |
#for each row (date) |
|
| 438 |
val <- df_tmp[[tile_coord]] |
|
| 439 |
if(val==1){
|
|
| 440 |
val<-0 #missing then not predicted |
|
| 441 |
}else{
|
|
| 442 |
val<-1 |
|
| 443 |
} |
|
| 444 |
|
|
| 445 |
tile_spdf$predicted <- val |
|
| 446 |
tile_spdf$tile_coord <- tile_coord |
|
| 447 |
tile_spdf$overlap <- 1 |
|
| 448 |
|
|
| 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 |
|
|
| 461 |
return(r) |
|
| 462 |
} |
|
| 463 | 406 |
|
| 464 | 407 |
|
| 465 | 408 |
### use rasterize |
| ... | ... | |
| 469 | 412 |
### Now use intersect to retain actual overlap |
| 470 | 413 |
|
| 471 | 414 |
for(i in 1:length(shps_tiles)){
|
| 472 |
overlap_intersect <- intersect(shps_tiles[[1]],shps_tiles[[i]]))
|
|
| 415 |
overlap_intersect <- intersect(shps_tiles[[1]],shps_tiles[[i]]) |
|
| 473 | 416 |
} |
| 474 | 417 |
|
| 475 | 418 |
overlap_intersect <- lapply(1:length(shps_tiles),FUN=function(i){intersect(shps_tiles[[1]],shps_tiles[[i]])})
|
Also available in: Unified diff
moving functions to function script global product assessment part0