대기환경학자를 위한 R (실전편 - openair 1)
이화여자대학교 2024
한국외국어대학교 최용주
2024-02-27
1차 정리된 파일 불러와서 자료 가공하기
기본 설정
setwd("~/Rclass/")설정된 작업 디렉토리를 확인하는 방법은 ‘getwd()’ 입니다.
getwd()## [1] "/Users/ychoi/Library/CloudStorage/OneDrive-hufs.ac.kr/Workspace_hufs/Rclass/Ehwa_240227"위와 같이 “D:/Rclass/”가 작업 디렉토리로 설정된 것을 볼 수 있습니다.
다음으로는 분석에 필요한 라이브러리를 로드해줍니다.
library("lubridate")##
## Attaching package: 'lubridate'## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, unionlibrary("plyr")aero.meteo.CSV 불러오기
DAT = read.csv(file = “파일명”, header = TRUE)aero.meteo = read.csv(file = "sample_data.csv", header = TRUE)
aero.meteo$Time = as.POSIXct(aero.meteo$Time, tz = "Asia/Seoul")
aero.meteo$Time[1]## [1] "2020-12-15 KST"그래프 그리기
- 기본적으로 그래프는 plot 명령어를 사용하여 그린다.
| argument | 설명 |
|---|---|
| x, y | 그래프에 사용될 x와 y값. 두 값의 갯수는 동일해야함. |
| type | 그래프의 형태 “p”는 점, “l”은 선, “b”는 점+선, “c” 결측치를 건너뛰는 선. |
| xlim, ylim | x축과 y축의 범위 (x1, x2) 혹은 (y1, y2). x1 > x2이면 감소하는 형태의 축. 기본설정은 자료 범위에 맞춰서 정해짐 |
| log | 로그축으로 변경, log = “x” 혹은 log = “y” 혹은 log = “xy” |
| main | 그림의 제목 |
| sub | 그림의 소제목 |
| xlab, ylab | x, y축의 라벨 |
| axes | 축을 표현하지 않음 “xaxt =” or “yaxt” to suppress just one of the axes. |
| col | 그래프 색깔 변경, “red”, “blue” |
| pch | 그래프 심볼, 1, 2, 3, …, 25 |
| cex | 크기, 실수 (real)로 입력. 숫자가 높을 수록 크기가 커짐 |
plot(aero.meteo$PM2.5 ~ aero.meteo$Time) # (1): 일반적인 점 그래프
plot(aero.meteo$PM2.5 ~ aero.meteo$Time, type = "l") # (2): 일반적인 선 그래프
plot(aero.meteo$PM2.5 ~ aero.meteo$Time, type = "l", ylim = c(0,100)) # (3): (2) + y축 범위 변경
plot(aero.meteo$PM2.5 ~ aero.meteo$Time, type = "l", ylim = c(0,100), col = "darkorange") # (4): (3) + 색깔 변경
plot(aero.meteo$PM2.5 ~ aero.meteo$Time, type = "l", ylim = c(0,100), col = "darkorange", xlab = "Time", ylab = "PM2.5") # (5): (4) + 이름 변경
plot(aero.meteo$PM2.5 ~ aero.meteo$Time, type = "l", ylim = c(0,100), col = "darkorange", xlab = "Time", ylab = expression("PM"[2.5]*" ("*mu*"g m"^-3*")")) # (6): (5) + 특수문자 및 아래첨자
{plot(aero.meteo$PM2.5 ~ aero.meteo$Time, type = "l", ylim = c(0,100), col = "darkorange", xlab = "Time", ylab = expression("PM"[2.5]*" ("*mu*"g m"^-3*")"))
points(aero.meteo$PM2.5 ~ aero.meteo$Time, col = "darkorange", pch = 22, cex = 0.7)
legend("topright", legend = "Organic", col = "darkorange", pch = 22, cex = 0.7)} # (7): (6) + 심볼 추가 + 레전드 추가 
{plot(aero.meteo$NO3. ~ aero.meteo$Time, type = "l", ylim = c(0,40), col = "darkorange", xlab = "Time", ylab = expression("PM"[1]*" Organic ("*mu*"g m"^-3*")"))
points(aero.meteo$NO3. ~ aero.meteo$Time, col = "darkorange", pch = 22, cex = 0.7)
points(aero.meteo$SO42. ~ aero.meteo$Time, col = "blue", pch = 21, cex = 0.7) # (8): (7) + SO4 자료 추가
lines(aero.meteo$SO42. ~ aero.meteo$Time, col = "blue")
legend("topright", legend = c("Nitrate", "Sulfate"), col = c("darkorange", "blue"), pch = c(22, 21), cex = 0.7)} 
- y축라벨과 숫자의 간격을 줄이는 방법
{plot(aero.meteo$NO3. ~ aero.meteo$Time, type = "l", ylim = c(0,40), col = "darkorange", xlab = "Time", ylab = "")
mtext(text = expression("PM"[2.5]*" Nitrate ("*mu*"g m"^-3*")"), line = 2, side = 2)} # (6): (5) + 특수문자 및 아래첨자 간격 줄이기기
* 그래프 저장하는 방법: jpeg, bmp, tiff, pdf등 여러 형태가 가능하나
png를 주로 사용함.
| argument | 설명 |
|---|---|
| filename | 저장할 파일 이름 |
| width | 그림의 폭 |
| height | 그림의 높이 |
| units | 그림의 폭과 높이에 사용된 단위; px (pixels, the default), in (inches), cm 또는 mm. |
| res | 해상도. 보통 200 이상의 값을 사용. 높을 수록 높은 품질 |
| pointsize | 심볼의 크기 |
| bg | 배경색 |
{
png(filename = "ts_no3.png", width = 5, height = 3, units = "in", res = 200) # 저장할 스케치북 열기
plot(aero.meteo$NO3. ~ aero.meteo$Time, type = "l", ylim = c(0,40), col = "darkorange", xlab = "Time", ylab = "")
mtext(text = expression("PM"[2.5]*" Nitrate ("*mu*"g m"^-3*")"), line = 2, side = 2)
dev.off() # 스케치북 닫기
} # (7): (6) 그래프 저장## quartz_off_screen
## 2산점도 및 추세선 추가
- lm (linear model): 선형 회귀직선
- abline: 선 추가 명령어
sl = lm(aero.meteo$NO3. ~ aero.meteo$NH4.) # y ~ x의 방정식 (y = bx + a)
(sl2 = summary(sl)) # sl 변수의 종합적 통계##
## Call:
## lm(formula = aero.meteo$NO3. ~ aero.meteo$NH4.)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.1004 -0.2690 0.0590 0.2844 2.1140
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.22913 0.03732 -6.14 1.33e-09 ***
## aero.meteo$NH4. 2.33457 0.01056 220.98 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7215 on 765 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.9846, Adjusted R-squared: 0.9846
## F-statistic: 4.883e+04 on 1 and 765 DF, p-value: < 2.2e-16{plot(aero.meteo$NO3. ~ aero.meteo$NH4., type = "p", col = "darkorange", xlab = "", ylab = "") #ORG와 SO4 산점도. 축 라벨은 없음
mtext(text = expression("PM"[2.5]*" Nitrate ("*mu*"g m"^-3*")"), line = 2, side = 2) # y축 라밸 추가
mtext(text = expression("PM"[2.5]*" Ammonium" *" ("*mu*"g m"^-3*")"), line = 2, side = 1) # x축 라밸 추가
abline(sl, col = "red", lwd = 2) # 추세선 추가
abline(a = 0, b = 1, col = "grey50", lty = 2) # 1:1 라인 추가 (lty는 선의 형태)
mtext(text = paste("R = ", round(sqrt(sl2$r.squared), 2)), line = 0, side = 3, adj = 1) # 상관계수 (R) 표기
mtext(text = paste("Y = ", round(sl$coefficients[2], 2), "X+", round(sl$coefficients[1], 2)), line = -1, side = 3, adj = 1) # 추세식 표기
} 
Openair
openair는 대기 자료 분석에 가장 많이 사용되는 패키지로써, 대기환경에서 가장 많이 사용되는 패키지중 하나이다. 특히 복잡한 통계 및 그래프를 짧은 명령어로 충분히 실행되기에 초보자에게도 진입장벽이 높지 않다.
openair 패키지 설치: openair가 처음이라면 아래 명령어를 사용하여 설치할 수 있다.
install.packages("openair")- 패키지 불러오기: 기존에 설치가 되었거나, 설치가 완료되었다면, 아래 명령어로 lubridate 패키지와 함께 불러온다.
library("lubridate")
library("openair")summary plot
- 시계열 자료의 전체적인 경향을 파악하기 좋은 그래프. 자료 및 결측치의 분포, 평균, 중간값 등에 대한 정보를 제공함함.
- 중요한 것은 시간은 date, 풍속은 ws, 풍향은 wd로 인식함
aero.meteo$date = aero.meteo$Time # date로 시간 컬럼을 변경
summaryPlot(aero.meteo) # 불광동의 summary plot
summaryPlot(aero.meteo, percentile = 0.95) # 상위 5% 값은 제거
summaryPlot(aero.meteo, na.len = 10) # 적어도 10개이상 연속적인 결측치를 보여줌
summaryPlot(aero.meteo, col.data = "green") # 자료 색깔을 녹색으로
summaryPlot(aero.meteo, col.mis = "yellow") # 결측치를 노란색으로
summaryPlot(aero.meteo, col.dens = "black") # density plot선을 검은색으로
summaryPlot(aero.meteo[, c(6:16, 26)]) # 컬럼 2번째와 5-7번째만 보여주기, date 컬럼은 있어야함!!
summaryPlot(subset(aero.meteo, select = c(date, PM2.5, OC, EC))) # 앞선 명령어의 대안 windrose
- 풍향, 풍속에 대한 정보를 제공
windRose(aero.meteo) # 기본 바랑장미
windRose(aero.meteo, type = "month", layout = c(2, 1)) # 월별 바랑장미
windRose(aero.meteo[is.na(aero.meteo$PM2.5) == F,], type = "PM2.5", layout = c(4, 1)) # PM2.5 구간에 따른 바람장비
### pollution rose
- 오염장미
pollutionRose(aero.meteo, pollutant = "PM2.5") # PM2.5 농도 구간별 풍향의 빈도 수
pollutionRose(aero.meteo, pollutant = "PM2.5", normalise = TRUE) # normalized된 PM2.5 농도 구간별 풍향의 빈도 수
Polar plot
- 풍향/풍속에 따른 오염물질의 평균 농도를 보여줌
polarPlot(aero.meteo, pollutant = "PM2.5") 
polarPlot(aero.meteo, pollutant = "PM2.5", statistic = "cpf", percentile = 75) # Conditional Probability Function
polarPlot(aero.meteo, pollutant = "PM2.5", statistic = "cpf", percentile = c(0,10)) # Conditional Probability Function
polarPlot(aero.meteo, pollutant = "PM2.5", statistic = "cpf", percentile = c(50,60)) # Conditional Probability Function
polarPlot(aero.meteo, pollutant = "PM2.5", type = "daylight") # ”year”, ”hour”, ”month”, ”season”, ”weekday”, ”weekend”, ”monthyear”, ”daylight”
Calendar plot
- 달력에 원하는 오염물질의 농도등을 표시하는 그림. 보다 직관적이기에 발표자료에 많이 사용됨.
calendarPlot(aero.meteo, pollutant = "PM2.5") # 불광 PM2.5 농도 기준, 표시
Trend 분석 (smoothTrend, TheilSen)
smoothTrend(aero.meteo, pollutant = "PM2.5", deseason = TRUE, simulate = F, ylab = "concentration (ppb)",
main = "monthly mean deseasonalised NO2 (bootstrap uncertainties)") # smooth trend를 계산. GAM 모델을 사용하여 smooth를 수행
smoothTrend(aero.meteo, pollutant = "PM2.5", deseason = TRUE, simulate = F, ylab = "concentration (ppb)",
main = "monthly mean deseasonalised NO2 (bootstrap uncertainties)", type = "wd") # 풍향을 기준으로 분류## Warning: 2 missing wind direction line(s) removed
TheilSen(aero.meteo, pollutant = "PM2.5", deseason = TRUE, ylab = "PM2.5", main = "Theil-Sen Trend") # Theil-sen 기울기 계산산## Taking bootstrap samples. Please wait.
TheilSen(aero.meteo, pollutant = "PM2.5", deseason = TRUE, ylab = "PM2.5", main = "Theil-Sen Trend", type = "season") # 계절절을 기준으로 분류## Taking bootstrap samples. Please wait.
### timeVariation 분석
- 월별, 시간별, 요일별 분석을 수행 가능
# 기본 그림
timeVariation(aero.meteo, pollutant = "PM2.5")
# 여러 항목들의 시간변화 비교 (절대값이 다르기에 nomarlized 적용)
timeVariation(aero.meteo, pollutant = c("PM2.5", "OC", "EC"), normalise = TRUE)
그 외 유용한 함수
scatterPlot(aero.meteo, x = "OC", y = "EC", method = "density", col = "jet", type = "season")## (loaded the KernSmooth namespace)
aero.meteo = rollingMean(aero.meteo, pollutant = "PM2.5", hours = 8, new.name = "PM2.5.8", data.thresh = 75) # 8시간 이동 평균trendLevel(aero.meteo, pollutant = "PM2.5", cols = "jet", y = "hour")