R core와 Rstudio 설치하기

> R cran 웹페이지에서 해당 OS에 해당하는 인스톨 파일 다운로드

Rstudio 다운로드 페이지

Rstudio 다운로드 페이지

> 텍스트 에디터 notepad++ 설치하기 https://notepad-plus-plus.org/downloads/ (맥이면, Sublime Text 추천)

notepad++ 다운로드 페이지

notepad++ 다운로드 페이지

> Rstudio를 실행하면 다음과 같은 화면이 나옴.

Rstudio 실행화면

Rstudio 실행화면

R 시작

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

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

> 기상청에서 운영하는 국내 기후변화감시소 정보 (http://www.climate.go.kr/home/04_watch/01.html)

> WMO에서 운영하는 GAW (global atmospheric watch) 측정자료 다운로드 사이트 (https://ebas-data.nilu.no/)

> 오늘은 안면도 (태안) 기후변화감시소에서 장기간 측정한 CO2와 CH4 자료 분석 (다른 이름으로 저장)

> https://gaw.kishou.go.jp/ 에서 조회 및 다운로드 가능

> CO2 자료 열어보기

co2 자료 구성

co2 자료 구성

> 밑으로 좀 더 내리면…

co2 자료 구성

co2 자료 구성

> R의 기본은 자료 읽기와 쓰기

텍스트 데이터

  • 텍스트 데이터는 데이터와 데이터를 구분하기 위해서 공백(blank), 콤마(comma), 탭(tab)와 같은 구분자(separator)를 사용한다. 텍스트 데이터를 R에서 읽어오기 위해서는 read.table() 함수를 사용하며, 그 사용 방법은 다음과 같다.
argument 설명
file R에서 읽어올 외부 데이터가 있는 파일의 위치(경로)와 파일명(확장자 포함)을 문자형으로 지정한다.
header 논리형으로 TRUE를 지정하면 외부 데이터에 있는 변수명을 R 데이터에서도 동일하게 사용하며, FALSE를 지정하면 외부 데이터의 변수명을 사용하지 않고 R에서 자체적으로 변수명을 지정한다. 그 이름은 “V1”, “V2” 식으로 된다.
sep 구분자로 문자형 형태로 지정한다. 공백이면 ” “, 콤마이면”,“, 탭이면”\t”
stringsAsFactors 문자형 데이터를 요인(factor)으로 자동으로 변경할 지를 설정한다. FALSE를 지정하면 문자형 데이터를 요인으로 변경하지 않고 문자형으로 읽어들인다.
na.strings 데이터에 결측치가 있거나 결측치로 지정할 내용을 문자형으로 지정한다
nrow 읽어들일 행의 갯수를 지정
skip 읽지 않고 넘길 행의 갯수
fileEncoding encoding 타입 윈도우는 WINDOWS-1252, 리눅스/맥은 UTF-8 
co2 = read.table(file = "co2_tap_surface-flask_2_3001-9999_monthly.txt", skip = 209, header = F)
head(co2)
##    V1   V2 V3 V4 V5 V6 V7   V8 V9 V10 V11 V12 V13    V14      V15 V16   V17
## 1 TAP 1990 11  1  0  0  0 -999 -9  -9  -9  -9  -9 360.79 -999.999  -9 36.73
## 2 TAP 1990 12  1  0  0  0 -999 -9  -9  -9  -9  -9 361.32 -999.999  -9 36.73
## 3 TAP 1991  1  1  0  0  0 -999 -9  -9  -9  -9  -9 362.87 -999.999  -9 36.73
## 4 TAP 1991  2  1  0  0  0 -999 -9  -9  -9  -9  -9 364.76 -999.999  -9 36.73
## 5 TAP 1991  3  1  0  0  0 -999 -9  -9  -9  -9  -9 364.85 -999.999  -9 36.73
## 6 TAP 1991  4  1  0  0  0 -999 -9  -9  -9  -9  -9 364.25 -999.999  -9 36.73
##      V18      V19 V20      V21      V22      V23 V24 V25 V26 V27
## 1 126.13 -999.999  20 -999.999 -999.999 -999.999   1  -9  -9 158
## 2 126.13 -999.999  20 -999.999 -999.999 -999.999   1  -9  -9 158
## 3 126.13 -999.999  20 -999.999 -999.999 -999.999   1  -9  -9 158
## 4 126.13 -999.999  20 -999.999 -999.999 -999.999   1  -9  -9 158
## 5 126.13 -999.999  20 -999.999 -999.999 -999.999   1  -9  -9 158
## 6 126.13 -999.999  20 -999.999 -999.999 -999.999   1  -9  -9 158

> 컬럼 이름 변경 colnames()

  • site_gaw_id, year, month, day, hour, minute, second, year, month, day, hour, minute, second, value, value_unc, nvalue, latitude, longitude, altitude, elevation, intake_height, flask_no, ORG_QCflag, QCflag, instrument, measurement_method scale
colnames(co2) =  c("site_gaw_id", "year", "month", "day", "hour", "minute", "second", "year1", "month1", "day1", "hour1", "minute1", "second1", "value", "value_unc", "nvalue", "latitude", "longitude", "altitude", "elevation", "intake_height", "flask_no", "ORG_QCflag", "QCflag", "instrument", "measurement_method", "scale")
head(co2)
##   site_gaw_id year month day hour minute second year1 month1 day1 hour1 minute1
## 1         TAP 1990    11   1    0      0      0  -999     -9   -9    -9      -9
## 2         TAP 1990    12   1    0      0      0  -999     -9   -9    -9      -9
## 3         TAP 1991     1   1    0      0      0  -999     -9   -9    -9      -9
## 4         TAP 1991     2   1    0      0      0  -999     -9   -9    -9      -9
## 5         TAP 1991     3   1    0      0      0  -999     -9   -9    -9      -9
## 6         TAP 1991     4   1    0      0      0  -999     -9   -9    -9      -9
##   second1  value value_unc nvalue latitude longitude altitude elevation
## 1      -9 360.79  -999.999     -9    36.73    126.13 -999.999        20
## 2      -9 361.32  -999.999     -9    36.73    126.13 -999.999        20
## 3      -9 362.87  -999.999     -9    36.73    126.13 -999.999        20
## 4      -9 364.76  -999.999     -9    36.73    126.13 -999.999        20
## 5      -9 364.85  -999.999     -9    36.73    126.13 -999.999        20
## 6      -9 364.25  -999.999     -9    36.73    126.13 -999.999        20
##   intake_height flask_no ORG_QCflag QCflag instrument measurement_method scale
## 1      -999.999 -999.999   -999.999      1         -9                 -9   158
## 2      -999.999 -999.999   -999.999      1         -9                 -9   158
## 3      -999.999 -999.999   -999.999      1         -9                 -9   158
## 4      -999.999 -999.999   -999.999      1         -9                 -9   158
## 5      -999.999 -999.999   -999.999      1         -9                 -9   158
## 6      -999.999 -999.999   -999.999      1         -9                 -9   158

> 필요한 행만 추출하기 (str 함수로 전체적인 자료 형태를 보기)

str(co2)
## 'data.frame':    372 obs. of  27 variables:
##  $ site_gaw_id       : chr  "TAP" "TAP" "TAP" "TAP" ...
##  $ year              : int  1990 1990 1991 1991 1991 1991 1991 1991 1991 1991 ...
##  $ month             : int  11 12 1 2 3 4 5 6 7 8 ...
##  $ day               : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ hour              : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ minute            : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ second            : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ year1             : int  -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 ...
##  $ month1            : int  -9 -9 -9 -9 -9 -9 -9 -9 -9 -9 ...
##  $ day1              : int  -9 -9 -9 -9 -9 -9 -9 -9 -9 -9 ...
##  $ hour1             : int  -9 -9 -9 -9 -9 -9 -9 -9 -9 -9 ...
##  $ minute1           : int  -9 -9 -9 -9 -9 -9 -9 -9 -9 -9 ...
##  $ second1           : int  -9 -9 -9 -9 -9 -9 -9 -9 -9 -9 ...
##  $ value             : num  361 361 363 365 365 ...
##  $ value_unc         : num  -1000 -1000 -1000 -1000 -1000 ...
##  $ nvalue            : int  -9 -9 -9 -9 -9 -9 -9 -9 -9 -9 ...
##  $ latitude          : num  36.7 36.7 36.7 36.7 36.7 ...
##  $ longitude         : num  126 126 126 126 126 ...
##  $ altitude          : num  -1000 -1000 -1000 -1000 -1000 ...
##  $ elevation         : int  20 20 20 20 20 20 20 20 20 20 ...
##  $ intake_height     : num  -1000 -1000 -1000 -1000 -1000 ...
##  $ flask_no          : num  -1000 -1000 -1000 -1000 -1000 ...
##  $ ORG_QCflag        : num  -1000 -1000 -1000 -1000 -1000 ...
##  $ QCflag            : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ instrument        : int  -9 -9 -9 -9 -9 -9 -9 -9 -9 -9 ...
##  $ measurement_method: int  -9 -9 -9 -9 -9 -9 -9 -9 -9 -9 ...
##  $ scale             : int  158 158 158 158 158 158 158 158 158 158 ...
co2 = co2[ , c("site_gaw_id", "year", "month", "day", "hour", "minute", "second","value", "QCflag", "scale")]
head(co2)
##   site_gaw_id year month day hour minute second  value QCflag scale
## 1         TAP 1990    11   1    0      0      0 360.79      1   158
## 2         TAP 1990    12   1    0      0      0 361.32      1   158
## 3         TAP 1991     1   1    0      0      0 362.87      1   158
## 4         TAP 1991     2   1    0      0      0 364.76      1   158
## 5         TAP 1991     3   1    0      0      0 364.85      1   158
## 6         TAP 1991     4   1    0      0      0 364.25      1   158

> CO2 농도 변화 그림 그려보기

plot(co2$value)

plot(co2$value, type = "l")

plot(co2$value, type = "b")

plot(co2$value, type = "l", col = "red")

plot(co2$value, type = "l", col = "red", lty = 2)

plot(co2$value, type = "l", col = "red", lty = 2, lwd = 2)

> 내부 데이터 저장

  • 텍스트 파일은 주로 txt (write.table) 또는 csv로 저장함 (write.csv 함수 사용)
argument 설명
x 저장할 데이터 프레임
file 저장할 파일명
append 기존 파일 뒤에 병합하여 저장하는 옵션, 논리형 (T 또는 F)
sep 각 열의 구분자 (seperator)
row.names 행 이름 (주로 NA로 지정)
col.names 열 이름 (주로 지정하지 않음)
fileEncoding 파일 encoding 형식 (주로 지정하지 않음) 
write.csv(file = "tap_co2.csv", co2, row.names = F) # 행번호를 없애기 위해 F를 씀.

시계열 (time-series) 자료 처리

> 분석에 필요한 라이브러리를 설치하기. 라이브러리의 종류는 수백가지라, 각각의 목적에 맞는 패키지를 검색 후 설치해야함. 시계열 자료 처리에 용이한 패키지는 ‘lubridate’

install.packages("lubridate")
library("lubridate")
## 
## Attaching package: 'lubridate'
## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union

> 지난 시간에 저장한 자료 불러오기

co2 = read.csv(file = "tap_co2.csv", header = T) # 행번호를 없애기 위해 F를 씀.

> ‘연-월-일’ 형식으로 날짜 변수를 만들어주기. (예, 2023-03-15)

co2 파일을 보면 year, month, day를 사용하면 ‘연-월-일’을 만들 수 있음. 각 변수를 문자열로 붙이는 함수는 ’paste’ 또는 ’paste0’임.

함수명 설명 입력내용 결과 내용
paste 지정해 준 구분자 (seperator; sep = “?”)를 이용하여 문자를 합침. 지정안하면 공백 paste(3, 4) “3 4”
paste0 구분자 없이 문자를 합침 paste0(3, 4) “34”

paste 함수에 구분자 ’-’를 사용하여 날짜 형식을 만들어줌. 그리고 ’Date’라는 변수로 저장

co2$Date = paste(co2$year, co2$month, co2$day, sep = "-")

’Date’를 1행부터 10행까지 찍어보면,

co2$Date[1:10]
##  [1] "1990-11-1" "1990-12-1" "1991-1-1"  "1991-2-1"  "1991-3-1"  "1991-4-1" 
##  [7] "1991-5-1"  "1991-6-1"  "1991-7-1"  "1991-8-1"

우리가 원하는 형식과 다른 것을 확인할 수 있음. 따라서 날짜 형식을 ’날짜’로 변경해줘야함. 명령어는 as.POSIXct. ’tz’는 timezone

argument 설명
x 날짜 형식으로 읽을 변수
tz 시간대, 지정안하면 GMT (UTC) 
co2$Date = as.POSIXct(co2$Date, tz = "Asia/Seoul")
co2$Date[1:10]
##  [1] "1990-11-01 KST" "1990-12-01 KST" "1991-01-01 KST" "1991-02-01 KST"
##  [5] "1991-03-01 KST" "1991-04-01 KST" "1991-05-01 KST" "1991-06-01 KST"
##  [9] "1991-07-01 KST" "1991-08-01 KST"

이제 시간에 따른 co2 농도 그림을 그리면,

plot(co2$value ~ co2$Date, ylab = "CO2 (ppm)", xlab = "Year", type = "l", col = "darkorange")
points(co2$value ~ co2$Date, col = "darkorange", pch = 15, cex = 0.5)

성공적으로 날짜 정보를 설정함. 날짜 정보를 이용하여 월별 또는 연도별 평균, 합, 표준편차등을 구할 수 있음. 엑셀에서는 부분합 기능. 함수는 aggregate

argument 설명
x 계산할 항목, 1개 혹은 그 이상도 가능
by 그룹화 할 항목
format 날짜 형식에서 원하는 형식 추출 (%Y, %m, %d, %H, %M, %S, %w
FUN 계산할 함수 (mean, sum, min, max, sd)
na.rm not available 자료 취급 방법 (계산시 제거하려면 T) 
DAT = aggregate(co2$value, by = list(format(co2$Date, "%Y")), FUN = mean, na.rm  = T)
DAT
##    Group.1        x
## 1     1990 361.0550
## 2     1991 360.1542
## 3     1992 360.9367
## 4     1993 360.9033
## 5     1994 361.4808
## 6     1995 364.3700
## 7     1996 366.7417
## 8     1997 367.3242
## 9     1998 370.9333
## 10    1999 372.0475
## 11    2000 373.5692
## 12    2001 375.8825
## 13    2002 377.9625
## 14    2003 379.6200
## 15    2004 381.6567
## 16    2005 383.5620
## 17    2006 387.0350
## 18    2007 389.9058
## 19    2008 392.2200
## 20    2009 390.4392
## 21    2010 396.3575
## 22    2011 398.4542
## 23    2012 401.4908
## 24    2013 405.0858
## 25    2014 405.9192
## 26    2015 409.5167
## 27    2016 412.4875
## 28    2017 414.6800
## 29    2018 414.6875
## 30    2019 418.8958
## 31    2020 421.4425
## 32    2021 424.4608
plot(DAT[,2] ~ DAT[,1], type = "b", lwd = 2, ylab = "CO2 (ppm)", xlab = "Year")

#### > 계절적인 변화를 배제하고 순수한 트렌드만 보고 싶다면, 아래와 같이 자료를 처리함. #### ts() 함수를 사용하여 시계열 자료로 변환해줌. | argument | 설명| |:—| :———————– | |data | 시계열화할 항목| |start | 데이터 수집이 시작된 날짜, c(1990, 11)| |end | 데이터 수집이 종료된 날짜 (필수는 아님)| |frequency | 1년 중 시료를 채취한 빈도, 메일 (365), 매월 (12), 혹은 매분기 (4)|

co2_ts = ts(co2$value, frequency = 12, start = c(1990,11))
plot.ts(co2_ts)

이 시계열의 추세, 계절 및 불규칙 구성 요소를 추정하려면 decompse 함수 사용.

de_co2_ts = decompose(co2_ts)

계절, 추세 및 불규칙 구성 요소의 추정값은 de_co2_ts\(seasonal, de_co2_ts\)trend 및 de_co2_ts$random에 저장됨.

de_co2_ts$seasonal
##             Jan        Feb        Mar        Apr        May        Jun
## 1990                                                                  
## 1991  4.2060451  4.9366146  4.6485451  2.4326701 -0.7606354 -4.4716910
## 1992  4.2060451  4.9366146  4.6485451  2.4326701 -0.7606354 -4.4716910
## 1993  4.2060451  4.9366146  4.6485451  2.4326701 -0.7606354 -4.4716910
## 1994  4.2060451  4.9366146  4.6485451  2.4326701 -0.7606354 -4.4716910
## 1995  4.2060451  4.9366146  4.6485451  2.4326701 -0.7606354 -4.4716910
## 1996  4.2060451  4.9366146  4.6485451  2.4326701 -0.7606354 -4.4716910
## 1997  4.2060451  4.9366146  4.6485451  2.4326701 -0.7606354 -4.4716910
## 1998  4.2060451  4.9366146  4.6485451  2.4326701 -0.7606354 -4.4716910
## 1999  4.2060451  4.9366146  4.6485451  2.4326701 -0.7606354 -4.4716910
## 2000  4.2060451  4.9366146  4.6485451  2.4326701 -0.7606354 -4.4716910
## 2001  4.2060451  4.9366146  4.6485451  2.4326701 -0.7606354 -4.4716910
## 2002  4.2060451  4.9366146  4.6485451  2.4326701 -0.7606354 -4.4716910
## 2003  4.2060451  4.9366146  4.6485451  2.4326701 -0.7606354 -4.4716910
## 2004  4.2060451  4.9366146  4.6485451  2.4326701 -0.7606354 -4.4716910
## 2005  4.2060451  4.9366146  4.6485451  2.4326701 -0.7606354 -4.4716910
## 2006  4.2060451  4.9366146  4.6485451  2.4326701 -0.7606354 -4.4716910
## 2007  4.2060451  4.9366146  4.6485451  2.4326701 -0.7606354 -4.4716910
## 2008  4.2060451  4.9366146  4.6485451  2.4326701 -0.7606354 -4.4716910
## 2009  4.2060451  4.9366146  4.6485451  2.4326701 -0.7606354 -4.4716910
## 2010  4.2060451  4.9366146  4.6485451  2.4326701 -0.7606354 -4.4716910
## 2011  4.2060451  4.9366146  4.6485451  2.4326701 -0.7606354 -4.4716910
## 2012  4.2060451  4.9366146  4.6485451  2.4326701 -0.7606354 -4.4716910
## 2013  4.2060451  4.9366146  4.6485451  2.4326701 -0.7606354 -4.4716910
## 2014  4.2060451  4.9366146  4.6485451  2.4326701 -0.7606354 -4.4716910
## 2015  4.2060451  4.9366146  4.6485451  2.4326701 -0.7606354 -4.4716910
## 2016  4.2060451  4.9366146  4.6485451  2.4326701 -0.7606354 -4.4716910
## 2017  4.2060451  4.9366146  4.6485451  2.4326701 -0.7606354 -4.4716910
## 2018  4.2060451  4.9366146  4.6485451  2.4326701 -0.7606354 -4.4716910
## 2019  4.2060451  4.9366146  4.6485451  2.4326701 -0.7606354 -4.4716910
## 2020  4.2060451  4.9366146  4.6485451  2.4326701 -0.7606354 -4.4716910
## 2021  4.2060451  4.9366146  4.6485451  2.4326701 -0.7606354 -4.4716910
##             Jul        Aug        Sep        Oct        Nov        Dec
## 1990                                              1.4102257  2.9826285
## 1991 -5.8867326 -5.2423021 -3.5610799 -0.6942882  1.4102257  2.9826285
## 1992 -5.8867326 -5.2423021 -3.5610799 -0.6942882  1.4102257  2.9826285
## 1993 -5.8867326 -5.2423021 -3.5610799 -0.6942882  1.4102257  2.9826285
## 1994 -5.8867326 -5.2423021 -3.5610799 -0.6942882  1.4102257  2.9826285
## 1995 -5.8867326 -5.2423021 -3.5610799 -0.6942882  1.4102257  2.9826285
## 1996 -5.8867326 -5.2423021 -3.5610799 -0.6942882  1.4102257  2.9826285
## 1997 -5.8867326 -5.2423021 -3.5610799 -0.6942882  1.4102257  2.9826285
## 1998 -5.8867326 -5.2423021 -3.5610799 -0.6942882  1.4102257  2.9826285
## 1999 -5.8867326 -5.2423021 -3.5610799 -0.6942882  1.4102257  2.9826285
## 2000 -5.8867326 -5.2423021 -3.5610799 -0.6942882  1.4102257  2.9826285
## 2001 -5.8867326 -5.2423021 -3.5610799 -0.6942882  1.4102257  2.9826285
## 2002 -5.8867326 -5.2423021 -3.5610799 -0.6942882  1.4102257  2.9826285
## 2003 -5.8867326 -5.2423021 -3.5610799 -0.6942882  1.4102257  2.9826285
## 2004 -5.8867326 -5.2423021 -3.5610799 -0.6942882  1.4102257  2.9826285
## 2005 -5.8867326 -5.2423021 -3.5610799 -0.6942882  1.4102257  2.9826285
## 2006 -5.8867326 -5.2423021 -3.5610799 -0.6942882  1.4102257  2.9826285
## 2007 -5.8867326 -5.2423021 -3.5610799 -0.6942882  1.4102257  2.9826285
## 2008 -5.8867326 -5.2423021 -3.5610799 -0.6942882  1.4102257  2.9826285
## 2009 -5.8867326 -5.2423021 -3.5610799 -0.6942882  1.4102257  2.9826285
## 2010 -5.8867326 -5.2423021 -3.5610799 -0.6942882  1.4102257  2.9826285
## 2011 -5.8867326 -5.2423021 -3.5610799 -0.6942882  1.4102257  2.9826285
## 2012 -5.8867326 -5.2423021 -3.5610799 -0.6942882  1.4102257  2.9826285
## 2013 -5.8867326 -5.2423021 -3.5610799 -0.6942882  1.4102257  2.9826285
## 2014 -5.8867326 -5.2423021 -3.5610799 -0.6942882  1.4102257  2.9826285
## 2015 -5.8867326 -5.2423021 -3.5610799 -0.6942882  1.4102257  2.9826285
## 2016 -5.8867326 -5.2423021 -3.5610799 -0.6942882  1.4102257  2.9826285
## 2017 -5.8867326 -5.2423021 -3.5610799 -0.6942882  1.4102257  2.9826285
## 2018 -5.8867326 -5.2423021 -3.5610799 -0.6942882  1.4102257  2.9826285
## 2019 -5.8867326 -5.2423021 -3.5610799 -0.6942882  1.4102257  2.9826285
## 2020 -5.8867326 -5.2423021 -3.5610799 -0.6942882  1.4102257  2.9826285
## 2021 -5.8867326 -5.2423021 -3.5610799 -0.6942882
  • 추정된 계절 요인은 1월~12월에 대해 제공되며 매년 동일함. 가장 큰 계절적 요인은 2월(약 4.94)이고 가장 낮은 월은 7월(약 -5.89)로 매년 2월에 CO2가 최고조에 달하고 7월에 최저점으로 나타남.

>“plot()” 함수를 사용하여 시계열의 추정된 추세, 계절 및 불규칙 구성 요소를 그릴 수 있음.

plot(de_co2_ts)

> 그래프 저장는 다양한 형식으로 저장할 수 있음, jpeg, bmp, tiff, pdf등 여러 형태가 가능하나 png()를 주로 사용함.

argument 설명
filename 저장할 파일 이름
height 그림의 폭
units 그림의 높이
frequency 1년 중 시료를 채취한 빈도, 메일 (365), 매월 (12), 혹은 매분기 (4)
res 해상도. 보통 200 이상의 값을 사용. 높을 수록 높은 품질
pointsize 심볼의 크기
bg 배경색 
{
  png(filename = "deseasonal_co2.png", width = 6, height = 8, units = "in", res = 200) # 저장할 스케치북 열기
  plot(de_co2_ts, col = "darkorange")
  dev.off() # 스케치북 닫기
} # (7): (6) 그래프 저장
## png 
##   2

CO2와 기상변수 시간자료 처리해보기

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

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

> 안면도에서 측정한 CO2와 기상변수 다운로드 하기

> CO2 자료 읽기

co2 = read.table(file = "co2_amy_surface-insitu_39_9999-9999_hourly.txt", skip = 214, header = F, na.strings = "-999.999")
head(co2)
##    V1   V2 V3 V4 V5 V6 V7   V8 V9 V10 V11 V12 V13 V14 V15 V16     V17    V18
## 1 AMY 1999  1  1  0  0  0 -999 -9  -9  -9  -9  -9  NA  NA  -9 36.5386 126.33
## 2 AMY 1999  1  1  1  0  0 -999 -9  -9  -9  -9  -9  NA  NA  -9 36.5386 126.33
## 3 AMY 1999  1  1  2  0  0 -999 -9  -9  -9  -9  -9  NA  NA  -9 36.5386 126.33
## 4 AMY 1999  1  1  3  0  0 -999 -9  -9  -9  -9  -9  NA  NA  -9 36.5386 126.33
## 5 AMY 1999  1  1  4  0  0 -999 -9  -9  -9  -9  -9  NA  NA  -9 36.5386 126.33
## 6 AMY 1999  1  1  5  0  0 -999 -9  -9  -9  -9  -9  NA  NA  -9 36.5386 126.33
##   V19 V20 V21 V22 V23 V24 V25 V26 V27
## 1 108  42  66  NA  NA   3   1   9  90
## 2 108  42  66  NA  NA   3   1   9  90
## 3 108  42  66  NA  NA   3   1   9  90
## 4 108  42  66  NA  NA   3   1   9  90
## 5 108  42  66  NA  NA   3   1   9  90
## 6 108  42  66  NA  NA   3   1   9  90

> 컬럼 이름 변경 colnames()

  • site_gaw_id, year, month, day, hour, minute, second, year, month, day, hour, minute, second, value, value_unc, nvalue, latitude, longitude, altitude, elevation, intake_height, flask_no, ORG_QCflag, QCflag, instrument, measurement_method scale
colnames(co2) =  c("site_gaw_id", "year", "month", "day", "hour", "minute", "second", "year1", "month1", "day1", "hour1", "minute1", "second1", "value", "value_unc", "nvalue", "latitude", "longitude", "altitude", "elevation", "intake_height", "flask_no", "ORG_QCflag", "QCflag", "instrument", "measurement_method", "scale")
head(co2)
##   site_gaw_id year month day hour minute second year1 month1 day1 hour1 minute1
## 1         AMY 1999     1   1    0      0      0  -999     -9   -9    -9      -9
## 2         AMY 1999     1   1    1      0      0  -999     -9   -9    -9      -9
## 3         AMY 1999     1   1    2      0      0  -999     -9   -9    -9      -9
## 4         AMY 1999     1   1    3      0      0  -999     -9   -9    -9      -9
## 5         AMY 1999     1   1    4      0      0  -999     -9   -9    -9      -9
## 6         AMY 1999     1   1    5      0      0  -999     -9   -9    -9      -9
##   second1 value value_unc nvalue latitude longitude altitude elevation
## 1      -9    NA        NA     -9  36.5386    126.33      108        42
## 2      -9    NA        NA     -9  36.5386    126.33      108        42
## 3      -9    NA        NA     -9  36.5386    126.33      108        42
## 4      -9    NA        NA     -9  36.5386    126.33      108        42
## 5      -9    NA        NA     -9  36.5386    126.33      108        42
## 6      -9    NA        NA     -9  36.5386    126.33      108        42
##   intake_height flask_no ORG_QCflag QCflag instrument measurement_method scale
## 1            66       NA         NA      3          1                  9    90
## 2            66       NA         NA      3          1                  9    90
## 3            66       NA         NA      3          1                  9    90
## 4            66       NA         NA      3          1                  9    90
## 5            66       NA         NA      3          1                  9    90
## 6            66       NA         NA      3          1                  9    90

> 필요한 행만 추출하기 (str 함수로 전체적인 자료 형태를 보기)

str(co2)
## 'data.frame':    204984 obs. of  27 variables:
##  $ site_gaw_id       : chr  "AMY" "AMY" "AMY" "AMY" ...
##  $ year              : int  1999 1999 1999 1999 1999 1999 1999 1999 1999 1999 ...
##  $ month             : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ day               : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ hour              : int  0 1 2 3 4 5 6 7 8 9 ...
##  $ minute            : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ second            : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ year1             : int  -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 ...
##  $ month1            : int  -9 -9 -9 -9 -9 -9 -9 -9 -9 -9 ...
##  $ day1              : int  -9 -9 -9 -9 -9 -9 -9 -9 -9 -9 ...
##  $ hour1             : int  -9 -9 -9 -9 -9 -9 -9 -9 -9 -9 ...
##  $ minute1           : int  -9 -9 -9 -9 -9 -9 -9 -9 -9 -9 ...
##  $ second1           : int  -9 -9 -9 -9 -9 -9 -9 -9 -9 -9 ...
##  $ value             : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ value_unc         : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ nvalue            : int  -9 -9 -9 -9 -9 -9 -9 -9 -9 -9 ...
##  $ latitude          : num  36.5 36.5 36.5 36.5 36.5 ...
##  $ longitude         : num  126 126 126 126 126 ...
##  $ altitude          : int  108 108 108 108 108 108 108 108 108 108 ...
##  $ elevation         : int  42 42 42 42 42 42 42 42 42 42 ...
##  $ intake_height     : int  66 66 66 66 66 66 66 66 66 66 ...
##  $ flask_no          : logi  NA NA NA NA NA NA ...
##  $ ORG_QCflag        : logi  NA NA NA NA NA NA ...
##  $ QCflag            : int  3 3 3 3 3 3 3 3 3 3 ...
##  $ instrument        : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ measurement_method: int  9 9 9 9 9 9 9 9 9 9 ...
##  $ scale             : int  90 90 90 90 90 90 90 90 90 90 ...
co2 = co2[ , c("site_gaw_id", "year", "month", "day", "hour", "minute", "second","value", "QCflag", "scale")]
head(co2)
##   site_gaw_id year month day hour minute second value QCflag scale
## 1         AMY 1999     1   1    0      0      0    NA      3    90
## 2         AMY 1999     1   1    1      0      0    NA      3    90
## 3         AMY 1999     1   1    2      0      0    NA      3    90
## 4         AMY 1999     1   1    3      0      0    NA      3    90
## 5         AMY 1999     1   1    4      0      0    NA      3    90
## 6         AMY 1999     1   1    5      0      0    NA      3    90

paste 함수에 구분자 ’-’를 사용하여 날짜 형식을 만들어줌. 그리고 ’Time’라는 변수로 저장

co2$Time = paste(paste(co2$year, co2$month, co2$day, sep = "-"), paste(co2$hour, co2$minute, co2$second, sep = ":"))

’Time’를 1행부터 10행까지 찍어보면,

co2$Time[1:10]
##  [1] "1999-1-1 0:0:0" "1999-1-1 1:0:0" "1999-1-1 2:0:0" "1999-1-1 3:0:0"
##  [5] "1999-1-1 4:0:0" "1999-1-1 5:0:0" "1999-1-1 6:0:0" "1999-1-1 7:0:0"
##  [9] "1999-1-1 8:0:0" "1999-1-1 9:0:0"

우리가 원하는 형식과 다른 것을 확인할 수 있음. 따라서 날짜 형식을 ’날짜’로 변경해줘야함. 명령어는 as.POSIXct. ’tz’는 timezone

co2$Time = as.POSIXct(co2$Time, tz = "Asia/Seoul")
co2$Time[1:10]
##  [1] "1999-01-01 00:00:00 KST" "1999-01-01 01:00:00 KST"
##  [3] "1999-01-01 02:00:00 KST" "1999-01-01 03:00:00 KST"
##  [5] "1999-01-01 04:00:00 KST" "1999-01-01 05:00:00 KST"
##  [7] "1999-01-01 06:00:00 KST" "1999-01-01 07:00:00 KST"
##  [9] "1999-01-01 08:00:00 KST" "1999-01-01 09:00:00 KST"

> 기상자료도 똑같이 처리

met = read.table(file = "co2_amy_surface-insitu_39_9999-9999_hourly_met.txt", skip = 140, header = F, na.strings = c("-9999.9", "-99.9", "-999.9", "-999.999"))

colnames(met) =  c("site_gaw_id", "year", "month", "day", "hour", "minute", "second", "wind_direction", "wind_speed", "relative_humidity", "precipitation_amount", "air_pressure", "air_temperature", "dew_point_temperature", "sea_water_temperature", "sea_surface_water_temperature", "sea_water_salinity", "sea_surface_water_salinity", "latitude", "longitude", "altitude", "elevation")
str(met)
## 'data.frame':    181810 obs. of  22 variables:
##  $ site_gaw_id                  : chr  "AMY" "AMY" "AMY" "AMY" ...
##  $ year                         : int  1999 1999 1999 1999 1999 1999 1999 1999 1999 1999 ...
##  $ month                        : int  3 3 3 3 3 3 3 3 3 3 ...
##  $ day                          : int  29 29 29 29 29 29 29 29 30 30 ...
##  $ hour                         : int  16 17 18 19 20 21 22 23 0 1 ...
##  $ minute                       : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ second                       : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ wind_direction               : num  NA NA 25.8 25.7 41.7 34.8 36.5 37 47.1 70.3 ...
##  $ wind_speed                   : num  NA NA 1.6 0.5 0.3 1.9 1.1 1.4 2 1.1 ...
##  $ relative_humidity            : logi  NA NA NA NA NA NA ...
##  $ precipitation_amount         : logi  NA NA NA NA NA NA ...
##  $ air_pressure                 : logi  NA NA NA NA NA NA ...
##  $ air_temperature              : logi  NA NA NA NA NA NA ...
##  $ dew_point_temperature        : logi  NA NA NA NA NA NA ...
##  $ sea_water_temperature        : logi  NA NA NA NA NA NA ...
##  $ sea_surface_water_temperature: logi  NA NA NA NA NA NA ...
##  $ sea_water_salinity           : logi  NA NA NA NA NA NA ...
##  $ sea_surface_water_salinity   : logi  NA NA NA NA NA NA ...
##  $ latitude                     : num  36.5 36.5 36.5 36.5 36.5 ...
##  $ longitude                    : num  126 126 126 126 126 ...
##  $ altitude                     : logi  NA NA NA NA NA NA ...
##  $ elevation                    : int  42 42 42 42 42 42 42 42 42 42 ...
met = met[ , c("site_gaw_id", "year", "month", "day", "hour", "minute", "second","wind_direction", "wind_speed")]
head(co2)
##   site_gaw_id year month day hour minute second value QCflag scale
## 1         AMY 1999     1   1    0      0      0    NA      3    90
## 2         AMY 1999     1   1    1      0      0    NA      3    90
## 3         AMY 1999     1   1    2      0      0    NA      3    90
## 4         AMY 1999     1   1    3      0      0    NA      3    90
## 5         AMY 1999     1   1    4      0      0    NA      3    90
## 6         AMY 1999     1   1    5      0      0    NA      3    90
##                  Time
## 1 1999-01-01 00:00:00
## 2 1999-01-01 01:00:00
## 3 1999-01-01 02:00:00
## 4 1999-01-01 03:00:00
## 5 1999-01-01 04:00:00
## 6 1999-01-01 05:00:00
met$Time = paste(paste(met$year, met$month, met$day, sep = "-"), paste(met$hour, met$minute, met$second, sep = ":"))
met$Time = as.POSIXct(met$Time, tz = "Asia/Seoul")
met$Time[1:10]
##  [1] "1999-03-29 16:00:00 KST" "1999-03-29 17:00:00 KST"
##  [3] "1999-03-29 18:00:00 KST" "1999-03-29 19:00:00 KST"
##  [5] "1999-03-29 20:00:00 KST" "1999-03-29 21:00:00 KST"
##  [7] "1999-03-29 22:00:00 KST" "1999-03-29 23:00:00 KST"
##  [9] "1999-03-30 00:00:00 KST" "1999-03-30 01:00:00 KST"

> 두 데이터 셋의 병합. (시간을 기준으로): co2와 met는 행의 갯수가 서로 다름

nrow(co2); nrow(met)
## [1] 204984
## [1] 181810

merge() 함수를 이용해서 두 데이터 셋을 합치기.

argument 설명
x, y 합칠 데이터 프레임 두 개
by 합칠 기준이되는 컬럼 이름. 지금이라면 “Time
all 모든 자료 다 합침
all.x x 데이터 기준 모든 자료 다 합침
all.y y 데이터 기준 모든 자료 다 합침 
co2.met = merge(co2, met, by = "Time")
head(co2.met)
##                  Time site_gaw_id.x year.x month.x day.x hour.x minute.x
## 1 1999-03-29 16:00:00           AMY   1999       3    29     16        0
## 2 1999-03-29 17:00:00           AMY   1999       3    29     17        0
## 3 1999-03-29 18:00:00           AMY   1999       3    29     18        0
## 4 1999-03-29 19:00:00           AMY   1999       3    29     19        0
## 5 1999-03-29 20:00:00           AMY   1999       3    29     20        0
## 6 1999-03-29 21:00:00           AMY   1999       3    29     21        0
##   second.x  value QCflag scale site_gaw_id.y year.y month.y day.y hour.y
## 1        0 382.01      2    90           AMY   1999       3    29     16
## 2        0 383.47      2    90           AMY   1999       3    29     17
## 3        0 386.09      2    90           AMY   1999       3    29     18
## 4        0     NA      3    90           AMY   1999       3    29     19
## 5        0 384.64      2    90           AMY   1999       3    29     20
## 6        0 385.10      2    90           AMY   1999       3    29     21
##   minute.y second.y wind_direction wind_speed
## 1        0        0             NA         NA
## 2        0        0             NA         NA
## 3        0        0           25.8        1.6
## 4        0        0           25.7        0.5
## 5        0        0           41.7        0.3
## 6        0        0           34.8        1.9

과제: co2.met 자료로 부터 매 시간평균 CO2를 구하고 (0시부터 23시), 그래프로 그려오시오.

필요한 열만 선택해서 붙여보기!

co2.met = merge(co2[,c("Time", "value")], met[,c("Time", "wind_direction", "wind_speed")], by = "Time")

co2.met를 csv 파일로 저장하기!

write.csv(file = "AMY_co2_met.csv", co2.met, row.names = F)

Openair

openair는 대기 자료 분석에 가장 많이 사용되는 패키지로써, 대기환경에서 가장 많이 사용되는 패키지중 하나이다. 특히 복잡한 통계 및 그래프를 짧은 명령어로 충분히 실행되기에 초보자에게도 진입장벽이 높지 않다.

openair 패키지 설치: openair가 처음이라면 아래 명령어를 사용하여 설치할 수 있다.

install.packages("openair")

패키지 불러오기: 기존에 설치가 되었거나, 설치가 완료되었다면, 아래 명령어로 lubridate 패키지와 함께 불러온다.

library("lubridate")
library("openair")

자료 불러오기: csv 파일이기에 read.csv 명령어를 이용하여 앞서 저장한 자료 (AMY_co2_met.csv)를 불러들인다. 우선 working directory 설정부터…

setwd("D:/HUFS/Workspace/Rclass/CC_2023/")
co2.met = read.csv(file = "AMY_co2_met.csv", header = T)

시간형태로 읽어들이기! as.POSIXct 명령어를 사용 (!! 주의할 점, openair는 date를 시간으로 인식함!)

co2.met$date = as.POSIXct(co2.met$Time, tz = "Asia/Seoul")  #연, 월, 일을 인식시켜줌으로써 손쉽게 시간형태로 변경

summary plot

  • 시계열 자료의 전체적인 경향을 파악하기 좋은 그래프. 자료 및 결측치의 분포, 평균, 중간값 등에 대한 정보를 제공함함
summaryPlot(co2.met) # 안면의 summary plot

summaryPlot(co2.met, percentile = 0.95) # 상위 5% 값은 제거
summaryPlot(co2.met, na.len = 10) # 적어도 10개이상 연속적인 결측치를 보여줌
summaryPlot(co2.met, col.data = "green") # 자료 색깔을 녹색으로
summaryPlot(co2.met, col.mis = "yellow") # 결측치를 노란색으로
summaryPlot(co2.met, col.dens = "black") # density plot선을 검은색으로
summaryPlot(co2.met[, c(5,2,4)]) # 컬럼 2번째와 5-7번째만 보여주기
summaryPlot(subset(co2.met, select = c(date, value, wind_speed))) # 앞선 명령어의 대안

windrose

  • 풍향, 풍속에 대한 정보를 제공 (!! 주의할 점, openair는 풍향, 풍속을 wd, ws로 인식함)
colnames(co2.met)[c(2, 3,4)] = c("co2","wd", "ws") # 풍향, 풍속의 컬럼 이름을 wd, ws로 교체. (교체 안할 시 일일이 지정해야하는 번거로움..)

windRose(co2.met) # 기본 바랑장미

windRose(co2.met, type = "month", layout = c(4, 3)) # 월별 바랑장미

windRose(co2.met[is.na(co2.met$co2) == F,], type = "co2", layout = c(4, 1)) # PM2.5 구간에 따른 바람장비

### pollution rose

  • 오염장미
pollutionRose(co2.met, pollutant = "co2") # PM2.5 농도 구간별 풍향의 빈도 수

pollutionRose(co2.met, pollutant = "co2", normalise = TRUE) # normalized된 PM2.5 농도 구간별 풍향의 빈도 수

Polar plot

  • 풍향/풍속에 따른 오염물질의 평균 농도를 보여줌
polarPlot(co2.met, pollutant = "co2") 

polarPlot(co2.met, pollutant = "co2", statistic = "cpf", percentile = 75) # Conditional Probability Function

polarPlot(co2.met, pollutant = "co2", statistic = "cpf", percentile = c(0,10)) # Conditional Probability Function

polarPlot(co2.met, pollutant = "co2", statistic = "cpf", percentile = c(50,60)) # Conditional Probability Function

polarPlot(co2.met, pollutant = "co2", type = "daylight") # ”year”, ”hour”, ”month”, ”season”, ”weekday”, ”weekend”, ”monthyear”, ”daylight”

Calendar plot

  • 달력에 원하는 오염물질의 농도등을 표시하는 그림. 보다 직관적이기에 발표자료에 많이 사용됨.
st = as.POSIXct("2018-12-01") # 분석 시작 날짜
et = as.POSIXct("2019-03-31") # 분석 종료 날짜

calendarPlot(co2.met[co2.met$date >= st & co2.met$date <= et,], pollutant = "co2") # 용산 PM2.5 농도 기준, 표시

calendarPlot(co2.met[co2.met$date >= st + years(1) & co2.met$date <= et + years(1),], pollutant = "co2") # 분석 기간 기준을 1년 후로 미룸. 용산 PM2.5 농도 기준, 표시

calendarPlot(co2.met[co2.met$date >= st & co2.met$date <= et,], pollutant = "co2", annotate = "ws") # 관악 PM2.5 기준 농도가 아닌, 풍향/풍속 표시

Trend 분석 (smoothTrend, TheilSen)

smoothTrend(co2.met, pollutant = "co2", deseason = TRUE, simulate = F,  ylab = "concentration (ppm)",
            main = "monthly mean deseasonalised CO2 (bootstrap uncertainties)") # smooth trend를 계산. GAM 모델을 사용하여 smooth를 수행
## Warning in StructTS(myts): possible convergence problem: 'optim' gave code = 52
## and message 'ERROR: ABNORMAL_TERMINATION_IN_LNSRCH'

smoothTrend(co2.met, pollutant = "co2", deseason = TRUE, simulate = F,  ylab = "concentration (ppm)",
            main = "monthly mean deseasonalised CO2 (bootstrap uncertainties)", type = "wd") # 풍향을 기준으로 분류
## Warning: 5744 missing wind direction line(s) removed

## Warning: possible convergence problem: 'optim' gave code = 52 and message
## 'ERROR: ABNORMAL_TERMINATION_IN_LNSRCH'

## Warning: possible convergence problem: 'optim' gave code = 52 and message
## 'ERROR: ABNORMAL_TERMINATION_IN_LNSRCH'

TheilSen(co2.met, pollutant = "co2", deseason = TRUE,  ylab = "CO2", main = "Theil-Sen Trend") # Theil-sen 기울기 계산산
## Warning in StructTS(myts): possible convergence problem: 'optim' gave code = 52
## and message 'ERROR: ABNORMAL_TERMINATION_IN_LNSRCH'
## [1] "Taking bootstrap samples. Please wait."

TheilSen(co2.met, pollutant = "co2", deseason = TRUE,  ylab = "CO2", main = "Theil-Sen Trend", type = "season") # 계절절을 기준으로 분류
## Warning in StructTS(myts): possible convergence problem: 'optim' gave code = 52
## and message 'ERROR: ABNORMAL_TERMINATION_IN_LNSRCH'
## [1] "Taking bootstrap samples. Please wait."
## [1] "Taking bootstrap samples. Please wait."
## Warning in StructTS(myts): possible convergence problem: 'optim' gave code = 52
## and message 'ERROR: ABNORMAL_TERMINATION_IN_LNSRCH'
## [1] "Taking bootstrap samples. Please wait."
## [1] "Taking bootstrap samples. Please wait."

### timeVariation 분석

  • 월별, 시간별, 요일별 분석을 수행 가능
# 기본 그림
timeVariation(co2.met, pollutant = "co2")

# 여러 항목들의 시간변화 비교 (절대값이 다르기에 nomarlized 적용)
timeVariation(co2.met, pollutant = c("co2", "wd", "ws"), normalise = TRUE)

# factor (예: 기준)을 적용하여 전,후 차이를 비교. 아래는 코로나 전, 후의 NO2의 변화
co2.met[co2.met$date < as.POSIXct("2010-01-01"), "millenium"] = "Before 2010" # 코로나 전 지정
co2.met[co2.met$date >= as.POSIXct("2010-01-01"), "millenium"] = "After 2010" # 코로나 후 지정
timeVariation(co2.met, pollutant = "co2", group = "millenium", difference = TRUE)

timeVariation(co2.met, pollutant = "co2", group = "millenium", difference = TRUE, ylim = c(350,450))

그 외 유용한 함수

co2.met = rollingMean(co2.met, pollutant = "co2", hours = 8, new.name = "co2.8", data.thresh = 75) # 8시간 이동 평균
trendLevel(co2.met, pollutant = "co2", cols = "jet", y = "hour")

airkorea 자료 처리

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

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

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

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

> 엑셀파일 자료 읽기

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

install.packages("readxl")

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

library(readxl)
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()
## [1] "LC_COLLATE=English_United States.utf8;LC_CTYPE=English_United States.utf8;LC_MONETARY=English_United States.utf8;LC_NUMERIC=C;LC_TIME=English_United States.utf8"
Sys.setlocale("LC_ALL","C") # 강제 언어 삭제
## [1] "C"
unique(nier_site$type)
##  [1] "<U+B3C4><U+C2DC><U+B300><U+AE30><U+CE21><U+C815><U+B9DD>" "<U+B3C4><U+B85C><U+BCC0><U+B300><U+AE30><U+CE21><U+C815><U+B9DD>"
##  [3] "<U+C720><U+D574><U+B300><U+AE30><U+BB3C><U+C9C8><U+CE21><U+C815><U+B9DD>" "PM2.5<U+C131><U+BD84><U+CE21><U+C815><U+B9DD>"
##  [5] "<U+B300><U+AE30><U+C911><U+AE08><U+C18D><U+CE21><U+C815><U+B9DD>" "<U+C0B0><U+C131><U+AC15><U+D558><U+BB3C><U+CE21><U+C815><U+B9DD>"
##  [7] "<U+C885><U+D569><U+B300><U+AE30><U+CE21><U+C815><U+C18C>" "<U+AD11><U+D654><U+D559><U+B300><U+AE30><U+C624><U+C5FC><U+BB3C><U+C9C8><U+CE21><U+C815><U+B9DD>"
##  [9] "<U+AD50><U+C678><U+B300><U+AE30><U+CE21><U+C815><U+B9DD>" "<U+D56D><U+B9CC><U+CE21><U+C815><U+B9DD>"
## [11] "<U+AD6D><U+AC00><U+BC30><U+ACBD><U+B18D><U+B3C4><U+CE21><U+C815><U+B9DD>" "<U+B300><U+AE30><U+D658><U+ACBD><U+C5F0><U+AD6C><U+C18C>"

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

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) # 지도 추가
}