Vẽ timeserie bằng ggplot

Một bài tutorial nhỏ về cách vẽ đồ thị về sự dao động của nồng độ chất X theo thời gian + vùng sáng tối tương ứng với ban ngày/ban đêm (mùa đông ban ngày 7h - 21h, mùa hè từ 9h đến 18h) Tần số đo là 5 phút trong vòng 1 năm. Nhưng hình vẽ sẽ chỉ giới hạn trong vòng 1 tuần.

Lồng vào cách vẽ, mình muốn giới thiệu sơ qua:

• Signal được xử lí bằng 3 cách: moving average, savitzky golay và median filter
• Cách tô màu sáng tối tương ứng với ban ngày và ban đêm

Mục đích của plot này là để so sánh 3 cách xử lí tín hiệu và để xem nồng độ chất X có pattern theo chu kì ban ngày/ban đêm không.

Load dữ liệu

Do không tiện công bố data nên, bước này mình bỏ trống.

Data có dạng sau

head(df0)

##                  Time        X tick
## 1 2015-08-19 00:03:07 112.9517  raw
## 2 2015-08-19 00:08:07 133.6279  raw
## 3 2015-08-19 00:13:07 137.6039  raw
## 4 2015-08-19 00:18:07 196.4246  raw
## 5 2015-08-19 00:23:07 102.3617  raw
## 6 2015-08-19 00:28:07 158.9293  raw

2. Tạo function để xử lí dữ liệu

library(prospectr)
  filtration <- function(fnc,tick,df10){
    df <- df10
    line <- dim(df)[1]
    if (fnc == "movav"){
      df00 <- movav(df$X,5)
      n <- line - length(df00)
      dat <- data.frame(Time = df$Time[(n/2+1):(nrow(df)-n/2)], X = df00, tick = tick, stringsAsFactors = F)
    }
    if (fnc =="savitzky"){
      df00 <- savitzkyGolay(df$X,p=3,w=5,m=0)
      n <- line - length(df00)
      dat <- data.frame(Time = df$Time[(n/2+1):(nrow(df)-n/2)], X = df00, tick = tick, stringsAsFactors = F)
    }
    if (fnc == "median"){
      df00 <- robfilter::med.filter(df$X, width = 6)
      df00 <- as.matrix(df00$level$MED)
      dat <- data.frame(Time = df$Time, X = df00, tick = tick, stringsAsFactors = F)
    }
    return(dat)
  }
  
  #moving average
  df3 <- filtration("movav","moving average",df0)
  head(df3)

##                  Time        X           tick
## 1 2015-08-19 00:13:07 136.5940 moving average
## 2 2015-08-19 00:18:07 145.7895 moving average
## 3 2015-08-19 00:23:07 142.8064 moving average
## 4 2015-08-19 00:28:07 146.5619 moving average
## 5 2015-08-19 00:33:07 143.6309 moving average
## 6 2015-08-19 00:38:07 163.6619 moving average

  #savitzky Golay
  df4 <- filtration("savitzky","savitzky-golay",df0)
  head(df4)

##                  Time        X           tick
## 1 2015-08-19 00:13:07 161.5416 savitzky-golay
## 2 2015-08-19 00:18:07 152.6038 savitzky-golay
## 3 2015-08-19 00:23:07 149.5842 savitzky-golay
## 4 2015-08-19 00:28:07 122.7507 savitzky-golay
## 5 2015-08-19 00:33:07 141.4129 savitzky-golay
## 6 2015-08-19 00:38:07 147.9980 savitzky-golay

  #median filter
  df5 <- filtration("median","median filter",df0)
  head(df5)

##                  Time        X          tick
## 1 2015-08-19 00:03:07 135.6159 median filter
## 2 2015-08-19 00:08:07 135.6159 median filter
## 3 2015-08-19 00:13:07 135.6159 median filter
## 4 2015-08-19 00:18:07 135.6159 median filter
## 5 2015-08-19 00:23:07 135.6159 median filter
## 6 2015-08-19 00:28:07 146.9926 median filter

3. Vẽ hình

Tạo matrix cho vùng sáng tối. Từ tháng 4 đến tháng 9, ban ngày từ 7h - 21h. Từ tháng 10 đén tháng 3, ban ngày từ 9h-18h

  df0 <- df0 %>% mutate(day = day(df0$Time), 
                        hour = hour(df0$Time),
                        min = minutes(df0$Time),
                        sec = seconds(df0$Time))
  mat_summer <- df0 %>% group_by(day) %>% filter(hour==7) %>% slice(1) %>% 
    ungroup()
  
  mat_winter <- df0 %>% group_by(day) %>% filter(hour==9) %>% slice(1) %>% 
    ungroup()
  
  if(debut_month %in% c("Apr,May,Jun,Jul,Aug,Sep")){
    xmin <- mat_summer$Time
    xmax <- mat_summer$Time + 60*60*12
  }else{
    xmin <- mat_winter$Time
    xmax <- mat_winter$Time + 60*60*12
  }
  mat <- data.frame(xmin = xmin, xmax=xmax)

Tạo function plot

  ploti <- function(dataraw,datasmooth){
    p <- ggplot() + geom_line(data = dataraw, aes(x= Time, y = X,color = tick), size = 0.25) +
      geom_point(data = datasmooth, aes(x= Time, y = X,color = tick), size = 0.5) + theme_bw()+
      xlab("") + ylab("X(µg/L)") + theme(legend.position="top") +
      scale_x_datetime(labels = scales::date_format("%Y-%m-%d"))+
      theme(legend.title=element_blank())+
      geom_rect(data=mat, aes(xmin=xmin, xmax=xmax, ymin=-Inf, ymax = Inf), 
                fill='yellow', alpha=0.2) 
    
  }
  
  p1 <- ploti(df0,df3)
  p2 <- ploti(df0,df4)
  p3 <- ploti(df0,df5)

Kết quả cuối cùng

  gridExtra::grid.arrange(p1,p2,p3, top = "Raw data")