Revision 7213ab40
Added by Benoit Parmentier about 9 years ago
| climate/research/oregon/interpolation/global_run_scalingup_mosaicing.R | ||
|---|---|---|
| 5 | 5 |
#Analyses, figures, tables and data are also produced in the script. |
| 6 | 6 |
#AUTHOR: Benoit Parmentier |
| 7 | 7 |
#CREATED ON: 04/14/2015 |
| 8 |
#MODIFIED ON: 10/26/2015
|
|
| 8 |
#MODIFIED ON: 10/28/2015
|
|
| 9 | 9 |
#Version: 5 |
| 10 | 10 |
#PROJECT: Environmental Layers project |
| 11 | 11 |
#COMMENTS: analyses run for reg4 1992 for test of mosaicing using 1500x4500km and other tiles |
| ... | ... | |
| 25 | 25 |
library(sp) # Spatial pacakge with class definition by Bivand et al. |
| 26 | 26 |
library(spdep) # Spatial pacakge with methods and spatial stat. by Bivand et al. |
| 27 | 27 |
library(rgdal) # GDAL wrapper for R, spatial utilities |
| 28 |
library(gstat) # Kriging and co-kriging by Pebesma et al. |
|
| 28 |
#library(gstat) # Kriging and co-kriging by Pebesma et al.
|
|
| 29 | 29 |
library(fields) # NCAR Spatial Interpolation methods such as kriging, splines |
| 30 | 30 |
library(raster) # Hijmans et al. package for raster processing |
| 31 | 31 |
library(gdata) # various tools with xls reading, cbindX |
| ... | ... | |
| 36 | 36 |
library(reshape) # Change shape of object, summarize results |
| 37 | 37 |
library(plotrix) # Additional plotting functions |
| 38 | 38 |
library(plyr) # Various tools including rbind.fill |
| 39 |
library(spgwr) # GWR method |
|
| 40 |
library(automap) # Kriging automatic fitting of variogram using gstat |
|
| 41 |
library(rgeos) # Geometric, topologic library of functions |
|
| 39 |
#library(spgwr) # GWR method
|
|
| 40 |
#library(automap) # Kriging automatic fitting of variogram using gstat
|
|
| 41 |
#library(rgeos) # Geometric, topologic library of functions
|
|
| 42 | 42 |
#RPostgreSQL # Interface R and Postgres, not used in this script |
| 43 | 43 |
library(gridExtra) |
| 44 | 44 |
#Additional libraries not used in workflow |
| ... | ... | |
| 49 | 49 |
|
| 50 | 50 |
#### FUNCTION USED IN SCRIPT |
| 51 | 51 |
|
| 52 |
function_mosaicing <-"global_run_scalingup_mosaicing_function_10232015.R"
|
|
| 52 |
function_mosaicing <-"global_run_scalingup_mosaicing_function_10282015.R"
|
|
| 53 | 53 |
|
| 54 |
in_dir_script <-"/home/parmentier/Data/IPLANT_project/env_layers_scripts" |
|
| 54 |
in_dir_script <-"/home/parmentier/Data/IPLANT_project/env_layers_scripts" #NCEAS UCSB |
|
| 55 |
#in_dir_script <- "/nobackupp8/bparmen1/env_layers_scripts" #NASA NEX |
|
| 55 | 56 |
source(file.path(in_dir_script,function_mosaicing)) |
| 56 | 57 |
|
| 57 | 58 |
############################################ |
| ... | ... | |
| 61 | 62 |
|
| 62 | 63 |
#in_dir <- "/data/project/layers/commons/NEX_data/mosaicing_data_test" #PARAM1 |
| 63 | 64 |
in_dir <- "/data/project/layers/commons/NEX_data/output_run10_1500x4500_global_analyses_pred_1992_10052015" #PARAM4 |
| 64 |
in_dir_tiles <- "/data/project/layers/commons/NEX_data/output_run10_1500x4500_global_analyses_pred_1992_10052015/tiles" |
|
| 65 |
#in_dir <- "/data/project/layers/commons/NEX_data/mosaicing_data_test/reg1" #North America |
|
| 65 |
#in_dir <- "/nobackupp8/bparmen1/output_run10_1500x4500_global_analyses_pred_1992_10052015" #NEX |
|
| 66 |
|
|
| 67 |
in_dir_tiles <- file.path(in_dir,"tiles") |
|
| 68 |
#in_dir_tiles <- "/nobackupp8/bparmen1/output_run10_1500x4500_global_analyses_pred_1992_10052015/tiles" #North America |
|
| 66 | 69 |
#in_dir <- "/data/project/layers/commons/NEX_data/mosaicing_data_test/reg2" #Europe |
| 67 | 70 |
#in_dir <- "/data/project/layers/commons/NEX_data/mosaicing_data_test/reg4" #South America |
| 68 | 71 |
#in_dir <- "/data/project/layers/commons/NEX_data/mosaicing_data_test/reg5" #Africa |
| ... | ... | |
| 81 | 84 |
create_out_dir_param <- FALSE #PARAM 8 |
| 82 | 85 |
|
| 83 | 86 |
#if daily mosaics NULL then mosaicas all days of the year |
| 84 |
day_to_mosaic <- c("19920101","19920102","19920103") #PARAM9
|
|
| 87 |
day_to_mosaic <- c("19920101","19920102","19920103","19920104","19920105") #PARAM9
|
|
| 85 | 88 |
|
| 86 | 89 |
#CRS_WGS84 <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +towgs84=0,0,0") #Station coords WGS84 #CONSTANT1
|
| 87 | 90 |
#CRS_locs_WGS84<-CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +towgs84=0,0,0") #Station coords WGS84
|
| ... | ... | |
| 94 | 97 |
num_cores <- 6 #PARAM 12 |
| 95 | 98 |
|
| 96 | 99 |
infile_mask <- "/data/project/layers/commons/NEX_data/output_run10_1500x4500_global_analyses_pred_1992_10052015/r_mask_reg4.tif" |
| 100 |
#infile_mask <- "/nobackupp8/bparmen1/output_run10_1500x4500_global_analyses_pred_1992_10052015/r_mask_reg4.tif" |
|
| 97 | 101 |
|
| 98 | 102 |
#tb_accuracy_name <- file.path(in_dir,paste("tb_diagnostic_v_NA","_",out_suffix_str,".txt",sep=""))
|
| 99 | 103 |
tb_accuracy_name <- "/data/project/layers/commons/NEX_data/output_run10_1500x4500_global_analyses_pred_1992_10052015/tb_diagnostic_v_NA_run10_1500x4500_global_analyses_pred_1992_10052015.txt" |
| 104 |
#tb_accuracy_name <- "/nobackupp8/bparmen1/output_run10_1500x4500_global_analyses_pred_1992_10052015/tb_diagnostic_v_NA_run10_1500x4500_global_analyses_pred_1992_10052015.txt" |
|
| 105 |
|
|
| 106 |
#python script and gdal on NEX NASA: |
|
| 107 |
#mosaic_python <- "/nobackupp6/aguzman4/climateLayers/sharedCode/" |
|
| 108 |
#python_bin <- "/nobackupp6/aguzman4/climateLayers/sharedModules/bin" |
|
| 109 |
#python script and gdal on Atlas NCEAS |
|
| 110 |
mosaic_python <- "/data/project/layers/commons/NEX_data/sharedCode" |
|
| 111 |
python_bin <- "/usr/bin" |
|
| 100 | 112 |
|
| 113 |
#algorithm <- "python" #if R use mosaic function for R, if python use modified gdalmerge script from Alberto Guzmann |
|
| 114 |
algorithm <- "R" #if R use mosaic function for R, if python use modified gdalmerge script from Alberto Guzmann |
|
| 115 |
|
|
| 101 | 116 |
########################## START SCRIPT ############################## |
| 102 | 117 |
|
| 103 | 118 |
|
| ... | ... | |
| 160 | 175 |
r1 <- raster(lf_mosaic[[1]][1]) |
| 161 | 176 |
#r2 <- raster(lf_mosaic2[2]) |
| 162 | 177 |
|
| 163 |
plot(r1) |
|
| 178 |
#plot(r1)
|
|
| 164 | 179 |
#plot(r2) |
| 165 |
r1_ac <- raster(lf_accuracy_raster[[1]][1]) |
|
| 166 |
plot(r1_ac) |
|
| 180 |
#r1_ac <- raster(lf_accuracy_raster[[1]][1])
|
|
| 181 |
#plot(r1_ac)
|
|
| 167 | 182 |
|
| 168 | 183 |
################################# |
| 169 | 184 |
#### Second mosaic tiles for the variable and accuracy metrid |
| ... | ... | |
| 178 | 193 |
out_suffix_tmp <- paste(interpolation_method,y_var_name,day_to_mosaic[i],out_suffix,sep="_") |
| 179 | 194 |
#debug(mosaicFiles) |
| 180 | 195 |
#can also loop through methods!!! |
| 181 |
python_bin <- "/usr/bin/" #python gdal bin, this may change on NEX! |
|
| 196 |
#python_bin <- "/usr/bin/" #python gdal bin, on Atlas NCEAS |
|
| 197 |
#python_bin <- "/nobackupp6/aguzman4/climateLayers/sharedModules/bin" #on NEX |
|
| 182 | 198 |
mosaic_edge_obj_prediction <- mosaicFiles(lf_mosaic[[i]], |
| 183 | 199 |
mosaic_method="use_edge_weights", |
| 184 | 200 |
num_cores=num_cores, |
| 185 | 201 |
r_mask_raster_name=infile_mask, |
| 186 | 202 |
python_bin=python_bin, |
| 187 |
df_points=NULL,NA_flag=NA_flag_val, |
|
| 188 |
file_format=file_format,out_suffix=out_suffix_tmp, |
|
| 203 |
mosaic_python=mosaic_python, |
|
| 204 |
algorithm=algorithm, |
|
| 205 |
df_points=NULL, |
|
| 206 |
NA_flag=NA_flag_val, |
|
| 207 |
file_format=file_format, |
|
| 208 |
out_suffix=out_suffix_tmp, |
|
| 189 | 209 |
out_dir=out_dir) |
| 190 | 210 |
|
| 191 | 211 |
mosaic_method <- "use_edge_weights" #this is distance from edge |
| ... | ... | |
| 198 | 218 |
num_cores=num_cores, |
| 199 | 219 |
r_mask_raster_name=infile_mask, |
| 200 | 220 |
python_bin=python_bin, |
| 201 |
df_points=NULL,NA_flag=NA_flag_val, |
|
| 202 |
file_format=file_format,out_suffix=out_suffix_tmp, |
|
| 221 |
mosaic_python=mosaic_python, |
|
| 222 |
algorithm=algorithm, |
|
| 223 |
df_points=NULL, |
|
| 224 |
NA_flag=NA_flag_val, |
|
| 225 |
file_format=file_format, |
|
| 226 |
out_suffix=out_suffix_tmp, |
|
| 203 | 227 |
out_dir=out_dir) |
| 204 | 228 |
|
| 205 | 229 |
list_mosaic_obj[[i]] <- list(prediction=mosaic_edge_obj_prediction,accuracy=mosaic_edge_obj_accuracy) |
| ... | ... | |
| 209 | 233 |
###### PART 2: Analysis and figures for the outputs of mosaic function ##### |
| 210 | 234 |
|
| 211 | 235 |
#### compute and aspect and slope with figures |
| 212 |
list_lf_mosaic_obj <- vector("list",length(day_to_mosaic))
|
|
| 213 |
lf_mean_mosaic <- vector("list",length(mosaicing_method))#2methods only
|
|
| 214 |
l_method_mosaic <- vector("list",length(mosaicing_method))
|
|
| 215 |
list_out_suffix <- vector("list",length(mosaicing_method))
|
|
| 236 |
#list_lf_mosaic_obj <- vector("list",length(day_to_mosaic))
|
|
| 237 |
#lf_mean_mosaic <- vector("list",length(mosaicing_method))#2methods only
|
|
| 238 |
#l_method_mosaic <- vector("list",length(mosaicing_method))
|
|
| 239 |
#list_out_suffix <- vector("list",length(mosaicing_method))
|
|
| 216 | 240 |
|
| 217 |
for(i in 1:length(day_to_mosaic)){
|
|
| 218 |
list_lf_mosaic_obj[[i]] <- list.files(path=out_dir,pattern=paste("*",day_to_mosaic[i],
|
|
| 219 |
"_.*.RData",sep="")) |
|
| 220 |
lf_mean_mosaic[[i]] <- unlist(lapply(list_lf_mosaic_obj[[i]],function(x){load_obj(x)[["mean_mosaic"]]}))
|
|
| 221 |
l_method_mosaic[[i]] <- paste(unlist(lapply(list_lf_mosaic_obj[[i]],function(x){load_obj(x)[["method"]]})),day_to_mosaic[i],sep="_")
|
|
| 222 |
list_out_suffix[[i]] <- unlist(paste(l_method_mosaic[[i]],day_to_mosaic[[i]],out_suffix,sep="_")) |
|
| 223 |
} |
|
| 241 |
#for(i in 1:length(day_to_mosaic)){
|
|
| 242 |
# list_lf_mosaic_obj[[i]] <- list.files(path=out_dir,pattern=paste("*",day_to_mosaic[i],
|
|
| 243 |
# "_.*.RData",sep=""))
|
|
| 244 |
# lf_mean_mosaic[[i]] <- unlist(lapply(list_lf_mosaic_obj[[i]],function(x){load_obj(x)[["mean_mosaic"]]}))
|
|
| 245 |
# l_method_mosaic[[i]] <- paste(unlist(lapply(list_lf_mosaic_obj[[i]],function(x){load_obj(x)[["method"]]})),day_to_mosaic[i],sep="_")
|
|
| 246 |
# list_out_suffix[[i]] <- unlist(paste(l_method_mosaic[[i]],day_to_mosaic[[i]],out_suffix,sep="_"))
|
|
| 247 |
#}
|
|
| 224 | 248 |
|
| 225 | 249 |
|
| 226 |
list_param_plot_mosaic <- list(lf_mosaic=unlist(lf_mean_mosaic), |
|
| 227 |
method=unlist(l_method_mosaic), |
|
| 228 |
out_suffix=unlist(list_out_suffix)) |
|
| 250 |
#list_param_plot_mosaic <- list(lf_mosaic=unlist(lf_mean_mosaic),
|
|
| 251 |
# method=unlist(l_method_mosaic),
|
|
| 252 |
# out_suffix=unlist(list_out_suffix))
|
|
| 229 | 253 |
|
| 230 | 254 |
#plot and produce png movie... |
| 231 | 255 |
#plot_mosaic(1,list_param=list_param_plot_mosaic) |
| 232 |
num_cores <- 4
|
|
| 256 |
num_cores <- 1
|
|
| 233 | 257 |
l_png_files <- mclapply(1:length(unlist(lf_mean_mosaic)),FUN=plot_mosaic, |
| 234 | 258 |
list_param= list_param_plot_mosaic, |
| 235 | 259 |
mc.preschedule=FALSE,mc.cores = num_cores) |
| 236 | 260 |
|
| 261 |
lf_plot<- list.files(pattern="r_m_use.*.mask.*.tif$") |
|
| 262 |
lf_mean_mosaic <- lf_plot |
|
| 263 |
|
|
| 264 |
|
|
| 265 |
list_param_plot_mosaic <- list(lf_raster_fname=unlist(lf_mean_mosaic[1:3]), |
|
| 266 |
screenRast=FALSE, |
|
| 267 |
l_dates=day_to_mosaic, |
|
| 268 |
out_dir_str=out_dir, |
|
| 269 |
out_prefix_str <- "dailyTmax_", |
|
| 270 |
out_suffix_str=out_suffix) |
|
| 271 |
#plot_screen_raster_val(3,list_param_plot_mosaic) |
|
| 272 |
#debug(plot_screen_raster_val) |
|
| 273 |
|
|
| 274 |
num_cores <- 3 |
|
| 275 |
l_png_files <- mclapply(1:length(unlist(lf_mean_mosaic)[1:3]),FUN=plot_screen_raster_val, |
|
| 276 |
list_param= list_param_plot_mosaic, |
|
| 277 |
mc.preschedule=FALSE,mc.cores = num_cores) |
|
| 278 |
|
|
| 279 |
|
|
| 280 |
list_param_plot_mosaic <- list(lf_raster_fname=unlist(lf_mean_mosaic[4:6]), |
|
| 281 |
screenRast=FALSE, |
|
| 282 |
l_dates=day_to_mosaic, |
|
| 283 |
out_dir_str=out_dir, |
|
| 284 |
out_prefix_str <- paste("rmse_",out_suffix,sep=""),
|
|
| 285 |
out_suffix_str=paste("rmse_",out_suffix,sep=""))
|
|
| 286 |
num_cores <- 3 |
|
| 287 |
l_png_files_rmse <- mclapply(1:length(unlist(lf_mean_mosaic)[4:6]),FUN=plot_screen_raster_val, |
|
| 288 |
list_param= list_param_plot_mosaic, |
|
| 289 |
mc.preschedule=FALSE,mc.cores = num_cores) |
|
| 237 | 290 |
|
| 238 | 291 |
############### |
| 239 | 292 |
|
Also available in: Unified diff
testing python option on NEX