airkorea 자료 처리

>여러 측정 지점의 자료를 처리하고 지도위에 표출해보자!

> 작업공간 (디렉토리)을 설정해줌.

setwd("D:/HUFS/Workspace/Rclass/CC_2023/")

> airkorea 사이트에서 제공하는 2022년 12월 자료 다운받기

> 엑셀파일 자료 읽기

‘readxl’ 패키지가 필요! 없다면 install.package 고고!

install.packages("readxl")

‘readxl’ 패키지가 있다면 불러오기!

library(readxl)
library(lubridate)
## 
## 다음의 패키지를 부착합니다: 'lubridate'
## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union
argument 설명
path R에서 읽어올 외부 데이터가 있는 파일의 위치(경로)와 파일명(확장자 포함)을 문자형으로 지정한다.
sheet 엑셀 파일의 있는 시트(sheet) 중에서 어떤 시트인지를 지정하는 것으로 시트명을 알고 있으면 문자형으로 지정하며, 시트의 위치(index)를 알면 수치형으로 지정하면 된다.
col_names 논리형으로 엑셀에 있는 변수명을 사용할 지를 의미한다. TRUE이면 엑셀에 있는 변수명을 사용한다 
DAT = read_excel(path = “파일명”, sheet = “Sheet1”, col_names = TRUE)
DAT = read_excel(path = “파일명”, sheet = “Sheet1”, col_names = TRUE)

파일 읽어오기

nier_aqm = read_excel(path = "2022_12.xlsx", sheet = "Sheet1", col_names = TRUE)
nier_site = read_excel(path = "nier_site.xlsx", sheet = "Sheet1", col_names = TRUE, na = "NA")

파일 형식 살펴보기

str(nier_aqm)
## tibble [480,622 × 12] (S3: tbl_df/tbl/data.frame)
##  $ sido   : chr [1:480622] "서울 중구" "서울 중구" "서울 중구" "서울 중구" ...
##  $ type   : chr [1:480622] "도시대기" "도시대기" "도시대기" "도시대기" ...
##  $ code   : num [1:480622] 111121 111121 111121 111121 111121 ...
##  $ name   : chr [1:480622] "중구" "중구" "중구" "중구" ...
##  $ date   : num [1:480622] 2.02e+09 2.02e+09 2.02e+09 2.02e+09 2.02e+09 ...
##  $ SO2    : num [1:480622] 0.003 0.003 0.003 0.003 0.003 0.003 0.003 0.003 0.003 0.003 ...
##  $ CO     : num [1:480622] 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.4 0.5 0.4 ...
##  $ O3     : num [1:480622] 0.029 0.03 0.028 0.03 0.029 0.028 0.02 0.009 0.01 0.013 ...
##  $ NO2    : num [1:480622] 0.006 0.006 0.006 0.005 0.006 0.008 0.015 0.026 0.029 0.026 ...
##  $ PM10   : num [1:480622] 15 12 14 15 15 16 14 21 25 24 ...
##  $ PM25   : num [1:480622] 6 3 4 6 7 6 7 7 8 9 ...
##  $ address: chr [1:480622] "서울 중구 덕수궁길 15" "서울 중구 덕수궁길 15" "서울 중구 덕수궁길 15" "서울 중구 덕수궁길 15" ...
str(nier_site)
## tibble [882 × 8] (S3: tbl_df/tbl/data.frame)
##  $ sido   : chr [1:882] "서울" "서울" "서울" "서울" ...
##  $ city   : chr [1:882] "서울" "서울" "서울" "서울" ...
##  $ code   : num [1:882] 111121 111122 111122 111123 111124 ...
##  $ name   : chr [1:882] "중구" "한강대로" "한강대로" "종로구" ...
##  $ address: chr [1:882] "중구 덕수궁길 15 (시청서소문별관 3동)" "용산구 한강대로 405 서울역앞" "용산구 한강대로 405 서울역앞" "종로구 종로35가길 19 (종로5.6가동 주민센터)" ...
##  $ lon    : num [1:882] 127 127 127 127 127 ...
##  $ lat    : num [1:882] 37.6 37.5 37.5 37.6 37.6 ...
##  $ type   : chr [1:882] "도시대기측정망" "도로변대기측정망" "유해대기물질측정망" "도시대기측정망" ...
unique(nier_site$type)
##  [1] "도시대기측정망"           "도로변대기측정망"        
##  [3] "유해대기물질측정망"       "PM2.5성분측정망"         
##  [5] "대기중금속측정망"         "산성강하물측정망"        
##  [7] "종합대기측정소"           "광화학대기오염물질측정망"
##  [9] "교외대기측정망"           "항만측정망"              
## [11] "국가배경농도측정망"       "대기환경연구소"

만약 한글이 깨진다면?

Sys.getlocale()
Sys.setlocale("LC_ALL","C") # 강제 언어 삭제
unique(nier_site$type)

측정소 코드를 기준으로 위경도 자료 합치기 (도시대기와 도로변 만)

nier_site = nier_site[nier_site$type %in% c("도시대기측정망", "도로변대기측정망"),]
nier_aqm = nier_aqm[is.na(nier_aqm$PM10) == F & is.na(nier_aqm$PM25) == F,]

nier_aqm_latlon = merge(nier_aqm, nier_site[,c("code", "lon", "lat")], by = "code")

> 위경도를 기준으로 측정소 위치 표시하기

지도를 그리려면 ‘maps’ 패키지가 필요! 없다면 install.package 고고!

install.packages("maps")

‘maps’ 패키지가 있다면 불러오기!

library(maps)

그림 그리기

{plot(nier_aqm_latlon$lat ~ nier_aqm_latlon$lon) # 위경도에 따른 심볼 그림
map("world", add = T) # 지도 추가
}

이제 꾸며보기!!

{plot(nier_aqm_latlon$lat ~ nier_aqm_latlon$lon, xlab = "", ylab = "") # 위경도에 따른 심볼 그림
mtext(text = expression("Latitude ("*degree*"E)"), line = 2, side = 2) # y축 이름
mtext(text = expression("Longtitude ("*degree*"E)"), line = 2, side = 1) # x축 이름
map("world", add = T) # 지도 추가
}

측정소별 일평균 구하기

nier_aqm_latlon$date = as.POSIXct(as.character(nier_aqm_latlon$date), format = "%Y%m%d%H") - hours(1)
d.nier_aqm_latlon = aggregate(nier_aqm_latlon[6:11], by = list(format(nier_aqm_latlon$date, "%Y-%m%-%d"), nier_aqm_latlon$code), mean, na.rm = T)
colnames(d.nier_aqm_latlon)[1:2] = c("date", "ID")
d.nier_aqm_latlon$date = as.POSIXct(d.nier_aqm_latlon$date)

서울 자료만 추출해서 일별 그래프 그리기

d.seoul = d.nier_aqm_latlon[substr(d.nier_aqm_latlon$ID, 1, 2) == 11, ]
station = unique(d.seoul$ID)
{
  plot(PM25 ~ date, data = d.seoul[d.seoul$ID == station[1],], type = "l", col = "red", xlab = "", ylab = "", ylim = c(0,80))
  points(PM25 ~ date, data = d.seoul[d.seoul$ID == station[1],], col = "red")
  mtext(text = expression("Date"), line = 2, side = 1) # x축 이름
  mtext(text = expression("PM"[2.5]*" ("*mu*"/m"^3*")"), line = 2, side = 2) # y축 이름
  
  for (i in 2:40){
    lines(PM25 ~ date, data = d.seoul[d.seoul$ID == station[i],], col = i)
    points(PM25 ~ date, data = d.seoul[d.seoul$ID == station[i],], col = i)
  }
  legend("topleft", legend = station, fill = 1:40, ncol = 7, cex = 0.7)
}

자동으로 색깔을 지정해서 그려보기!

my_col = rainbow(40)
{
  plot(PM25 ~ date, data = d.seoul[d.seoul$ID == station[1],], type = "l", col = my_col[1], xlab = "", ylab = "", ylim = c(0,80))
  points(PM25 ~ date, data = d.seoul[d.seoul$ID == station[1],], col = my_col[1])
  mtext(text = expression("Date"), line = 2, side = 1) # x축 이름
  mtext(text = expression("PM"[2.5]*" ("*mu*"/m"^3*")"), line = 2, side = 2) # y축 이름
  
  for (i in 2:32){
    lines(PM25 ~ date, data = d.seoul[d.seoul$ID == station[i],], col = my_col[i])
    points(PM25 ~ date, data = d.seoul[d.seoul$ID == station[i],], col = my_col[i])
  }
  legend("topleft", legend = station, fill = my_col[1:40], ncol = 7, cex = 0.7)
}

서울지역 측정소별 boxplot

boxplot(d.seoul$PM25 ~ d.seoul$ID)

서울지역 측정소별 boxplot 중간값 순위로 정렬

group_ordered = with(d.seoul, reorder(ID, PM25, median, na.rm = T))
# group_ordered
boxplot(d.seoul$PM25 ~ group_ordered)