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)

서울지역 측정소별 boxplot 꾸미기

levels(group_ordered)

##  [1] "111181" "111122" "111282" "111191" "111142" "111221" "111125" "111262"
##  [9] "111171" "111261" "111212" "111231" "111213" "111152" "111251" "111141"
## [17] "111291" "111131" "111312" "111263" "111281" "111301" "111275" "111242"
## [25] "111161" "111273" "111232" "111264" "111124" "111201" "111151" "111154"
## [33] "111241" "111311" "111274" "111202" "111162" "111121" "111143" "111123"

x.label = as.numeric(levels(group_ordered))
i=1
for (i in 1:length(x.label)){ 
  loc = which(x.label[i] == nier_site$code)
  x.label[i] = nier_site$name[loc]
}

##### 한줄로 표현
x.label = as.numeric(levels(group_ordered))
for (i in 1:length(x.label)) x.label[i] = unlist(nier_site[which(x.label[i] == nier_site$code), "name"])

par(mai=c(0.7,0.8,0.1,0.1)) #b, l, t, r
boxplot(d.seoul$PM25 ~ group_ordered, xlab = "", ylab = expression("PM"[2.5]) , col="#69b3a2", boxwex=0.4 , main="", xaxt="n")
text(cex=1, x=(1:40)+0.5, y= 2, x.label, xpd=TRUE, srt=90, pos=2)

서울지역 측정소별 월평균 계산

m.seoul = aggregate(d.seoul[,3:8], by = list(d.seoul$ID, format(d.seoul$date, "%Y-%m-01")), median, na.rm = T)
colnames(m.seoul)[1:2] = c("code", "date")
m.seoul$date = as.POSIXct(m.seoul$date, tz = "Asia/Seoul")

서울지역 측정소 월평균과 측정소 위경도 병합

m.seoul_latlon = merge(m.seoul, nier_site[,c("code", "lon", "lat")], by = "code")

서울지역 측정소 그리기

{plot(m.seoul_latlon$lat ~ m.seoul_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) # 지도 추가
}

서울 shapefile 가져오기 서울 shp 다운로드

install.packages("rgdal")

library("rgdal")

## 필요한 패키지를 로딩중입니다: sp

## Please note that rgdal will be retired during 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
## See https://r-spatial.org/r/2022/04/12/evolution.html and https://github.com/r-spatial/evolution
## rgdal: version: 1.6-5, (SVN revision 1199)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.5.2, released 2022/09/02
## Path to GDAL shared files: C:/Users/YChoi/AppData/Local/R/win-library/4.2/rgdal/gdal
## GDAL binary built with GEOS: TRUE 
## Loaded PROJ runtime: Rel. 8.2.1, January 1st, 2022, [PJ_VERSION: 821]
## Path to PROJ shared files: C:/Users/YChoi/AppData/Local/R/win-library/4.2/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.6-0
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.

Seoul = readOGR("Seoul.shp")

## Warning: OGR support is provided by the sf and terra packages among others

## Warning: OGR support is provided by the sf and terra packages among others

## Warning: OGR support is provided by the sf and terra packages among others

## Warning: OGR support is provided by the sf and terra packages among others

## Warning: OGR support is provided by the sf and terra packages among others

## Warning: OGR support is provided by the sf and terra packages among others

## Warning: OGR support is provided by the sf and terra packages among others

## OGR data source with driver: ESRI Shapefile 
## Source: "D:\HUFS\Workspace\Rclass\CC_2023\Seoul.shp", layer: "Seoul"
## with 25 features
## It has 3 fields

Seoul = spTransform(Seoul, CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"))

plot(Seoul)

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

## Warning in SpatialPolygons2map(database, namefield = namefield): database does
## not (uniquely) contain the field 'name'.

R을 활용한 기후변화 자료처리 -심화과정-

2023년 1학기 기후변화론

한국외국어대학교 최용주

2023-05-18

airkorea 자료 처리

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

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

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

대기환경자료 다운로드 / 측정소 위경도 자료

> 엑셀파일 자료 읽기

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

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

파일 읽어오기

파일 형식 살펴보기

만약 한글이 깨진다면?

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

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

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

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

그림 그리기

이제 꾸며보기!!

측정소별 일평균 구하기

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

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

서울지역 측정소별 boxplot

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

서울지역 측정소별 boxplot 꾸미기

서울지역 측정소별 월평균 계산

서울지역 측정소 월평균과 측정소 위경도 병합

서울지역 측정소 그리기

서울 shapefile 가져오기 서울 shp 다운로드