Revision 3ccd131a
Added by Benoit Parmentier over 8 years ago
| climate/research/oregon/interpolation/global_run_scalingup_mosaicing_function.R | ||
|---|---|---|
| 650 | 650 |
#plot(r_ref) |
| 651 | 651 |
|
| 652 | 652 |
if(mosaic_method%in%c("use_linear_weights","use_sine_weights","use_edge_weights")){
|
| 653 |
|
|
| 654 |
|
|
| 655 |
##################### |
|
| 656 |
###### PART 4: compute the weighted mean with the mosaic function ##### |
|
| 657 | 653 |
|
| 658 | 654 |
if(algorithm=="python"){
|
| 659 | 655 |
|
| ... | ... | |
| 835 | 831 |
|
| 836 | 832 |
##Add here the int and scaling? |
| 837 | 833 |
#note that the nodata was fixed... |
| 838 |
|
|
| 839 |
|
|
| 840 | 834 |
#if not null use the value specificied in the parameters |
| 841 | 835 |
|
| 842 | 836 |
python_cmd <- file.path(python_bin,"gdal_calc.py") |
| ... | ... | |
| 851 | 845 |
"--overwrite",sep=" ") #division by zero is problematic... |
| 852 | 846 |
system(cmd_str3) |
| 853 | 847 |
|
| 854 |
###Starting rescaling |
|
| 848 |
#writeRaster(r_m_weighted_mean, NAflag=NA_flag_val,filename=raster_name,overwrite=TRUE) |
|
| 849 |
|
|
| 850 |
###Starting rescaling, switched by masking first then screening to avoid potential problems |
|
| 855 | 851 |
##Can merge one and 2 with parentheses operations!!, make this a function? |
| 856 | 852 |
##Reclassify with valid range: -100,100 |
| 857 |
|
|
| 853 |
raster_name <- r_m_weighted_mean_raster_name |
|
| 858 | 854 |
max_val <- valid_range[2]*scaling #set min_valid |
| 859 |
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=""))
|
|
| 855 |
raster_name_rec1 <- file.path(out_dir_str,paste("r_m_",mosaic_method,"_weighted_mean_rec1_",out_suffix,"_tmp",".tif",sep=""))
|
|
| 860 | 856 |
#rec_tmp1 <- file.path(out_dir_str,paste("r_m_",mosaic_method,"_weighted_mean_rec_",out_suffix,".tif",sep=""))
|
| 861 | 857 |
cmd_str4 <- paste(python_cmd, |
| 862 |
paste("-A ", r_m_weighted_mean_raster_name,sep=""),
|
|
| 863 |
paste("--outfile=",r_m_weighted_mean_raster_name_rec1,sep=""),
|
|
| 858 |
paste("-A ", raster_name,sep=""),
|
|
| 859 |
paste("--outfile=",raster_name_rec1,sep=""),
|
|
| 864 | 860 |
paste("--type=",data_type,sep=""),
|
| 865 | 861 |
paste("--NoDataValue=",NA_flag_val,sep=""),
|
| 866 | 862 |
"--co='COMPRESS=LZW'", |
| ... | ... | |
| 868 | 864 |
"--overwrite",sep=" ") #division by zero is problematic... |
| 869 | 865 |
system(cmd_str4) |
| 870 | 866 |
|
| 871 |
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=""))
|
|
| 867 |
raster_name_rec2 <- file.path(out_dir_str,paste("r_m_",mosaic_method,"_weighted_mean_rec2_",out_suffix,"_tmp",".tif",sep=""))
|
|
| 872 | 868 |
min_val <- valid_range[1]*scaling #set min_valid as a input |
| 873 | 869 |
cmd_str5 <- paste(python_cmd, |
| 874 |
paste("-A ", r_m_weighted_mean_raster_name,sep=""),
|
|
| 875 |
paste("--outfile=",r_m_weighted_mean_raster_name_rec2,sep=""),
|
|
| 870 |
paste("-A ", raster_name,sep=""),
|
|
| 871 |
paste("--outfile=",raster_name_rec2,sep=""),
|
|
| 876 | 872 |
paste("--type=",data_type,sep=""),
|
| 877 | 873 |
"--co='COMPRESS=LZW'", |
| 878 | 874 |
paste("--NoDataValue=",NA_flag_val,sep=""),
|
| ... | ... | |
| 881 | 877 |
system(cmd_str5) |
| 882 | 878 |
|
| 883 | 879 |
##Now combine A and B and multiply by mask?? |
| 884 |
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=""))
|
|
| 880 |
r_m_weighted_mean_raster_name_rec <- file.path(out_dir_str,paste("r_m_",mosaic_method,"_weighted_mean_masked_",out_suffix,".tif",sep=""))
|
|
| 885 | 881 |
cmd_str6 <- paste(python_cmd, |
| 886 |
paste("-A ", r_m_weighted_mean_raster_name_rec1,sep=""),
|
|
| 887 |
paste("-B ", r_m_weighted_mean_raster_name_rec2,sep=""),
|
|
| 888 |
paste("-C ", r_m_weighted_mean_raster_name,sep=""), #if mask exists
|
|
| 882 |
paste("-A ", raster_name_rec1,sep=""),
|
|
| 883 |
paste("-B ", raster_name_rec2,sep=""),
|
|
| 884 |
paste("-C ", raster_name,sep=""), #if mask exists
|
|
| 889 | 885 |
paste("--outfile=",r_m_weighted_mean_raster_name_rec,sep=""),
|
| 890 | 886 |
paste("--type=",data_type,sep=""),
|
| 891 | 887 |
"--co='COMPRESS=LZW'", |
| ... | ... | |
| 894 | 890 |
"--overwrite",sep=" ") #division by zero is problematic... |
| 895 | 891 |
system(cmd_str6) |
| 896 | 892 |
|
| 897 |
cmd_mosaic_logfile <- file.path(out_dir,paste("cmd_mosaic_",out_suffix,".txt",sep=""))
|
|
| 898 | 893 |
#check if file exists first... |
| 899 | 894 |
|
| 900 | 895 |
### End of rescaling section |
| ... | ... | |
| 914 | 909 |
|
| 915 | 910 |
#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" |
| 916 | 911 |
|
| 917 |
#writeRaster(r_m_weighted_mean, NAflag=NA_flag_val,filename=raster_name,overwrite=TRUE) |
|
| 918 | 912 |
#now use the mask |
| 919 | 913 |
if(!is.null(r_mask_raster_name)){
|
| 920 | 914 |
#different extent between mask and output if match extent is false!! |
| 921 | 915 |
#match resolution and extent first |
| 922 | 916 |
|
| 923 | 917 |
#lf_files <- c(r_m_weighted_mean_raster_name) #file(s) to be matched |
| 924 |
lf_files <- c(r_m_weighted_mean_raster_name_rec) |
|
| 918 |
lf_files <- c(r_m_weighted_mean_raster_name_rec) #match to mask
|
|
| 925 | 919 |
rast_ref <- r_mask_raster_name |
| 926 | 920 |
##Maching resolution is probably only necessary for the r mosaic function |
| 927 | 921 |
#Modify later to take into account option R or python... |
| ... | ... | |
| 939 | 933 |
##Now combine A and B and multiply by mask?? This must be faster than the mask option in R |
| 940 | 934 |
#The mask must be in 1,NA format with 1 being valid values being considered in roi. |
| 941 | 935 |
#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 |
|
|
| 936 |
#browser() |
|
| 937 |
cmd_mosaic_logfile <- file.path(out_dir,paste("cmd_mosaic_",out_suffix,".txt",sep=""))
|
|
| 938 |
|
|
| 944 | 939 |
cmd_str7 <- paste(python_cmd, |
| 945 | 940 |
paste("-A ", r_m_weighted_mean_raster_name_matched,sep=""),
|
| 946 | 941 |
paste("-B ", r_mask_raster_name,sep=""),
|
| ... | ... | |
| 958 | 953 |
r_mosaiced <- setMinMax(r_mosaiced) |
| 959 | 954 |
NAvalue(r_mosaiced) <- NA_flag_val |
| 960 | 955 |
raster_name <- r_m_weighted_mean_mask_raster_name |
| 961 |
|
|
| 956 |
browser() |
|
| 957 |
#-32,768 is NA |
|
| 962 | 958 |
}else{
|
| 963 | 959 |
raster_name <- r_m_weighted_mean_raster_name |
| 964 | 960 |
} |
| 965 |
} |
|
| 966 | 961 |
|
| 962 |
} #end of weighted |
|
| 963 |
|
|
| 964 |
### |
|
| 967 | 965 |
if(mosaic_method=="unweighted"){
|
| 968 | 966 |
#### Fourth use original images |
| 969 | 967 |
#macth file to mosaic extent using the original predictions |
Also available in: Unified diff
more debugging of mask and testing of code