Revision 03a2bc79
Added by Benoit Parmentier about 8 years ago
climate/research/oregon/interpolation/global_product_assessment_part2.R | ||
---|---|---|
149 | 149 |
python_bin <- "/usr/bin" #PARAM 30 |
150 | 150 |
|
151 | 151 |
day_start <- "1984101" #PARAM 12 arg 12 |
152 |
day_end <- "19991231" #PARAM 13 arg 13
|
|
152 |
day_end <- "20141231" #PARAM 13 arg 13
|
|
153 | 153 |
#date_start <- day_start |
154 | 154 |
#date_end <- day_end |
155 | 155 |
|
... | ... | |
180 | 180 |
in_dir_list_filename <- NULL #if NULL, use the in_dir directory to search for info |
181 | 181 |
countries_shp <-"/data/project/layers/commons/NEX_data/countries.shp" #Atlas |
182 | 182 |
lf_raster <- NULL #list of raster to consider |
183 |
item_no <- 13 |
|
183 | 184 |
|
184 | 185 |
##################### START SCRIPT ################# |
185 | 186 |
|
... | ... | |
210 | 211 |
r_stack <- stack(lf_raster,quick=T) #this is very fast now with the quick option! |
211 | 212 |
} |
212 | 213 |
|
214 |
#### check for missing dates from list of tif |
|
215 |
|
|
216 |
###This should be a function!! |
|
217 |
|
|
218 |
##### This could be moved in a separate file!! |
|
219 |
############### PART4: Checking for mosaic produced for given region ############## |
|
220 |
## From list of mosaic files predicted extract dates |
|
221 |
## Check dates predicted against date range for a given date range |
|
222 |
## Join file information to centroids of tiles data.frame |
|
223 |
#list_dates_produced <- unlist(mclapply(1:length(lf_mosaic_list),FUN=extract_date,x=lf_mosaic_list,item_no=13,mc.preschedule=FALSE,mc.cores = num_cores)) |
|
224 |
#list_dates_produced <- mclapply(1:2,FUN=extract_date,x=lf_mosaic_list,item_no=13,mc.preschedule=FALSE,mc.cores = 2) |
|
225 |
item_no <- 13 |
|
226 |
date_start <- day_start |
|
227 |
date_end <- day_end |
|
228 |
#day_start <- "1984101" #PARAM 12 arg 12 |
|
229 |
#day_end <- "20141231" #PARAM 13 arg 13 |
|
230 |
|
|
231 |
check_missing <- function(lf, pattern_str=NULL,in_dir=".",date_start="1984101",date_end="20141231",item_no=13,out_suffix=""){ |
|
232 |
#Function to check for missing files such as mosaics or predictions for tiles etc. |
|
233 |
#The function assumes the name of the files contain "_". |
|
234 |
#INPUTS: |
|
235 |
#lf |
|
236 |
#pattern_str |
|
237 |
#in_dir |
|
238 |
#date_start |
|
239 |
#date_end |
|
240 |
#item_no |
|
241 |
#out_suffix |
|
242 |
#OUTPUTS |
|
243 |
# |
|
244 |
# |
|
245 |
|
|
246 |
##### Start script ##### |
|
247 |
|
|
248 |
out_dir <- in_dir |
|
249 |
|
|
250 |
list_dates_produced <- unlist(mclapply(1:length(lf), |
|
251 |
FUN = extract_date, |
|
252 |
x = lf, |
|
253 |
item_no = item_no, |
|
254 |
mc.preschedule = FALSE, |
|
255 |
mc.cores = num_cores)) |
|
256 |
|
|
257 |
list_dates_produced_date_val <- as.Date(strptime(list_dates_produced, "%Y%m%d")) |
|
258 |
month_str <- format(list_dates_produced_date_val, "%b") ## Month, char, abbreviated |
|
259 |
year_str <- format(list_dates_produced_date_val, "%Y") ## Year with century |
|
260 |
day_str <- as.numeric(format(list_dates_produced_date_val, "%d")) ## numeric month |
|
261 |
df_files <- data.frame(lf = basename(lf), |
|
262 |
date = list_dates_produced_date_val, |
|
263 |
month_str = month_str, |
|
264 |
year = year_str, |
|
265 |
day = day_str, |
|
266 |
dir = dirname(lf)) |
|
267 |
|
|
268 |
df_files_fname <- file.path(out_dir, paste0("df_files_", out_suffix, ".txt")) |
|
269 |
write.table(df_files,file = df_files_fname,sep = ",",row.names = F) |
|
270 |
|
|
271 |
#undebug(finding_missing_dates ) |
|
272 |
missing_dates_obj <- finding_missing_dates (date_start,date_end,list_dates_produced_date_val) |
|
273 |
|
|
274 |
df_time_series <- missing_dates_obj$df_dates |
|
275 |
df_time_series$date <- as.character(df_time_series$date) |
|
276 |
df_files$date <- as.character(df_files$date) |
|
277 |
|
|
278 |
df_time_series <- merge(df_time_series,df_files,by="date",all=T) #outer join to keep missing dates |
|
279 |
|
|
280 |
df_time_series$month_str <- format(as.Date(df_time_series$date), "%b") ## Month, char, abbreviated |
|
281 |
df_time_series$year_str <- format(as.Date(df_time_series$date), "%Y") ## Year with century |
|
282 |
df_time_series$day <- as.numeric(format(as.Date(df_time_series$date), "%d")) ## numeric month |
|
283 |
|
|
284 |
df_time_series_fname <- file.path(out_dir,paste0("df_time_series_",out_suffix,".txt")) #add the name of var later (tmax) |
|
285 |
write.table(df_time_series,file= df_time_series_fname,sep=",",row.names = F) |
|
286 |
|
|
287 |
df_time_series_obj <- list(df_raster_fname,df_time_series_fname,df_time_series) |
|
288 |
names(df_time_series_obj) <- c("df_raster_fname","df_time_series_fname","df_time_series") |
|
289 |
|
|
290 |
## report in text file missing by year and list of dates missing in separate textfile!! |
|
291 |
return(df_time_series_obj) |
|
292 |
} |
|
293 |
|
|
294 |
|
|
295 |
##### |
|
213 | 296 |
NAvalue(r_stack) |
214 | 297 |
plot(r_stack,y=6,zlim=c(-10000,10000)) #this is not rescaled |
215 | 298 |
#plot(r_stack,zlim=c(-50,50),col=matlab.like(255)) |
216 | 299 |
var_name <- "dailyTmax" |
217 |
debug(plot_and_animate_raster_time_series) |
|
300 |
|
|
301 |
|
|
302 |
|
|
303 |
|
|
304 |
|
|
305 |
|
|
306 |
#debug(plot_and_animate_raster_time_series) |
|
218 | 307 |
|
219 | 308 |
metric_name <- "var_pred" #use RMSE if accuracy |
220 | 309 |
#df_raster <- read.table("df_raster_global_assessment_reg6_10102016.txt",sep=",",header=T) |
... | ... | |
234 | 323 |
out_suffix=out_suffix, |
235 | 324 |
out_dir=out_dir) |
236 | 325 |
|
237 |
## Create function here: |
|
238 |
|
|
239 |
plot_and_animate_raster_time_series <- function(lf_raster, item_no,region_name,var_name,metric_name,NA_flag_val,filenames_figures=NULL,frame_speed=60,animation_format=".gif",zlim_val=NULL,plot_figure=T,generate_animation=T,num_cores=2,out_suffix="",out_dir="."){ |
|
240 |
#Function to generate figures and animation for a list of raster |
|
241 |
# |
|
242 |
# |
|
243 |
#INPUTS |
|
244 |
#1) lf_raster |
|
245 |
#2) filenames_figures |
|
246 |
#2) NAvalue |
|
247 |
#3) item_no |
|
248 |
#4) region_name, |
|
249 |
#5) var_name |
|
250 |
#6) metric_name |
|
251 |
#7) frame_speed |
|
252 |
#8) animation_format |
|
253 |
#9) zlim_val |
|
254 |
#10) plot_figure |
|
255 |
#11) generate_animation |
|
256 |
#12) num_cores |
|
257 |
#13) out_suffix |
|
258 |
#14) out_dir |
|
259 |
#OUTPUTS |
|
260 |
# |
|
261 |
# |
|
262 |
|
|
263 |
|
|
264 |
|
|
265 |
#lf_mosaic_list <- lf_raster |
|
266 |
variable_name <- var_name |
|
267 |
|
|
268 |
if(!is.null(plot_figure)){ |
|
269 |
#item_no <- 13 |
|
270 |
list_dates_produced <- unlist(mclapply(1:length(lf_raster), |
|
271 |
FUN = extract_date, |
|
272 |
x = lf_raster, |
|
273 |
item_no = item_no, |
|
274 |
mc.preschedule = FALSE, |
|
275 |
mc.cores = num_cores)) |
|
276 |
|
|
277 |
list_dates_produced_date_val <- as.Date(strptime(list_dates_produced, "%Y%m%d")) |
|
278 |
month_str <- format(list_dates_produced_date_val, "%b") ## Month, char, abbreviated |
|
279 |
year_str <- format(list_dates_produced_date_val, "%Y") ## Year with century |
|
280 |
day_str <- as.numeric(format(list_dates_produced_date_val, "%d")) ## numeric month |
|
281 |
df_raster <- data.frame(lf = basename(lf_raster), |
|
282 |
date = list_dates_produced_date_val, |
|
283 |
month_str = month_str, |
|
284 |
year = year_str, |
|
285 |
day = day_str, |
|
286 |
dir = dirname(lf_raster)) |
|
287 |
|
|
288 |
df_raster_fname <- file.path(out_dir, paste0("df_raster_", out_suffix, ".txt")) |
|
289 |
write.table(df_raster,file = df_raster_fname,sep = ",",row.names = F) |
|
290 |
|
|
291 |
############### PART5: Make raster stack and display maps ############# |
|
292 |
#### Extract corresponding raster for given dates and plot |
|
293 |
|
|
294 |
r_stack <- stack(lf_raster,quick=T) |
|
295 |
l_dates <- list_dates_produced_date_val #[1:11] |
|
296 |
|
|
297 |
#undebug(plot_raster_mosaic) |
|
298 |
zlim_val <- zlim_val |
|
299 |
|
|
300 |
### Now run for the full time series |
|
301 |
#13.26 Western time: start |
|
302 |
#l_dates <- list_dates_produced_date_val |
|
303 |
#r_stack_subset <- r_stack |
|
304 |
#zlim_val <- NULL |
|
305 |
out_suffix_str <- paste0(var_name,"_",metric_name,"_",out_suffix) |
|
306 |
list_param_plot_raster_mosaic <- list(l_dates,r_stack,NA_flag_val,out_dir, |
|
307 |
out_suffix_str,region_name,variable_name,zlim_val) |
|
308 |
names(list_param_plot_raster_mosaic) <- c("l_dates","r_mosaiced_scaled","NA_flag_val_mosaic","out_dir", |
|
309 |
"out_suffix", "region_name","variable_name","zlim_val") |
|
310 |
|
|
311 |
#lf_mosaic_plot_fig <- mclapply(1:length(l_dates[1:11]), |
|
312 |
# FUN = plot_raster_mosaic, |
|
313 |
# list_param = list_param_plot_raster_mosaic, |
|
314 |
# mc.preschedule = FALSE, |
|
315 |
# mc.cores = num_cores) |
|
316 |
|
|
317 |
##start at 12.29 |
|
318 |
##finished at 15.23 (for reg 6 with 2,991 figures) |
|
319 |
lf_mosaic_plot_fig <- mclapply(1:length(l_dates), |
|
320 |
FUN = plot_raster_mosaic, |
|
321 |
list_param = list_param_plot_raster_mosaic, |
|
322 |
mc.preschedule = FALSE, |
|
323 |
mc.cores = num_cores) |
|
324 |
|
|
325 |
if (is.null(zlim_val)) { |
|
326 |
out_suffix_movie <- paste("min_max_", out_suffix_str, sep = "") |
|
327 |
} else{ |
|
328 |
zlim_val_str <- paste(zlim_val, sep = "_", collapse = "_") |
|
329 |
out_suffix_movie <- paste(zlim_val_str, "_", out_suffix, sep = "") |
|
330 |
} |
|
331 |
filenames_figures_mosaic <- paste0("list_figures_animation_", out_suffix_movie, ".txt") |
|
332 |
|
|
333 |
write.table(unlist(lf_mosaic_plot_fig),filenames_figures_mosaic,row.names = F,col.names = F,quote = F) |
|
334 |
|
|
335 |
} |
|
336 |
|
|
337 |
### Part 2 generate movie |
|
338 |
|
|
339 |
if(generate_animation==TRUE){ |
|
340 |
|
|
341 |
out_suffix_str <- paste0(var_name,"_",metric_name,"_",out_suffix) |
|
342 |
|
|
343 |
if (is.null(zlim_val)) { |
|
344 |
out_suffix_movie <- paste("min_max_", out_suffix, sep = "") |
|
345 |
} else{ |
|
346 |
zlim_val_str <- paste(zlim_val, sep = "_", collapse = "_") |
|
347 |
out_suffix_movie <- paste(zlim_val_str, "_", out_suffix, sep = "") |
|
348 |
} |
|
349 |
|
|
350 |
#already provided as a parameter |
|
351 |
#filenames_figures_mosaic <- paste0("list_figures_animation_", out_suffix_movie, ".txt") |
|
352 |
#write.table(unlist(lf_mosaic_plot_fig),filenames_figures_mosaic,row.names = F,col.names = F,quote = F) |
|
353 |
|
|
354 |
#now generate movie with imageMagick |
|
355 |
#frame_speed <- 60 |
|
356 |
#animation_format <- ".gif" |
|
357 |
#out_suffix_str <- out_suffix |
|
358 |
#started |
|
359 |
#debug(generate_animation_from_figures_fun) |
|
360 |
#lf_mosaic_plot_fig <- read.table(filenames_figure,sep=",") |
|
361 |
|
|
362 |
#out_filename_figure_animation <- generate_animation_from_figures_fun(filenames_figures = unlist(lf_mosaic_plot_fig[1:11]), |
|
363 |
# frame_speed = frame_speed, |
|
364 |
# format_file = animation_format, |
|
365 |
# out_suffix = out_suffix_str, |
|
366 |
# out_dir = out_dir, |
|
367 |
# out_filename_figure_animation = "test2_reg6_animation.gif") |
|
368 |
|
|
369 |
#started 17.36 Western time on Oct 10 and 18.18 |
|
370 |
# 15.58 oct 11 for 16.38 for reg6 pred (about 2991) |
|
371 |
out_filename_figure_animation <- generate_animation_from_figures_fun(filenames_figures = filenames_figures_mosaic, |
|
372 |
frame_speed = frame_speed, |
|
373 |
format_file = animation_format, |
|
374 |
out_suffix = out_suffix_movie, |
|
375 |
out_dir = out_dir, |
|
376 |
out_filename_figure_animation = NULL) |
|
377 |
} |
|
378 |
|
|
379 |
## prepare object to return |
|
380 |
|
|
381 |
figure_animation_obj <- list(filenames_figures_mosaic,out_filename_figure_animation) |
|
382 |
names(figure_animation_obj) <- c("filenames_figures_mosaic","out_filename_figure_animation") |
|
383 |
return(figure_animation_obj) |
|
384 |
|
|
385 |
} |
|
386 | 326 |
|
387 | 327 |
############ Now accuracy |
388 | 328 |
#### PLOT ACCURACY METRICS: First test #### |
Also available in: Unified diff
adding checking missing tiles and moving plotting and animate raster time series to function script