plot_observation_predictions_time_series <- function(df_pix_time_series,var_name1,var_name2,station_id,scaling=NULL,out_dir=".",out_suffix=""){

  #This function plots time series for observed and predicted points.

  #This assumes that time series with residuals were extracted from raster stacks.

  #

  #Inputs:

  #1) df_pix_time_series

  #2) var_name1: name of variables use

  #3) var_name2:

  #4) station_id

  #5) scaling: if null no scaling

  #6) out_dir

  #7) out_suffix

  #

  ########## Start script #########

  ## Screen of NaN and make sure we have NA

  df_pix_ts[is.na(df_pix_ts)] <- NA

  #Scale if necessary

  if(!is.null(scaling)){

    df_pix_ts[[var_pred_mosaic]] <-  df_pix_ts[[var_pred_mosaic]]*scaling 

  }

  #### var_name 1

  d_z_obs <- zoo(as.numeric(df_pix_ts[[y_var_name]]),as.Date(df_pix_ts$date))

  plot(d_z_obs)  

  #### var_name 2

  d_z_var <- zoo(as.numeric(df_pix_ts[[var_pred_mosaic]]),as.Date(df_pix_ts$date)) #make sure date is a date object !!!

  #names(d_z_var) <- var_pred_mosaic

  plot(d_z_var)  

  d_z_diff <- d_z_var - d_z_obs 

  #d_z_res <- d_z_var - d_z_obs 

  #d_z_res <- zoo(as.numeric(df_pix_ts[[paste0("res_",var_pred_mosaic)]]),as.Date(df_pix_ts$date))

  plot(d_z_res)  

  d_z_all <- merge(d_z_var,d_z_obs,d_z_res)

  names(d_z_all) <- c("pred","obs","res")

  d_z <-  merge(d_z_var,d_z_obs)

  names(d_z) <- c("pred","obs")

  range(d_z_res,na.rm=T)

  mean(d_z_res,na.rm=T)

  quantile(d_z_res,na.rm=T)

  range_dates <- range(as.Date(df_pix_time_series$date))

  #day_start <- "1984-01-01" #PARAM 12 arg 12

  #day_end <- "2014-12-31" #PARAM 13 arg 13

  #start_date <- as.Date(day_start)

  #end_date <- as.Date(day_end)

  start_year <- year(range_dates[1])

  end_year <- year(range_dates[2])

  #range_year <- range(df_pix_ts$year)

  #start_year <- range_year[1]

  #end_year <- range_year[2]

  var="tmax"

  res_pix <- 1000

  #res_pix <- 480

  col_mfrow <- 2

  row_mfrow <- 1

  png_filename <-  file.path(out_dir,paste("Figure_5a_time_series_profile_",region_name,"_",id_name,y_var_name,out_suffix,".png",sep =""))

  title_str <- paste("Daily", var_name1,"and", var_name2,"for", station_id," for the ", start_year,"-",end_year," time period",sep=" ")

  png(filename=png_filename,width = col_mfrow * res_pix,height = row_mfrow * res_pix)

  #this is the whole time series

  #plot(d_z_var,ylab="tmax in deg C",xlab="Daily time steps",

  #     main=title_str,cex=3,font=2,

  #     cex.main=1.5,cex.lab=1.5,font.lab=2,

  #     lty=3)

  #range(df_pix_ts[[var_pred_mosaic]],na.rm=T)

  #plot(d_z,plot.type="single",col=c("blue","red"))

  col_str <- c("blue","red")

  plot(d_z,plot.type="single",col=col_str,

       ylab="tmax in deg C",xlab="Daily time steps",

       main=title_str,cex=1,font=2,type="l",

       cex.main=1.5,cex.lab=1.5,font.lab=2)

  #plot(df_pix_ts[[var_pred_mosaic]] ~ as.Date(df_pix_ts$date),ylab="tmax in deg C",xlab="Daily time steps",

  #   main=title_str,cex=1,font=2,type="l",

  #   cex.main=1.5,cex.lab=1.5,font.lab=2,

  #   lty=3)

  #lines(as.numeric(df_pix_ts[[y_var_name]]) ~ as.Date(df_pix_ts$date),ylab="tmax in deg C",xlab="Daily time steps",

  #   main=title_str,cex=1,font=2,col="red",type="l",

  #   cex.main=1.5,cex.lab=1.5,font.lab=2,

  #   lty=3)

  dev.off()

  #for list_windows go in a loop

  if(!is.null(list_windows)){

    for(i in 1:length(list_windows)){

      ### get smaller window

      window_range <- list_windows[i]

      #day_start <- "1994-01-01" #PARAM 12 arg 12

      #day_end <- "1999-12-31" #PARAM 13 arg 13

      day_start <- window_range[1]

      day_end <- windwow_range[2]

      start_date <- as.Date(day_start)

      end_date <- as.Date(day_end)

      start_year <- year(start_date)

      end_year <- year(end_date)

      ### now select from the original time series

      #

      d_z_w <- window(d_z_tmp,start=start_date,end=end_date)

      names(d_z) <- var_pred_mosaic

      df_var <- as.data.frame(d_z)

      d_z_tmp <- d_z_obs

      d_z_tmp <- d_z

      d_z_tmp <- window(d_z_tmp,start=start_date,end=end_date)

      plot(d_z_tmp,col=c("blue","red"),plot.type="single")

      df_z <- as.data.frame(d_z)

      names(d_z) <- y_var_name

      df_obs <- as.data.frame(d_z)

      names(d_z) <- y_var_name

      res_pix <- 1000

      #res_pix <- 480

      col_mfrow <- 2

      row_mfrow <- 1

      png_filename <-  file.path(out_dir,paste("Figure5b_time_series_profile_",region_name,"_",out_suffix,".png",sep =""))

      title_str <- paste("Predicted daily ", station_id," ",var," for the ", start_year,"-",end_year," time period",sep="")

      png(filename=png_filename,width = col_mfrow * res_pix,height = row_mfrow * res_pix)

      plot(df_pix_ts[[var_pred_mosaic]] ~ as.Date(df_pix_ts$date),ylab="tmax in deg C",xlab="Daily time steps",

      main=title_str,cex=1,font=2,type="l",

      cex.main=1.5,cex.lab=1.5,font.lab=2,

      lty=3)

      lines(as.numeric(df_pix_ts[[y_var_name]]) ~ as.Date(df_pix_ts$date),ylab="tmax in deg C",xlab="Daily time steps",

      main=title_str,cex=1,font=2,col="red",type="l",

      cex.main=1.5,cex.lab=1.5,font.lab=2,

      lty=3)

      dev.off()

    }

  }

  #return()

}