Revision 3b4e827a
Added by Benoit Parmentier over 8 years ago
| climate/research/oregon/interpolation/global_run_scalingup_mosaicing_function.R | ||
|---|---|---|
| 508 | 508 |
out_suffix_str_tmp <- paste0(out_suffix,"_tmp") |
| 509 | 509 |
#} |
| 510 | 510 |
|
| 511 |
if(use_int==TRUE){
|
|
| 512 |
data_type <- "Int16" #should be a parameter!! |
|
| 513 |
}else{
|
|
| 514 |
data_type <- "Float32" |
|
| 515 |
} |
|
| 516 |
if(is.null(scaling)){
|
|
| 517 |
scaling <- 1 |
|
| 518 |
} |
|
| 519 |
valid_range <- c(-100,100) #pass this as parameter!! (in the next update) |
|
| 520 |
|
|
| 511 | 521 |
lf_r_weights <- vector("list",length=length(lf_mosaic))
|
| 512 | 522 |
|
| 513 | 523 |
############### |
| ... | ... | |
| 823 | 833 |
|
| 824 | 834 |
##Add here the int and scaling? |
| 825 | 835 |
#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 |
} |
|
| 836 |
|
|
| 837 |
|
|
| 835 | 838 |
#if not null use the value specificied in the parameters |
| 836 | 839 |
|
| 837 | 840 |
python_cmd <- file.path(python_bin,"gdal_calc.py") |
| ... | ... | |
| 840 | 843 |
paste("-B ", r_weights_sum_raster_name,sep=""),
|
| 841 | 844 |
paste("--outfile=",r_m_weighted_mean_raster_name,sep=""),
|
| 842 | 845 |
paste("--type=",data_type,sep=""),
|
| 846 |
"--co='COMPRESS=LZW'", |
|
| 843 | 847 |
paste("--NoDataValue=",NA_flag_val,sep=""),
|
| 844 |
#"--calc='A/B'","--overwrite",sep=" ") #division by zero is problematic... |
|
| 845 | 848 |
paste("--calc='(A/B)*",scaling,"'",sep=""),
|
| 846 | 849 |
"--overwrite",sep=" ") #division by zero is problematic... |
| 847 | 850 |
system(cmd_str) |
| 848 | 851 |
|
| 849 |
##Can merge one and 2 with parentheses operations!! |
|
| 852 |
###Starting rescaling |
|
| 853 |
##Can merge one and 2 with parentheses operations!!, make this a function? |
|
| 850 | 854 |
##Reclassify with valid range: -100,100 |
| 851 |
max_val <- 100*scaling #set min_valid |
|
| 855 |
|
|
| 856 |
max_val <- valid_range[2]*scaling #set min_valid |
|
| 852 | 857 |
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 | 858 |
#rec_tmp1 <- file.path(out_dir_str,paste("r_m_",mosaic_method,"_weighted_mean_rec_",out_suffix,".tif",sep=""))
|
| 854 | 859 |
cmd_str4 <- paste(python_cmd, |
| ... | ... | |
| 856 | 861 |
paste("--outfile=",r_m_weighted_mean_raster_name_rec1,sep=""),
|
| 857 | 862 |
paste("--type=",data_type,sep=""),
|
| 858 | 863 |
paste("--NoDataValue=",NA_flag_val,sep=""),
|
| 864 |
"--co='COMPRESS=LZW'", |
|
| 859 | 865 |
paste("--calc='(A>",max_val,")*",NA_flag_val,"'",sep=""),
|
| 860 | 866 |
"--overwrite",sep=" ") #division by zero is problematic... |
| 861 | 867 |
system(cmd_str4) |
| 862 | 868 |
|
| 863 | 869 |
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
|
|
| 870 |
min_val <- valid_range[1]*scaling #set min_valid as a input
|
|
| 865 | 871 |
cmd_str5 <- paste(python_cmd, |
| 866 | 872 |
paste("-A ", r_m_weighted_mean_raster_name,sep=""),
|
| 867 | 873 |
paste("--outfile=",r_m_weighted_mean_raster_name_rec2,sep=""),
|
| 868 | 874 |
paste("--type=",data_type,sep=""),
|
| 875 |
"--co='COMPRESS=LZW'", |
|
| 869 | 876 |
paste("--NoDataValue=",NA_flag_val,sep=""),
|
| 870 | 877 |
paste("--calc='(A<",min_val,")*",NA_flag_val,"'",sep=""),
|
| 871 | 878 |
"--overwrite",sep=" ") #division by zero is problematic... |
| 872 | 879 |
system(cmd_str5) |
| 873 | 880 |
|
| 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 | 881 |
##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 |
|
| 882 |
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=""))
|
|
| 888 | 883 |
cmd_str6 <- paste(python_cmd, |
| 889 | 884 |
paste("-A ", r_m_weighted_mean_raster_name_rec1,sep=""),
|
| 890 | 885 |
paste("-B ", r_m_weighted_mean_raster_name_rec2,sep=""),
|
| 891 | 886 |
paste("-C ", r_m_weighted_mean_raster_name,sep=""), #if mask exists
|
| 892 |
paste("--outfile=",r_m_weighted_mean_raster_name_rec3,sep=""),
|
|
| 887 |
paste("--outfile=",r_m_weighted_mean_raster_name_rec,sep=""),
|
|
| 893 | 888 |
paste("--type=",data_type,sep=""),
|
| 889 |
"--co='COMPRESS=LZW'", |
|
| 894 | 890 |
paste("--NoDataValue=",NA_flag_val,sep=""),
|
| 895 | 891 |
paste("--calc='(((((A+B))<",min_val,")*",NA_flag_val,")+1)*C'",sep=""),
|
| 896 | 892 |
"--overwrite",sep=" ") #division by zero is problematic... |
| 897 | 893 |
system(cmd_str6) |
| 898 | 894 |
|
| 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) |
|
| 895 |
### End of rescaling section |
|
| 896 |
#r_m_use_edge_weights_weighted_mean_rec_gam_CAI_dailyTmax_19910101_reg4_tmp.tif |
|
| 897 |
|
|
| 898 |
#browser() |
|
| 910 | 899 |
|
| 911 | 900 |
cmd_mosaic_logfile <- file.path(out_dir,paste("cmd_mosaic_",out_suffix,".txt",sep=""))
|
| 912 | 901 |
|
| ... | ... | |
| 922 | 911 |
#different extent between mask and output if match extent is false!! |
| 923 | 912 |
#match resolution and extent first |
| 924 | 913 |
|
| 925 |
lf_files <- c(r_m_weighted_mean_raster_name) #file(s) to be matched |
|
| 914 |
#lf_files <- c(r_m_weighted_mean_raster_name) #file(s) to be matched |
|
| 915 |
lf_files <- c(r_m_weighted_mean_raster_name_rec) |
|
| 926 | 916 |
rast_ref <- r_mask_raster_name |
| 927 | 917 |
##Maching resolution is probably only necessary for the r mosaic function |
| 928 | 918 |
#Modify later to take into account option R or python... |
| ... | ... | |
| 998 | 988 |
mask(raster( r_m_mean_raster_name_matched),mask=raster(r_mask_raster_name), |
| 999 | 989 |
filename=r_m_mean_mask_raster_name,overwrite=TRUE) |
| 1000 | 990 |
raster_name <- r_m_mean_mask_raster_name |
| 991 |
r_raster <- raster(raster_name) |
|
| 992 |
r_raster <- setMinMax(r_raster) #set correct min and max in the file |
|
| 1001 | 993 |
}else{
|
| 1002 | 994 |
raster_name <- r_m_mean_raster_name |
| 995 |
r_raster <- raster(raster_name) |
|
| 996 |
r_raster <- setMinMax(r_raster) #set correct min and max in the file |
|
| 1003 | 997 |
} |
| 1004 | 998 |
|
| 1005 | 999 |
} |
| ... | ... | |
| 1139 | 1133 |
m <- c(-Inf, -100, NA, |
| 1140 | 1134 |
-100, 100, 1, |
| 1141 | 1135 |
100, Inf, NA) #can change the thresholds later |
| 1136 |
|
|
| 1142 | 1137 |
rclmat <- matrix(m, ncol=3, byrow=TRUE) |
| 1143 | 1138 |
rc <- reclassify(r_test, rclmat,filename=paste("rc_tmp_",i,".tif",sep=""),dataType="FLT4S",overwrite=T)
|
| 1144 | 1139 |
file_name <- unlist(strsplit(basename(lf_m[i]),"_")) |
Also available in: Unified diff
adding compression and debugging for valid range and rescaling of mosaicing