Revision 2b15a9a3
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: 12/17/2015
|
|
| 7 |
#MODIFIED ON: 12/19/2015
|
|
| 8 | 8 |
#Version: 2 |
| 9 | 9 |
#PROJECT: Environmental Layers project |
| 10 | 10 |
#COMMENTS: first commit of function script to test mosaicing using 1500x4500km and other tiles |
| ... | ... | |
| 12 | 12 |
#1) Make this is a script/function callable from the shell/bash |
| 13 | 13 |
#2) Improve performance: there will be a need to improve efficiency for the workflow. |
| 14 | 14 |
|
| 15 |
#Functions: available: |
|
| 15 |
#Error message for gdal_proximity: |
|
| 16 |
#ERROR 1: Source and proximity bands are not the same size. |
|
| 17 |
#gdal_proximity.py give an error when tries to replace an existent output file with different band. |
|
| 18 |
#If the output already exists and was created by a different band input file, this error is displayed: ERROR 1: Source and proximity bands are not the same size. |
|
| 19 |
#If you remove the existent file, it works fine. |
|
| 20 |
#available: |
|
| 16 | 21 |
#See below |
| 17 | 22 |
|
| 18 | 23 |
################################################################################################# |
| ... | ... | |
| 35 | 40 |
library(reshape) # Change shape of object, summarize results |
| 36 | 41 |
library(plotrix) # Additional plotting functions |
| 37 | 42 |
library(plyr) # Various tools including rbind.fill |
| 38 |
#library(spgwr) # GWR method, not on NEX
|
|
| 43 |
library(spgwr) # GWR method, not on NEX |
|
| 39 | 44 |
library(automap) # Kriging automatic fitting of variogram using gstat, not on NEX |
| 40 | 45 |
library(rgeos) # Geometric, topologic library of functions |
| 41 |
#RPostgreSQL # Interface R and Postgres, not used in this script
|
|
| 46 |
library(RPostgreSQL) # Interface R and Postgres, not used in this script
|
|
| 42 | 47 |
library(gridExtra) |
| 43 | 48 |
#Additional libraries not used in workflow |
| 44 | 49 |
library(pgirmess) # Krusall Wallis test with mulitple options, Kruskalmc {pgirmess}
|
| ... | ... | |
| 318 | 323 |
r_init[1:n_row,1] <- 1 |
| 319 | 324 |
r_init[1:n_row,n_col] <- 1 |
| 320 | 325 |
#r_dist <- distance(r_init) |
| 321 |
srcfile <- file.path(out_dir,paste("feature_target_",tile_no,".tif",sep=""))
|
|
| 326 |
srcfile <- file.path(out_dir,paste("feature_target_",tile_no,"_",Sys.getpid(),".tif",sep=""))
|
|
| 322 | 327 |
|
| 323 | 328 |
writeRaster(r_init,filename=srcfile,overwrite=T) |
| 324 |
dstfile <- file.path(out_dir,paste("feature_target_edge_distance",tile_no,".tif",sep=""))
|
|
| 329 |
#Sys.getpid |
|
| 330 |
dstfile <- file.path(out_dir,paste("feature_target_edge_distance",tile_no,"_",Sys.getpid(),".tif",sep=""))
|
|
| 325 | 331 |
n_values <- "1" |
| 326 | 332 |
|
| 327 | 333 |
cmd_str <- paste("gdal_proximity.py", srcfile, dstfile,"-values",n_values,sep=" ")
|
| ... | ... | |
| 548 | 554 |
names(list_param_create_weights) <- c("lf","df_points","r_feature","method","out_dir_str")
|
| 549 | 555 |
#num_cores <- 11 |
| 550 | 556 |
#debug(create_weights_fun) |
| 551 |
#weights_obj <- create_weights_fun(1,list_param=list_param_create_weights)
|
|
| 557 |
#weights_obj <- create_weights_fun(3,list_param=list_param_create_weights)
|
|
| 552 | 558 |
|
| 553 | 559 |
#This is the function creating the weights by tile. Distance from the centroids needs to be change from distance to |
| 554 | 560 |
#the edges...can use rows and columsn to set edges to 1 and 0 for the others. |
| ... | ... | |
| 796 | 802 |
#writeRaster(r_m_weighted_mean, NAflag=NA_flag_val,filename=raster_name,overwrite=TRUE) |
| 797 | 803 |
#now use the mask |
| 798 | 804 |
if(!is.null(r_mask_raster_name)){
|
| 805 |
#different extent between mask and output if match extent is false!! |
|
| 806 |
#match resolution and extent first |
|
| 807 |
|
|
| 808 |
lf_files <- c(r_m_weighted_mean_raster_name) #file(s) to be matched |
|
| 809 |
rast_ref <- r_mask_raster_name |
|
| 810 |
##Maching resolution is probably only necessary for the r mosaic function |
|
| 811 |
#Modify later to take into account option R or python... |
|
| 812 |
list_param_raster_match <- list(lf_files,rast_ref,file_format,python_bin,out_suffix,out_dir) |
|
| 813 |
names(list_param_raster_match) <- c("lf_files","rast_ref","file_format","python_bin","out_suffix","out_dir_str")
|
|
| 814 |
|
|
| 815 |
#undebug(raster_match) |
|
| 816 |
r_m_weighted_mean_raster_name_matched <- raster_match(1,list_param_raster_match) |
|
| 817 |
|
|
| 799 | 818 |
r_m_weighted_mean_mask_raster_name <- file.path(out_dir,paste("r_m_",mosaic_method,"_weighted_mean_mask_",out_suffix,".tif",sep=""))
|
| 800 |
mask(raster(r_m_weighted_mean_raster_name),mask=raster(r_mask_raster_name), |
|
| 819 |
mask(raster(r_m_weighted_mean_raster_name_matched),mask=raster(r_mask_raster_name),
|
|
| 801 | 820 |
filename=r_m_weighted_mean_mask_raster_name,overwrite=TRUE) |
| 802 | 821 |
raster_name <- r_m_weighted_mean_mask_raster_name |
| 803 | 822 |
}else{
|
| ... | ... | |
| 842 | 861 |
|
| 843 | 862 |
|
| 844 | 863 |
if(!is.null(r_mask_raster_name)){
|
| 864 |
#different extent between mask and output if match extent is false!! |
|
| 865 |
#match resolution and extent first |
|
| 866 |
|
|
| 867 |
lf_files <- c(r_m_mean_raster_name) #file(s) to be matched |
|
| 868 |
rast_ref <- r_mask_raster_name |
|
| 869 |
##Maching resolution is probably only necessary for the r mosaic function |
|
| 870 |
#Modify later to take into account option R or python... |
|
| 871 |
list_param_raster_match <- list(lf_files,rast_ref,file_format,python_bin,out_suffix,out_dir) |
|
| 872 |
names(list_param_raster_match) <- c("lf_files","rast_ref","file_format","python_bin","out_suffix","out_dir_str")
|
|
| 873 |
|
|
| 874 |
#undebug(raster_match) |
|
| 875 |
r_m_mean_raster_name_matched <- raster_match(1,list_param_raster_match) |
|
| 876 |
|
|
| 845 | 877 |
r_m_mean_mask_raster_name <- file.path(out_dir,paste("r_m_",method_str,"_unweighted_mean_mask_",out_suffix,".tif",sep=""))
|
| 846 |
mask(raster(r_m_mean_raster_name),mask=raster(r_mask_raster_name),
|
|
| 878 |
mask(raster( r_m_mean_raster_name_matched),mask=raster(r_mask_raster_name),
|
|
| 847 | 879 |
filename=r_m_mean_mask_raster_name,overwrite=TRUE) |
| 848 | 880 |
raster_name <- r_m_mean_mask_raster_name |
| 849 | 881 |
}else{
|
| ... | ... | |
| 1424 | 1456 |
#test_lf <- lapply(1,FUN=generate_residuals_raster,list_param=list_param_generate_residuals_raster) |
| 1425 | 1457 |
|
| 1426 | 1458 |
list_pred_res_obj <- mclapply(1:length(lf),FUN=generate_residuals_raster,list_param=list_param_generate_residuals_raster,mc.preschedule=FALSE,mc.cores = num_cores) |
| 1427 |
|
|
| 1459 |
## Add to df_raster_pred_tiles |
|
| 1460 |
|
|
| 1461 |
|
|
| 1428 | 1462 |
#write output |
| 1429 | 1463 |
accuracy_residuals_obj <-list(list_pred_res_obj,data_df,df_raster_pred_tiles) |
| 1430 | 1464 |
names(accuracy_residuals_obj)<-c("list_pred_res_obj","data_df","df_raster_pred_tiles")
|
| ... | ... | |
| 1434 | 1468 |
return(accuracy_residuals_obj) |
| 1435 | 1469 |
} |
| 1436 | 1470 |
|
| 1471 |
##Also found in accuracy assessment script:global_run_scalingup_assessment_part1_functions_02112015.R |
|
| 1472 |
#This extract a data.frame object from raster prediction obj and combine them in one data.frame |
|
| 1473 |
extract_from_list_obj<-function(obj_list,list_name){
|
|
| 1474 |
#extract object from list of list. This useful for raster_prediction_obj |
|
| 1475 |
library(plyr) |
|
| 1476 |
list_tmp<-vector("list",length(obj_list))
|
|
| 1477 |
for (i in 1:length(obj_list)){
|
|
| 1478 |
tmp <- obj_list[[i]] |
|
| 1479 |
if(inherits(tmp,"try-error")){
|
|
| 1480 |
print(paste("no model generated or error in list",sep=" ")) #change message for any model type...
|
|
| 1481 |
list_tmp[[i]] <- NULL #double bracket to return data.frame |
|
| 1482 |
}else{
|
|
| 1483 |
#tmp<-obj_list[[i]][[list_name]] #double bracket to return data.frame |
|
| 1484 |
list_tmp[[i]] <- as.data.frame(tmp[[list_name]]) |
|
| 1485 |
} |
|
| 1486 |
# |
|
| 1487 |
#tmp<-obj_list[[i]][[list_name]] #double bracket to return data.frame |
|
| 1488 |
#list_tmp[[i]]<- as.data.frame(tmp) #if spdf |
|
| 1489 |
} |
|
| 1490 |
# |
|
| 1491 |
#list_tmp <-list_tmp[!is.null(list_tmp)] |
|
| 1492 |
list_tmp <- list_tmp[unlist(lapply(list_tmp,FUN=function(x){!is.null(x)}))]
|
|
| 1493 |
|
|
| 1494 |
tb_list_tmp<-do.call(rbind.fill,list_tmp) #long rownames |
|
| 1495 |
#tb_list_tmp<-do.call(rbind,list_tmp) #long rownames |
|
| 1496 |
return(tb_list_tmp) #this is a data.frame |
|
| 1497 |
} |
|
| 1498 |
|
|
| 1437 | 1499 |
##################### END OF SCRIPT ###################### |
Also available in: Unified diff
adding functions for residuals kriging options and clean up function mosaicing script