대기환경 측정자료 시각화를 위한 R

국립환경과학원 대기공학과

한국외국어대학교 최용주

2022-09-16

Time series 그래프

필요한 라이브러리 설치 (일회성)

install.packages("readxl") # 엑셀파일을 읽는 라이브러리 설치
install.packages("lubridate") # 엑셀파일을 읽는 라이브러리 설치

라이브러리 불러오기

library("readxl") # 엑셀파일을 읽는 라이브러리 설치
library("lubridate") # 엑셀파일을 읽는 라이브러리 설치
## 
## 다음의 패키지를 부착합니다: 'lubridate'
## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union
library("openair") # 엑셀파일을 읽는 라이브러리 설치

작업 디렉토리 지정

setwd("D:/HUFS/Workspace/Rclass/NIER_220918/") # 설정
getwd() #확인 
## [1] "D:/HUFS/Workspace/Rclass/NIER_220918"

측정자료 (excel 파일) 불러오기

SH = read_xlsx(path = "시화산단AQMS.xlsx", sheet = 1, skip = 5)   # 시화산단 자료
## New names:
## * `` -> ...1
## * 등급 -> 등급...2
## * `1시간` -> `1시간...3`
## * 등급 -> 등급...4
## * `1시간` -> `1시간...5`
## * ...
WS = read_xlsx(path = "원시동AQMS.xlsx", sheet = 1, skip = 5)     # 원시동 자료
## New names:
## * `` -> ...1
## * 등급 -> 등급...2
## * `1시간` -> `1시간...3`
## * 등급 -> 등급...4
## * `1시간` -> `1시간...5`
## * ...
AWS = read_xlsx(path = "안산(545) AWS.xlsx", sheet = 1, skip = 1) # AWS 자료

불러온 자료의 컬럼 이름 확인 및 필요한 컬럼만 추출

colnames(SH)
##  [1] "...1"       "등급...2"   "1시간...3"  "등급...4"   "1시간...5" 
##  [6] "등급...6"   "1시간...7"  "등급...8"   "1시간...9"  "등급...10" 
## [11] "1시간...11" "등급...12"  "1시간...13"
colnames(SH[,c(seq(1,13,2))])
## [1] "...1"       "1시간...3"  "1시간...5"  "1시간...7"  "1시간...9" 
## [6] "1시간...11" "1시간...13"
SH = SH[,c(seq(1,13,2))]
WS = WS[,c(seq(1,13,2))]
colnames(SH) = colnames(WS) = c("date", "PM10", "PM2.5", "O3", "NO2", "CO", "SO2")

자료의 유형 파악

str(SH)
## tibble [600 x 7] (S3: tbl_df/tbl/data.frame)
##  $ date : POSIXct[1:600], format: "2022-07-26 01:00:00" "2022-07-26 02:00:00" ...
##  $ PM10 : chr [1:600] "22" "16" "27" "25" ...
##  $ PM2.5: chr [1:600] "10" "9" "16" "17" ...
##  $ O3   : chr [1:600] "0.020" "0.021" "0.010" "0.002" ...
##  $ NO2  : chr [1:600] "0.017" "0.017" "0.028" "0.033" ...
##  $ CO   : chr [1:600] "0.4" "0.4" "0.4" "0.5" ...
##  $ SO2  : chr [1:600] "0.001" "0.001" "0.001" "0.002" ...

모두 ’문자형’이기에 숫자형으로 바꿔줌

for (i in 2:ncol(SH)){
  SH[,i] = as.numeric(unlist(SH[,i]))
  WS[,i] = as.numeric(unlist(WS[,i]))
}

AWS 자료도 컬럼 이름을 수정

colnames(AWS)
##  [1] "지점"          "지점명"        "일시"          "기온(°C)"    
##  [5] "풍향(deg)"     "풍속(m/s)"     "강수량(mm)"    "현지기압(hPa)"
##  [9] "해면기압(hPa)" "습도(%)"
colnames(AWS) = c("site", "name", "date", "temp", "wd", "ws", "prec", "pres", "s.pres", "rh")

AWS 자료도 컬럼 이름을 수정ㅜ

SH = merge(SH, AWS, by = "date", all.x = T)
WS = merge(WS, AWS, by = "date", all.x = T)
summaryPlot(SH[,1:7])
##     date1     date2      PM10     PM2.5        O3       NO2        CO       SO2 
## "POSIXct"  "POSIXt" "numeric" "numeric" "numeric" "numeric" "numeric" "numeric"

plot(SH$PM10 ~ SH$date) # (1): 일반적인 점 그래프

plot(SH$PM10 ~ SH$date, type = "l") # (2): 일반적인 선 그래프

plot(SH$PM10 ~ SH$date, type = "l", ylim = c(0,50)) # (3): (2) + y축 범위 변경

plot(SH$PM10 ~ SH$date, type = "l", ylim = c(0,50), col = "darkorange") # (4): (3) + 색깔 변경

plot(SH$PM10 ~ SH$date, type = "l", ylim = c(0,50), col = "darkorange", xlab = "Time", ylab = "PM10") # (5): (4) + 이름 변경

plot(SH$PM10 ~ SH$date, type = "l", ylim = c(0,50), col = "darkorange", xlab = "Time", ylab = expression("PM"[10]*"  ("*mu*"g m"^-3*")")) # (6): (5) + 특수문자 및 아래첨자 

{
  plot(SH$PM10 ~ SH$date, type = "l", ylim = c(0, 70), col = "darkorange", xlab = "Time", ylab ="") 
  points(SH$PM10 ~ SH$date, col = "darkorange", pch = 22, cex = 0.7)
  mtext(text = expression("PM"[10]*" ("*mu*"g m"^-3*")"), line = 2, side = 2) # (6): (5) + 특수문자 및 아래첨자 간격 줄이기기
  
  points(WS$PM10 ~ WS$date, col = "blue", pch = 21, cex = 0.7) # (8): (7) + SO4 자료 추가
  lines(WS$PM10 ~ WS$date, col = "blue")

  legend("topright", legend = c("시화공단", "원시"), col = c("darkorange", "blue"), pch = 22, cex = 0.7) # (7): (6) + 심볼 추가 + 레전드 추가 
} 

기타 그래프

Bar 그래프

DAT = SH[SH$date >= as.POSIXct("2022-08-15", tz = "UTC") & SH$date < as.POSIXct("2022-08-18", tz = "UTC"),]
DAT$PM10_25 = DAT$PM10 - DAT$PM2.5
barplot(t(DAT[,c("PM10_25", "PM2.5")]), col = c("darkorange", "skyblue"), beside = F, names.arg = as.Date(DAT$date))

바람장미 그래프

polarplot

Scatterplot

기본 상호작용하는 지도 생성하기

패키지 설치

install.packages("leaflet")

패키지 불러오기

library("leaflet")
library("geosphere")

측정소 위경도 및 AWS와 거리 계산

lon_nier = c(126.7272389, 126.7885028 , 126.8418193)
lat_nier = c(37.3351039, 37.3052549, 37.2802273)

loc = data.frame(group = c("A", "B", "A", "B"),
                    long = c(lon_nier[c(1,3)], lon_nier[c(2,3)]),
                    lat = c(lat_nier[c(1,3)], lat_nier[c(2,3)]))

dist.sh = distm(x = c(lon_nier[1], lat_nier[1]), y = c(lon_nier[3], lat_nier[3]))/1000
dist.ws = distm(x = c(lon_nier[2], lat_nier[2]), y = c(lon_nier[3], lat_nier[3]))/1000
{
  leaflet() %>%
  addTiles() %>%
  setView( lng= mean(lon_nier), lat = mean(lat_nier), zoom = 13) %>%
  addProviderTiles("Esri.WorldImagery") %>%
  addMarkers(lng = 126.7272389, 
             lat = 37.3351039, 
             popup = paste("시화공단 </br>도시대기측정망</br> AWS 기준", round(dist.sh, 1), "km")) %>%
  addMarkers(lng = 126.7885028, 
             lat = 37.3052549, 
             popup = paste("원시동 </br>도시대기측정망 </br> AWS 기준", round(dist.ws, 1), "km")) %>%
  addMarkers(lng = 126.8418193, 
             lat = 37.2802273, 
             popup = "AWS 측정소</br>") %>%
  addPolylines(data = loc, lng = ~long[1:2], lat = ~lat[1:2], col = "red") %>%
  addPolylines(data = loc, lng = ~long[3:4], lat = ~lat[3:4], col = "blue")
}

지리학적 분석

무인드론에 설치된 OPC (optical particle couneter) 자료를 분석해보자

필요한 라이브러리 설치 (일회성)

install.packages("readxl") # 엑셀파일을 읽는 라이브러리 설치
install.packages("maps")   # google map 을 불러오는 라이브러리 설치
install.packages("squash") # color bar를 줄 수 있는 라이브러리 설치
install.packages("plyr")   # 다소 복잡한 통계 함수를 위한 라이브러리 설치

필요한 라이브러리 호출

library("readxl")
library("ggmap")
## 필요한 패키지를 로딩중입니다: ggplot2
## Google's Terms of Service: https://cloud.google.com/maps-platform/terms/.
## Please cite ggmap if you use it! See citation("ggmap") for details.
library("squash")
library("plyr")

작업 디렉토리 지정

setwd("D:/HUFS/Workspace/Rclass/NIER_220918/") # 설정
getwd() #확인 
## [1] "D:/HUFS/Workspace/Rclass/NIER_220918"

측정자료 (excel 파일) 불러오기

DAT = read_xlsx(path = "opc_data.xlsx", sheet = 1, skip = 1)

불러온 자료 확인

head(DAT)
## # A tibble: 6 x 14
##   일자_텍스트 시간_텍스트 코드   지역  OPC_자료사용 일자_오전오후_~    TSP  PM10
##   <chr>       <chr>       <chr>  <chr> <chr>        <chr>            <dbl> <dbl>
## 1 2022-07-26  9:55:22     2022-~ 1구역 X            2022-07-26 오전~  50.6  49.9
## 2 2022-07-26  9:55:28     2022-~ 1구역 X            2022-07-26 오전~  68.5  61.3
## 3 2022-07-26  9:55:34     2022-~ 1구역 X            2022-07-26 오전~  53.1  52.7
## 4 2022-07-26  9:55:40     2022-~ 1구역 X            2022-07-26 오전~  66.6  48.7
## 5 2022-07-26  9:55:46     2022-~ 1구역 X            2022-07-26 오전~  86.3  67.1
## 6 2022-07-26  9:55:52     2022-~ 1구역 X            2022-07-26 오전~  65.3  56.7
## # ... with 6 more variables: PM4 <dbl>, PM25 <dbl>, PM1 <dbl>,
## #   OPC_Latitude <dbl>, OPC_Longitude <dbl>, OPC_Height <dbl>

자료의 컬럼별 통계값 확인

summary(DAT)
##  일자_텍스트        시간_텍스트            코드               지역          
##  Length:45402       Length:45402       Length:45402       Length:45402      
##  Class :character   Class :character   Class :character   Class :character  
##  Mode  :character   Mode  :character   Mode  :character   Mode  :character  
##                                                                             
##                                                                             
##                                                                             
##                                                                             
##  OPC_자료사용       일자_오전오후_시간      TSP               PM10        
##  Length:45402       Length:45402       Min.   :   5.60   Min.   :   5.60  
##  Class :character   Class :character   1st Qu.:  19.20   1st Qu.:  19.20  
##  Mode  :character   Mode  :character   Median :  25.95   Median :  25.70  
##                                        Mean   :  36.97   Mean   :  28.38  
##                                        3rd Qu.:  37.23   3rd Qu.:  36.10  
##                                        Max.   :4477.70   Max.   :1135.10  
##                                        NA's   :36862     NA's   :36862    
##       PM4               PM25             PM1          OPC_Latitude  
##  Min.   :   5.60   Min.   :  5.50   Min.   :  4.60   Min.   :37.29  
##  1st Qu.:  18.80   1st Qu.: 17.70   1st Qu.: 15.00   1st Qu.:37.31  
##  Median :  25.10   Median : 24.00   Median : 21.90   Median :37.33  
##  Mean   :  27.12   Mean   : 25.36   Mean   : 22.43   Mean   :37.33  
##  3rd Qu.:  34.40   3rd Qu.: 31.90   3rd Qu.: 28.40   3rd Qu.:37.34  
##  Max.   :1123.20   Max.   :963.60   Max.   :308.70   Max.   :37.35  
##  NA's   :36862     NA's   :36862    NA's   :36862    NA's   :37258  
##  OPC_Longitude     OPC_Height    
##  Min.   :126.7   Min.   :-19.00  
##  1st Qu.:126.7   1st Qu.: 37.00  
##  Median :126.7   Median :105.00  
##  Mean   :126.7   Mean   : 90.92  
##  3rd Qu.:126.8   3rd Qu.:126.00  
##  Max.   :126.8   Max.   :217.00  
##  NA's   :37258   NA's   :37258

PM2.5 자료릐 NA값이 아닌 자료만을 취합함 (엑셀 필터 기능)

DAT = DAT[is.na(DAT$PM25) == F, ]

단순한 그래프화

plot(DAT$OPC_Latitude ~ DAT$OPC_Longitude) # (1) y ~ x  형태의 일반적인 점 그래프

plot(DAT$OPC_Latitude ~ DAT$OPC_Longitude, xlab = "Longitude", ylab = "Latitude") # (2): (1) + 축 추가

plot(DAT$OPC_Latitude ~ DAT$OPC_Longitude, xlab = "Longitude", ylab = "Latitude", ylim = c(37.29, 37.36)) # (3): (2) + y축 범위 변경

plot(DAT$OPC_Latitude ~ DAT$OPC_Longitude, xlab = "Longitude", ylab = "Latitude", ylim = c(37.29, 37.36), col = "darkorange") # (4): (3) + 색깔 변경

plot(DAT$OPC_Latitude ~ DAT$OPC_Longitude, ylim = c(37.29, 37.36), col = "darkorange", ylab = expression("Latitude ("*degree*"N)"), xlab = expression("Longitude ("*degree*"E)")) # (5): (4) + 특수부호

단순한 그래프화 (축 이름 별도 지정)

{
  plot(DAT$OPC_Latitude ~ DAT$OPC_Longitude, ylim = c(37.29, 37.36), col = "darkorange", ylab = "", xlab = "") 
  mtext(text = expression("Latitude ("*degree*"N)"), line = 2, side = 2) # y축 추가
  mtext(text = expression("Longitude ("*degree*"E)"), line = 2, side = 1) # x축 추가
}

PM2.5 농도에 따른 color bar

PM2.5의 농도 범위 파악

(zlim = range(DAT$PM25, na.rm = T)) # 결과를 zlim에 저장
## [1]   5.5 963.6

최고농도의 범위를 90퍼센타일 값으로 변경

(zlim[2] = round_any(quantile(DAT$PM25, 0.90, na.rm  = T), accuracy = 1, f = trunc))
## 90% 
##  39
col.map = makecmap(zlim,  n = 100, colFn = jet, outlier = 'red') # 팔레트 제작 jet, heat, bluered, blueorange, rainbow, greyscale, coolheat
DAT$PM25.1 = DAT$PM25
DAT[DAT$PM25 > zlim[2] & is.na(DAT$PM25) == F, "PM25.1"] = zlim[2]

{
  plot(DAT$OPC_Latitude ~ DAT$OPC_Longitude, type = "p", xlab = "", ylab = "",  # 위경도 그림
      col = cmap(DAT$PM25.1, col.map), main = expression("PM"[2.5]*" ("*mu*"g m"^-3*")")) # ㅒORG 농도에 따른 심볼 색 추가
  mtext(text = expression("Latitude ("*degree*"E)"), line = 2, side = 2) # y축 이름
  mtext(text = expression("Longtitude ("*degree*"E)"), line = 2, side = 1) # x축 이름
  vkey(col.map, expression("PM"[2.5]),y = 37.3, x = 126.81, stretch = 0.2) # 레전드 추가
} 

구글 지도 위에 표시

key = "AIzaSyDGUtMAIyCV7ejSian3e3M3NgFTK4rr8IU" #본인이 할당 받은 구글 키를 입력 (참고: https://duopix.co.kr/google-map-key/)
register_google(key)
loc = c(mean(DAT$OPC_Longitude, na.rm = T), mean(DAT$OPC_Latitude, na.rm = T)) # 비행 경로의 평균 위경도를 구글 지도의 중심점으로 잡음. MA가 있기에 na.rm = T라는 명령어를 넣어줌
qmap( loc , zoom = 12, maptype = "satellite") + # 중심점을 기준으로 zoom 값으로 지도를 확대 (큰 수) 혹은 축소 (작은 수). maptype은 satellite, terrain, hybrid, roadmap가 있음
  geom_point(aes(x = OPC_Longitude, y = OPC_Latitude), data = DAT, col = cmap(DAT$PM25.1, col.map))
Google Map을 사용한 위경도 그림

Google Map을 사용한 위경도 그림

qmap( loc , zoom = 12, maptype = "terrain") + 
  geom_point(aes(x = OPC_Longitude, y = OPC_Latitude), data = DAT, col = cmap(DAT$PM25.1, col.map))
Google Map을 사용한 위경도 그림

Google Map을 사용한 위경도 그림

qmap( loc , zoom = 12, maptype = "hybrid") + # hybrid
  geom_point(aes(x = OPC_Longitude, y = OPC_Latitude), data = DAT, col = cmap(DAT$PM25.1, col.map))
Google Map을 사용한 위경도 그림

Google Map을 사용한 위경도 그림

qmap( loc , zoom = 12, maptype = "roadmap") + # roadmap
  geom_point(aes(x = OPC_Longitude, y = OPC_Latitude), data = DAT, col = cmap(DAT$PM25.1, col.map))
Google Map을 사용한 위경도 그림

Google Map을 사용한 위경도 그림

3차원 그래프 완성

library("rgl") # OPEN GL (open graphics library) 2D, 3D 렌더링을 해주는 라이브러리
par3d(windowRect = c(10, 10, 800, 800)) # 3D 윈도우 생성
plot3d(DAT$OPC_Longitude, DAT$OPC_Latitude, DAT$OPC_Height, col = cmap(DAT$PM25.1, col.map), 
       size = 2, type = 'p', xlab = "Lon", ylab = "Lat", zlab = "Alt", zlim = c(0,300))
movie3d(spin3d(axis = c(0, 0, 0.5), rpm = 5), duration = 20, dir = getwd())
Google Map을 사용한 위경도 그림

Google Map을 사용한 위경도 그림

DAT$alt_bin2 = factor(DAT$OPC_Height, levels = as.character(seq(0,250,10)))

par(mai=c(0.7,0.7,0.5,0.1))
boxplot(DAT$PM25 ~ DAT$alt_bin2, horizontal = T, border = "firebrick3", xlab = "", col = "white", 
        pch=21, ylab="", cex = 0.7, ylim = c(0, 80), main = "OPC @시화공단")
mtext(expression("Altitude (m)"), side = 2, line = 2)
mtext(expression("PM"[2.5]*" ("*mu*"g/m"^3*")"), side = 1, line = 2)

abline(h = 3.5, lty = 2, col = "forestgreen")
abline(h = 6.5, lty = 2, col = "forestgreen")