R에 들어가기에 앞서..

R 소프트웨어는 통계 처리를 목적으로 만들어진 오픈 소스 프로그램이다. 기존의 상용화된 통계 프로그램인 ‘S’ 보다 가격면 뿐만 아닌 성능면에서 능가하는 프로그램이 되길 기원하는 마음으로 알파벳 순서가 한글자 앞인 ’R’이라고 명명했다. R을 사용하기에 앞서 R core와 R studio가 설치되어야한다.

R studio는 보다 편하게 R를 쓸 수 있게 해주는 GUI 프로그램이다.

문법

기본 문법

RStudio 프로그램을 이용하여 R의 기본적인 문법과 데이터의 유형을 학습한다. R에서는 콘솔(console)에도 명령어를 입력하여 실행할 수 있지만 가능하면 스크립트 (script)를 이용하여 명령어를 실행하고 저장하는 것이 좋다. 여기서는 RStudio에서 제공하는 스크립트를 기준으로 설명한다.

  • #

    주석(comment)의 기능으로 프로그램 전반적인 내용, 명령어의 내용 등이 무엇을 의미하는지 알 수 있도록 사용자가 설명을 달아주는 기능이다. # 뒤에 있는 한 줄이 주석으로 처리되며 R에서 지정된 문법을 검사하지 않는다. 단 한 줄만 주석으로 처리되므로 다른 줄도 하고 싶으면 해당 줄 앞에 #를 써 주어야 한다.

  • ;

    세미콜론은 하나의 명령어가 끝났음을 알려주는 기능이다. 하지만 한 줄에 하나의 명령밖에 없으면 세미콜론을 해 주지 않아도 명령어가 끝났음을 인식한다.

  • Enter

    다음 줄로 이동할 때 사용한다.

  • Ctrl + Enter

    R의 명령어를 실행하는 기능이다. 명령어가 한 줄인 경우는 마우스의 위치는 해당 줄의 아무 곳에 있어도 상관이 없다. 명령어가 두 줄 이상인 경우에는 반드시 해당 명령어가 있는 곳을 블록잡고 실행시켜야 한다.

  • Shift + Enter

    함수를 사용할 때에 함수에 들어가는 값을 argument라고 한다. 함수 안에 들어가는 argument가 많아지면 하나의 명령어가 길게 표현된다. 그러면 명령어를 이해하는 데에 불편함이 있다. Shift + Enter를 하면 동일한 위치에 다른 argument를 올 수 있도록 해 준다.

  • 대소문자

    R은 대소문자를 구별한다. 이것은 “case sensitive하다” 라고 표현한다. 소문자 x와 대문자 X는 전혀 다른 것을 의미하기 때문에 주의하기 바란다.

  • 산술연산자 (+, -, *, /, …)

기호 의미
+ 더하기
- 빼기
* 곱하기
/ 나누기
** 거듭제곱
^ 거듭제곱
%/% 몫
%% 나머지
  • 비교 연산자(Relational Operator)
연산자 설명 입력내용 결과 내용
> 크다 3 > 4 FALSE
>= 크거나 같다 3 >= 4 TRUE
< 작다 3 < 4 TRUE
<= 작거나 같다 3 <= 4 TRUE
== 같다 3 == 4 FALSE
!= 같지 않다 3 != 4 TRUE
! 아니다 !(3 == 4) TRUE
  • 할당 연산자(allocation operator)
연산자 설명 입력내용
<- 오른쪽의 값을 왼쪽의 이름에 저장함 x <- 3
= ’’ y = 4
-> 왼쪽의 값을 오른쪽의 이름에 저장함 5 -> z
  • 논리 연산자(logical operator)
연산자 설명 입력내용 결과 내용
& 또는 && AND (조건1) & (조건2)
| 또는 || OR (조건1) | (조건2)
! NOT !(조건)
  • 수학함수 (sin, cos, tan, …), 지수/로그(^, log, log10) 등등 숫자형 변수간의 수학연산이 가능하다
함수명 설명 입력내용 결과 내용
abs() 절대값 abs(-3) 3
sqrt() 제곱근 sqrt(16) 4
pi 원주율 pi 3.141593
sign() 부호 sign(-3) -1
round() 반올림 round(2.345, digits=2) 2.35
ceiling() 무조건 올림 ceiling(2.3) 3
floor() 무조건 내림 floor(2.7) 2
exp() 지수 exp(10) 22026.47
log() 자연로그 log(10) 2.302585
log10() 사용로그 log10(10) 1
log2() 로그(2) log2(10) 3.321928
logb() 일반화 로그 logb(10, base=3) 2.095903
factorial() 계승 factorial(4) 24
choose() 조합 choose(4, 2) 6
prod() 곱 prod(1:4) 24
sin() 사인(sine) sin(0.5) 0.4794255
cos() 코사인(cosine) cos(0.5) 0.8775826
tan() 탄젠트(tangent) tan(0.5) 0.5463025

데이터 유형

  • 수치형(numeric): 숫자로 되어 있으며, 정수형(integer)과 실수형(double)이 있다.
  • 문자형(character): 하나의 문자 또는 문자열로 되어 있으며, “” 또는 ’’로 묶여 있다.
  • 논리형(logical): 참과 거짓의 논리값으로 TRUE(or T)이나 FALSE(or F)를 가진다.
  • 복소수형(complex): 실수와 허수로 이루어진 복소수이다
  • NULL: 객체(object)로서 존재하지 않는 객체로 지정할 때 사용함.
  • NA: Not Available의 약자로 결측치(missing value)를 의미함
  • NaN: Not available Number의 약자로 수학적으로 계산이 불가능한 수를 의미한다. 예를 들면 sqrt(-3)로 음수에 대한 제곱근은 구할 수 없다.
  • Inf: Infinite의 양자로 양의 무한대이다.
  • -Inf: 음의 무한대이다.
  • 데이터의 유형을 알려주는 함수: mode() 혹은 is.~ 함수
x1 = 3; x2 = "NIER"; x3 = FALSE; x4 = 3-2i
mode(x1); mode(x2); mode(x3); mode(x4)
## [1] "numeric"
## [1] "character"
## [1] "logical"
## [1] "complex"
함수명 설명 입력내용 결과 내용
is.numeric() 수치형 여부 is.numeric(데이터) TRUE or FALSE
is.integer() 정수형 여부 is.integer(데이터) TRUE or FALSE
is.double() 실수형 여부 is.double(데이터) TRUE or FALSE
is.character() 문자형 여부 is.character(데이터) TRUE or FALSE
is.logical() 논리형 여부 is.logical(데이터) TRUE or FALSE
is.complex() 복소수형 여부 is.complex(데이터) TRUE or FALSE
is.null() NULL 여부 is.null(데이터) TRUE or FALSE
is.na() NA 여부 is.na(데이터) TRUE or FALSE
is.finite() 유한 수치 여부 is.finite(데이터) TRUE or FALSE
is.infinite() 무한 수치 여부 is.infinite(데이터) TRUE or FALSE
  • 데이터의 유형의 우선순위: 문자형(character) > 복소수형(complex) > 수치형(numeric) > 논리형(logical)
x2 = c(x1, x2, x3, x4)
x2
## [1] "3"     "NIER"  "FALSE" "3-2i"
  • 데이터의 유형을 강제적으로 변경하는 함수
함수명 설명 입력내용 결과 내용
as.numeric() 수치형으로 변환 as.numeric(데이터) 변환되거나 NA
as.integer() 정수형으로 변환 as.integer(데이터) 변환되거나 NA
as.double() 실수형으로 변환 as.double(데이터) 변환되거나 NA
as.character() 문자형으로 변환 as.character(데이터) 변환되거나 NA
as.logical() 논리형으로 변환 as.logical(데이터) 변환되거나 NA
as.complex() 복소수형으로 변환 as.complex(데이터) 변환되거나 NA

데이터 형태

벡터(vector): 스칼라의 확장으로 한 개 이상의 데이터를 갖는다

벡터의 생성
  • c() : combine 또는 concatenate. 기본적으로는 스칼라가 들어가며, 스칼라와 스칼라 사이는 콤마로 구분한다. 또한 벡터도 들어갈 수 있다
v1 = c(3, 10, 12)
v2 = c("Kim", "Lee", "Park")
v3 = c(TRUE, FALSE, FALSE, FALSE)
v4 = c(v1, v2, v3)
v1; v2; v3; v4
## [1]  3 10 12
## [1] "Kim"  "Lee"  "Park"
## [1]  TRUE FALSE FALSE FALSE
##  [1] "3"     "10"    "12"    "Kim"   "Lee"   "Park"  "TRUE"  "FALSE" "FALSE"
## [10] "FALSE"
  • : : 수치형에만 적용되며, 1씩 증가되거나 1씩 감소되는 규칙이 있는 벡터를 생성. start:end 구조로 사용하며, start, end는 숫자이고, start > end이면 1씩 감소되고, start < end이면 1씩 증가되며, start=end 이면 start 또는 end가 된다. 시작은 무조건 start에서 시작해서 end를 넘지 않는다
(v1 = 1:5)
## [1] 1 2 3 4 5
(v2 = 5:1)
## [1] 5 4 3 2 1
(v3 = -3.3:5)
## [1] -3.3 -2.3 -1.3 -0.3  0.7  1.7  2.7  3.7  4.7
(v4 = 5:-3.3)
## [1]  5  4  3  2  1  0 -1 -2 -3
  • seq() : sequence의 약자, 1 이외의 증감이 되는 규칙 있는 수치형 벡터를 생성
argument 설명
from 시작값
to 끝값
by 얼마씩 증가 또는 감소시킬지를 정하는 단계값. 단, 감소시킬 경우에는 부호가 음수이어야 한다. 그렇지 않으면 에러(error)가 발생. 
(v1 = seq(from = 1, to = 5, by = 1))
## [1] 1 2 3 4 5
(v2 = seq(from = 1, to = 5, by = 0.5))
## [1] 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0
(v3 = seq(from = 5, to = 1, by = -0.5))
## [1] 5.0 4.5 4.0 3.5 3.0 2.5 2.0 1.5 1.0
(v4 = seq(5, 1, -0.5))
## [1] 5.0 4.5 4.0 3.5 3.0 2.5 2.0 1.5 1.0
벡터의 속성
  • 데이터의 유형: mode(), is.numeric(), is.character(), is.logical(),
mode(v1)
## [1] "numeric"
is.character(v1)
## [1] FALSE
is.numeric(v1)
## [1] TRUE
  • 데이터의 이름: names(벡터)
names(v1) = c("A", "B", "C", "D", "E")
v1
## A B C D E 
## 1 2 3 4 5
  • 데이터의 추출: 대괄호([])를 사용. 양의 정수는 색인(index)으로 벡터가 가지는 데이터의 위치를 의미
weight = c(57, 81, 65, 49, 72)
weight[1]
## [1] 57
weight[2]
## [1] 81
weight[2:4]
## [1] 81 65 49
weight[c(1, 4, 5)]
## [1] 57 49 72
weight[-c(1, 4, 5)]
## [1] 81 65
벡터의 연산
  • 벡터의 길이가 동일한 경우: 연산에 사용되는 벡터들의 길이가 동일한 경우로서 벡터들 간의 사칙연산을 할 수 있고, 최종적인 결과는 벡터가 된다. 벡터들 간의 연산이 될 때에는 각 벡터에 있는 동일한 위치의 값들 간이 연산이 된다. 예를 들어, 수치형 벡터인 v1과 v2가 있고 v1 + v2의 연산을 하면, v1[1] + v2[1]의 연산이 된다.
v1 = 1:3
v2 = 4:6
v1 + v2
## [1] 5 7 9
v1 - v2
## [1] -3 -3 -3
v1 * v2
## [1]  4 10 18
v1 / v2
## [1] 0.25 0.40 0.50
v1 ** v2
## [1]   1  32 729
  • 벡터의 길이가 동일하지 않은 경우: 연산과정에서 데이터의 개수가 적은 쪽의 벡터는 데이터 개수가 많은 쪽의 벡터와 동일하게 데이터의 개수를 맞춘다.
v1 = 1:3
v2 = 1:6
v1 + v2 
## [1] 2 4 6 5 7 9
v2 = 1:10
v1 + v2
## Warning in v1 + v2: 두 객체의 길이가 서로 배수관계에 있지 않습니다
##  [1]  2  4  6  5  7  9  8 10 12 11

요인(factor)

  • 데이터를 질적 자료(또는 범주형 자료)로 변환. 질적자료로 변경되면 집단별로 통계분석을 할 수 있다. factor() 함수를 사용하며, 그 사용 방법은 다음과 같다
argument 설명
x 벡터를 지정한다.
levels 그룹으로 지정할 문자형 벡터를 지정하며, levels를 쓰지 않으면 오름차순으로 구분하여 자체적으로 그룹을 지정한다.
labels levels에 대한 문자형 벡터를 지정한다.
ordered levels에 대해 특정한 순서를 정하고 싶으면 TRUE를 지정한다. 
bt = c("A", "O", "A", "AB", "O", "B")
bt
## [1] "A"  "O"  "A"  "AB" "O"  "B"
bt_factor = factor(bt)
bt_factor
## [1] A  O  A  AB O  B 
## Levels: A AB B O
levels(bt_factor)
## [1] "A"  "AB" "B"  "O"
levels(bt_factor) = c("A형", "AB형", "B형", "O형")
bt_factor
## [1] A형  O형  A형  AB형 O형  B형 
## Levels: A형 AB형 B형 O형

리스트

  • 주로 통계분석의 결과를 저장할 때 사용하는 데이터 형태이며, 리스트의 원소로 스칼라, 벡터, 행렬, 배열, 데이터 프레임, 요인 그리고 리스트 자신을 가질 수 있다.
s1 = 10
v1 = 1:10
m1 = matrix(1:4, nrow = 2, ncol = 2)
a1 = array(1:8, dim = c(2, 2, 2))
df1 = data.frame(id = 1:3, age = c(20, 30, 40))
f1 = factor(c("A", "O", "A", "AB", "O", "B"))
result = list(s1, v1, m1, a1, df1, f1)
result
## [[1]]
## [1] 10
## 
## [[2]]
##  [1]  1  2  3  4  5  6  7  8  9 10
## 
## [[3]]
##      [,1] [,2]
## [1,]    1    3
## [2,]    2    4
## 
## [[4]]
## , , 1
## 
##      [,1] [,2]
## [1,]    1    3
## [2,]    2    4
## 
## , , 2
## 
##      [,1] [,2]
## [1,]    5    7
## [2,]    6    8
## 
## 
## [[5]]
##   id age
## 1  1  20
## 2  2  30
## 3  3  40
## 
## [[6]]
## [1] A  O  A  AB O  B 
## Levels: A AB B O
  • 리스트(list)의 원소를 추출하는 방법
    • 대괄호([index]) 하나를 쓰는 것: 대괄호를 하나를 쓰면 최종적인 결과가 리스트가 되며,
    • 대괄호를 두 개([[index]]: 대괄호를 두 개를 쓰면 최종적인 결과가 해당 원소의 데이터 형태가 된다
result[2]
## [[1]]
##  [1]  1  2  3  4  5  6  7  8  9 10
result[[2]]
##  [1]  1  2  3  4  5  6  7  8  9 10

생초보 시작

기본 설정

R을 시작하게 되면, 처음으로 작업 디렉토리를 설정해줘야합니다. 먼저, D 드라이브에 ‘Rclass’라는 폴더를 만들어줍니다. 설정 방법은 ’setwd(“디렉토리”)’ 입니다.

setwd("D:/Rclass/")

설정된 작업 디렉토리를 확인하는 방법은 ‘getwd()’ 입니다.

getwd()
## [1] "D:/Rclass"

위와 같이 “D:/Rclass/”가 작업 디렉토리로 설정된 것을 볼 수 있습니다.

다음으로는 분석에 필요한 라이브러리를 설치해줍니다. 라이브러리의 종류는 수백가지라, 각각의 목적에 맞는 패키지를 검색 후 설치해야합니다. 대기과학에서 주로 사용한다고 말하기는 힘들지만, 제가 기본으로 불러오는 패키지는 lubridate와 ‘plyr’ 입니다.

install.packages("lubridate")
install.packages("plyr")

’lubridate’는 시계열 자료를 다루는데 매우 유용하며, ’plyr’은 자료의 반올림, 내림, 버림을 하는 데 유용합니다.

앞서 ‘install.packages’를 실행하여 해당 패키지를 설치했다면, 앞으로는 ’library’ 명령어를 이용해서 패키지를 불러오면 됩니다.

library("lubridate")
## 
## 다음의 패키지를 부착합니다: 'lubridate'
## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union
library("plyr")

오류 메세지 없이 불러왔다면, 이제 문제 없이 해당 패키지를 사용할 수 있습니다.

외부 데이터 읽어오기

텍스트 데이터

  • 텍스트 데이터는 데이터와 데이터를 구분하기 위해서 공백(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
  • CSV 데이터
DAT = read.csv(file = “파일명”, header = TRUE)

엑셀 데이터

install.packages("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)

본 수업에서는 과학원 항공측정 자료 일부 (flight_aerosol.csv, flight_gas.csv)를 이용하여 실습해본다. 두 파일 모두 csv 파일이기에 read.csv 명령어를 사용하여 불러온다. 각각 aero와 gas라는 이름에 할당한다.

aero = read.csv(file =  "flight_aerosol.csv", header = TRUE)#, fileEncoding = "WINDOWS-1252")
gas = read.csv(file = "flight_gas.csv", header = TRUE)#, fileEncoding = "WINDOWS-1252")

데이터 핸들링

설명 방법 사용법
행의 개수 nrow(데이터명) nrow(aero)
열의 개수 ncol(데이터명) ncol(aero)
행의 이름 row.names(데이터명) row.names(aero)
열의 이름 colnames(데이터명) colnames(aero)
행과 열이 갯수 dim(데이터명) dim(aero)
데이터의 구조 str(데이터명) str(aero)
  • Tip
paste("R", 1:53940, sep = "")

데이터 추출하기

  • 데이터의 일부 보기
  • head(데이터명), tail(데이터명)
head(aero)
##                G_Time        ORG          SO4          NO3          NH4
## 1 2020-12-09 12:06:42 0.13565736 -0.005863614 -0.004942043  0.001469853
## 2 2020-12-09 12:06:43 0.11113925 -0.010245030  0.004178098  0.000242704
## 3 2020-12-09 12:06:44 0.02192263 -0.010380049 -0.006855341 -0.000740776
## 4 2020-12-09 12:06:45 0.06655706 -0.001110576  0.007406887  0.001082767
## 5 2020-12-09 12:06:46 0.09853335 -0.003626055  0.002752427  0.001341103
## 6 2020-12-09 12:06:47 0.07565946  0.002965167  0.003072852  0.001166844
##            CHl 풍상풍하 rBC  SC WD_marker
## 1  0.014723084        n NaN NaN        NA
## 2 -0.028423049        n  NA  NA        NA
## 3  0.041331194        n  NA  NA        NA
## 4 -0.005300041        n  NA  NA        NA
## 5 -0.013113281        n  NA  NA        NA
## 6  0.009039521        n  NA  NA        NA
head(aero, n = 3)
##                G_Time        ORG          SO4          NO3          NH4
## 1 2020-12-09 12:06:42 0.13565736 -0.005863614 -0.004942043  0.001469853
## 2 2020-12-09 12:06:43 0.11113925 -0.010245030  0.004178098  0.000242704
## 3 2020-12-09 12:06:44 0.02192263 -0.010380049 -0.006855341 -0.000740776
##           CHl 풍상풍하 rBC  SC WD_marker
## 1  0.01472308        n NaN NaN        NA
## 2 -0.02842305        n  NA  NA        NA
## 3  0.04133119        n  NA  NA        NA
tail(aero)
##                   G_Time ORG SO4 NO3 NH4 CHl 풍상풍하 rBC SC WD_marker
## 8327 2020-12-09 14:25:28  NA  NA  NA  NA  NA        n  NA NA        NA
## 8328 2020-12-09 14:25:29  NA  NA  NA  NA  NA        n  NA NA        NA
## 8329 2020-12-09 14:25:30  NA  NA  NA  NA  NA        n  NA NA        NA
## 8330 2020-12-09 14:25:31  NA  NA  NA  NA  NA        n  NA NA        NA
## 8331 2020-12-09 14:25:32  NA  NA  NA  NA  NA        n  NA NA        NA
## 8332 2020-12-09 14:25:33  NA  NA  NA  NA  NA        n  NA NA        NA
tail(aero, n = 3)
##                   G_Time ORG SO4 NO3 NH4 CHl 풍상풍하 rBC SC WD_marker
## 8330 2020-12-09 14:25:31  NA  NA  NA  NA  NA        n  NA NA        NA
## 8331 2020-12-09 14:25:32  NA  NA  NA  NA  NA        n  NA NA        NA
## 8332 2020-12-09 14:25:33  NA  NA  NA  NA  NA        n  NA NA        NA
  • 열 추출하기
aero[1:10, 1] # aero 첫번째 열 중 1-10행까지만 벡터 형태로 출력
##  [1] "2020-12-09 12:06:42" "2020-12-09 12:06:43" "2020-12-09 12:06:44"
##  [4] "2020-12-09 12:06:45" "2020-12-09 12:06:46" "2020-12-09 12:06:47"
##  [7] "2020-12-09 12:06:48" "2020-12-09 12:06:49" "2020-12-09 12:06:50"
## [10] "2020-12-09 12:06:51"
aero[1:10, 1, drop=FALSE] # 데이터 프레임 형태로 출력
##                 G_Time
## 1  2020-12-09 12:06:42
## 2  2020-12-09 12:06:43
## 3  2020-12-09 12:06:44
## 4  2020-12-09 12:06:45
## 5  2020-12-09 12:06:46
## 6  2020-12-09 12:06:47
## 7  2020-12-09 12:06:48
## 8  2020-12-09 12:06:49
## 9  2020-12-09 12:06:50
## 10 2020-12-09 12:06:51
aero[1:10 , c(1, 3, 4)] # 1, 3, 4번째 열만 출력
##                 G_Time          SO4          NO3
## 1  2020-12-09 12:06:42 -0.005863614 -0.004942043
## 2  2020-12-09 12:06:43 -0.010245030  0.004178098
## 3  2020-12-09 12:06:44 -0.010380049 -0.006855341
## 4  2020-12-09 12:06:45 -0.001110576  0.007406887
## 5  2020-12-09 12:06:46 -0.003626055  0.002752427
## 6  2020-12-09 12:06:47  0.002965167  0.003072852
## 7  2020-12-09 12:06:48  0.003797006 -0.003239110
## 8  2020-12-09 12:06:49  0.000301463 -0.005191548
## 9  2020-12-09 12:06:50  0.001864705  0.003566221
## 10 2020-12-09 12:06:51  0.012906581 -0.000038600
aero[1:10 , 2:5] # 2부터 다섯번째 열까지 출력
##           ORG          SO4          NO3          NH4
## 1  0.13565736 -0.005863614 -0.004942043  0.001469853
## 2  0.11113925 -0.010245030  0.004178098  0.000242704
## 3  0.02192263 -0.010380049 -0.006855341 -0.000740776
## 4  0.06655706 -0.001110576  0.007406887  0.001082767
## 5  0.09853335 -0.003626055  0.002752427  0.001341103
## 6  0.07565946  0.002965167  0.003072852  0.001166844
## 7  0.01100334  0.003797006 -0.003239110  0.001834071
## 8  0.04356024  0.000301463 -0.005191548  0.000703719
## 9  0.13631469  0.001864705  0.003566221  0.001714683
## 10 0.03757255  0.012906581 -0.000038600  0.003916938
aero[1:10 , seq(from = 1, to = ncol(aero), by = 2)] #첫번째 열부터 마지막 열까지 홀수만 출력
##                 G_Time          SO4          NH4 풍상풍하  SC
## 1  2020-12-09 12:06:42 -0.005863614  0.001469853        n NaN
## 2  2020-12-09 12:06:43 -0.010245030  0.000242704        n  NA
## 3  2020-12-09 12:06:44 -0.010380049 -0.000740776        n  NA
## 4  2020-12-09 12:06:45 -0.001110576  0.001082767        n  NA
## 5  2020-12-09 12:06:46 -0.003626055  0.001341103        n  NA
## 6  2020-12-09 12:06:47  0.002965167  0.001166844        n  NA
## 7  2020-12-09 12:06:48  0.003797006  0.001834071        n  NA
## 8  2020-12-09 12:06:49  0.000301463  0.000703719        n  NA
## 9  2020-12-09 12:06:50  0.001864705  0.001714683        n  NA
## 10 2020-12-09 12:06:51  0.012906581  0.003916938        n  NA
aero[1:10 , c("SO4", "NO3", "CHl")] # 다음과 같은 열 이름만 출력력
##             SO4          NO3          CHl
## 1  -0.005863614 -0.004942043  0.014723084
## 2  -0.010245030  0.004178098 -0.028423049
## 3  -0.010380049 -0.006855341  0.041331194
## 4  -0.001110576  0.007406887 -0.005300041
## 5  -0.003626055  0.002752427 -0.013113281
## 6   0.002965167  0.003072852  0.009039521
## 7   0.003797006 -0.003239110  0.014645658
## 8   0.000301463 -0.005191548  0.009893852
## 9   0.001864705  0.003566221 -0.004059999
## 10  0.012906581 -0.000038600  0.003152022
  • grep: 변수명 중에서 특정 문자로 시작하거나 끝나거나 포함하고 있는 것을 추출할 때
argument 설명
pattern 특정한 문자로 시작하거나 끝나거나 포함하고 있는지를 문자형으로 표현한다. a로 시작하는 것은 “^a”, a로 끝나는 것은 “a$”, a를 포함하는 것은 “a”로 지정한다.x 문자형 벡터를 지정한다.
value 논리형으로 TRUE으로 pattern를 만족하는 벡터의 값을 반환하고,FALSE를 지정하면 pattern를 만족하는 벡터의 위치(index)를 반환한다. 
grep("O", colnames(aero), value = TRUE) # 헤더이름 중 O가 들어가는 이름
## [1] "ORG" "SO4" "NO3"
grep("O", colnames(aero), value = FALSE) # 헤더이름 중 O가 들어가는 순서 
## [1] 2 3 4
aero[1:10, grep("O", colnames(aero), value = TRUE)] # 헤더이름 중 O가 들어가는 이름 중 10열까지 출력
##           ORG          SO4          NO3
## 1  0.13565736 -0.005863614 -0.004942043
## 2  0.11113925 -0.010245030  0.004178098
## 3  0.02192263 -0.010380049 -0.006855341
## 4  0.06655706 -0.001110576  0.007406887
## 5  0.09853335 -0.003626055  0.002752427
## 6  0.07565946  0.002965167  0.003072852
## 7  0.01100334  0.003797006 -0.003239110
## 8  0.04356024  0.000301463 -0.005191548
## 9  0.13631469  0.001864705  0.003566221
## 10 0.03757255  0.012906581 -0.000038600
  • substr: 문자 중에서 일부를 추출. 엑셀의 mid 함수와 동일
argument 설명
x 문자형 벡터를 지정한다.
start 추출할 문자의 첫 번째 위치를 정수형으로 지정.
stop 추출한 문자의 마지막 위치를 정수형으로 지정. 
substr(colnames(aero), start = 1, stop = 2) # 변수명 중에서 처음 두 개의 문자를 추출
##  [1] "G_"   "OR"   "SO"   "NO"   "NH"   "CH"   "풍상" "rB"   "SC"   "WD"
substr(colnames(aero), 1, 2) # 변수명 중에서 처음 두 개의 문자를 추출
##  [1] "G_"   "OR"   "SO"   "NO"   "NH"   "CH"   "풍상" "rB"   "SC"   "WD"
aero[1:10 , substr(colnames(aero), start = 1, stop = 2) == "SO"] # 변수명 중에서 처음 두 개의 문자를 추출한 것이 “SO”와같은 열을 추출 (1-10행)
##  [1] -0.005863614 -0.010245030 -0.010380049 -0.001110576 -0.003626055
##  [6]  0.002965167  0.003797006  0.000301463  0.001864705  0.012906581

원하는 케이스의 행 추출

엑셀의 필터링 기능
  • aero 데이터에서 배출원 기준 풍상 지역과 풍하 지역을 나타내는 ‘풍상풍하’ 변수가 있다. 풍상풍하에는 “u (upwind)”, “d (downwind)”, “n (not relevant)”가 있다. 이중에서 각각을 추출하여 ‘upwind’, ‘downwind’, ‘nr’ 이라는 변수로 할당
upwind = aero[aero$풍상풍하 == "u" , ] # 풍상풍하열 값 중 'u'만 추출하여 upwind에 할당
downwind = aero[aero$풍상풍하 == "d" , ]
nr = aero[aero$풍상풍하 == "n" , ]
  • factor를 이용하면 그룹에 따른 boxplot을 그릴 수 있다.
boxplot(aero$ORG ~ aero$풍상풍하) #factor를 이용한 boxplot

* 이와 마찬가지로 원하는 조건에 따라 추출 가능

aero[aero$ORG >= 5, ] # Organic이 5 ug m-3 이상인 자료
aero[(aero$풍상풍하 == 'u') & (aero$ORG >= 5), ] # Organic이 5 ug m-3 이상이고 풍상지역인 자료
aero[(aero$풍상풍하 == 'u') | (aero$ORG >= 5), ] # Organic이 5 ug m-3 이상이거나 풍상지역인 자료
aero[(aero$풍상풍하 == 'u') | (aero$ORG >= 5), 2:3] # Organic이 5 ug m-3 이상이거나 풍상지역인 자료 중 2번째, 3번째 컬럼
aero[(aero$풍상풍하 == 'u') | (aero$ORG >= 5), c("ORG", "SO4")] # Organic이 5 ug m-3 이상이거나 풍상지역인 자료 중 "ORG", "SO4" 컬럼

시계열 자료 처리

  • aero와 gas의 첫 번째 행 (G_Time)의 형식은 모두 YYYY-mm-dd HH:MI:SS의 형태임
aero$G_Time[1]
gas$G_Time[1]
  • 위 형태는 R에서 시간 형태가 아닌 텍스트 형태로 인식했기 때문에, ‘lubridate’ 라이브러리의 ‘as.POSIXct()’ 한수를 이용하여 시간 형태로 변경시켜줘야 함.
argument 설명
x 변경할 변수
tz 사용자가 원하는 time zone. 미지정 시 “GMT”. 한국은 “Asia/Seoul”. 
aero$G_Time = as.POSIXct(aero$G_Time, tz = "Asia/Seoul")
gas$G_Time = as.POSIXct(gas$G_Time, tz = "Asia/Seoul")
aero$G_Time[1]
## [1] "2020-12-09 12:06:42 KST"

  • 뒤에 KST가 붙은 것을 알 수 있다. 이제 R은 해당 텍스트를 시간형식으로 인식한다.

  • 해당 항공 측정자료는 초단위로 측정되었다. 만약 일정한 시간단위로 평균을 내고 싶다면, ‘plyr’ 패키지 중 round_any 명령어를 사용하여 새로운 ’시간 자료를 생성한다. 만약 분 단위면 60 (초), 시간 단위면 60(초)*60(분)이다.

as.numeric(aero$G_Time[1])
## [1] 1607483202
aero$G_Time2 = round_any(aero$G_Time, 60, trunc) # trunc = 0분 0-59초 = 0분, ceiling = 0분 0-59초 = 1분
aero$G_Time2[1]
## Warning in as.POSIXlt.POSIXct(x, tz): unknown timezone 'KST'
## [1] "2020-12-09 03:06:00 GMT"
  • 변경된 G_Time2는 KST가 아닌 GMT로 출력되어 이를 as.POSIXlt()를 사용하여 KST로 바꿔준다.
aero$G_Time2 = as.POSIXlt(aero$G_Time2, tz = "Asia/Seoul")
aero$G_Time2[1]
## [1] "2020-12-09 12:06:00 KST"
gas$G_Time2 = round_any(gas$G_Time, 60, trunc) # gas 역시 동일하게 처리
gas$G_Time2 = as.POSIXlt(gas$G_Time2, tz = "Asia/Seoul")
  • 1분 평균 계산 (엑셀의 부분합): aggregate 함수를 사용하여 1분 평균의 부분합을 계산함.
argument 설명
x 계산할 항목
by 그룹화할 항목
FUN 계산할 함수 (mean, median, max, min, …)
na.action 자료 중에 NA가 포함되었을 시 처리 방법 
mi.aero = aggregate(aero[,2:6], by = list(format(aero$G_Time2, "%Y-%m-%d %H:%M:00")), FUN = mean, na.rm = T) # 부분합 결과를 mi.aero에 할당
mi.aero$Group.1[1]
## [1] "2020-12-09 12:06:00"
mi.aero$Group.1 = as.POSIXct(mi.aero$Group.1, tz= "Asia/Seoul") # mi.aero$Group.1을 시간형식으로 변환 as.POSIXlt 함수 사용용
mi.aero$Group.1[1]
## [1] "2020-12-09 12:06:00 KST"
mi.gas = aggregate(gas[,2:36], by = list(format(gas$G_Time2, "%Y-%m-%d %H:%M:00")), FUN = mean, na.rm = T)
mi.gas$Group.1 = as.POSIXct(mi.gas$Group.1)
mi.gas$Group.1[1]
## [1] "2020-12-09 12:06:00 KST"
  • aero와 gas 자료 합치기 (merge 함수 사용; 엑셀의 vlookup): 기준이 되는 열을 지정후 두 자료를 병합. 엑셀과 달리 기준이 되는 열을 여러개 지정가능함.
argument 설명
x, y 합칠 데이터 셋
by 기준이 될 열 이름
all 겹치지 않는 행이 존재할 시 살려둘지, 삭제할지를 정함
all.x x 데이터 셋의 겹치지 않는 행을 살려둠 (TRUE일 경우)
all.y y 데이터 셋의 겹치지 않는 행을 살려둠 (TRUE일 경우) 
flight = merge(mi.gas, mi.aero, by = "Group.1", all = T) #분단위 평균의 gas와 aero 자료를 병합
flight.all = merge(gas, aero, by = "G_Time", all = T) #초단위 모든 gas와 aero 자료를 병합

내부 데이터 저장

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

그래프 그리기

  • 기본적으로 그래프는 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(flight.all$ORG ~ flight.all$G_Time) # (1): 일반적인 점 그래프

plot(flight.all$ORG ~ flight.all$G_Time, type = "l") # (2): 일반적인 선 그래프

plot(flight.all$ORG ~ flight.all$G_Time, type = "l", ylim = c(0,12)) # (3): (2) + y축 범위 변경

plot(flight.all$ORG ~ flight.all$G_Time, type = "l", ylim = c(0,12), col = "darkorange") # (4): (3) + 색깔 변경

plot(flight.all$ORG ~ flight.all$G_Time, type = "l", ylim = c(0,12), col = "darkorange", xlab = "Time", ylab = "Organic") # (5): (4) + 이름 변경

plot(flight.all$ORG ~ flight.all$G_Time, type = "l", ylim = c(0,12), col = "darkorange", xlab = "Time", ylab = expression("PM"[1]*" Organic ("*mu*"g m"^-3*")")) # (6): (5) + 특수문자 및 아래첨자 
{plot(flight.all$ORG ~ flight.all$G_Time, type = "l", ylim = c(0,12), col = "darkorange", xlab = "Time", ylab = expression("PM"[1]*" Organic ("*mu*"g m"^-3*")")) 
points(flight.all$ORG ~ flight.all$G_Time, col = "darkorange", pch = 22, cex = 0.7)
legend("topright", legend = "Organic", col = "darkorange", pch = 22, cex = 0.7)} # (7): (6) + 심볼 추가 + 레전드 추가 

{plot(flight.all$ORG ~ flight.all$G_Time, type = "l", ylim = c(0,12), col = "darkorange", xlab = "Time", ylab = expression("PM"[1]*" Organic ("*mu*"g m"^-3*")")) 
points(flight.all$ORG ~ flight.all$G_Time, col = "darkorange", pch = 22, cex = 0.7)

points(flight.all$SO4 ~ flight.all$G_Time, col = "blue", pch = 21, cex = 0.7) # (8): (7) + SO4 자료 추가
lines(flight.all$SO4 ~ flight.all$G_Time, col = "blue")
legend("topright", legend = c("Organic", "Sulfate"), col = c("darkorange", "blue"), pch = c(22, 21), cex = 0.7)}

  • y축라벨과 숫자의 간격을 줄이는 방법
{plot(flight.all$ORG ~ flight.all$G_Time, type = "l", ylim = c(0,12), col = "darkorange", xlab = "Time", ylab = "")
mtext(text = expression("PM"[1]*" Organic ("*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_org.png", width = 6, height = 3, units = "in", res = 200) # 저장할 스케치북 열기
  plot(flight.all$ORG ~ flight.all$G_Time, type = "l", ylim = c(0,12), col = "darkorange", xlab = "Time", ylab = "")
  mtext(text = expression("PM"[1]*" Organic ("*mu*"g m"^-3*")"), line = 2, side = 2)
  dev.off() # 스케치북 닫기
} # (7): (6) 그래프 저장
## png 
##   2

산점도 및 추세선 추가

  • lm (linear model): 선형 회귀직선
  • abline: 선 추가 명령어
sl = lm(flight.all$ORG ~ flight.all$SO4) # y ~ x의 방정식 (y = bx + a)
(sl2 = summary(sl)) # sl 변수의 종합적 통계
## 
## Call:
## lm(formula = flight.all$ORG ~ flight.all$SO4)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.1148 -0.4811 -0.2413  0.2963  9.0917 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     0.34381    0.02482   13.85   <2e-16 ***
## flight.all$SO4  0.62229    0.01235   50.41   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6837 on 5519 degrees of freedom
##   (결측으로 인하여 2811개의 관측치가 삭제되었습니다.)
## Multiple R-squared:  0.3153, Adjusted R-squared:  0.3151 
## F-statistic:  2541 on 1 and 5519 DF,  p-value: < 2.2e-16
{plot(flight.all$ORG ~ flight.all$SO4, type = "p", col = "darkorange", xlab = "", ylab = "") #ORG와 SO4 산점도. 축 라벨은 없음
mtext(text = expression("PM"[1]*" Organic ("*mu*"g m"^-3*")"), line = 2, side = 2) # y축 라밸 추가
mtext(text = expression("PM"[1]*" SO"[4]*" ("*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) # 추세식 표기
}

지도위에 자료 표출하기

위, 경도 자료를 이용하여 지도위에 표시하기

  • maps 라이브러리 호출
library("maps") # 단순한 지도를 그릴 수 있는 라이브러리 불러오기
## 
## 다음의 패키지를 부착합니다: 'maps'
## The following object is masked from 'package:plyr':
## 
##     ozone
{plot(flight.all$Lat ~ flight.all$Long, type = "p", col = "darkorange", xlab = "", ylab = "", xlim = c(126, 127)) # 위경도에 따른 심볼 그림
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) # 지도 추가
}

농도에 따른 색깔을 표시

  • squash 라이브러리 호출
library("squash") # color bar를 줄 수 있는 라이브러리 호출
(zlim = range(flight.all$ORG, na.rm = T)) # ORG의 농도 범위 파악
## [1] -0.1826231 11.0010180
(zlim[2] = round_any(zlim[2]/3, accuracy = 1, f = trunc))
## [1] 3
col.map = makecmap(zlim,  n = 100, colFn = jet, outlier = 'red') # 팔레트 제작 jet, heat, bluered, blueorange, rainbow, greyscale, coolheat
flight.all$ORG1 = flight.all$ORG
flight.all[flight.all$ORG > zlim[2] & is.na(flight.all$ORG) == F, "ORG1"] = zlim[2]

{plot(flight.all$Lat ~ flight.all$Long, type = "p", xlab = "", ylab = "", xlim = c(126, 127), # 위경도 그림
      col = cmap(flight.all$ORG1, col.map), main = expression("PM"[1]*" Organic ("*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축 이름
map("world", add = T) # 지도 추가
vkey(col.map, paste("Conc."),y = 36.6, x = 127, stretch = 0.05) # 레전드 추가
}

구글 지도 위에 표시

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

Google Map을 사용한 위경도 그림 

qmap( loc , zoom = 9, maptype = "terrain") + # terrain
  geom_point(aes(x = Long, y = Lat), data = flight.all, col = cmap(flight.all$ORG1, map = col.map))

qmap( loc , zoom = 9, maptype = "hybrid") + # hybrid
  geom_point(aes(x = Long, y = Lat), data = flight.all, col = cmap(flight.all$ORG1, map = col.map))

qmap( loc , zoom = 9, maptype = "roadmap") + # roadmap
  geom_point(aes(x = Long, y = Lat), data = flight.all, col = cmap(flight.all$ORG1, map = col.map))

3차원 그래프

library("rgl") # OPEN GL (open graphics library) 2D, 3D 렌더링을 해주는 라이브러리
par3d(windowRect = c(10, 10, 800, 800)) # 3D 윈도우 생성
plot3d(flight.all$Long, flight.all$Lat, flight.all$Z, col = cmap(flight.all$ORG1, map = col.map), 
       size = 2, type = 'p', xlab = "Lon", ylab = "Lat", zlab = "Alt", zlim = c(0,1000))
# movie3d(spin3d(axis = c(0, 0, 0.5), rpm = 5), duration = 20, dir = getwd())
# system("cmd.exe", input = paste("cd", getwd(), "&& convert -loop 0 movie.gif philae.gif")) # change GIF animation cycles to inf
# rgl.snapshot("Fig/vertical_3d.png", fmt="png", top = TRUE)

RNetCDF, OceanView 패키지

netcdf 자료 읽기와 시각화 하기

  • netcdf는 3차원 이상 형태의 자료 구조. 일반적으로 시간+위도+경도의 차원으로 구성됨
NetCDF

NetCDF 

install.packages("RNetCDF")
install.packages("OceanView")
  • ‘RNetCDF’ 및 ‘OceanView’ 라이브러리 불러오기
library("RNetCDF")
library("OceanView")
## 필요한 패키지를 로딩중입니다: plot3D
## 필요한 패키지를 로딩중입니다: plot3Drgl
library("maps")
  • 폴더 안에 있는 입력 파일 읽기
nc.path = "nc_file/" # 입력 파일의 위치
omi_list = list.files(nc.path, "omi") # omi라는 단어가 들어간 파일명
era_list = list.files(nc.path, "era") # era라는 단어가 들어간 파일명

N2019 = open.nc(paste0(nc.path, era_list[1])) # 2019년 ERA5 파일 읽기
N2020 = open.nc(paste0(nc.path, era_list[2])) # 2020년 ERA5 파일 읽기

o2019 = open.nc(paste0(nc.path, omi_list[1])) # 2019년 OMI 파일 읽기
o2020 = open.nc(paste0(nc.path, omi_list[2])) # 2020년 OMI 파일 읽기

# print.nc(N2020) # N2020 내 변수명 확인 
# print.nc(o2019) # o2019 내 변수명 확인
  • 변수 가져오기
# netcdf 파일의 경우, 자료 용량을 줄이기 위해, 실수가 아닌 정수로 자료를 저장 후, 기울기 (scale factor)와 편차 (add_offset)을 사용하여 변환해야함. 따라서 해당 식을 함수로 만들어 사용함.
getData = function(DATA, varname){
  sclt = var.get.nc(DATA, varname)
  sclt[sclt == -32768] = NA
  sclt = sclt * att.get.nc(DATA, varname, "scale_factor") + att.get.nc(DATA, varname, "add_offset")
}

lat   = var.get.nc(N2019, "latitude") # 변수 중 위도 가져오기
lon   = var.get.nc(N2019, "longitude") # 경도 가져오기
time = var.get.nc(N2019, "time") # 시간 가져오기

xwind1 = getData(N2019, "u10") # 2019년 10 m u 벡터 
xwind2 = getData(N2020, "u10") # 2020년 10 m u 벡터 

ywind1 = getData(N2019, "v10") # 2019년 10 m v 벡터 
ywind2 = getData(N2020, "v10") # 2020년 10 m v 벡터 

no2.1 = var.get.nc(o2019, "OMNO2d_003_ColumnAmountNO2TropCloudScreened") # 2019년대류권 NO2 컬럼 VCD
no2.2 = var.get.nc(o2020, "OMNO2d_003_ColumnAmountNO2TropCloudScreened") # 2020년대류권 NO2 컬럼 VCD
lat.omi   = var.get.nc(o2019, "lat")
lon.omi   = var.get.nc(o2019, "lon")
  
no2 = (no2.2 - no2.1)/1e15 # 2020-2019년 NO2 차이. 
dim(no2) # no2의 크기 확인
## [1] 200 120
xwind.19 = apply(xwind1 , c(1,2), mean, na.rm=T) # ERA5는 월별 자료이기에 연평균 계산
ywind.19 = apply(ywind1 , c(1,2), mean, na.rm=T)
xwind.20 = apply(xwind2 , c(1,2), mean, na.rm=T)
ywind.20 = apply(ywind2 , c(1,2), mean, na.rm=T)

F2019 = quiver2D(xwind.19, ywind.19, x = lon, y = lat, speed.max = 5, # 2019년 바람장 
                 by = 10, scale = 1, arr.max = 0.1, col="black", 
                 xlim=c(110,145), ylim=c(25,45))

F2020 = quiver2D(xwind.20, ywind.20, x = lon, y = lat, speed.max = 5, # 2020년 바람장 (표시 안함)
                 by = 10, scale = 1, arr.max = 0.1, col="black", 
                 xlim=c(110,145), ylim=c(25,45), plot = FALSE )

ct1 = F2020$length - F2019$length # 2020년과 2019년 풍속의 차이
ct1[ct1 >= 0] = "red" # 풍속이 빨라졌으면 적색
ct1[ct1 < 0] = "blue" # 풍속이 느려졌으면 청색

rng = c(-2, 2) # colormap 설정
no2[no2 > rng[2]] = rng[2] # colormap 최대값보다 크면 치환
no2[no2 < rng[1]] = rng[1] # colormap 최소값보다 작으면 치환
par(mai = c(0.8,0.8,0.5,0.7)) # 그림의 여백 조정

F1 = quiver2D(xwind.20, ywind.20, x = lon, y = lat, speed.max = 5, # 2020년 바람장을 그리고 풍속의 증감은 색깔로 표시
              by = 5, scale = 2, arr.max = 0.1, col = ct1, 
              image = list(z = no2, col = jet.col(100), x = lon.omi, y = lat.omi, clim = rng), # 2020년과 2019년 omi NO2 농도 차이를 이미지로 투영
              main = "", xlim=c(100, 145), ylim=c(20, 50), xlab = "", ylab = "")
mtext(expression(paste("Longitude ",degree*E)), side = 1, line = 2) # x축 라벨 추가
mtext(expression(paste("Latitude ",degree*N)), side = 2, line = 2) # y축 라벨 추가
map("world", add=T, lwd = 2, col = "black") # 지도 추가

colkey(clim = rng, add = TRUE, width = 0.5, length = 0.5,  # legend 추가
       cex.axis = 0.75, clab = expression(Delta*"NO"[2]), cex.clab = 0.75, line.clab = 0.5, 
       col = jet.col(100))
  
legend("bottomright", legend = paste("max ws= ", formatC(F1$speed.max, digits = 2), "m/s"), bg = "white", cex = 0.7) # 풍속 최대값 추가
mtext(outer = TRUE, text = expression("OMI "*Delta*"NO"[2]*" TropVCD between 2020 and 2019"), side = 3, line = -2, cex = 1) # 제목 추가가

Openair

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

  • 본 수업에 사용되는 장기간 측정 자료는 아래 링크에서 다운받을 수 있다. 1_gangbuk.csv / 2_gwanak.csv / 3_yongsan.csv

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

install.packages("openair")
  • 패키지 불러오기: 기존에 설치가 되었거나, 설치가 완료되었다면, 아래 명령어로 lubridate 패키지와 함께 불러온다.
library("lubridate")
library("openair")
  • 자료 불러오기: csv 파일이기에 read.csv 명령어를 이용하여 자료를 불러들인다. 세 지점 자료를 각각 GB, GA, YS에 할당한다.
GB = read.csv(file = "1_gangbuk.csv", header = T)
GA = read.csv(file = "2_gwanak.csv", header = T)
YS = read.csv(file = "3_yongsan.csv", header = T)
  • 세 지점의 데이터를 하나의 데이터로 합치기: 단, 합치게 되면 지점에 대한 정보를 모르기에, 합치기전 지점의 정보를 입력해주고 AQM이라는 데이터 셋에 합친다. 합칠 때애는 rbind라는 명령어를 사용 (주의할점은 합치려는 데이터셋의 column 이름이 모두 동일해야함)
GB$site = "gangbuk"
GA$site = "gwanak"
YS$site = "yongsan"
# 합치기
AQM = rbind(GB, GA, YS)
  • 시간형태로 변경: YYYYMMDDHH인 시간의 형태를 R이 인식할 수 있도록 작업함. (1) 연, 월, 일을 따로 구분하거나, (2) simple하게 as.POSIXct 명령어를 사용
# 복잡한 방법: substr 함수를 이용하여 텍스트을 연, 월, 일로 분리
yy = substr(AQM$YYYYMMDDHH, 1,4)
mm = substr(AQM$YYYYMMDDHH, 5,6)
dd = substr(AQM$YYYYMMDDHH, 7,8)
hh = substr(AQM$YYYYMMDDHH, 9,10)

AQM$date = paste(yy,mm,dd, sep ="-") # paste 함수를 이용하여 연-월-일 형태로 변경
AQM$date = as.POSIXct(AQM$date) + hours(hh) - hours(1) # as.POSIXct 함수를 사용하여 시간 형태로 변경하고, lubridate hour 함수를 사용하여 시간을 더해줌.

# 단순한 방법
AQM$date = as.POSIXct(as.character(AQM$YYYYMMDDHH), format = "%Y%m%d%H") - hours(1) #연, 월, 일을 인식시켜줌으로써 손쉽게 시간형태로 변경
  • 필요한 컬럼만 정리
colnames(AQM) #컬럼 이름 확인
##  [1] "YYYYMMDDHH" "TYPE"       "시도"       "도시"       "시군구"    
##  [6] "측정소명"   "측정소코드" "정렬코드"   "SO2"        "NO2"       
## [11] "O3"         "CO"         "PM10"       "PM2.5"      "Nox"       
## [16] "NO"         "풍향"       "풍향빈도"   "풍속"       "온도"      
## [21] "습도"       "현지기압"   "일사"       "강수량"     "해면기압"  
## [26] "site"       "date"
AQM = AQM[,c(27,9:26)] # 필요한 컬럼만 추출
for (nc in 2:ncol(AQM)){ # for 반복문을 활용하여 -90이하의 값을 NA로 치환함. ncol(AQM)은 AQM 컬럼 수임.
  ID =  which(is.na(AQM[,nc]) == T | AQM[,nc] < -90) # which는 조건에 부합하는 자료를 파악.
  AQM[ID,nc] = NA # 해당 자료를 NA로 변경경
}

summary plot

  • 시계열 자료의 전체적인 경향을 파악하기 좋은 그래프. 자료 및 결측치의 분포, 평균, 중간값 등에 대한 정보를 제공함함
YS = AQM[AQM$site == "yongsan",] # 용산의 자료만을 추출
summaryPlot(YS) # 용산의  summary plot
##       date1       date2         SO2         NO2          O3          CO 
##   "POSIXct"    "POSIXt"   "numeric"   "numeric"   "numeric"   "numeric" 
##        PM10       PM2.5         Nox          NO        풍향    풍향빈도 
##   "integer"   "integer"   "numeric"   "numeric"   "integer" "character" 
##        풍속        온도        습도    현지기압        일사      강수량 
##   "numeric"   "numeric"   "numeric"   "numeric"   "logical"   "numeric" 
##    해면기압        site 
##   "numeric" "character"

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

windrose

  • 풍향, 풍속에 대한 정보를 제공
colnames(AQM)[c(10,12)] = c("wd", "ws") # 풍향, 풍속의 컬럼 이름을 wd, ws로 교체. (교체 안할 시 일일이 지정해야하는 번거로움..)
YS = AQM[AQM$site == "yongsan",] # 용산의 자료만을 추출
GA = AQM[AQM$site == "gwanak",] # 관악의 자료만을 추출
GB = AQM[AQM$site == "gangbuk",] # 강북의 자료만을 추출

windRose(YS) # 기본 바랑장미

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

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

### pollution rose

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

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

Polar plot

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

polarPlot(YS, pollutant = "PM2.5", statistic = "cpf", percentile = 75) # Conditional Probability Function

polarPlot(YS, pollutant = "PM2.5", statistic = "cpf", percentile = c(0,10)) # Conditional Probability Function

polarPlot(YS, pollutant = "PM2.5", statistic = "cpf", percentile = c(50,60)) # Conditional Probability Function

polarPlot(YS, pollutant = "PM2.5", type = "daylight") # ”year”, ”hour”, ”month”, ”season”, ”weekday”, ”weekend”, ”monthyear”, ”daylight”

polarPlot(AQM[AQM$site == "gangbuk",], pollutant = "PM2.5") # 강북지점의 polarPlot

polarPlot(AQM[AQM$site == "gwanak",], pollutant = "PM2.5") # 관악지점의 polarPlot

Calendar plot

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

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

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

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

Trend 분석 (smoothTrend, TheilSen)

smoothTrend(YS, pollutant = "NO2", deseason = TRUE, simulate = F,  ylab = "concentration (ppb)",
            main = "monthly mean deseasonalised NO2 (bootstrap uncertainties)") # smooth trend를 계산. GAM 모델을 사용하여 smooth를 수행

smoothTrend(YS, pollutant = "NO2", deseason = TRUE, simulate = F,  ylab = "concentration (ppb)",
            main = "monthly mean deseasonalised NO2 (bootstrap uncertainties)", type = "wd") # 풍향을 기준으로 분류
## Warning: 2603 missing wind direction line(s) removed

TheilSen(YS, pollutant = "PM2.5", deseason = TRUE,  ylab = "PM2.5", main = "Theil-Sen Trend") # Theil-sen 기울기 계산산
## [1] "Taking bootstrap samples. Please wait."

TheilSen(YS, pollutant = "PM2.5", deseason = TRUE,  ylab = "PM2.5", main = "Theil-Sen Trend", type = "season") # 계절절을 기준으로 분류
## [1] "Taking bootstrap samples. Please wait."
## [1] "Taking bootstrap samples. Please wait."
## [1] "Taking bootstrap samples. Please wait."
## [1] "Taking bootstrap samples. Please wait."

### timeVariation 분석

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

# 여러 항목들의 시간변화 비교 (절대값이 다르기에 nomarlized 적용)
timeVariation(YS, pollutant = c("PM2.5", "CO", "NO2", "O3"), normalise = TRUE)

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

그 외 유용한 함수

scatterPlot(YS, x = "PM10", y = "PM2.5", method = "density", col = "jet", type = "season")
## (loaded the KernSmooth namespace)

YS = rollingMean(YS, pollutant = "O3", hours = 8, new.name = "O3.8", data.thresh = 75) # 8시간 이동 평균
trendLevel(GB, pollutant = "PM2.5", cols = "jet", y = "hour")