서울시의 미세/초미세먼지(1/n)
서울특별시 시간별 (초)미세먼지
<R로 하는 빅데이터 분석: 데이터 전처리와 시각화>를 활용하는 예시입니다.
압축된 데이터 읽기
서울특별시 시간별 (초)미세먼지 데이터를 살펴보자.
먼저 데이터 파일은 서울특별시_시간별 (초)미세먼지_20210531.zip으로 다음의 파일이 압축되어 있다. * 서울시 대기질 자료 제공_2020-2021.csv( 9MB) * 서울시 대기질 자료 제공_2016-2019.csv(27MB) * 서울시 대기질 자료 제공_2012-2015.csv(27MB) * 서울시 대기질 자료 제공_2008-2011.csv(27MB)
zip 파일을 압축 해제하여 4개의 csv 파일에서 자료를 읽어올 수도 있지만, 책의 <압축 파일에서 읽어오기>(p.133)을 적용할 수 있는 적절한 데이터 구성이다. 이 방법은 압축을 해제 하지 않기때문에 총 90MB의 디스크를 절약할 수 있다. 저자의 16GB 메모리 컴퓨터에서 90MB는 그다지 크지 않은 용량이기 때문에 압축 파일을 굳이 해제할 필요가 없어 보인다.
다음은 책의 진행 과정을 따라 간다. 먼저 압축된 파일 목록을 확인하자.
library(readr)
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionfn_zip = 'data/서울특별시_시간별 (초)미세먼지_20210531.zip'
fns = as.character(unzip(fn_zip, list=TRUE)$Name)
cat("* 압축된 파일 목록\n")## * 압축된 파일 목록cat(paste0(fns, collapse='\n'))## 서울시 대기질 자료 제공_2008-2011.csv
## 서울시 대기질 자료 제공_2012-2015.csv
## 서울시 대기질 자료 제공_2016-2019.csv
## 서울시 대기질 자료 제공_2020-2021.csv이제 압축된 파일을 읽어서 하나의 데이터프레임으로 합친다.
lst = vector("list", length(fns))
for (i in seq_along(fns)) {
lst[[i]] = read_csv(unz(fn_zip, fns[i]),
col_types = cols())
#lst[[i]]$src = fns[i]
}
dat <- data.table::rbindlist(lst)
head(dat)## Warning in do.call("cbind", lapply(x, function(col, ...) {: unable to translate
## '<c0><U+03FD><c3>' to native encoding## Warning in do.call("cbind", lapply(x, function(col, ...) {: unable to translate
## '<U+00B1><U+00B8><U+00BA><d0>' to native encoding## Warning in do.call("cbind", lapply(x, function(col, ...) {: unable to translate
## '<U+00B9><U+033C><U+00BC><U+00B8><d5><c1><f6>(PM10)' to native encoding## Warning in do.call("cbind", lapply(x, function(col, ...) {: unable to translate
## '<c3><U+02B9><U+033C><U+00BC><U+00B8><d5><c1><f6>(PM25)' to native encoding## <c0><U+03FD><c3> <U+00B1><U+00B8><U+00BA><d0>
## 1: 2011-12-31 23:00 <c6><f2><U+00B1><d5>
## 2: 2011-12-31 23:00 <U+00B0><U+00AD><U+00B3><U+00B2><U+00B1><U+00B8>
## 3: 2011-12-31 23:00 <U+00B0><U+00AD><U+00B5><U+00BF><U+00B1><U+00B8>
## 4: 2011-12-31 23:00 <U+00B0><U+00AD><U+00BA><U+03F1><U+00B8>
## 5: 2011-12-31 23:00 <U+00B0><U+00AD><U+00BC><U+00AD><U+00B1><U+00B8>
## 6: 2011-12-31 23:00 <U+00B0><fc><U+00BE><U+01F1><U+00B8>
## <U+00B9><U+033C><U+00BC><U+00B8><d5><c1><f6>(PM10)
## 1: 89
## 2: 91
## 3: 89
## 4: 68
## 5: 79
## 6: 84
## <c3><U+02B9><U+033C><U+00BC><U+00B8><d5><c1><f6>(PM25)
## 1: 61
## 2: 58
## 3: 59
## 4: 49
## 5: 69
## 6: 55읽긴 했는데, 소위 문자가 깨진다. 인코딩 문제인 듯 하니 인코딩을 새롭게 설정해서 읽어본다. (인코딩에 대한 자세한 내용은 책의 11장 문자열에서 자세히 다룬다.) 많이 쓰는 인코딩은 UTF-8과 CP949이고, readr::read_csv()의 기본값은 UTF-8이므로 CP949로 바꿔본다. (이때 CP949라고 하면 Microsoft의 CodePage를 의미할 수도 있고, IBM의 CodePage를 의미할 수도 있으니 주의할 필요가 있다. p205)
lst = vector("list", length(fns))
for (i in seq_along(fns)) {
lst[[i]] = read_csv(unz(fn_zip, fns[i]),
col_types = cols(),
locale = locale(encoding = 'CP949'))
#lst[[i]]$src = fns[i]
}
dat <- data.table::rbindlist(lst)
head(dat)## 일시 구분 미세먼지(PM10) 초미세먼지(PM25)
## 1: 2011-12-31 23:00 평균 89 61
## 2: 2011-12-31 23:00 강남구 91 58
## 3: 2011-12-31 23:00 강동구 89 59
## 4: 2011-12-31 23:00 강북구 68 49
## 5: 2011-12-31 23:00 강서구 79 69
## 6: 2011-12-31 23:00 관악구 84 55제대로 읽었다! 그런데 위의 head(dat) 결과를 보면, 일시가 <chr>이다. 우선 이를 수정해주자.
dat$일시 = as.POSIXct(dat$일시, format = "%Y-%m-%d %H:%M")
# 이 부분은 시간이 조금 걸린다.
# 위의 read_csv(, col_types=)로 해결할 수도 있다.
head(dat)## 일시 구분 미세먼지(PM10) 초미세먼지(PM25)
## 1: 2011-12-31 23:00:00 평균 89 61
## 2: 2011-12-31 23:00:00 강남구 91 58
## 3: 2011-12-31 23:00:00 강동구 89 59
## 4: 2011-12-31 23:00:00 강북구 68 49
## 5: 2011-12-31 23:00:00 강서구 79 69
## 6: 2011-12-31 23:00:00 관악구 84 55데이터 살펴보기
<14.2 데이터 프레임의 모든 변수에 대해 요약통계치 구하기>에는 데이터를 요약하는 여러가지 방법이 소개되어 있다. 가장 첫 걸음은 아마도 summary()가 아닐까?
summary(dat)## 일시 구분 미세먼지(PM10)
## Min. :2008-01-01 01:00:00 Length:3044017 Min. : 1.00
## 1st Qu.:2011-05-18 02:00:00 Class :character 1st Qu.: 25.00
## Median :2014-09-22 02:00:00 Mode :character Median : 38.00
## Mean :2014-09-21 09:07:08 Mean : 45.56
## 3rd Qu.:2018-01-28 10:00:00 3rd Qu.: 57.00
## Max. :2021-05-31 23:00:00 Max. :1354.00
## NA's :77449
## 초미세먼지(PM25)
## Min. :-273.00
## 1st Qu.: 13.00
## Median : 20.00
## Mean : 24.59
## 3rd Qu.: 31.00
## Max. :9940.00
## NA's :62650일시를 보면 데이터는 2008-01-01에서 시작하여 2021-05-31에 끝난다. 구분에는 별다른 정보가 드러나지 않기 때문에 table()을 해본다. 미세먼지, 초미세먼지를 NA도 상당히 많고, 초미세먼지(PM25)는 최솟값이 음수이다. 일단 미세먼지와 초미세먼지의 측정 단위를 확인해 볼 필요가 있다!
table(dat$구분)##
## 강남구 강동구 강북구 강서구 관악구 광진구 구로구 금천구
## 115348 116463 116911 116192 117403 117434 117480 116742
## 노원구 도봉구 동대문구 동작구 마포구 서대문구 서초구 성동구
## 117489 117203 116717 117168 116992 117024 117569 117570
## 성북구 송파구 양천구 영등포구 용산구 은평구 종로구 중구
## 117316 117538 116649 117216 117531 116951 116919 117491
## 중랑구 평균
## 117104 117597검색 결과 미세먼지(PM-10) 농도와 초미세먼지(PM-2.5) 농도의 단위는 모두 \(\mu g/m^3\) 또는 ppm(백만분율; parts per million)인 듯 합니다.1 (책의 부록 B.측정단위에는 단위를 관리하는 방법을 소개합니다.) 그리고 최솟값은 0이므로 초미세먼지 농도가 음수라면 이는 잘못된 데이터로 보입니다.
sum(dat$`초미세먼지(PM25)`<0, na.rm=TRUE)## [1] 48잘못된 정보가 많지는 않습니다. 언제 어디서 이런 잘못된 정보가 생겼는지 한번 살펴봅니다.
dat[`초미세먼지(PM25)`<0] %>% head## 일시 구분 미세먼지(PM10) 초미세먼지(PM25)
## 1: 2009-12-04 00:00:00 양천구 NA -273
## 2: 2009-12-03 23:00:00 양천구 NA -266
## 3: 2009-12-03 22:00:00 양천구 NA -270
## 4: 2009-12-03 21:00:00 양천구 NA -266
## 5: 2009-12-03 20:00:00 양천구 NA -269
## 6: 2009-12-03 19:00:00 양천구 NA -2722009년 12월 2일과 12월 3일의 양천구에 문제가 좀 있었나보네요. 다른 지역도 함께 확인해보죠.
dat[as.Date(`일시`)==as.Date('2009-12-02')] %>% head## 일시 구분 미세먼지(PM10) 초미세먼지(PM25)
## 1: 2009-12-03 08:00:00 평균 80 47
## 2: 2009-12-03 08:00:00 강남구 123 72
## 3: 2009-12-03 08:00:00 강동구 99 81
## 4: 2009-12-03 08:00:00 강북구 100 64
## 5: 2009-12-03 08:00:00 강서구 63 60
## 6: 2009-12-03 08:00:00 관악구 61 38다른 지역은 크게 문제가 없는 듯 하니 모두 NA로 바꾸는 게 나을 듯 합니다.
dat[`초미세먼지(PM25)`<0, `초미세먼지(PM25)`:=NA]이제 서울 평균 미세먼지와 초미세먼지 농도의 분포를 히스토그램으로 확인해보죠. 데이터가 300여만 건으로 적지 않은 수이므로 데이터테이블의 함수를 사용하겠습니다. 데이터 프레임이나 티블보다 훨씬 빠르게 결과를 얻을 수 있습니다(책 10장).
library(ggplot2)
dat[`구분` == '평균',
ggplot(.SD, aes(x=`미세먼지(PM10)`)) +
geom_histogram() +
labs(title='미세먼지(PM10) 농도(평균)')
]## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
이 히스토그램을 보니 1000 이상의 값이 진짜일까란 생각이 드네요. 아마도 시계열로 시각화하면 시간의 연속성에 따라 300 이상의 극단값이 실제값인지 아니면 오류인지 확인이 가능할 것이란 생각이 듭니다. 초미세먼지의 분포도 그려보죠.
dat[`구분` == '평균',
ggplot(.SD, aes(x=`초미세먼지(PM25)`)) +
geom_histogram() +
labs(title='초미세먼지(PM2.5) 농도(평균)')
]## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
이제 시계열로 그려보죠.
dat[`구분` == '평균',
ggplot(.SD, aes(x=`일시`, y=`미세먼지(PM10)`)) +
geom_line() +
labs(title='미세먼지(PM10) 농도(평균)')
]
아무래도 정보가 너무 많습니다. 년도별, 월별, 일별 평균을 구해서 다시 그려보죠. 먼저 열이름에 한글이 거추장스러우니 열이름을 바꾼 후 시작합니다.
library(lubridate)##
## Attaching package: 'lubridate'## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, unionlibrary(data.table)##
## Attaching package: 'data.table'## The following objects are masked from 'package:lubridate':
##
## hour, isoweek, mday, minute, month, quarter, second, wday, week,
## yday, year## The following objects are masked from 'package:dplyr':
##
## between, first, lastdat[ , byyear := year(`일시`)]
setnames(dat, '미세먼지(PM10)', 'pm10')
setnames(dat, '초미세먼지(PM25)', 'pm25')
setnames(dat, '일시', 'datetime')
setnames(dat, '구분', 'place')이제 년도별 평균을 그려보면 다음과 같습니다. (아래 그림과 비교를 위해 수평선 두 개도 첨가합니다.)
dat[, date_year:=as.Date(as.character(byyear), "%Y")]
dat2 <-
dat[place=='평균',.(mean_pm10=mean(pm10, na.rm=TRUE)),
by=date_year]
min_mean_pm10 = dat2[, min(mean_pm10, na.rm=TRUE)]
max_mean_pm10 = dat2[, max(mean_pm10, na.rm=TRUE)]
dat2[, ggplot(.SD, aes(x=date_year, y=mean_pm10)) +
geom_line() + geom_point() +
scale_x_date(date_labels="%y", date_breaks="1 year") +
geom_hline(yintercept = min_mean_pm10, linetype='dashed') +
geom_hline(yintercept = max_mean_pm10, linetype='dashed')
]
요약 통계치를 쓰지 않은 위의 그래프와 비교를 해보죠. 평균은 극단값에 영향을 많이 받습니다. 만약 .1-절사평균을 쓴다면 어떨까요? .1-절사평균은 mean(, trim=0.05)로 구할 수 있습니다.
dat2 <-
dat[,.(trmean_pm10=base::mean(pm10, na.rm=TRUE, trim=0.05)),
by=date_year]
dat2[, ggplot(.SD, aes(x=date_year, y=trmean_pm10)) +
geom_line() + geom_point() +
scale_x_date(date_labels="%y", date_breaks="1 year") +
geom_hline(yintercept = min_mean_pm10, linetype='dashed') +
geom_hline(yintercept = max_mean_pm10, linetype='dashed')
]
그래프의 모양이 비슷한 걸 보니 단순히 극단값의 문제는 아닌 듯 합니다.
다시 한번 상자그림으로 년도별로 농도의 분포를 살펴보죠.
source('utils.R')
dat[place=='평균',
ggplot(.SD, aes(x=as.character(byyear), y=pm10)) +
geom_boxplot() + x_label_vertical() +
geom_hline(yintercept = min_mean_pm10, linetype='dashed', color = 'blue') +
geom_hline(yintercept = max_mean_pm10, linetype='dashed', color = 'blue')
#+
#scale_x_date(date_labels="%y", date_breaks="1 year")
]
outlier가 너무 많아서 서로 식별하기가 쉽지 않네요. 바이올린은 괜찮을까요?
dat[place=='평균',
ggplot(.SD, aes(x=as.character(byyear), y=pm10)) +
geom_violin() + x_label_vertical() +
geom_hline(yintercept = min_mean_pm10, linetype='dashed', color = 'blue') +
geom_hline(yintercept = max_mean_pm10, linetype='dashed', color = 'blue')
]
조금 나아졌지만 아직 충분하진 않아 보입니다. ggplot의 확장패키지 ggridges를 활용해봅시다.
dat3 <-
dat[place=='평균']
dat3[ , pm10cut:=cut(pm10,
breaks=seq(round(min(pm10),-1),
round(max(pm10),-1),
10),
include.lowest=TRUE)]
dat4 <- dat3[ , .(freq = nrow(.SD)) ,by=.(byyear, pm10cut)]library(ggridges)
dat4[, ggplot(.SD,
aes(x=as.numeric(pm10cut),
y=byyear,
height = freq,
group=byyear)) +
geom_ridgeline()
]
흠. 뭔가 이상하죠? 년도를 나타내는 byyear와 빈도수를 나타내는 freq의 스케일 문제입니다. freq를 적당하게 조절합니다.
dat4 <- dat3[ , .(freq = nrow(.SD)) ,by=.(byyear, pm10cut)]
dat4[, freq:=freq/1100]dat4[, ggplot(.SD,
aes(x=as.numeric(pm10cut),
y=byyear,
height = freq,
group=byyear)) +
geom_ridgeline() +
scale_y_reverse()
]
2021년 분포가 지나치게 낮은데 2021년은 데이터가 5월까지만 존재하기 때문입니다. 이를 적절하게 수정합니다.
dat3[ , ny:=nrow(.SD), by=byyear]
dat4 <- dat3[ , .(density = nrow(.SD)/ny) ,by=.(byyear, pm10cut)]
dat4[, ggplot(.SD,
aes(x=as.numeric(pm10cut),
y=byyear,
height = density*10,
group=byyear)) +
geom_ridgeline() +
scale_y_reverse()
]
PM10 농도 100까지만 자세히 보고 싶어집니다.
dat4[, ggplot(.SD,
aes(x=as.numeric(pm10cut),
y=byyear,
height = density*10,
group=byyear)) +
geom_ridgeline() +
scale_y_reverse() +
coord_cartesian(xlim=c(0,20))
]
좀더 자세히 살펴보죠.
dat3 <-
dat[place=='평균']
dat3[ , pm10cut:=cut(pm10,
breaks=seq(round(min(pm10),0),
round(max(pm10),0),
1),
include.lowest=TRUE)]
dat4 <- dat3[ , .(freq = nrow(.SD)) ,by=.(byyear, pm10cut)]
dat3[ , ny:=nrow(.SD), by=byyear]
dat4 <- dat3[ , .(density = nrow(.SD)/ny) ,by=.(byyear, pm10cut)]
dat4[, ggplot(.SD,
aes(x=as.numeric(pm10cut),
y=byyear,
fill=factor(byyear),
height = density*200,
group=byyear)) +
geom_ridgeline(col='black') +
scale_y_reverse() +
coord_cartesian(xlim=c(0,100))
]
dat3 <-
dat[place=='평균']
dat3[, fac_year:=factor(byyear, levels=seq(2021,2008))]
dat3[,
ggplot(.SD, aes(x=pm10,
y=fac_year,
fill = fac_year)) +
stat_density_ridges(alpha=0.6) +
scale_fill_cyclical(values = c('lightskyblue', 'royalblue3')) +
geom_vline(xintercept=min_mean_pm10, linetype='dashed') +
geom_vline(xintercept=max_mean_pm10, linetype='dashed')
] ## Picking joint bandwidth of 3.36
stat_density_rdiges()를 활용하여 분포를 추정해볼 수도 있습니다. 위의 그래프와 다른 점은 무엇일까요? 분포를 부르러운 곡선으로 추정하지만, 데이터가 존재하지 않는 지점에서도 그래프가 그려지고 있다는 점입니다.
일단 오늘은 여기까지 하겠습니다.
(2부에서 계속)
참고 자료
Leave a comment