Revision 02113802
Added by Benoit Parmentier over 9 years ago
| climate/research/oregon/interpolation/global_run_scalingup_mosaicing.R | ||
|---|---|---|
| 313 | 313 |
out_dir_str <- out_dir |
| 314 | 314 |
|
| 315 | 315 |
lf_r_weights <- vector("list",length=length(lf_mosaic))
|
| 316 |
|
|
| 316 |
|
|
| 317 | 317 |
############### |
| 318 | 318 |
### PART 2: prepare weights using tile rasters ############ |
| 319 | 319 |
#methods availbable:use_sine_weights,use_edge,use_linear_weights |
| ... | ... | |
| 400 | 400 |
## Rasters tiles vary slightly in resolution, they need to be matched for the mosaic. Resolve issue in the |
| 401 | 401 |
#mosaic funciton using gdal_merge to compute a reference image to mach. |
| 402 | 402 |
|
| 403 |
cmd_str <- paste("python","/usr/bin/gdal_merge.py","-o avg.tif",paste(lf_mosaic,collapse=" "))
|
|
| 403 |
rast_ref <- file.path(out_dir,paste("avg_",out_suffix,file_format,sep="")) #this is a the ref
|
|
| 404 |
|
|
| 405 |
cmd_str <- paste("python","/usr/bin/gdal_merge.py","-o ",rast_ref,paste(lf_mosaic,collapse=" "))
|
|
| 404 | 406 |
system(cmd_str) |
| 405 | 407 |
|
| 406 | 408 |
## Create raster image for original predicted images with matching resolution and extent to the mosaic (reference image) |
| 407 | 409 |
|
| 408 |
rast_ref <- file.path(out_dir,"avg",out_suffix,file_format) #this is a the ref
|
|
| 410 |
#rast_ref <- file.path(out_dir,"avg.tif")
|
|
| 409 | 411 |
r_ref <- raster(rast_ref) |
| 410 | 412 |
#plot(r_ref) |
| 411 | 413 |
|
| 412 |
if(method_str%in%c("use_linear_weights","use_sine_weights","use_edge_weights")){
|
|
| 414 |
if(mosaic_method%in%c("use_linear_weights","use_sine_weights","use_edge_weights")){
|
|
| 413 | 415 |
|
| 414 |
lf_files <- unlist(list_r_weights)
|
|
| 416 |
lf_files <- unlist(list_weights) |
|
| 415 | 417 |
|
| 416 | 418 |
list_param_raster_match <- list(lf_files,rast_ref,file_format,out_suffix,out_dir) |
| 417 | 419 |
names(list_param_raster_match) <- c("lf_files","rast_ref","file_format","out_suffix","out_dir_str")
|
| ... | ... | |
| 419 | 421 |
#debug(raster_match) |
| 420 | 422 |
#r_test <- raster(raster_match(1,list_param_raster_match)) |
| 421 | 423 |
|
| 422 |
list_linear_weights_m <- mclapply(1:length(lf_files),FUN=raster_match,list_param=list_param_raster_match,mc.preschedule=FALSE,mc.cores = num_cores)
|
|
| 424 |
list_weights_m <- mclapply(1:length(lf_files),FUN=raster_match,list_param=list_param_raster_match,mc.preschedule=FALSE,mc.cores = num_cores) |
|
| 423 | 425 |
|
| 424 |
lf_files <- unlist(list_r_weights_prod)
|
|
| 426 |
lf_files <- unlist(list_weights_prod) |
|
| 425 | 427 |
list_param_raster_match <- list(lf_files,rast_ref,file_format,out_suffix,out_dir) |
| 426 | 428 |
names(list_param_raster_match) <- c("lf_files","rast_ref","file_format","out_suffix","out_dir_str")
|
| 427 | 429 |
|
| 428 | 430 |
#num_cores <-11 |
| 429 |
list_linear_weights_prod_m <- mclapply(1:length(lf_files),FUN=raster_match,list_param=list_param_raster_match,mc.preschedule=FALSE,mc.cores = num_cores) |
|
| 430 |
lf_files <- unlist(list_linear_r_weights) |
|
| 431 |
|
|
| 432 |
list_param_raster_match <- list(lf_files,rast_ref,file_format,out_suffix,out_dir) |
|
| 433 |
names(list_param_raster_match) <- c("lf_files","rast_ref","file_format","out_suffix","out_dir_str")
|
|
| 434 |
|
|
| 435 |
#debug(raster_match) |
|
| 436 |
#r_test <- raster(raster_match(1,list_param_raster_match)) |
|
| 437 |
|
|
| 438 |
list_linear_weights_m <- mclapply(1:length(lf_files),FUN=raster_match,list_param=list_param_raster_match,mc.preschedule=FALSE,mc.cores = num_cores) |
|
| 439 |
|
|
| 440 |
lf_files <- unlist(list_linear_r_weights_prod) |
|
| 441 |
list_param_raster_match <- list(lf_files,rast_ref,file_format,out_suffix,out_dir) |
|
| 442 |
names(list_param_raster_match) <- c("lf_files","rast_ref","file_format","out_suffix","out_dir_str")
|
|
| 443 |
|
|
| 444 |
num_cores <-11 |
|
| 445 |
list_linear_weights_prod_m <- mclapply(1:length(lf_files),FUN=raster_match,list_param=list_param_raster_match,mc.preschedule=FALSE,mc.cores = num_cores) |
|
| 446 |
|
|
| 431 |
list_weights_prod_m <- mclapply(1:length(lf_files),FUN=raster_match,list_param=list_param_raster_match,mc.preschedule=FALSE,mc.cores = num_cores) |
|
| 447 | 432 |
|
| 448 | 433 |
##################### |
| 449 | 434 |
###### PART 4: compute the weighted mean with the mosaic function ##### |
| ... | ... | |
| 456 | 441 |
|
| 457 | 442 |
list_args_weights <- list_weights_m |
| 458 | 443 |
list_args_weights_prod <- list_weights_prod_m |
| 459 |
method_str <- list_methods
|
|
| 444 |
method_str <- method
|
|
| 460 | 445 |
|
| 461 | 446 |
#list_args_weights <- (mixedsort(list.files(pattern="r_weights_m_.*.tif"))) |
| 462 | 447 |
list_args_weights <- lapply(1:length(list_args_weights), FUN=function(i,x){raster(x[[i]])},x=list_args_weights)
|
| ... | ... | |
| 482 | 467 |
|
| 483 | 468 |
} |
| 484 | 469 |
|
| 485 |
if(method_mosaic=="unweighted"){
|
|
| 470 |
if(mosaic_method=="unweighted"){
|
|
| 486 | 471 |
#### Fourth use original images |
| 487 | 472 |
#macth file to mosaic extent using the original predictions |
| 488 | 473 |
lf_files <- lf_mosaic |
| ... | ... | |
| 511 | 496 |
raster_name <- file.path(out_dir,paste("r_m_mean_",out_suffix,".tif",sep=""))
|
| 512 | 497 |
writeRaster(r_m_mean, NAflag=NA_flag_val,filename=raster_name,overwrite=TRUE) #unweighted mean |
| 513 | 498 |
|
| 514 |
r_m_mean_unweighted <- paste("r_m_mean_",out_suffix,".tif",sep="")
|
|
| 499 |
#r_m_mean_unweighted <- paste("r_m_mean_",out_suffix,".tif",sep="")
|
|
| 500 |
list_weights <- NULL |
|
| 501 |
list_weights_prod <- NULL |
|
| 515 | 502 |
|
| 516 | 503 |
} |
| 517 | 504 |
|
| 518 | 505 |
#Create return object |
| 519 |
mosaic_obj <- list(raster_name,list_r_weights,list_r_weights_prod,method)
|
|
| 506 |
mosaic_obj <- list(raster_name,list_weights,list_weights_prod,mosaic_method)
|
|
| 520 | 507 |
names(mosaic_obj) <- c("mean_mosaic","r_weights","r_weigths_prod","method")
|
| 521 | 508 |
return(mosaic_obj) |
| 522 | 509 |
} |
| ... | ... | |
| 535 | 522 |
|
| 536 | 523 |
y_var_name <- "dailyTmax" #PARAM1 |
| 537 | 524 |
interpolation_method <- c("gam_CAI") #PARAM2
|
| 538 |
region_name <- "reg4" #PARAM 13 #reg4 is South America, Africa Region 5
|
|
| 525 |
region_name <- "reg2" #PARAM 13 #reg4 South America, Africa reg5,Europe reg2, North America reg1, Asia reg3
|
|
| 539 | 526 |
|
| 540 | 527 |
out_suffix <- paste(region_name,"_","mosaic_run10_1500x4500_global_analyses_06152015",sep="") |
| 541 | 528 |
#PARAM3 |
| ... | ... | |
| 556 | 543 |
file_format <- ".tif" #PARAM 9 |
| 557 | 544 |
NA_value <- -9999 #PARAM10 |
| 558 | 545 |
NA_flag_val <- NA_value |
| 559 |
|
|
| 546 |
|
|
| 547 |
num_cores <- 11 |
|
| 560 | 548 |
tile_size <- "1500x4500" #PARAM 11 |
| 561 | 549 |
mulitple_region <- TRUE #PARAM 12 |
| 562 | 550 |
|
| ... | ... | |
| 567 | 555 |
|
| 568 | 556 |
####### PART 1: Read in data and process data ######## |
| 569 | 557 |
|
| 558 |
#make this a loop?, fistt use sept 1, 2010 data |
|
| 559 |
#out_suffix <- paste(day_to_mosaic[2],out_suffix,sep="_") |
|
| 560 |
|
|
| 570 | 561 |
in_dir <- file.path(in_dir,region_name) |
| 571 | 562 |
out_dir <- in_dir |
| 572 | 563 |
if(create_out_dir_param==TRUE){
|
| ... | ... | |
| 578 | 569 |
|
| 579 | 570 |
setwd(out_dir) |
| 580 | 571 |
|
| 581 |
lf_mosaic <-list.files(path=file.path(in_dir), |
|
| 572 |
lf_mosaic1 <-list.files(path=file.path(in_dir), |
|
| 573 |
pattern=paste(".*.",day_to_mosaic[1],".*.tif$",sep=""),full.names=T) #choosing date 2...20100901
|
|
| 574 |
|
|
| 575 |
lf_mosaic2 <-list.files(path=file.path(in_dir), |
|
| 582 | 576 |
pattern=paste(".*.",day_to_mosaic[2],".*.tif$",sep=""),full.names=T) #choosing date 2...20100901
|
| 583 | 577 |
#lf_mosaic <- lf_mosaic[1:20] |
| 584 |
r1 <- raster(lf_mosaic[1]) |
|
| 585 |
r2 <- raster(lf_mosaic[2]) |
|
| 578 |
r1 <- raster(lf_mosaic1[1])
|
|
| 579 |
r2 <- raster(lf_mosaic2[2])
|
|
| 586 | 580 |
|
| 587 | 581 |
plot(r1) |
| 588 | 582 |
plot(r2) |
| 589 | 583 |
|
| 590 |
lf <- sub(".tif","",lf_mosaic)
|
|
| 584 |
lf <- sub(".tif","",lf_mosaic2)
|
|
| 591 | 585 |
tx<-strsplit(as.character(lf),"_") |
| 592 | 586 |
|
| 593 | 587 |
lat<- as.character(lapply(1:length(tx),function(i,x){x[[i]][13]},x=tx))
|
| ... | ... | |
| 602 | 596 |
proj4string(df_centroids) <- projection(r1) |
| 603 | 597 |
df_points <- df_centroids |
| 604 | 598 |
#methods availbable:use_sine_weights,use_edge,use_linear_weights |
| 599 |
out_suffix_str <- paste(day_to_mosaic[2],out_suffix,sep="_") |
|
| 600 |
|
|
| 601 |
#debug(mosaicFiles) |
|
| 602 |
mosaic_edge_20100901_obj <- mosaicFiles(lf_mosaic2,mosaic_method="use_edge_weights", |
|
| 603 |
num_cores=num_cores, |
|
| 604 |
python_bin=NULL, |
|
| 605 |
df_points=NULL,NA_flag=NA_flag_val, |
|
| 606 |
file_format=file_format,out_suffix=out_suffix_str, |
|
| 607 |
out_dir=out_dir) |
|
| 608 |
#debug(mosaicFiles) |
|
| 609 |
|
|
| 610 |
mosaic_unweighted_20100901_obj <- mosaicFiles(lf_mosaic2,mosaic_method="unweighted", |
|
| 611 |
num_cores=num_cores, |
|
| 612 |
python_bin=NULL, |
|
| 613 |
df_points=NULL,NA_flag=NA_flag_val, |
|
| 614 |
file_format=file_format,out_suffix=out_suffix_str, |
|
| 615 |
out_dir=out_dir) |
|
| 616 |
|
|
| 617 |
r2_unweighted <-raster(mosaic_unweighted_20100901_obj$mean_mosaic) |
|
| 618 |
r2_edge <-raster(mosaic_edge_20100901_obj$mean_mosaic) |
|
| 619 |
|
|
| 620 |
#plot(r2_edge) |
|
| 621 |
|
|
| 622 |
mosaic_edge_20100831_obj <- mosaicFiles(lf_mosaic1,mosaic_method="use_edge_weights", |
|
| 623 |
num_cores=num_cores, |
|
| 624 |
python_bin=NULL, |
|
| 625 |
df_points=NULL,NA_flag=NA_flag_val, |
|
| 626 |
file_format=file_format,out_suffix=out_suffix_str, |
|
| 627 |
out_dir=out_dir) |
|
| 628 |
|
|
| 629 |
mosaic_edge_20100831_obj <- mosaicFiles(lf_mosaic1,mosaic_method="unweighted", |
|
| 630 |
num_cores=num_cores, |
|
| 631 |
python_bin=NULL, |
|
| 632 |
df_points=NULL,NA_flag=NA_flag_val, |
|
| 633 |
file_format=file_format,out_suffix=out_suffix_str, |
|
| 634 |
out_dir=out_dir) |
|
| 605 | 635 |
|
| 606 |
debug(mosaicFiles) |
|
| 607 |
mosaic_edge_20100901_obj <- mosaicFiles(lf_mosaic,mosaic_method="use_edge_weights",python_bin=NULL,NA_flag,file_format,num_cores,out_suffix=out_suffix) |
|
| 608 | 636 |
|
| 609 | 637 |
## |
| 610 | 638 |
|
| ... | ... | |
| 616 | 644 |
#m_val= sum_val/weight_sum #mean value |
| 617 | 645 |
# |
| 618 | 646 |
|
| 619 |
#### Now run unweighted |
|
| 620 |
|
|
| 621 |
#.... |
|
| 622 | 647 |
|
| 623 | 648 |
##################### |
| 624 |
###### PART 5: Now plot of the weighted mean and unweighted mean with the mosaic function #####
|
|
| 649 |
###### PART 2: Analysis and figures for the outputs of mosaic function #####
|
|
| 625 | 650 |
|
| 626 | 651 |
#### compute and aspect and slope with figures |
| 627 | 652 |
|
Also available in: Unified diff
debugging and testing mosaicFiles function