Revision 35290bf2
Added by Benoit Parmentier over 8 years ago
| climate/research/oregon/interpolation/NASA2016_conference_temperature_predictions.R | ||
|---|---|---|
| 22 | 22 |
|
| 23 | 23 |
################################################################################################# |
| 24 | 24 |
|
| 25 |
#### FUNCTION USED IN SCRIPT |
|
| 26 |
|
|
| 27 |
#function_analyses_paper1 <-"contribution_of_covariates_paper_interpolation_functions_07182014.R" #first interp paper |
|
| 28 |
#function_analyses_paper2 <-"multi_timescales_paper_interpolation_functions_08132014.R" |
|
| 29 |
#function_global_run_assessment_part2 <- "global_run_scalingup_assessment_part2_functions_0923015.R" |
|
| 30 |
|
|
| 31 |
############################################ |
|
| 32 |
#### Parameters and constants |
|
| 33 |
|
|
| 34 | 25 |
|
| 35 | 26 |
### Loading R library and packages |
| 36 | 27 |
#library used in the workflow production: |
| ... | ... | |
| 85 | 76 |
source(file.path(script_path,function_assessment_part2_functions)) #source all functions used in this script |
| 86 | 77 |
source(file.path(script_path,function_assessment_part3)) #source all functions used in this script |
| 87 | 78 |
|
| 88 |
### Parameters, constants and arguments ### |
|
| 79 |
pre_process_raster_mosaic_fun <- function(i,list_param){
|
|
| 80 |
|
|
| 81 |
|
|
| 82 |
## Extract parameters |
|
| 83 |
|
|
| 84 |
lf <- list_param$lf |
|
| 85 |
python_bin <- list_param$python_bin |
|
| 86 |
infile_mask <- list_param$infile_mask |
|
| 87 |
scaling <- list_param$scaling |
|
| 88 |
mask_pred <- list_param$mask_pred |
|
| 89 |
NA_flag_val <- list_param$NA_flag_val |
|
| 90 |
out_suffix <- list_param$out_suffix |
|
| 91 |
out_dir <- list_param$out_dir |
|
| 92 |
|
|
| 93 |
raster_name_in <- lf[i] |
|
| 94 |
|
|
| 95 |
#Step 1: match extent and resolution |
|
| 96 |
|
|
| 97 |
lf_files <- c(raster_name_in) #match to mask |
|
| 98 |
rast_ref <- infile_mask |
|
| 99 |
##Maching resolution is probably only necessary for the r mosaic function |
|
| 100 |
#Modify later to take into account option R or python... |
|
| 101 |
list_param_raster_match <- list(lf_files,rast_ref,file_format,python_bin,out_suffix,out_dir) |
|
| 102 |
names(list_param_raster_match) <- c("lf_files","rast_ref","file_format","python_bin","out_suffix","out_dir_str")
|
|
| 103 |
r_pred_matched <- raster_match(1,list_param_raster_match) |
|
| 104 |
raster_name_in <- c(r_pred_matched) |
|
| 105 |
|
|
| 106 |
#Step 2: mask |
|
| 107 |
if(mask_pred==TRUE){
|
|
| 108 |
r_mask <- raster(infile_mask) |
|
| 109 |
extension_str <- extension(raster_name_in) |
|
| 110 |
raster_name_tmp <- gsub(extension_str,"",basename(raster_name_in)) |
|
| 111 |
raster_name <- file.path(out_dir,paste(raster_name_tmp,"_masked.tif",sep = "")) |
|
| 112 |
r_pred <- mask(raster(raster_name_in),r_mask,filename = raster_name,overwrite = TRUE) |
|
| 113 |
} |
|
| 114 |
|
|
| 115 |
NAvalue(r_pred) <- NA_flag_val |
|
| 116 |
#r_pred <- setMinMax(r_pred) |
|
| 117 |
|
|
| 118 |
#Step 3: remove scaling factor |
|
| 119 |
raster_name_in <- filename(r_pred) |
|
| 120 |
extension_str <- extension(raster_name_in) |
|
| 121 |
raster_name_tmp <- gsub(extension_str,"",basename(filename(r_pred))) |
|
| 122 |
raster_name_out <- file.path(out_dir,paste(raster_name_tmp,"_rescaled.tif",sep = "")) |
|
| 123 |
#r_pred <- overlay(r_pred, fun=function(x){x*1/scaling},filename=raster_name,overwrite=T)
|
|
| 124 |
|
|
| 125 |
#raster_name_in <- "comp_r_m_use_edge_weights_weighted_mean_gam_CAI_dailyTmax_19990102_reg4_1999_m_gam_CAI_dailyTmax_19990102_reg4_1999_m__meeting_NASA_reg4_04292016_masked.tif" |
|
| 126 |
python_cmd <- file.path(python_bin,"gdal_calc.py") |
|
| 127 |
cmd_str3 <- paste(python_cmd, |
|
| 128 |
paste("-A ", raster_name_in,sep=""),
|
|
| 129 |
paste("--outfile=",raster_name_out,sep=""),
|
|
| 130 |
#paste("--type=","Int32",sep=""),
|
|
| 131 |
"--co='COMPRESS=LZW'", |
|
| 132 |
paste("--NoDataValue=",NA_flag_val,sep=""),
|
|
| 133 |
paste("--calc='(A)*",scaling,"'",sep=""),
|
|
| 134 |
"--overwrite",sep=" ") #division by zero is problematic... |
|
| 135 |
#system(cmd_str3) |
|
| 136 |
system(cmd_str3) |
|
| 137 |
#NAvalue(r_pred) <- NA_flag_val |
|
| 138 |
#r_pred <- |
|
| 139 |
r_pred <- setMinMax(raster(raster_name_out)) |
|
| 140 |
|
|
| 141 |
return(raster_name_out) |
|
| 142 |
} |
|
| 143 |
|
|
| 144 |
############################### |
|
| 145 |
####### Parameters, constants and arguments ### |
|
| 89 | 146 |
|
| 90 | 147 |
CRS_locs_WGS84<-CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +towgs84=0,0,0") #constant 1
|
| 91 | 148 |
|
| ... | ... | |
| 101 | 158 |
y_var_month <- "TMin" |
| 102 | 159 |
} |
| 103 | 160 |
|
| 104 |
|
|
| 105 | 161 |
##Add for precip later... |
| 106 | 162 |
if (var == "TMAX") {
|
| 107 | 163 |
variable_name <- "maximum temperature" |
| ... | ... | |
| 157 | 213 |
|
| 158 | 214 |
df_centroids_fname <- "/data/project/layers/commons/NEX_data/climateLayers/out/reg4/mosaic/output_reg4_1999/df_centroids_19990701_reg4_1999.txt" |
| 159 | 215 |
|
| 160 |
raster_name_lf <- c("/data/project/layers/commons/NEX_data/climateLayers/out/reg4/mosaic/int_mosaics/comp_r_m_use_edge_weights_weighted_mean_gam_CAI_dailyTmax_19990101_reg4_1999_m_gam_CAI_dailyTmax_19990101_reg4_1999.tif",
|
|
| 161 |
"/data/project/layers/commons/NEX_data/climateLayers/out/reg4/mosaic/int_mosaics/comp_r_m_use_edge_weights_weighted_mean_gam_CAI_dailyTmax_19990102_reg4_1999_m_gam_CAI_dailyTmax_19990102_reg4_1999.tif", |
|
| 162 |
"/data/project/layers/commons/NEX_data/climateLayers/out/reg4/mosaic/int_mosaics/comp_r_m_use_edge_weights_weighted_mean_gam_CAI_dailyTmax_19990103_reg4_1999_m_gam_CAI_dailyTmax_19990103_reg4_1999.tif", |
|
| 163 |
"/data/project/layers/commons/NEX_data/climateLayers/out/reg4/mosaic/int_mosaics/comp_r_m_use_edge_weights_weighted_mean_gam_CAI_dailyTmax_19990701_reg4_1999_m_gam_CAI_dailyTmax_19990701_reg4_1999.tif", |
|
| 164 |
"/data/project/layers/commons/NEX_data/climateLayers/out/reg4/mosaic/int_mosaics/comp_r_m_use_edge_weights_weighted_mean_gam_CAI_dailyTmax_19990702_reg4_1999_m_gam_CAI_dailyTmax_19990702_reg4_1999.tif", |
|
| 165 |
"/data/project/layers/commons/NEX_data/climateLayers/out/reg4/mosaic/int_mosaics/comp_r_m_use_edge_weights_weighted_mean_gam_CAI_dailyTmax_19990703_reg4_1999_m_gam_CAI_dailyTmax_19990703_reg4_1999.tif") |
|
| 216 |
#raster_name_lf <- c("/data/project/layers/commons/NEX_data/climateLayers/out/reg4/mosaic/int_mosaics/comp_r_m_use_edge_weights_weighted_mean_gam_CAI_dailyTmax_19990101_reg4_1999_m_gam_CAI_dailyTmax_19990101_reg4_1999.tif",
|
|
| 217 |
# "/data/project/layers/commons/NEX_data/climateLayers/out/reg4/mosaic/int_mosaics/comp_r_m_use_edge_weights_weighted_mean_gam_CAI_dailyTmax_19990102_reg4_1999_m_gam_CAI_dailyTmax_19990102_reg4_1999.tif", |
|
| 218 |
# "/data/project/layers/commons/NEX_data/climateLayers/out/reg4/mosaic/int_mosaics/comp_r_m_use_edge_weights_weighted_mean_gam_CAI_dailyTmax_19990103_reg4_1999_m_gam_CAI_dailyTmax_19990103_reg4_1999.tif", |
|
| 219 |
# "/data/project/layers/commons/NEX_data/climateLayers/out/reg4/mosaic/int_mosaics/comp_r_m_use_edge_weights_weighted_mean_gam_CAI_dailyTmax_19990701_reg4_1999_m_gam_CAI_dailyTmax_19990701_reg4_1999.tif", |
|
| 220 |
# "/data/project/layers/commons/NEX_data/climateLayers/out/reg4/mosaic/int_mosaics/comp_r_m_use_edge_weights_weighted_mean_gam_CAI_dailyTmax_19990702_reg4_1999_m_gam_CAI_dailyTmax_19990702_reg4_1999.tif", |
|
| 221 |
# "/data/project/layers/commons/NEX_data/climateLayers/out/reg4/mosaic/int_mosaics/comp_r_m_use_edge_weights_weighted_mean_gam_CAI_dailyTmax_19990703_reg4_1999_m_gam_CAI_dailyTmax_19990703_reg4_1999.tif") |
|
| 222 |
|
|
| 223 |
raster_name_lf <- c("/data/project/layers/commons/NEX_data/climateLayers/out/reg4/mosaic/int_mosaics/comp_r_m_use_edge_weights_weighted_mean_gam_CAI_dailyTmax_19920101_reg4_1992_m_gam_CAI_dailyTmax_19920101_reg4_1992.tif",
|
|
| 224 |
"/data/project/layers/commons/NEX_data/climateLayers/out/reg4/mosaic/int_mosaics/comp_r_m_use_edge_weights_weighted_mean_gam_CAI_dailyTmax_19920102_reg4_1992_m_gam_CAI_dailyTmax_19920102_reg4_1992.tif", |
|
| 225 |
"/data/project/layers/commons/NEX_data/climateLayers/out/reg4/mosaic/int_mosaics/comp_r_m_use_edge_weights_weighted_mean_gam_CAI_dailyTmax_19920103_reg4_1992_m_gam_CAI_dailyTmax_19920103_reg4_1992.tif", |
|
| 226 |
"/data/project/layers/commons/NEX_data/climateLayers/out/reg4/mosaic/int_mosaics/comp_r_m_use_edge_weights_weighted_mean_gam_CAI_dailyTmax_19920701_reg4_1992_m_gam_CAI_dailyTmax_19920701_reg4_1992.tif", |
|
| 227 |
"/data/project/layers/commons/NEX_data/climateLayers/out/reg4/mosaic/int_mosaics/comp_r_m_use_edge_weights_weighted_mean_gam_CAI_dailyTmax_19920702_reg4_1992_m_gam_CAI_dailyTmax_19920702_reg4_1992.tif", |
|
| 228 |
"/data/project/layers/commons/NEX_data/climateLayers/out/reg4/mosaic/int_mosaics/comp_r_m_use_edge_weights_weighted_mean_gam_CAI_dailyTmax_19920703_reg4_1992_m_gam_CAI_dailyTmax_19920703_reg4_1992.tif") |
|
| 229 |
|
|
| 230 |
#l_dates <- c("19990101","19990102","19990103","19990701","19990702","19990703")
|
|
| 231 |
l_dates <- c("19920101","19920102","19920103","19920701","19920702","19990703")
|
|
| 166 | 232 |
|
| 167 |
l_dates <- c("19990101","19990102","19990103","19990701","19990702","19990703")
|
|
| 168 | 233 |
|
| 169 | 234 |
##################### START SCRIPT ################# |
| 170 | 235 |
|
| ... | ... | |
| 189 | 254 |
|
| 190 | 255 |
#date_to_plot <- seq(as.Date(strptime(start_date,"%Y%m%d")), as.Date(strptime(end_date,"%Y%m%d")), 'day') |
| 191 | 256 |
#l_dates <- format(date_to_plot,"%Y%m%d") #format back to the relevant date format for files |
| 192 |
|
|
| 193 |
|
|
| 194 |
|
|
| 195 | 257 |
mask_pred <- TRUE |
| 196 | 258 |
list_param_pre_process <- list(raster_name_lf,python_bin,infile_mask,scaling,mask_pred,NA_flag_val,out_suffix,out_dir) |
| 197 | 259 |
names(list_param_pre_process) <- c("lf","python_bin","infile_mask","scaling","mask_pred","NA_flag_val","out_suffix","out_dir")
|
| ... | ... | |
| 201 | 263 |
#lf_mosaic_scaled <- mclapply(1:length(raster_name_lf),FUN=pre_process_raster_mosaic_fun,list_param=list_param_pre_process,mc.preschedule=FALSE,mc.cores = num_cores) |
| 202 | 264 |
lf_mosaic_scaled <- mclapply(1:length(raster_name_lf),FUN=pre_process_raster_mosaic_fun,list_param=list_param_pre_process,mc.preschedule=FALSE,mc.cores = num_cores) |
| 203 | 265 |
|
| 204 |
test <- pre_process_raster_mosaic_fun(2,list_param_pre_process) |
|
| 266 |
#test <- pre_process_raster_mosaic_fun(2,list_param_pre_process) |
|
| 267 |
lf_mosaic_scaled <- unlist(lf_mosaic_scaled) |
|
| 205 | 268 |
|
| 206 |
pre_process_raster_mosaic_fun <- function(i,list_param){
|
|
| 207 |
|
|
| 208 |
|
|
| 209 |
## Extract parameters |
|
| 210 |
|
|
| 211 |
lf <- list_param$lf |
|
| 212 |
python_bin <- list_param$python_bin |
|
| 213 |
infile_mask <- list_param$infile_mask |
|
| 214 |
scaling <- list_param$scaling |
|
| 215 |
mask_pred <- list_param$mask_pred |
|
| 216 |
NA_flag_val <- list_param$NA_flag_val |
|
| 217 |
out_suffix <- list_param$out_suffix |
|
| 218 |
out_dir <- list_param$out_dir |
|
| 219 |
|
|
| 220 |
raster_name_in <- lf[i] |
|
| 221 |
|
|
| 222 |
#Step 1: match extent and resolution |
|
| 223 |
|
|
| 224 |
lf_files <- c(raster_name_in) #match to mask |
|
| 225 |
rast_ref <- infile_mask |
|
| 226 |
##Maching resolution is probably only necessary for the r mosaic function |
|
| 227 |
#Modify later to take into account option R or python... |
|
| 228 |
list_param_raster_match <- list(lf_files,rast_ref,file_format,python_bin,out_suffix,out_dir) |
|
| 229 |
names(list_param_raster_match) <- c("lf_files","rast_ref","file_format","python_bin","out_suffix","out_dir_str")
|
|
| 230 |
r_pred_matched <- raster_match(1,list_param_raster_match) |
|
| 231 |
raster_name_in <- c(r_pred_matched) |
|
| 232 |
|
|
| 233 |
#Step 2: mask |
|
| 234 |
if(mask_pred==TRUE){
|
|
| 235 |
r_mask <- raster(infile_mask) |
|
| 236 |
extension_str <- extension(raster_name_in) |
|
| 237 |
raster_name_tmp <- gsub(extension_str,"",basename(raster_name_in)) |
|
| 238 |
raster_name <- file.path(out_dir,paste(raster_name_tmp,"_masked.tif",sep = "")) |
|
| 239 |
r_pred <- mask(raster(raster_name_in),r_mask,filename = raster_name,overwrite = TRUE) |
|
| 240 |
} |
|
| 241 |
|
|
| 242 |
NAvalue(r_pred) <- NA_flag_val |
|
| 243 |
#r_pred <- setMinMax(r_pred) |
|
| 244 |
|
|
| 245 |
#Step 3: remove scaling factor |
|
| 246 |
raster_name_in <- filename(r_pred) |
|
| 247 |
extension_str <- extension(raster_name_in) |
|
| 248 |
raster_name_tmp <- gsub(extension_str,"",basename(filename(r_pred))) |
|
| 249 |
raster_name_out <- file.path(out_dir,paste(raster_name_tmp,"_rescaled.tif",sep = "")) |
|
| 250 |
#r_pred <- overlay(r_pred, fun=function(x){x*1/scaling},filename=raster_name,overwrite=T)
|
|
| 251 |
|
|
| 252 |
#raster_name_in <- "comp_r_m_use_edge_weights_weighted_mean_gam_CAI_dailyTmax_19990102_reg4_1999_m_gam_CAI_dailyTmax_19990102_reg4_1999_m__meeting_NASA_reg4_04292016_masked.tif" |
|
| 253 |
python_cmd <- file.path(python_bin,"gdal_calc.py") |
|
| 254 |
cmd_str3 <- paste(python_cmd, |
|
| 255 |
paste("-A ", raster_name_in,sep=""),
|
|
| 256 |
paste("--outfile=",raster_name_out,sep=""),
|
|
| 257 |
#paste("--type=","Int32",sep=""),
|
|
| 258 |
"--co='COMPRESS=LZW'", |
|
| 259 |
paste("--NoDataValue=",NA_flag_val,sep=""),
|
|
| 260 |
paste("--calc='(A)*",scaling,"'",sep=""),
|
|
| 261 |
"--overwrite",sep=" ") #division by zero is problematic... |
|
| 262 |
#system(cmd_str3) |
|
| 263 |
system(cmd_str3) |
|
| 264 |
#NAvalue(r_pred) <- NA_flag_val |
|
| 265 |
#r_pred <- |
|
| 266 |
r_pred <- setMinMax(raster(raster_name_out)) |
|
| 267 |
|
|
| 268 |
return(raster_name_out) |
|
| 269 |
} |
|
| 269 |
r_mosaic_scaled <- stack(lf_mosaic_scaled) |
|
| 270 |
NAvalue(r_mosaic_scaled)<- -3399999901438340239948148078125514752.000 |
|
| 271 |
plot(r_mosaic_scaled,y=6,zlim=c(-50,50)) |
|
| 272 |
plot(r_mosaic_scaled,zlim=c(-50,50),col=matlab.like(255)) |
|
| 270 | 273 |
|
| 271 | 274 |
month_name <- month.name() |
| 272 | 275 |
l_dates <- as.Date(strptime(date_proc,"%Y%m%d")) |
| ... | ... | |
| 275 | 278 |
#s.range <- s.range+c(5,-5) |
| 276 | 279 |
#col.breaks <- pretty(s.range, n=200) |
| 277 | 280 |
#lab.breaks <- pretty(s.range, n=100) |
| 278 |
temp.colors <- colorRampPalette(c('blue', 'white', 'red'))
|
|
| 279 |
max_val<-s.range[2] |
|
| 280 |
min_val <-s.range[1] |
|
| 281 |
#temp.colors <- colorRampPalette(c('blue', 'white', 'red'))
|
|
| 282 |
#max_val<-s.range[2]
|
|
| 283 |
#min_val <-s.range[1]
|
|
| 281 | 284 |
#max_val<- -10 |
| 282 |
min_val <- 0 |
|
| 285 |
#min_val <- 0
|
|
| 283 | 286 |
|
| 284 | 287 |
#layout_m<-c(1,3) #one row two columns |
| 285 | 288 |
date_proc <- l_dates[i] |
| 289 |
levelplot(r_mosaic_scaled,zlim=c(-50,50),col.regions=matlab.like(255)) |
|
| 290 |
#levelplot(r_mosaic_scaled,zlim=c(-50,50),col.regions=matlab.like(255)) |
|
| 286 | 291 |
|
| 287 | 292 |
#png(paste("Figure7a__spatial_pattern_tmax_prediction_levelplot_",date_selected,out_prefix,".png", sep=""),
|
| 288 | 293 |
# height=480*layout_m[1],width=480*layout_m[2]) |
| 289 | 294 |
#plot(r_pred,col=temp.colors(255),zlim=c(-3500,4500)) |
| 290 |
plot(r_pred,col=matlab.like(255),zlim=c(-35,45))
|
|
| 295 |
#plot(r_pred,col=matlab.like(255),zlim=c(-40,50))
|
|
| 291 | 296 |
#paste(raster_name[1:7],collapse="_") |
| 292 | 297 |
#add filename option later |
| 293 | 298 |
|
| 299 |
for (i in 1:length(nlayers(r_mosaic_scaled))){
|
|
| 300 |
|
|
| 301 |
date_proc <- l_dates[i] |
|
| 302 |
r_pred <- subset(r_mosaic_scaled,i) |
|
| 303 |
|
|
| 304 |
res_pix <- 1200 |
|
| 305 |
#res_pix <- 480 |
|
| 306 |
col_mfrow <- 1 |
|
| 307 |
row_mfrow <- 1 |
|
| 308 |
|
|
| 309 |
png_filename <- file.path(out_dir_str,paste("Figure4_clim_mosaics_day_","_",date_proc,"_",reg_name,"_",out_suffix,".png",sep =""))
|
|
| 310 |
title_str <- paste("Predicted ",variable_name, " on ",date_proc , " ", sep = "")
|
|
| 311 |
|
|
| 312 |
png(filename=png_filename,width = col_mfrow * res_pix,height = row_mfrow * res_pix) |
|
| 313 |
plot(r_pred,main =title_str,cex.main =1.5,col=matlab.like(255)) |
|
| 314 |
dev.off() |
|
| 294 | 315 |
|
| 295 |
|
|
| 296 |
res_pix <- 1200 |
|
| 297 |
#res_pix <- 480 |
|
| 298 |
|
|
| 299 |
col_mfrow <- 1 |
|
| 300 |
row_mfrow <- 1 |
|
| 301 |
date_proc <- |
|
| 302 |
png_filename <- file.path(out_dir_str,paste("Figure4_clim_mosaics_day_","_",date_proc,"_",reg_name,"_",out_suffix,".png",sep =""))
|
|
| 303 |
|
|
| 304 |
png(filename=png_filename,width = col_mfrow * res_pix,height = row_mfrow * res_pix) |
|
| 305 |
|
|
| 306 |
plot(r_pred,main = paste("Predicted ",variable_name, " on ",date_proc , " ", ,sep = ""),cex.main =1.5)
|
|
| 307 |
|
|
| 308 |
dev.off() |
|
| 309 |
|
|
| 310 |
#col.regions=temp.colors(25)) |
|
| 311 |
#dev.off() |
|
| 316 |
} |
|
| 312 | 317 |
|
| 313 | 318 |
#### Extract time series |
| 314 | 319 |
|
| 315 | 320 |
#-65,-22 |
| 316 | 321 |
|
| 317 | 322 |
df_centroids <- read.table(df_centroids_fname,sep=",") |
| 323 |
coordinates(df_centroids)<- c("long","lat")
|
|
| 324 |
proj4string(df_centroids) <- CRS_locs_WGS84 |
|
| 325 |
|
|
| 326 |
extract(df_centroids,) |
|
| 318 | 327 |
|
| 319 | 328 |
raster() |
| 329 |
|
|
| 320 | 330 |
############################ END OF SCRIPT ################################## |
Also available in: Unified diff
moving pre process mosaic function and generating 6 mosaic figures for poster NASA conference