Revision 0b599be0
Added by Benoit Parmentier over 8 years ago
climate/research/oregon/interpolation/NASA2016_conference_temperature_predictions.R | ||
---|---|---|
5 | 5 |
#Figures and data for the AAG conference are also produced. |
6 | 6 |
#AUTHOR: Benoit Parmentier |
7 | 7 |
#CREATED ON: 05/01/2016 |
8 |
#MODIFIED ON: 05/02/2016
|
|
8 |
#MODIFIED ON: 05/03/2016
|
|
9 | 9 |
#Version: 1 |
10 | 10 |
#PROJECT: Environmental Layers project |
11 | 11 |
#COMMENTS: Initial commit, script based on part 2 of assessment, will be modified further for overall assessment |
... | ... | |
51 | 51 |
library(colorRamps) |
52 | 52 |
library(zoo) |
53 | 53 |
library(xts) |
54 |
library(lubridate) |
|
54 | 55 |
|
55 | 56 |
###### Function used in the script ####### |
56 | 57 |
|
57 | 58 |
#script_path <- "/nobackupp8/bparmen1/env_layers_scripts" #path to script |
58 | 59 |
script_path <- "/home/parmentier/Data/IPLANT_project/env_layers_scripts" #path to script |
59 | 60 |
|
60 |
#Mosaic related |
|
61 |
## NASA poster and paper related |
|
62 |
source(file.path(script_path,"NASA2016_conference_temperature_predictions_function_05032016b.R")) |
|
63 |
|
|
64 |
#Mosaic related on NEX |
|
61 | 65 |
#script_path <- "/home/parmentier/Data/IPLANT_project/env_layers_scripts" |
62 | 66 |
function_mosaicing_functions <- "global_run_scalingup_mosaicing_function_04232016.R" #PARAM12 |
63 | 67 |
function_mosaicing <-"global_run_scalingup_mosaicing_05012016.R" |
64 | 68 |
source(file.path(script_path,function_mosaicing)) #source all functions used in this script |
65 | 69 |
source(file.path(script_path,function_mosaicing_functions)) #source all functions used in this script |
66 | 70 |
|
67 |
#Assessment |
|
71 |
#Assessment on NEX
|
|
68 | 72 |
function_assessment_part1_functions <- "global_run_scalingup_assessment_part1_functions_02112015.R" #PARAM12 |
69 | 73 |
function_assessment_part1a <-"global_run_scalingup_assessment_part1a_01042016.R" |
70 | 74 |
function_assessment_part2 <- "global_run_scalingup_assessment_part2_02092016.R" |
... | ... | |
76 | 80 |
source(file.path(script_path,function_assessment_part2_functions)) #source all functions used in this script |
77 | 81 |
source(file.path(script_path,function_assessment_part3)) #source all functions used in this script |
78 | 82 |
|
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 |
plot_raster_mosaic <- function(i,list_param){ |
|
145 |
#Function to plot mosaic for poster |
|
146 |
# |
|
147 |
l_dates <- list_param$l_dates |
|
148 |
r_mosaiced_scaled <- list_param$r_mosaiced_scaled |
|
149 |
NA_flag_val <- list_param$NA_flag_val |
|
150 |
out_dir <- list_param$out_dir |
|
151 |
out_suffix <- list_param$out_suffix |
|
152 |
region_name <- list_param$region_name |
|
153 |
variable_name <- list_param$variable_name |
|
154 |
|
|
155 |
#for (i in 1:length(nlayers(r_mosaic_scaled))){ |
|
156 |
|
|
157 |
date_proc <- l_dates[i] |
|
158 |
r_pred <- subset(r_mosaic_scaled,i) |
|
159 |
NAvalue(r_pred)<- NA_flag_val |
|
160 |
|
|
161 |
date_proc <- l_dates[i] |
|
162 |
date_val <- as.Date(strptime(date_proc,"%Y%m%d")) |
|
163 |
#month_name <- month.name(date_val) |
|
164 |
month_str <- format(date_val, "%b") ## Month, char, abbreviated |
|
165 |
year_str <- format(date_val, "%Y") ## Year with century |
|
166 |
day_str <- as.numeric(format(date_val, "%d")) ## numeric month |
|
167 |
date_str <- paste(month_str," ",day_str,", ",year_str,sep="") |
|
168 |
|
|
169 |
res_pix <- 1200 |
|
170 |
#res_pix <- 480 |
|
171 |
col_mfrow <- 1 |
|
172 |
row_mfrow <- 1 |
|
173 |
|
|
174 |
png_filename <- file.path(out_dir,paste("Figure4_clim_mosaics_day_","_",date_proc,"_",region_name,"_",out_suffix,".png",sep ="")) |
|
175 |
title_str <- paste("Predicted ",variable_name, " on ",date_str , " ", sep = "") |
|
176 |
|
|
177 |
png(filename=png_filename,width = col_mfrow * res_pix,height = row_mfrow * res_pix) |
|
178 |
plot(r_pred,main =title_str,cex.main =1.5,col=matlab.like(255),zlim=c(-50,50), |
|
179 |
legend.shrink=0.8,legend.width=0.8) |
|
180 |
#axis.args = list(cex.axis = 1.6), #control size of legend z |
|
181 |
#legend.args=list(text='dNBR', side=4, line=2.5, cex=2.2)) |
|
182 |
#legend.args=list(text='dNBR', side=4, line=2.49, cex=1.6)) |
|
183 |
dev.off() |
|
184 |
|
|
185 |
return(png_filename) |
|
186 |
} |
|
187 |
|
|
188 |
extract_date <- function(i,x,item_no=NULL){ |
|
189 |
y <- unlist(strsplit(x[[i]],"_")) |
|
190 |
if(is.null(item_no)){ |
|
191 |
date_str <- y[length(y)-2] #count from end |
|
192 |
}else{ |
|
193 |
date_str <- y[item_no] |
|
194 |
} |
|
195 |
return(date_str) |
|
196 |
} |
|
197 |
|
|
198 | 83 |
############################### |
199 | 84 |
####### Parameters, constants and arguments ### |
200 | 85 |
|
... | ... | |
252 | 137 |
#python_bin <- "/nobackupp6/aguzman4/climateLayers/sharedModules2/bin" #PARAM 30 |
253 | 138 |
python_bin <- "/usr/bin" #PARAM 30 |
254 | 139 |
|
255 |
day_start <- "19840101" #PARAM 12 arg 12
|
|
256 |
day_end <- "20021231" #PARAM 13 arg 13
|
|
140 |
day_start <- "1986101" #PARAM 12 arg 12
|
|
141 |
day_end <- "19981231" #PARAM 13 arg 13
|
|
257 | 142 |
|
258 | 143 |
#infile_mask <- "/nobackupp8/bparmen1/NEX_data/regions_input_files/r_mask_LST_reg4.tif" |
259 | 144 |
infile_mask <- "/data/project/layers/commons/NEX_data/regions_input_files/r_mask_LST_reg4.tif" |
... | ... | |
275 | 160 |
l_dates <- c("19920101","19920102","19920103","19920701","19920702","19990703") |
276 | 161 |
|
277 | 162 |
df_points_extracted_fname <- "/data/project/layers/commons/NEX_data/climateLayers/out/reg4/mosaic/int_mosaics/data_points_extracted.txt" |
163 |
NA_flag_val_mosaic <- -3399999901438340239948148078125514752.000 |
|
278 | 164 |
|
279 | 165 |
##################### START SCRIPT ################# |
280 | 166 |
|
... | ... | |
326 | 212 |
#paste(raster_name[1:7],collapse="_") |
327 | 213 |
#add filename option later |
328 | 214 |
|
329 |
NA_flag_val_mosaic <- -3399999901438340239948148078125514752.000 |
|
215 |
#NA_flag_val_mosaic <- -3399999901438340239948148078125514752.000
|
|
330 | 216 |
|
331 | 217 |
list_param_plot_raster_mosaic <- list(l_dates,r_mosaic_scaled,NA_flag_val_mosaic,out_dir,out_suffix, |
332 | 218 |
region_name,variable_name) |
... | ... | |
335 | 221 |
|
336 | 222 |
lf_mosaic_plot_fig <- mclapply(1:length(lf_mosaic_scaled),FUN=plot_raster_mosaic,list_param=list_param_plot_raster_mosaic,mc.preschedule=FALSE,mc.cores = num_cores) |
337 | 223 |
|
338 |
|
|
339 |
|
|
340 | 224 |
############### PART2: temporal profile ############# |
341 | 225 |
#### Extract time series |
342 | 226 |
### |
... | ... | |
379 | 263 |
|
380 | 264 |
station_id <- 8 |
381 | 265 |
var_name <-paste0("ID_",station_id) |
382 |
#aggregate(sdf_tmp |
|
266 |
|
|
267 |
##Screen for unique date values |
|
383 | 268 |
if(max(unique_date_tb)>1){ |
384 | 269 |
# formula_str <- paste(var_name," ~ ","TRIP_START_DATE_f",sep="") |
385 | 270 |
var_pix <- aggregate(ID_8 ~ date, data = df_points, mean) #aggregate by date |
386 | 271 |
} |
387 | 272 |
|
388 |
|
|
389 |
#start_date <-input$dates[1] |
|
390 |
#end_date <-input$dates[2] |
|
273 |
var_pix$ID_8 <- var_pix$ID_8*scaling |
|
391 | 274 |
|
392 |
d_z_tmp <- zoo(df_points[,station_id],df_points$date) |
|
275 |
d_z_tmp <- zoo(var_pix$ID_8,var_pix$date) |
|
276 |
names(d_z_tmp)<-"ID_8" |
|
277 |
min(d_z_tmp$ID_8) |
|
278 |
max(d_z_tmp$ID_8) |
|
393 | 279 |
|
280 |
plot(d_z_tmp) |
|
281 |
|
|
282 |
day_start <- "1986-01-01" #PARAM 12 arg 12 |
|
283 |
day_end <- "1998-12-31" #PARAM 13 arg 13 |
|
284 |
|
|
285 |
start_date <- as.Date(day_start) |
|
286 |
end_date <- as.Date(day_end) |
|
287 |
start_year <- year(start_date) |
|
288 |
end_year <- year(end_date) |
|
289 |
|
|
290 |
d_z <- window(d_z_tmp,start=start_date,end=end_date) |
|
291 |
#d_z2 <- window(d_z_tmp2,start=start_date,end=end_date) |
|
292 |
|
|
293 |
title_str <- paste("Predicted daily ",variable_name," for the ", start_year,"-",end_year," time period",sep="") |
|
294 |
plot(d_z,ylab="tmax in deg C",xlab="daily time steps", |
|
295 |
main=title_str, |
|
296 |
lty=3) |
|
394 | 297 |
|
395 |
d_z <- window(d_z_tmp,start=start_date,end=end_date) |
|
396 | 298 |
|
397 |
|
|
398 | 299 |
#data_pixel <- data_df[id_selected,] |
399 | 300 |
#data_pixel$rainfall <- as.numeric(data_pixel$rainfall) |
400 | 301 |
#d_z_tmp <-zoo(data_pixel$rainfall,as.Date(data_pixel$date)) |
... | ... | |
402 | 303 |
#data_pixel <- as.data.frame(data_pixel) |
403 | 304 |
#d_z_tmp2 <- zoo(data_pixel[[var_name]],as.Date(data_pixel$date)) |
404 | 305 |
|
405 |
#start_date <-input$dates[1] |
|
406 |
#end_date <-input$dates[2] |
|
407 |
|
|
408 |
#d_z <- window(d_z_tmp,start=start_date,end=end_date) |
|
409 |
#d_z2 <- window(d_z_tmp2,start=start_date,end=end_date) |
|
410 |
#df2 <- as.data.frame(d_z2) |
|
411 |
#names(df2)<- var_name |
|
412 |
|
|
413 | 306 |
#df_tmp <- subset(data_var,data_var$ID_stat==id_name) |
414 | 307 |
#if(da) |
415 | 308 |
|
Also available in: Unified diff
moving function to script and fixing temporal profiles extracted