Revision 2907e644
Added by Benoit Parmentier about 9 years ago
| climate/research/oregon/interpolation/global_run_scalingup_mosaicing_function.R | ||
|---|---|---|
| 4 | 4 |
#Different options to explore mosaicing are tested. This script only contains functions. |
| 5 | 5 |
#AUTHOR: Benoit Parmentier |
| 6 | 6 |
#CREATED ON: 04/14/2015 |
| 7 |
#MODIFIED ON: 10/21/2015
|
|
| 7 |
#MODIFIED ON: 10/26/2015
|
|
| 8 | 8 |
#Version: 1 |
| 9 | 9 |
#PROJECT: Environmental Layers project |
| 10 | 10 |
#COMMENTS: first commit of function script to test mosaicing using 1500x4500km and other tiles |
| ... | ... | |
| 113 | 113 |
# |
| 114 | 114 |
############ |
| 115 | 115 |
|
| 116 |
lf <- list_param$lf #list of files to mosaic |
|
| 116 |
lf <- list_param$lf[[i]] #list of files to mosaic
|
|
| 117 | 117 |
tb <- list_param$tb #fitting or validation table with all days |
| 118 | 118 |
metric_name <- list_param$metric_name #RMSE, MAE etc. |
| 119 | 119 |
pred_mod_name <- list_param$pred_mod_name #mod1, mod_kr etc. |
| ... | ... | |
| 184 | 184 |
list_raster_name[[j]] <- outFilename |
| 185 | 185 |
} |
| 186 | 186 |
|
| 187 |
|
|
| 188 | 187 |
raster_created_obj <- list(list_raster_name,df_centroids) |
| 189 | 188 |
names(raster_created_obj) <- c("list_raster_name","df_centroids")
|
| 190 | 189 |
return(raster_created_obj) |
| ... | ... | |
| 538 | 537 |
|
| 539 | 538 |
## Rasters tiles vary slightly in resolution, they need to be matched for the mosaic. Resolve issue in the |
| 540 | 539 |
#mosaic function using gdal_merge to compute a reference image to mach. |
| 541 |
|
|
| 542 |
rast_ref <- file.path(out_dir,paste("avg_",out_suffix,file_format,sep="")) #this is a the ref
|
|
| 543 |
if(is.null(python_bin)){
|
|
| 540 |
#This step of creating a merged raster can be avoided if a reference maks image is given |
|
| 541 |
#this needs to be changed to avoid further bugs!!! |
|
| 542 |
|
|
| 543 |
if(!is.null(r_mask_raster_name)){
|
|
| 544 |
rast_ref <- r_mask_raster_name #the mask file is used as ref. |
|
| 545 |
#mask(raster(r_m_weighted_mean_raster_name),mask=r_mask_raster_name,filename=r_m_weighted_mean_mask_raster_name) |
|
| 546 |
#raster_name <- r_m_weighted_mean_mask_raster_name |
|
| 547 |
}else{
|
|
| 548 |
rast_ref <- file.path(out_dir,paste("avg_",out_suffix,file_format,sep="")) #this is a the ref
|
|
| 549 |
if(is.null(python_bin)){
|
|
| 544 | 550 |
python_bin="" |
| 545 |
} |
|
| 551 |
}
|
|
| 546 | 552 |
|
| 547 |
python_cmd <- file.path(python_bin,"gdal_merge.py") |
|
| 548 |
cmd_str <- paste("python",python_cmd,"-o ",rast_ref,paste(lf_mosaic,collapse=" "))
|
|
| 549 |
system(cmd_str) |
|
| 550 |
|
|
| 553 |
python_cmd <- file.path(python_bin,"gdal_merge.py") |
|
| 554 |
cmd_str <- paste("python",python_cmd,"-o ",rast_ref,paste(lf_mosaic,collapse=" "))
|
|
| 555 |
system(cmd_str) |
|
| 556 |
} |
|
| 557 |
|
|
| 551 | 558 |
## Create raster image for original predicted images with matching resolution and extent to the mosaic (reference image) |
| 552 | 559 |
|
| 553 | 560 |
#rast_ref <- file.path(out_dir,"avg.tif") |
| ... | ... | |
| 558 | 565 |
|
| 559 | 566 |
lf_files <- unlist(list_weights) |
| 560 | 567 |
|
| 568 |
|
|
| 569 |
##Maching resolution is probably only necessary for the r mosaic function |
|
| 570 |
#MOdify later to take into account option R or python... |
|
| 561 | 571 |
list_param_raster_match <- list(lf_files,rast_ref,file_format,python_bin,out_suffix,out_dir) |
| 562 | 572 |
names(list_param_raster_match) <- c("lf_files","rast_ref","file_format","python_bin","out_suffix","out_dir_str")
|
| 563 | 573 |
|
| ... | ... | |
| 592 | 602 |
#get the list of weights product into raster objects |
| 593 | 603 |
|
| 594 | 604 |
list_args_weights_prod <- lapply(1:length(list_args_weights_prod), FUN=function(i,x){raster(x[[i]])},x=list_args_weights_prod)
|
| 595 |
list_args_weights_prod$fun <- "sum" |
|
| 596 |
list_args_weights_prod$na.rm <- TRUE |
|
| 605 |
list_args_weights_prod$fun <- "sum" #use sum while mosaicing
|
|
| 606 |
list_args_weights_prod$na.rm <- TRUE #deal with NA by removal
|
|
| 597 | 607 |
r_weights_sum_raster_name <- file.path(out_dir,paste("r_weights_sum_m_",method_str,"_weighted_mean_",out_suffix,".tif",sep=""))
|
| 598 | 608 |
list_args_weights$filename <- r_weights_sum_raster_name |
| 599 | 609 |
list_args_weights$overwrite<- TRUE |
| ... | ... | |
| 608 | 618 |
#sum in addition ot the current option for mean. |
| 609 | 619 |
#This took 23 minutes! |
| 610 | 620 |
r_weights_sum <- do.call(mosaic,list_args_weights) #weights sum image mosaiced |
| 621 |
#This took 23 minutes with the R function |
|
| 611 | 622 |
r_prod_sum <- do.call(mosaic,list_args_weights_prod) #weights sum product image mosacied |
| 612 | 623 |
|
| 613 | 624 |
#r_m_weighted_mean <- r_prod_sum/r_weights_sum #this is the mosaic using weighted mean... |
| 614 |
#"gdal_calc.py -A input1.tif -B input2.tif --outfile=result.tif --calc="A+B"" |
|
| 615 |
|
|
| 625 |
|
|
| 616 | 626 |
r_m_weighted_mean_raster_name <- file.path(out_dir,paste("r_m_",method_str,"_weighted_mean_",out_suffix,".tif",sep=""))
|
| 617 | 627 |
|
| 618 |
#cmd_str <- paste("gdal_calc.py"
|
|
| 619 |
# " "-A input1.tif", |
|
| 620 |
# "-B input2.tif", |
|
| 621 |
# paste("--outfile=",r_m_weighted_mean_name,sep=""),
|
|
| 622 |
# --calc="A+B" |
|
| 623 |
#--NoDataValue=0) |
|
| 624 | 628 |
if(is.null(python_bin)){
|
| 625 | 629 |
python_bin="" |
| 626 | 630 |
} |
| ... | ... | |
| 630 | 634 |
paste("-A ", r_prod_sum_raster_name,sep=""),
|
| 631 | 635 |
paste("-B ", r_weights_sum_raster_name,sep=""),
|
| 632 | 636 |
paste("--outfile=",r_m_weighted_mean_raster_name,sep=""),
|
| 633 |
"--calc='A/B'","--overwrite",sep=" ") |
|
| 637 |
"--calc='A/B'","--overwrite",sep=" ") #division by zero is problematic...
|
|
| 634 | 638 |
system(cmd_str) |
| 635 | 639 |
|
| 636 | 640 |
#writeRaster(r_m_weighted_mean, NAflag=NA_flag_val,filename=raster_name,overwrite=TRUE) |
| 637 | 641 |
#now use the mask |
| 638 | 642 |
if(!is.null(r_mask_raster_name)){
|
| 639 | 643 |
r_m_weighted_mean_mask_raster_name <- file.path(out_dir,paste("r_m_",method_str,"_weighted_mean_mask_",out_suffix,".tif",sep=""))
|
| 640 |
mask(raster(r_m_weighted_mean_raster_name),mask=r_mask_raster_name,filename=r_m_weighted_mean_mask_raster_name) |
|
| 644 |
mask(raster(r_m_weighted_mean_raster_name),mask=raster(r_mask_raster_name), |
|
| 645 |
filename=r_m_weighted_mean_mask_raster_name,overwrite=TRUE) |
|
| 641 | 646 |
raster_name <- r_m_weighted_mean_mask_raster_name |
| 642 | 647 |
}else{
|
| 643 | 648 |
raster_name <- r_m_weighted_mean_raster_name |
| ... | ... | |
| 682 | 687 |
|
| 683 | 688 |
if(!is.null(r_mask_raster_name)){
|
| 684 | 689 |
r_m_mean_mask_raster_name <- file.path(out_dir,paste("r_m_",method_str,"_unweighted_mean_mask_",out_suffix,".tif",sep=""))
|
| 685 |
mask(raster(r_m_mean_raster_name),mask=r_mask_raster_name,filename=r_m_mean_mask_raster_name) |
|
| 690 |
mask(raster(r_m_mean_raster_name),mask=raster(r_mask_raster_name), |
|
| 691 |
filename=r_m_mean_mask_raster_name,overwrite=TRUE) |
|
| 686 | 692 |
raster_name <- r_m_mean_mask_raster_name |
| 687 | 693 |
}else{
|
| 688 | 694 |
raster_name <- r_m_mean_raster_name |
Also available in: Unified diff
adding mask option debugging