NA_value <- NA_flag_val 

  #metric_name <- "rmse" #to be added to the code later...

  #data_type <- "Int16" #, param 19, use int32 for output layers mosaiced

  ##### Select relevant day and create stack of missing tiles

  missing_tiles <- df_missing_tiles_day[i,]

  col_mfrow <- 1

  row_mfrow <- 1

  title_str <-  paste("Number of predicted pixels for ",variable_name," on ",date_str, sep = "")

  #my_col=c('blue','red','green')

  my_col <- rainbow(length(tb_freq$value))

  plot(r_day_predicted,col=my_col,legend=F,box=F,axes=F,main=title_str)

  col_mfrow <- 1

  row_mfrow <- 1

  title_str <-  paste("Number of predicted pixels for ",variable_name," on ",date_str, sep = "")

  #my_col=c('blue','red','green')

  my_col <- c("black","red")

  plot(r_missing_day,col=my_col,legend=F,box=F,axes=F,main=title_str)

  ### generate return object

  return(obj_number_day_predicted)

}

  setwd(out_dir)

  #list_outfiles <- vector("list", length=35) #collect names of output files, this should be dynamic?

  #list_outfiles_names <- vector("list", length=35) #collect names of output files

  r_overlap_m <- mask(r_overlap,

                  mask=r_mask,

                  filename=out_mosaic_name_overlap_masked,

                  options=c("COMPRESS=LZW"),#compress tif

                  overwrite=TRUE,

                  NAflag=NA_flag_val)

  col_mfrow <- 1

  row_mfrow <- 1

  title_str <-  paste("Maximum overlap: Number of predicted pixels for ",variable_name, sep = "")

    #my_col=c('blue','red','green')

  my_col <- rainbow(length(tb_freq_overlap$value))

                                                 FUN=function(i){mask(raster(list_tiles_predicted_m[i]),

                                                                      r_mask,filename=list_mask_out_file_name[i],

                                                                      overwrite=T,

                                                                      options=c("COMPRESS=LZW"))},

                                                 mc.preschedule=FALSE,

                                                 mc.cores = num_cores))                         

  list_param_generate_raster_number_pred <- list(list_tiles_predicted_masked,df_missing_tiles_day,r_overlap_m,

                                                 NA_flag_val,out_suffix,out_dir)

  names(list_param_generate_raster_number_pred) <- c("list_tiles_predicted_masked","df_missing_tiles_day","r_overlap_m",

                                                      "NA_flag_val","out_suffix","out_dir")

  #function_product_assessment_part0_functions <- "global_product_assessment_part0_functions_11152016b.R"

  #undebug(generate_raster_number_of_prediction_by_day)

  #4.51pm

  #5.10pm

  #test_number_pix_predictions <- generate_raster_number_of_prediction_by_day(1,list_param=list_param_generate_raster_number_pred)

  if(nrow(df_missing_tiles_day)>0){

  }

  predictions_tiles_missing_obj <- list(df_lf_tiles_time_series,df_missing_tiles_day,out_mosaic_name_overlap_masked,

  names(predictions_tiles_missing_obj) <- c("df_lf_tiles_time_series","df_missing_tiles_day","raster_name_overlap",

  return(predictions_tiles_missing_obj)

}