Revision 1eba4f7c
Added by Benoit Parmentier over 8 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: 04/11/2016
|
|
| 7 |
#MODIFIED ON: 04/20/2016
|
|
| 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 |
| ... | ... | |
| 463 | 463 |
return(raster_name) |
| 464 | 464 |
} |
| 465 | 465 |
|
| 466 |
mosaicFiles <- function(lf_mosaic,mosaic_method="unweighted",num_cores=1,r_mask_raster_name=NULL,python_bin=NULL,mosaic_python="/nobackupp6/aguzman4/climateLayers/sharedCode/gdal_merge_sum_noDataTest.py",algorithm="R",match_extent=TRUE,df_points=NULL,NA_flag_val=-9999,file_format=".tif",out_suffix=NULL,out_dir=NULL,tmp_files=FALSE){
|
|
| 466 |
mosaicFiles <- function(lf_mosaic,mosaic_method="unweighted",num_cores=1,r_mask_raster_name=NULL,python_bin=NULL,mosaic_python="/nobackupp6/aguzman4/climateLayers/sharedCode/gdal_merge_sum_noDataTest.py",algorithm="R",match_extent=TRUE,df_points=NULL,NA_flag_val=-9999,file_format=".tif",out_suffix=NULL,out_dir=NULL,tmp_files=FALSE,use_int=TRUE,scaling=100){
|
|
| 467 | 467 |
#This functions mosaics tiles/files give a list of files. |
| 468 | 468 |
#There are four options to mosaic: use_sine_weights,use_edge,use_linear_weights, unweighted |
| 469 | 469 |
#Sine weights fits sine fuctions across rows and column producing elliptical/spherical patterns from center |
| ... | ... | |
| 488 | 488 |
#12)algorithm: use R or python function |
| 489 | 489 |
#13)match extent: if TRUE match extent before mosaicing |
| 490 | 490 |
#14)tmp_files: if TRUE then keep temporary files |
| 491 |
#15)use_int: if TRUE use int values for mosaicing |
|
| 492 |
#16)scaling: numeric value to multiply the values before conversion to integer |
|
| 491 | 493 |
# |
| 492 | 494 |
#OUTPUT: |
| 493 | 495 |
# Object is produced with 3 components: |
| ... | ... | |
| 814 | 816 |
#r_m_weighted_mean <- r_prod_sum/r_weights_sum #this is the mosaic using weighted mean... |
| 815 | 817 |
|
| 816 | 818 |
r_m_weighted_mean_raster_name <- file.path(out_dir_str,paste("r_m_",mosaic_method,"_weighted_mean_",out_suffix,".tif",sep=""))
|
| 817 |
|
|
| 819 |
#r_m_weighted_mean_raster_name <- "test_tmp.tif" |
|
| 818 | 820 |
if(is.null(python_bin)){
|
| 819 | 821 |
python_bin="" |
| 820 | 822 |
} |
| 821 | 823 |
|
| 824 |
##Add here the int and scaling? |
|
| 825 |
#note that the nodata was fixed... |
|
| 826 |
#browser() |
|
| 827 |
if(use_int==TRUE){
|
|
| 828 |
data_type <- "int16" |
|
| 829 |
}else{
|
|
| 830 |
data_type <- "Float32" |
|
| 831 |
} |
|
| 832 |
if(is.null(scaling)){
|
|
| 833 |
scaling <- 1 |
|
| 834 |
} |
|
| 835 |
#if not null use the value specificied in the parameters |
|
| 836 |
|
|
| 822 | 837 |
python_cmd <- file.path(python_bin,"gdal_calc.py") |
| 823 | 838 |
cmd_str <- paste(python_cmd, |
| 824 | 839 |
paste("-A ", r_prod_sum_raster_name,sep=""),
|
| 825 | 840 |
paste("-B ", r_weights_sum_raster_name,sep=""),
|
| 826 | 841 |
paste("--outfile=",r_m_weighted_mean_raster_name,sep=""),
|
| 827 |
paste("NoDataValue=",NA_flag_val,sep=""),
|
|
| 828 |
"--calc='A/B'","--overwrite",sep=" ") #division by zero is problematic... |
|
| 842 |
paste("--type=",data_type,sep=""),
|
|
| 843 |
paste("--NoDataValue=",NA_flag_val,sep=""),
|
|
| 844 |
#"--calc='A/B'","--overwrite",sep=" ") #division by zero is problematic... |
|
| 845 |
paste("--calc='(A/B)*",scaling,"'",sep=""),
|
|
| 846 |
"--overwrite",sep=" ") #division by zero is problematic... |
|
| 829 | 847 |
system(cmd_str) |
| 848 |
|
|
| 849 |
##Can merge one and 2 with parentheses operations!! |
|
| 850 |
##Reclassify with valid range: -100,100 |
|
| 851 |
max_val <- 100*scaling #set min_valid |
|
| 852 |
r_m_weighted_mean_raster_name_rec1 <- file.path(out_dir_str,paste("r_m_",mosaic_method,"_weighted_mean_rec1_",out_suffix,"_tmp",".tif",sep=""))
|
|
| 853 |
#rec_tmp1 <- file.path(out_dir_str,paste("r_m_",mosaic_method,"_weighted_mean_rec_",out_suffix,".tif",sep=""))
|
|
| 854 |
cmd_str4 <- paste(python_cmd, |
|
| 855 |
paste("-A ", r_m_weighted_mean_raster_name,sep=""),
|
|
| 856 |
paste("--outfile=",r_m_weighted_mean_raster_name_rec1,sep=""),
|
|
| 857 |
paste("--type=",data_type,sep=""),
|
|
| 858 |
paste("--NoDataValue=",NA_flag_val,sep=""),
|
|
| 859 |
paste("--calc='(A>",max_val,")*",NA_flag_val,"'",sep=""),
|
|
| 860 |
"--overwrite",sep=" ") #division by zero is problematic... |
|
| 861 |
system(cmd_str4) |
|
| 862 |
|
|
| 863 |
r_m_weighted_mean_raster_name_rec2 <- file.path(out_dir_str,paste("r_m_",mosaic_method,"_weighted_mean_rec2_",out_suffix,"_tmp",".tif",sep=""))
|
|
| 864 |
min_val <- -100*scaling #set min_valid as a input |
|
| 865 |
cmd_str5 <- paste(python_cmd, |
|
| 866 |
paste("-A ", r_m_weighted_mean_raster_name,sep=""),
|
|
| 867 |
paste("--outfile=",r_m_weighted_mean_raster_name_rec2,sep=""),
|
|
| 868 |
paste("--type=",data_type,sep=""),
|
|
| 869 |
paste("--NoDataValue=",NA_flag_val,sep=""),
|
|
| 870 |
paste("--calc='(A<",min_val,")*",NA_flag_val,"'",sep=""),
|
|
| 871 |
"--overwrite",sep=" ") #division by zero is problematic... |
|
| 872 |
system(cmd_str5) |
|
| 873 |
|
|
| 874 |
#r_m_weighted_mean_raster_name_rec3 <- file.path(out_dir_str,paste("r_m_",mosaic_method,"_weighted_mean_rec3_",out_suffix,"_tmp",".tif",sep=""))
|
|
| 875 |
#min_valid <- -100*scaling #set min_valid as a input |
|
| 876 |
#cmd_str6 <- paste(python_cmd, |
|
| 877 |
# paste("-A ", r_m_weighted_mean_raster_name,sep=""),
|
|
| 878 |
# paste("--outfile=",r_m_weighted_mean_raster_name_rec2,sep=""),
|
|
| 879 |
# paste("--type=",data_type,sep=""),
|
|
| 880 |
# paste("--NoDataValue=",NA_flag_val,sep=""),
|
|
| 881 |
# paste("--calc='(",min_val,"<A<",max_val,")","'",sep=""),
|
|
| 882 |
# "--overwrite",sep=" ") #division by zero is problematic... |
|
| 883 |
#system(cmd_str6) |
|
| 884 |
|
|
| 885 |
##Now combine A and B and multiply by mask?? |
|
| 886 |
r_m_weighted_mean_raster_name_rec3 <- file.path(out_dir_str,paste("r_m_",mosaic_method,"_weighted_mean_rec3_",out_suffix,"_tmp",".tif",sep=""))
|
|
| 887 |
min_valid <- -100*scaling #set min_valid as a input |
|
| 888 |
cmd_str6 <- paste(python_cmd, |
|
| 889 |
paste("-A ", r_m_weighted_mean_raster_name_rec1,sep=""),
|
|
| 890 |
paste("-B ", r_m_weighted_mean_raster_name_rec2,sep=""),
|
|
| 891 |
paste("-C ", r_m_weighted_mean_raster_name,sep=""), #if mask exists
|
|
| 892 |
paste("--outfile=",r_m_weighted_mean_raster_name_rec3,sep=""),
|
|
| 893 |
paste("--type=",data_type,sep=""),
|
|
| 894 |
paste("--NoDataValue=",NA_flag_val,sep=""),
|
|
| 895 |
paste("--calc='(((((A+B))<",min_val,")*",NA_flag_val,")+1)*C'",sep=""),
|
|
| 896 |
"--overwrite",sep=" ") #division by zero is problematic... |
|
| 897 |
system(cmd_str6) |
|
| 898 |
|
|
| 899 |
#r_m_weighted_mean_raster_name_rec4 <- file.path(out_dir_str,paste("r_m_",mosaic_method,"_weighted_mean_rec4_",out_suffix,"_tmp",".tif",sep=""))
|
|
| 900 |
#min_valid <- -100*scaling #set min_valid as a input |
|
| 901 |
#cmd_str6 <- paste(python_cmd, |
|
| 902 |
# paste("-A ", r_m_weighted_mean_raster_name_rec3,sep=""),
|
|
| 903 |
# #paste("-C ", r_m_weighted_mean_raster_name_rec2,sep=""), #if mask exists
|
|
| 904 |
# paste("--outfile=",r_m_weighted_mean_raster_name_rec4,sep=""),
|
|
| 905 |
# paste("--type=",data_type,sep=""),
|
|
| 906 |
# paste("--NoDataValue=",NA_flag_val,sep=""),
|
|
| 907 |
# paste("--calc='(A+1)'",sep=""),
|
|
| 908 |
# "--overwrite",sep=" ") #division by zero is problematic... |
|
| 909 |
#system(cmd_str6) |
|
| 910 |
|
|
| 830 | 911 |
cmd_mosaic_logfile <- file.path(out_dir,paste("cmd_mosaic_",out_suffix,".txt",sep=""))
|
| 912 |
|
|
| 831 | 913 |
writeLines(cmd_str1,con=cmd_mosaic_logfile) #weights files to mosaic |
| 832 | 914 |
#writeLines(cmd_str2,con=file.path(out_dir,paste("cmd_mosaic_",out_suffix,".txt",sep=""))) #weights files to mosaic
|
| 833 | 915 |
cat(cmd_str2, file=cmd_mosaic_logfile, append=TRUE, sep = "\n") |
Also available in: Unified diff
valid value range with gdal calc in mosaic function