Project

General

Profile

« Previous | Next » 

Revision 03a2bc79

Added by Benoit Parmentier about 8 years ago

adding checking missing tiles and moving plotting and animate raster time series to function script

View differences:

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