서울특별시 시간별 (초)미세먼지 03: 극단값 확인
서울특별시 시간별 (초)미세먼지 03: 극단값 확인
<R로 하는 빅데이터 분석: 데이터 전처리와 시각화>를 활용하는 예시입니다.
압축된 데이터 읽기
전처리 과정은 첫 번째 포스팅를 참조하세요.
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.csvlst = 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'))
# col_types = 'Tfdd') # 각 열의 타입까지 지정 p129
# 이것도 가능하지만 파일에 약간 문제가 있음
# 예: 2021-05-31 9:00 (2021-05-31 09:00이어야 함)
}
dat <- data.table::rbindlist(lst)
dat[`초미세먼지(PM25)`<0, `초미세먼지(PM25)`:=NA]
library(data.table)##
## Attaching package: 'data.table'## The following objects are masked from 'package:dplyr':
##
## between, first, lastsetnames(dat, '미세먼지(PM10)', 'pm10')
setnames(dat, '초미세먼지(PM25)', 'pm25')
setnames(dat, '일시', 'datetime')
setnames(dat, '구분', 'place')
#dat[, datetime:=as.POSIXct(datetime, format='%Y-%m-%d %H:%M')]
dat$datetime2 = as.POSIXct(dat$datetime, format = '%Y-%m-%d %H:%M')
#dat <- dat[!is.na(datetime)]
# dat[datetime == "2021-05-31 9:00", ] 2021-05-31 9:00과 같은 문자열도 정확히 인식됨
dat[, datetime:=datetime2]
dat[, datetime2:=NULL]
dat[, date:=as.Date(datetime, tz='Asia/Seoul')] # 주의 tz를 확실히 정하자.
#dat[, date:=as.Date(format(datetime, "%Y-%m-%d"))]
dat[, time:=as.POSIXct(format(datetime, "%H:%M:%S"), format='%H:%M:%S')]
head(dat)## datetime place pm10 pm25 date time
## 1: 2011-12-31 23:00:00 평균 89 61 2011-12-31 2022-02-07 23:00:00
## 2: 2011-12-31 23:00:00 강남구 91 58 2011-12-31 2022-02-07 23:00:00
## 3: 2011-12-31 23:00:00 강동구 89 59 2011-12-31 2022-02-07 23:00:00
## 4: 2011-12-31 23:00:00 강북구 68 49 2011-12-31 2022-02-07 23:00:00
## 5: 2011-12-31 23:00:00 강서구 79 69 2011-12-31 2022-02-07 23:00:00
## 6: 2011-12-31 23:00:00 관악구 84 55 2011-12-31 2022-02-07 23:00:00그런데 말입니다. 시간이 아주 많이 걸리는 것은 아니지만, 어쨋든 항상 데이터를 zip에서 불러와서 전처리를 하는 과정을 거쳐야 합니다. 그냥 전처리가 완료된 데이터를 저장하면 되지 않을까요? 용량이 큰 데이터의 경우 바이너리, 특히 .feather로 저장하면 입출력이 상대적으로 빠릅니다(p132).
library(feather)
write_feather(dat, 'data/airpollution_seoul.feather')다시 읽어봅시다.
dat2 <- read_feather('data/airpollution_seoul.feather')제 컴퓨터에서는 1초도 안 걸리네요. 🙂 데이터가 동일한지 확인해 볼까요?
all.equal(dat, as.data.table(dat2))## [1] TRUE다음부터는 전처리 부분을 dat <- as.data.table(read_feather('data/airpollution_seoul.feather'))로 대신하겠습니다.
단위 설정
이번엔 pm10의 단위도 설정해 줍시다.(p340) 이 사이트에 따르면 PM10과 PM2.5의 단위는 모두 \(\mu g / m^3\) 입니다.
먼저 적절한 단위가 있는지 검색해봅니다. gram과 meter가 동시에 들어간 단위가 있을까요?
library(units)## udunits database from C:/Users/Home/Documents/R/win-library/4.1/units/share/udunits/udunits2.xmllibrary(stringr)
datUnits = valid_udunits()## udunits database from C:/Users/Home/Documents/R/win-library/4.1/units/share/udunits/udunits2.xmldatUnits %>% filter_all(any_vars(str_detect(., "meter") & str_detect(., "kilo")))## # A tibble: 3 x 11
## symbol symbol_aliases name_singular name_singular_a~ name_plural
## <chr> <chr> <chr> <chr> <chr>
## 1 "ua" "" astronomical_unit_BIPM_2006 "" ""
## 2 "at" "" technical_atmosphere "" ""
## 3 "" "" metric_horsepower "" ""
## # ... with 6 more variables: name_plural_aliases <chr>, def <chr>,
## # definition <chr>, comment <chr>, dimensionless <lgl>, source_xml <chr>없는 것 같습니다. 일단 단위를 microgram/meter^3로 지정합니다.
units(dat$pm10) = 'microgram/meter^3'단위는 부록 B에서 다루고 있습니다. 예를 들어 다음과 같은 계산이 가능합니다.
단위 활용
미세먼지 ‘보통’ 속 산책, ‘나쁨’ 속 휴식보다 해롭다라는 한국일보 기사를 보면, 성인 남성의 안정시 호흡량은 분당 6 L 정도입니다. 그렇다면 주어진 PM10 농도에서 시간 당 호흡하는 미세먼지(PM10)의 양은 얼마나 될까요?
persp_adult_per_hour = 6 * 60
units(persp_adult_per_hour) = "L / h"
dat$pm10[1] * persp_adult_per_hour## 32040 [L*microgram/h/m^3]그런데 단위가 조금 이상합니다. 우리가 바라던 결과와 조금 다른 것 같은데요. dat$pm10의 단위가 microgram/meter^3이므로 pers_adult_per_hour의 단위도 비슷하게 바꿔줍니다.
units(persp_adult_per_hour) = "meter^3 / h"
persp_adult_per_hour## 0.36 [m^3/h]그리고 다시 계산을 해보면,
print(dat$pm10[1])## 89 [microgram/m^3]dat$pm10[1] * persp_adult_per_hour## 32.04 [microgram/h]따라서 \(\textrm{PM}_{10}\) 이 89 \(\mu g / m^3\) 일 때, 성인의 안정시 호흡하는 미세먼지의 양은 \(32.04 \mu g /h\) 라고 추정할 수 있습니다.
극단값 확인
논의를 더 진행하기 전에, 우리가 일별 또는 주별로 평균을 구한 것은 시간별 데이터가 들쭉날쭉으로 잡음이 많이 포함되어 있다고 생각하기 때문입니다. 그런데 300 이상의 PM10 농도는 사실일까? 300이상 전후의 상황을 살펴봅시다.
library(ggplot2)
source('utils.R')
#dat[, time:=as.POSIXct(format(datetime, "%H:%M:%S"), format='%H:%M:%S')]
dat2 = dat[place == "평균", ]
dat2 = dat2[order(datetime),]
dat2 = dat2[pm10 > 300 |
shift(pm10) > 300 |
shift(pm10, n=2) > 300 |
shift(pm10, n=3) > 300 |
shift(pm10, type='lead') > 300 |
shift(pm10, type='lead', n=2) > 300 |
shift(pm10, type='lead', n=3) > 300,]## Error: both operands of the expression should be "units" objects단위를 설정하면 단위가 적절하지 않은 연산은 실행이 되지 않습니다. 위의 예도 마찬가지인데요. pm10의 단위는 \(\mu g / m^3\) 이므로 단위가 없는 300 또는 1000과 비교를 하는 것은 적절하지 않습니다. 300은 set_units(300, "microgram/m^3")으로 바꿔야 하는 것이죠.
set_units(300, "microgram/m^3")## 300 [microgram/m^3]dat$pm10[1] > set_units(300, "microgram/m^3") ## [1] FALSE편의를 위해 dust()란 함수를 만들어줍시다.
dust = function(x) {
return(set_units(x, "microgram/m^3"))
}
dust(300)## 300 [microgram/m^3]이제 set_units(300, "microgram/m^3") 대신 dust(300)으로 쓸 수 있습니다. 그리고 그래프를 그립니다.
dat2 = dat2[pm10 > dust(300) |
shift(pm10) > dust(300) |
shift(pm10, n=2) > dust(300) |
shift(pm10, n=3) > dust(300) |
shift(pm10, type='lead') > dust(300) |
shift(pm10, type='lead', n=2) > dust(300) |
shift(pm10, type='lead', n=3) > dust(300),]
#dat2[, date:=format(datetime, "%Y-%m-%d")]
dat2[, date:=as.Date(datetime, tz='Asia/Seoul')]
source('utils.R')
library(ggplot2)
library(ggforce) ## Registered S3 method overwritten by 'ggforce':
## method from
## scale_type.units units# p343 : ggplot에서 unit을 사용하기 위해
ggplot(dat2, aes(x=time, y=pm10)) +
geom_line() +
geom_point(data=dat2[pm10>dust(300),], aes(x=time, y=pm10)) +
facet_wrap(~date) +
scale_x_datetime(date_labels='%H') +
x_label_vertical() # p326## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
그런데 같은 날 PM10 농도가 300이 넘는 경우가 2번 이상일 경우에는 약간 이상하네요.
dat2 = dat[place == "평균", ]
dat2 = dat2[order(datetime),]
datDates = dat2[pm10 > dust(300) |
shift(pm10) > dust(300) |
shift(pm10, n=2) > dust(300) |
shift(pm10, n=3) > dust(300) |
shift(pm10, type='lead') > dust(300) |
shift(pm10, type='lead', n=2) > dust(300) |
shift(pm10, type='lead', n=3) > dust(300), .(date=unique(date))]
dat3 <- dat2[datDates,on='date']
#dat2[, date:=format(datetime, "%Y-%m-%d")]
dat3[, date:=as.Date(datetime, tz='Asia/Seoul')]
ggplot(dat3, aes(x=time, y=pm10)) +
geom_line() +
geom_point(data=dat2[pm10>dust(300),], aes(x=time, y=pm10)) +
facet_wrap(~date) +
scale_x_datetime(date_labels='%H') +
x_label_vertical() # p326
ggplot(dat3, aes(x=time, y=pm10)) +
geom_line() +
geom_point(data=dat2[pm10>dust(300),], aes(x=time, y=pm10)) +
facet_wrap(~date) +
scale_x_datetime(date_labels='%H') +
x_label_vertical() # p326
시간의 연속성을 보아 특별히 데이터의 오류라고 보긴 힘들 것 같습니다. 지역별 데이터의 경우도 한번 확인해 볼까요? 이번에 기준을 높여서 800으로 합시다.
dat2 = dat[place != "평균", ]
dat3 = dat2[,
.SD[pm10 > dust(800) |
shift(pm10) > dust(800) |
shift(pm10, n=2) > dust(800) |
shift(pm10, n=3) > dust(800) |
shift(pm10, type='lead') > dust(800) |
shift(pm10, type='lead', n=2) > dust(800) |
shift(pm10, type='lead', n=3) > dust(800),],by=place]
source('utils.R')
library(ggplot2)
library(ggforce)
# p343 : ggplot에서 unit을 사용하기 위해
ggplot(dat3, aes(x=time, y=pm10, col=place)) +
geom_line() +
geom_point(data=dat3[pm10>dust(800),], aes(x=time, y=pm10, col=place)) +
facet_wrap(~date) +
scale_x_datetime(date_labels='%H') + # p284
x_label_vertical() # p326
여러 관측소에서 관측이 되었고 데이터의 시간적 구성으로 보아 역시 데이터의 오류라고 보긴 힘들 것 같습니다. 그런데 2021년 5월 8일, 2015년 2월 23일, 2010년 11월 12일, 그리고 2009년 12월 26일에는 어떤 일이 있었을까요? 관련 신문 기사를 살펴보았다.
성탄절 덮친 최악의 황사(한국일보, 2009-12-25/빅카인즈 검색 결과)
성탄절인 25일 서울을 비롯한 수도권과 강원, 중부 지역이 올 겨울 들어 처음 발생한 짙은 황사로 뒤덮였다. 저녁에는 서울과 중부, 호남 지역에 눈이 내려 4년 만에 화이트 크리스마스를 연출하기도 했다. 강추위가 다시 닥치는 26일에도 황사가 전국에 영향을 줄 것으로 보인다.
기상청은 이날 오후 서울 경기남부 인천 충남 등에 황사 경보를, 경기북부 강원 대전 광주 경북 등에 황사주의보를 각각 발령했다. 미세먼지 농도 최고치(1시간 평균)는 오후 11시까지 서울 963㎍/㎥, 강화 1,045㎍/㎥, 수원 1,132㎍/㎥, 천안 826㎍/㎥를 기록했다. 미세먼지 농도가 ㎥당 400㎍ 미만일 때는 '옅은 황사', 400㎍ 이상 800㎍ 미만은 '짙은 황사'로 분류된다. 800㎍ 이상은 '매우 짙은 황사'로 불린다.
최악의 가을 황사(MBC, 2010-11-12/빅카인즈 검색 결과)
최악의 가을 황사 ● 앵커： 기상청은 전국을 강타한 이번 황사가 관측을 시작한 이후 가장 강력한 가을황사로 기록됐다고 밝혔습니다.백령도는 미세먼지의 농도가 ３제곱미터당 １６６４마이크로그램을 기록했고 춘천 １１１３，서울 １１９２， 광주 １０９４마이크로그램 등으로 평상시보다 높았습니다.기상청은 황사의 발원지인 중국 북부지역이 예년보다 기온이 높고 건조해 올가을과 겨울철에 황사가 다시 날아올 가능성이 있다고 밝혔습니다.
마스크 쓴 한반도..."우리 애 어린이집 보내도 될까"(동아일보, 2015-02-24/빅카인즈 검색 결과)
23일 한반도 전체가 황사가 불러온 ‘미세먼지 포비아(공포증)’에 휩싸였다. 기상청에 따르면 이날 오전 4시 서울의 1시간 평균 미세먼지 농도는 m³당 1044μg(마이크로그램·1μg은 100만분의 1g)까지 치솟았다. 서울의 경우 관측 역사상 다섯 번째로 높았다.
서울 송파구 삼전 제2경로당은 한산했다. 박정분 씨(78·여)는 “평소 경로당에 30명 가까이 오는데 오늘은 절반도 안 왔다”며 “‘외출하지 말라’는 자녀들의 성화에 오지 않은 것 같다”고 말했다. 점심시간 서울 광화문과 여의도, 강남 등지에는 마스크를 쓴 채 종종걸음을 치는 직장인을 쉽게 볼 수 있었다. 경기 용인시와 화성시의 삼성전자 반도체 공장은 이날 생산라인에 들어갈 때 거치는 ‘에어샤워’ 시간을 평소 15초에서 2배인 30초로 늘렸다. 또 라인 내 발먼지 제거 패드의 교체 주기를 평소 3시간에서 1시간으로 단축했다.신문기사를 검색하니 기록적인 황사는 고비 사막, 편서풍, 서풍과 관련이 깊어 보입니다.
그 밖의 이상값 확인 방법
이상값을 확인하는 방법은 다음과 같이 구분해볼 수 있습니다.
- 이론적인 배경
- 시계열의 경우 값의 연속성
- 2-변수 이상에서는 변수 간 관계
앞에서 \(PM_{10}\) 이 음수가 될 수 없다는 점을 이용하여 몇몇 이상값을 제거했습니다.
시계열 1변수의 경우, 위에서와 같이 시간적 연속성을 이용할 수 있겠죠.
만약 다른 변수와의 관계를 알고 있다면 이를 이용할 수도 있죠. 예를 들어 체중 196 m와 체중 50 kg는 그 자체로 이상값이라고 할 수 없지만, 한 사람이 키 196 m에 제충 50 kg라면 조금 이상하다고 생각할 수 있을 것입니다.
만약 다른 변수들이 있다면, 그리고 그 변수와의 관계를 알고 있다면 다른 이상값을 찾아낼 수도 있을 것입니다만, 현재 \(PM_{10}\) 와 \(PM_{2.5}\) 의 관계에 아는 바가 없으니 일단 넘어갑니다. 이후에 데이터 모델링에서 이 관계를 확인한 후에 다시 이상값을 확인해 볼 수 있겠죠.
(4부에서 계속 됩니다.)
Leave a comment