서울특별시 시간별 (초)미세먼지 02

x축 단위 : 일

일별로 평균을 구한 후 연도별로 시각화해봅니다. 먼저 일별 평균을 구합시다.

library(data.table)
dat[, year:=year(datetime)]
dat[, md:=as.Date(format(datetime, "2000-%m-%d"))]
# 평균에 대, year, date으로 
dat2 <- dat[place == "평균", 
    .(m_pm10=mean(pm10, na.rm=TRUE), 
      m_pm25=mean(pm25, na.rm=TRUE)), 
    by=.(year,md)]

그리고 시각화를 합니다.

library(ggplot2)
dat2[, ggplot(.SD, aes(x=md, y=m_pm10)) + 
       geom_line() + 
       facet_grid(year~.)  
       ]

극단값이 있어서 년도별로 비교하기가 쉽지 않네요. y축의 범위를 조절하고, gghighlight() 사용해봅니다.

library(gghighlight)
dat2[, ggplot(.SD, aes(x=md, y=m_pm10, col=factor(year))) + 
       geom_line() + 
       facet_grid(year~.) +
       coord_cartesian(ylim=c(0,quantile(m_pm10, 0.95))) + 
       scale_x_date(date_labels="%m") + #p284
       gghighlight() # p327
       ]

확실히 7-10월에는 진폭이 잦아지면서 PM10 농도 일평균이 감소하는 경향이 있어보입니다. 하지만 노이즈가 너무 많은 것 같습니다. 주별로 평균을 구해서 그래프를 그려보죠. ISO8601의 년-주-일 형식을 활용하여 년-주로 바꿔봅니다.

x축 단위 : 주

library(data.table)
dat[, year:=as.numeric(format(datetime, "%G"))] # p87
dat[, week:=as.numeric(format(datetime, "%V"))] # p87
# 평균에 대, year, date으로 
dat2 <- dat[place == "평균", 
    .(m_pm10=mean(pm10, na.rm=TRUE), 
      m_pm25=mean(pm25, na.rm=TRUE),
      md= min(md)), 
    by=.(year,week)]

최근의 추세만 한 번 비교해 봅시다.

library(gghighlight)
m_pm10_q = quantile(dat2$m_pm10, c(0.1, 0.3, 0.5, 0.7, 0.9))
dat_pm10_q = data.frame(y = m_pm10_q)
dat2[, year2:=ifelse(year<2015, "<=2014", year)]
dat2[, 
     ggplot(.SD, aes(x=week, y=m_pm10)) + 
       geom_line(aes(group=year, col=factor(year2)),
                 size = 1.3) + 
       geom_point(aes(col=factor(year2))) + 
       facet_grid(year2~.) +
       coord_cartesian(ylim=c(0,quantile(m_pm10, 0.98))) + 
       geom_hline(mapping = aes(yintercept=y, 
                                linetype=factor(y)), 
                  alpha = 0.2, 
                  #linetype = 'dotted', 
                  data=dat_pm10_q) + 
       gghighlight() + # p327
       theme(legend.position = 'bottom') # p314
       ]

## label_key: year

## Warning: ggrepel: 5 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

보조 수평선은 10,30,50,70,90-분위수로 했습니다만 좋음/나쁨 등의 기준으로 바꿔도 좋을 듯 합니다. 대기오염도 등급 기준은 다음과 같습니다.[^2]

# p139
library(xml2)
library(rvest)

## 
## Attaching package: 'rvest'

## The following object is masked from 'package:readr':
## 
##     guess_encoding

url <- "https://www.gyeongnam.go.kr/knhe/index.gyeong?menuCd=DOM_000001702001004000"
html = read_html(url, encoding = 'UTF-8')
tables = html %>% html_nodes("table")
df = html_table(tables[[2]], convert=FALSE) 

library(knitr)
df %>% kable

물론 주단위 평균과 농도 기준을 바로 비교하긴 힘들 것입니다. 왜냐하면 평균이란 극단값에 많이 좌우되기 때문이죠. 중앙값을 그려본다면 전체 시간별 농도 중에 50%는 해당값보다 낮았다는 것을 알 수 있습니다(다른 백분위수도 동일한 해석이 가능합니다).

반면에 평균이 좋은 면도 있습니다. 표본 크기가 클수록(데이터가 많을 수록) 평균의 분포는 정규분포에 가까워지죠. 이산 분포에서도 평균의 분포는 연속적이지만 중앙값의 분포는 그렇지 않습니다.

그리고 이렇게 평균 또는 중앙값을 구하는 이유는 노이즈를 제거하는 의미도 있습니다. 주 단위로 평균을 구하는 것이 노이즈를 충분히 제거하지 못하는 듯 합니다. ggplot의 geom_smooth()를 이용해서 알아서 노이즈를 제거하도록 해봅시다.

dat2[, ggplot(.SD, aes(x=week, y=m_pm10, col=factor(year))) + 
       geom_smooth(se=FALSE) + 
       facet_grid(factor(year2)~.)+
       gghighlight()]

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

전체 데이터를 통해 SE도 확인해 봅니다.

library(ggrepel)
dat3 = dat2[, .(md = md[week==max(week)],
                m_pm10=m_pm10[week==max(week)]) , by=year]
#dat2[, date:=as.Date(paste0(year,"-", week), format="%Y-%")]
dat2[, ggplot(.SD, aes(x=md, y=m_pm10)) + 
       geom_smooth(aes(group=year, col=year), se=FALSE, alpha=0.3, col='grey') + 
       geom_smooth() + 
       scale_x_date(date_labels="%m")
       #geom_label(aes(col=year)) + 
       #geom_text_repel(aes(label=year)) #+ 
       #geom_point(dat3, aes(x=md, y=m_pm10))
     ]

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

확실히 7-10월의 PM10 농도가 가장 낮고, 3-4월에 PM10 농도가 가장 높은 경향이 있습니다. 아마도 어떤 계절적 영향이 있는 것 같습니다.