Revision 82d6f44b
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/20/2016
|
|
| 7 |
#MODIFIED ON: 04/22/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,use_int=TRUE,scaling=100){
|
|
| 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,data_type="Float32",scaling=NULL,values_range=NULL){
|
|
| 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 |
#15)data_type: Float32 is default values for mosaicing |
|
| 492 |
#16)scaling: if NULL, use scaling 1, numeric value to multiply the values before conversion to integer |
|
| 493 |
#17)values_range: if NULL, don't screen |
|
| 493 | 494 |
# |
| 494 | 495 |
#OUTPUT: |
| 495 | 496 |
# Object is produced with 3 components: |
| ... | ... | |
| 508 | 509 |
out_suffix_str_tmp <- paste0(out_suffix,"_tmp") |
| 509 | 510 |
#} |
| 510 | 511 |
|
| 511 |
if(use_int==TRUE){
|
|
| 512 |
data_type <- "Int16" #should be a parameter!! |
|
| 513 |
}else{
|
|
| 514 |
data_type <- "Float32" |
|
| 515 |
} |
|
| 512 |
#if(data_type==NULL){
|
|
| 513 |
# data_type <- "Int16" #should be a parameter!!
|
|
| 514 |
#}else{
|
|
| 515 |
# data_type <- "Float32"
|
|
| 516 |
#}
|
|
| 516 | 517 |
if(is.null(scaling)){
|
| 517 | 518 |
scaling <- 1 |
| 518 | 519 |
} |
| 519 |
valid_range <- c(-100,100) #pass this as parameter!! (in the next update) |
|
| 520 |
valid_range <- values_range #if NULL don't screen values!! |
|
| 521 |
#valid_range <- c(-100,100) #pass this as parameter!! (in the next update) |
|
| 520 | 522 |
|
| 521 | 523 |
lf_r_weights <- vector("list",length=length(lf_mosaic))
|
| 522 | 524 |
|
| ... | ... | |
| 838 | 840 |
#if not null use the value specificied in the parameters |
| 839 | 841 |
|
| 840 | 842 |
python_cmd <- file.path(python_bin,"gdal_calc.py") |
| 841 |
cmd_str <- paste(python_cmd, |
|
| 843 |
cmd_str3 <- paste(python_cmd,
|
|
| 842 | 844 |
paste("-A ", r_prod_sum_raster_name,sep=""),
|
| 843 | 845 |
paste("-B ", r_weights_sum_raster_name,sep=""),
|
| 844 | 846 |
paste("--outfile=",r_m_weighted_mean_raster_name,sep=""),
|
| ... | ... | |
| 847 | 849 |
paste("--NoDataValue=",NA_flag_val,sep=""),
|
| 848 | 850 |
paste("--calc='(A/B)*",scaling,"'",sep=""),
|
| 849 | 851 |
"--overwrite",sep=" ") #division by zero is problematic... |
| 850 |
system(cmd_str) |
|
| 852 |
system(cmd_str3)
|
|
| 851 | 853 |
|
| 852 | 854 |
###Starting rescaling |
| 853 | 855 |
##Can merge one and 2 with parentheses operations!!, make this a function? |
| ... | ... | |
| 892 | 894 |
"--overwrite",sep=" ") #division by zero is problematic... |
| 893 | 895 |
system(cmd_str6) |
| 894 | 896 |
|
| 897 |
cmd_mosaic_logfile <- file.path(out_dir,paste("cmd_mosaic_",out_suffix,".txt",sep=""))
|
|
| 898 |
#check if file exists first... |
|
| 899 |
|
|
| 895 | 900 |
### End of rescaling section |
| 896 | 901 |
#r_m_use_edge_weights_weighted_mean_rec_gam_CAI_dailyTmax_19910101_reg4_tmp.tif |
| 897 | 902 |
|
| 898 |
#browser() |
|
| 903 |
#browser() #22minutes for one mosaic
|
|
| 899 | 904 |
|
| 900 |
cmd_mosaic_logfile <- file.path(out_dir,paste("cmd_mosaic_",out_suffix,".txt",sep=""))
|
|
| 905 |
#cmd_mosaic_logfile <- file.path(out_dir,paste("cmd_mosaic_",out_suffix,".txt",sep=""))
|
|
| 901 | 906 |
|
| 902 | 907 |
writeLines(cmd_str1,con=cmd_mosaic_logfile) #weights files to mosaic |
| 903 | 908 |
#writeLines(cmd_str2,con=file.path(out_dir,paste("cmd_mosaic_",out_suffix,".txt",sep=""))) #weights files to mosaic
|
| 904 | 909 |
cat(cmd_str2, file=cmd_mosaic_logfile, append=TRUE, sep = "\n") |
| 905 |
|
|
| 910 |
cat(cmd_str3, file=cmd_mosaic_logfile, append=TRUE, sep = "\n") |
|
| 911 |
cat(cmd_str4, file=cmd_mosaic_logfile, append=TRUE, sep = "\n") |
|
| 912 |
cat(cmd_str5, file=cmd_mosaic_logfile, append=TRUE, sep = "\n") |
|
| 913 |
cat(cmd_str6, file=cmd_mosaic_logfile, append=TRUE, sep = "\n") |
|
| 914 |
|
|
| 906 | 915 |
#cmd_str <- "/nobackupp6/aguzman4/climateLayers/sharedModules/bin/gdal_calc.py -A r_prod_weights_sum_m_use_edge_weights_weighted_mean_gam_CAI_dailyTmax_19920101_reg4_run10_1500x4500_global_analyses_pred_1992_10052015.tif -B r_weights_sum_m_use_edge_weights_weighted_mean_gam_CAI_dailyTmax_19920101_reg4_run10_1500x4500_global_analyses_pred_1992_10052015.tif --outfile='test2.tif' --calc='A/B' --overwrite" |
| 907 | 916 |
|
| 908 | 917 |
#writeRaster(r_m_weighted_mean, NAflag=NA_flag_val,filename=raster_name,overwrite=TRUE) |
| ... | ... | |
| 921 | 930 |
|
| 922 | 931 |
#undebug(raster_match) |
| 923 | 932 |
r_m_weighted_mean_raster_name_matched <- raster_match(1,list_param_raster_match) |
| 924 |
|
|
| 933 |
#output |
|
| 925 | 934 |
r_m_weighted_mean_mask_raster_name <- file.path(out_dir_str,paste("r_m_",mosaic_method,"_weighted_mean_mask_",out_suffix,".tif",sep=""))
|
| 926 |
mask(raster(r_m_weighted_mean_raster_name_matched),mask=raster(r_mask_raster_name), |
|
| 927 |
filename=r_m_weighted_mean_mask_raster_name,overwrite=TRUE) |
|
| 935 |
|
|
| 936 |
#mask(raster(r_m_weighted_mean_raster_name_matched),mask=raster(r_mask_raster_name), |
|
| 937 |
# filename=r_m_weighted_mean_mask_raster_name,overwrite=TRUE) |
|
| 938 |
|
|
| 939 |
##Now combine A and B and multiply by mask?? This must be faster than the mask option in R |
|
| 940 |
#The mask must be in 1,NA format with 1 being valid values being considered in roi. |
|
| 941 |
#r_m_weighted_mean_raster_name_rec <- file.path(out_dir_str,paste("r_m_",mosaic_method,"_weighted_mean_rec_",out_suffix,"_tmp",".tif",sep=""))
|
|
| 942 |
browser() |
|
| 943 |
|
|
| 944 |
cmd_str7 <- paste(python_cmd, |
|
| 945 |
paste("-A ", r_m_weighted_mean_raster_name_matched,sep=""),
|
|
| 946 |
paste("-B ", r_mask_raster_name,sep=""),
|
|
| 947 |
paste("--outfile=",r_m_weighted_mean_mask_raster_name,sep=""),
|
|
| 948 |
paste("--type=",data_type,sep=""),
|
|
| 949 |
"--co='COMPRESS=LZW'", |
|
| 950 |
paste("--NoDataValue=",NA_flag_val,sep=""),
|
|
| 951 |
paste("--calc='(A*B)'",sep=""),
|
|
| 952 |
"--overwrite",sep=" ") #division by zero is problematic... |
|
| 953 |
system(cmd_str7) |
|
| 954 |
cat(cmd_str7, file=cmd_mosaic_logfile, append=TRUE, sep = "\n") |
|
| 955 |
|
|
| 956 |
##Set min max and NA value |
|
| 957 |
r_mosaiced <- raster(r_m_weighted_mean_mask_raster_name) |
|
| 958 |
r_mosaiced <- setMinMax(r_mosaiced) |
|
| 959 |
NAvalue(r_mosaiced) <- NA_flag_val |
|
| 928 | 960 |
raster_name <- r_m_weighted_mean_mask_raster_name |
| 961 |
|
|
| 929 | 962 |
}else{
|
| 930 | 963 |
raster_name <- r_m_weighted_mean_raster_name |
| 931 | 964 |
} |
Also available in: Unified diff
fixing mask used and no data values