대기환경학자를 위한 R의 기본

국립환경과학원 2024

R에 들어가기에 앞서..

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

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

R core와 Rstudio 설치하기

R core 설치

https://cran.yu.ac.kr/에 접속

해당 OS에 해당하는 인스톨 파일 다운로드

Rstudio 다운로드 페이지

Rstudio 설치

https://posit.co/download/rstudio-desktop/

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

Rstudio 다운로드 페이지

텍스트 에디터 notepad++ 설치하기

https://notepad-plus-plus.org/downloads/

맥이면, Sublime Text 추천

notepad++ 다운로드 페이지

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

Rstudio 실행화면

문법

기본 문법

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

• 데이터 프레임(Data Frame)
• 데이터 프레임은 R에서 가장 중요한 자료형이다. 데이터 프레임은 행렬과 마찬가지의 모습을 하고 있지만 행렬과 달리 다양한 변수, 관측치(observations), 범주 등을 표현하기 위해 특화되어있다.

data.frame() 사용

dat = data.frame ( x = c (1 , 2 , 3 , 4 , 5) , y = c (2 , 4 , 6 , 8 , 10) )
dat

##   x  y
## 1 1  2
## 2 2  4
## 3 3  6
## 4 4  8
## 5 5 10

dat = data.frame ( x = c (1 , 2 , 3 , 4 , 5) ,
                  y = c (2 , 4 , 6 , 8 , 10) ,
                  z = c ('M', 'F', 'M', 'F', 'M') ) # Factor 형태의 z 열이 추가된 예

dat

##   x  y z
## 1 1  2 M
## 2 2  4 F
## 3 3  6 M
## 4 4  8 F
## 5 5 10 M

dat$v = c (3 , 6 , 9 , 12 , 15)
dat

##   x  y z  v
## 1 1  2 M  3
## 2 2  4 F  6
## 3 3  6 M  9
## 4 4  8 F 12
## 5 5 10 M 15

데이터 프레임 접근

• 데이터 프레임의 각 열은 ’$변수명’으로 접근할 수 있으며, 행이나 열의 인덱스를 사용해서도 데이터에 접근할 수 있다.

dat$x

## [1] 1 2 3 4 5

dat$y

## [1]  2  4  6  8 10

dat$z

## [1] "M" "F" "M" "F" "M"

dat$v

## [1]  3  6  9 12 15

dat[1,2]

## [1] 2

dat[c(1,3),2]

## [1] 2 6

dat[-1,-c(2,3)] # 벡터로 인덱스를 지정하거나, 또는 제외할 행 또는 열을 - 로 표시할 수 있다.

##   x  v
## 2 2  6
## 3 3  9
## 4 4 12
## 5 5 15

dat[,c("x","y")] # 또는 컬럼명을 지정할 수도 있다.

##   x  y
## 1 1  2
## 2 2  4
## 3 3  6
## 4 4  8
## 5 5 10

• str()과 head()를 사용하여 간략하게 살펴볼 수 있다.
• str() 함수는 R 객체의 내부 구조 (데이터 유형)를 보는데 사용
• head() 함수는 R 객체의 상부 일부분을 보여줌
• tail() 함수는 R 객체의 하부 일부분을 보여줌

str(dat)

## 'data.frame':    5 obs. of  4 variables:
##  $ x: num  1 2 3 4 5
##  $ y: num  2 4 6 8 10
##  $ z: chr  "M" "F" "M" "F" ...
##  $ v: num  3 6 9 12 15

head(dat)

##   x  y z  v
## 1 1  2 M  3
## 2 2  4 F  6
## 3 3  6 M  9
## 4 4  8 F 12
## 5 5 10 M 15

head(dat, n = 3)

##   x y z v
## 1 1 2 M 3
## 2 2 4 F 6
## 3 3 6 M 9

tail(dat)

##   x  y z  v
## 1 1  2 M  3
## 2 2  4 F  6
## 3 3  6 M  9
## 4 4  8 F 12
## 5 5 10 M 15

tail(dat, n = 3)

##   x  y z  v
## 3 3  6 M  9
## 4 4  8 F 12
## 5 5 10 M 15

• 데이터 프레임의 행 이름, 열 이름은 각각 rownames(), colnames() 함수로 지정할 수 있다.
• head(데이터명), tail(데이터명)

rownames(dat)

## [1] "1" "2" "3" "4" "5"

rownames(dat) = c("a", "b", "c", "d", "e")
colnames(dat)

## [1] "x" "y" "z" "v"

colnames(dat) = c("ab", "bc", "cd", "da")

• 주어진 값이 벡터에 존재하는지를 판별하는 %in% 연산자를 이용하면 특정 열만 선택하는 작업을 보다 손쉽게 할 수 있다. 예를들어 다음과 같이 a, b, c 열이 있는 데이터 프레임에서 b, c 열만 선택하는 경우를 보자.

dat[,colnames(dat) %in% c("ab")]

## [1] 1 2 3 4 5

• 반대로 ! 연산자를 사용해 특정 값들만 제외해서 열을 선택할 수도 있다.

dat[,!colnames(dat) %in% c("ab")]

##   bc cd da
## a  2  M  3
## b  4  F  6
## c  6  M  9
## d  8  F 12
## e 10  M 15

• grep: 변수명 중에서 특정 문자로 시작하거나 끝나거나 포함하고 있는 것을 추출할 때

argument	설명
pattern	특정한 문자로 시작하거나 끝나거나 포함하고 있는지를 문자형으로 표현한다. a로 시작하는 것은 “^a”, a로 끝나는 것은 “a$”, a를 포함하는 것은 “a”로 지정한다.x 문자형 벡터를 지정한다.
value	논리형으로 TRUE으로 pattern를 만족하는 벡터의 값을 반환하고,FALSE를 지정하면 pattern를 만족하는 벡터의 위치(index)를 반환한다.

grep("a", colnames(dat), value = TRUE) # 헤더이름 중 O가 들어가는 이름

## [1] "ab" "da"

grep("a", colnames(dat), value = FALSE) # 헤더이름 중 O가 들어가는 순서 

## [1] 1 4

dat[, grep("a", colnames(dat), value = TRUE)] # 헤더이름 중 O가 들어가는 이름 중 10열까지 출력

##   ab da
## a  1  3
## b  2  6
## c  3  9
## d  4 12
## e  5 15

• substr: 문자 중에서 일부를 추출. 엑셀의 mid 함수와 동일

argument	설명
x	문자형 벡터를 지정한다.
start	추출할 문자의 첫 번째 위치를 정수형으로 지정.
stop	추출한 문자의 마지막 위치를 정수형으로 지정.

substr(colnames(dat), start = 1, stop = 1) # 변수명 중에서 처음 두 개의 문자를 추출

## [1] "a" "b" "c" "d"

substr(colnames(dat), 1, 1) # 변수명 중에서 처음 두 개의 문자를 추출

## [1] "a" "b" "c" "d"

dat[, substr(colnames(dat), start = 1, stop = 1) == "b"] # 변수명 중에서 처음 1개의 문자를 추출한 것이 “a”와같은 열을 추출

## [1]  2  4  6  8 10

원하는 케이스의 행 추출

엑셀의 필터링 기능

• dat 데이터에서 배출원 셩별을 나타내는 ‘cd’ 변수가 있다. 각각을 추출하여 ‘male’, ’female’이라는 데이터셋으로 할당

male = dat[dat$cd == "M" , ] 
female = dat[dat$cd == "F" , ]

• factor를 이용하면 그룹에 따른 boxplot을 그릴 수 있다.

boxplot(dat$ab ~ dat$cd) #factor를 이용한 boxplot

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

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

조건

R의 기본적인 제어문, 연산, 함수는 기본적으로 다른 언어와 유사하지만 저장된 데이터를 한번에 다루는 벡터연산을 수행한다는 점이 차별화임

if: if (조건) 참일 때 표현

x = 5
if (x > 0) "Positive number"

## [1] "Positive number"

if else: if(조건) 참일 때 표현 else 거짓일 때 표현

x = - 5
if (x > 0) "Positive number" else "Negative number"

## [1] "Negative number"

• 더 단순하게 아래 처럼 ifelse 함수로 표현도 가능함.

x = - 5
ifelse (x > 0, "Positive number", "Negative number")

## [1] "Negative number"

• 표현이 한 줄로 정리되지 않을 시, 아래와 같이 중괄호 { }에 묶어서 사용할 수 있음.

x = - 5
if (x > 0){
  "Positive number"
  } else {
    "Negative number"
}

## [1] "Negative number"

• 조건이 두개 이상일 시는 else if 함수로 추가 가능함.

x = 0
if (x > 0){
  "Positive number"
  } else if (x == 0) {
    "zero"
  } else {
    "Negative number"
  }

## [1] "zero"

반복문

for: for (반복할 변수명 in 반복할 내용) 표현

for (i in 1:5) print(1:i) # i를 1에서 5까지 반복하면서 1부터 1씩 증가하는 배열을 프린트 하라

## [1] 1
## [1] 1 2
## [1] 1 2 3
## [1] 1 2 3 4
## [1] 1 2 3 4 5

• 마찬가지로 중괄호를 이용하여 표현을 묶을 수 있다. 반복문의 경우, 조건문과 달리 print 혹은 cat 명령어를 쓰지 않으면 background에서 실행됨. 따라서 표출하고 싶다면 두 명령어를 사용해야함.

for (i in 1:5){
  print(1:i) # i를 1에서 5까지 반복하면서 1부터 1씩 증가하는 배열을 프린트 하라
}

## [1] 1
## [1] 1 2
## [1] 1 2 3
## [1] 1 2 3 4
## [1] 1 2 3 4 5

for(n in c(2,5,10,20,50)){
  x = stats::rnorm(n) # stats라는 패키지 안의 rnorm이라는 함수 사용 (Normal Distribution)
  cat(n, ": ", sum(x^2), "\n", sep = "") # n, :, x 제곱의 합, 줄 바꿈, 구분자는 "")
}

## 2: 4.929371
## 5: 6.591888
## 10: 9.960744
## 20: 19.29739
## 50: 61.77835

f = factor(sample(letters[1:5], 10, replace = TRUE)) #letters (a부터 z까지) 중 처음 다섯글자를 대상으로 무작위 샘플링 10개를 해서 factor로 변환
for(i in unique(f)) print(i) # f 중 겹치지 않는 unique한 값만 하나씩 출력해라

## [1] "a"
## [1] "c"
## [1] "b"
## [1] "d"
## [1] "e"

while: while (조건) 표현

i = 0
while ( i < 10) {
  print ( i )
  i = i + 1
  }

## [1] 0
## [1] 1
## [1] 2
## [1] 3
## [1] 4
## [1] 5
## [1] 6
## [1] 7
## [1] 8
## [1] 9

함수

function: ‘함수명 <- function(인자, 인자, ….) { 함수 본문 }’

fibo = function(n){
  if ( n == 1 | n == 2){
    x = n*2
  } else {
    x = n^2
  }
  return(x)
}
fibo(5)

## [1] 25

생초보 시작

기본 설정

R을 시작하게 되면, 처음으로 작업 디렉토리를 설정해줘야합니다. 먼저, D 드라이브에 ‘Rclass’라는 폴더를 만들어줍니다. 되도록이면, 한글이름을 사용하지 말고, 단순한 경로를 사용합니다. (바탕화면 노노..) 설정 방법은 ’setwd(“디렉토리”)’ 입니다.

setwd("~/Rclass/")

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

getwd()

## [1] "D:/OneDrive - hufs.ac.kr/Workspace_hufs/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)

본 강좌에서는 sample_data.xlsx를 사용할 예정입니다 (파일명이 빈 공백등은 되도록 노노..).

• 본 수업에 사용되는 입력자료는 아래 링크에서 다운받을 수 있습니다. sample_data.xlsx

library("readxl")
aero = read_excel(path =  "sample_data.xlsx", sheet = "Sheet1", col_names = 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)

##                  Date Temperature Wind direction Wind speed PM2.5       OC
## 1 2020-12-15 00:00:00       -12.3           21.5        0.2    10 1.944157
## 2 2020-12-15 01:00:00       -13.1          118.1        0.1    12 2.048968
## 3 2020-12-15 02:00:00       -13.6          125.1        0.3    13 1.805130
## 4 2020-12-15 03:00:00       -12.4           30.6        0.6    14 2.083666
## 5 2020-12-15 04:00:00       -13.7           95.3        0.5    13 1.766200
## 6 2020-12-15 05:00:00       -12.2          323.0        0.7    11 1.400701
##          EC   SO42-    NO3-     Cl-     Na+    NH4+      K+    Mg2+    Ca2+
## 1 0.6255070 1.40232 1.39549 0.23159 0.03549 0.85731 0.03115 0.07053 0.01698
## 2 0.6197119 1.45922 1.49809 0.22669 0.03399 0.89731 0.02605 0.05833 0.01718
## 3 0.5682204 1.48532 1.54819 0.22779 0.03159 0.92621 0.03265 0.04803 0.01488
## 4 0.5837869 1.54722 1.53099 0.24299 0.03029 0.91341 0.03495 0.04203 0.01028
## 5 0.5339275 1.45522 1.34259 0.25099 0.02999 0.78641 0.03495 0.02513 0.01728
## 6 0.4581114 1.39652 1.17579 0.28119 0.02659 0.74611 0.03425 0.01333 0.01998

head(aero, n = 3)

##                  Date Temperature Wind direction Wind speed PM2.5       OC
## 1 2020-12-15 00:00:00       -12.3           21.5        0.2    10 1.944157
## 2 2020-12-15 01:00:00       -13.1          118.1        0.1    12 2.048968
## 3 2020-12-15 02:00:00       -13.6          125.1        0.3    13 1.805130
##          EC   SO42-    NO3-     Cl-     Na+    NH4+      K+    Mg2+    Ca2+
## 1 0.6255070 1.40232 1.39549 0.23159 0.03549 0.85731 0.03115 0.07053 0.01698
## 2 0.6197119 1.45922 1.49809 0.22669 0.03399 0.89731 0.02605 0.05833 0.01718
## 3 0.5682204 1.48532 1.54819 0.22779 0.03159 0.92621 0.03265 0.04803 0.01488

tail(aero)

##                    Date Temperature Wind direction Wind speed PM2.5       OC
## 763 2021-01-15 18:00:00         5.1          297.0        2.1    59 7.341742
## 764 2021-01-15 19:00:00         3.4          249.2        0.9    47 6.241189
## 765 2021-01-15 20:00:00         2.5          353.6        0.6    44 5.638757
## 766 2021-01-15 21:00:00         1.6          238.3        0.9    42 5.481069
## 767 2021-01-15 22:00:00        -0.3          336.9        0.6    35 4.791733
## 768 2021-01-15 23:00:00        -0.5          317.7        0.8    39 5.795961
##           EC   SO42-     NO3-     Cl-     Na+    NH4+      K+    Mg2+    Ca2+
## 763 2.089880 3.66652 19.95149 1.93789 0.47189 8.24251 0.15435 0.09103 0.12778
## 764 1.609157 4.93662 18.84049 1.71229 0.32849 8.42511 0.16545 0.07553 0.11468
## 765 1.279726 5.04612 13.01529 1.15939 0.24319 6.20201 0.12595 0.07133 0.12818
## 766 1.285806 4.56892  9.08269 0.67019 0.11739 4.66061 0.10895 0.05663 0.14198
## 767 1.184745 4.14502  7.94549 0.73569 0.14109 3.97071 0.08135 0.06813 0.13818
## 768 1.371843 4.38362  5.56269 0.64329 0.12409 3.31881 0.11415 0.07043 0.10958

tail(aero, n = 3)

##                    Date Temperature Wind direction Wind speed PM2.5       OC
## 766 2021-01-15 21:00:00         1.6          238.3        0.9    42 5.481069
## 767 2021-01-15 22:00:00        -0.3          336.9        0.6    35 4.791733
## 768 2021-01-15 23:00:00        -0.5          317.7        0.8    39 5.795961
##           EC   SO42-    NO3-     Cl-     Na+    NH4+      K+    Mg2+    Ca2+
## 766 1.285806 4.56892 9.08269 0.67019 0.11739 4.66061 0.10895 0.05663 0.14198
## 767 1.184745 4.14502 7.94549 0.73569 0.14109 3.97071 0.08135 0.06813 0.13818
## 768 1.371843 4.38362 5.56269 0.64329 0.12409 3.31881 0.11415 0.07043 0.10958

• 열 추출하기

aero[1:10, 1] # aero 첫번째 열 중 1-10행까지만 벡터 형태로 출력

##  [1] "2020-12-15 00:00:00 UTC" "2020-12-15 01:00:00 UTC"
##  [3] "2020-12-15 02:00:00 UTC" "2020-12-15 03:00:00 UTC"
##  [5] "2020-12-15 04:00:00 UTC" "2020-12-15 05:00:00 UTC"
##  [7] "2020-12-15 06:00:00 UTC" "2020-12-15 07:00:00 UTC"
##  [9] "2020-12-15 08:00:00 UTC" "2020-12-15 09:00:00 UTC"

aero[1:10, 1, drop=FALSE] # 데이터 프레임 형태로 출력

##                   Date
## 1  2020-12-15 00:00:00
## 2  2020-12-15 01:00:00
## 3  2020-12-15 02:00:00
## 4  2020-12-15 03:00:00
## 5  2020-12-15 04:00:00
## 6  2020-12-15 05:00:00
## 7  2020-12-15 06:00:00
## 8  2020-12-15 07:00:00
## 9  2020-12-15 08:00:00
## 10 2020-12-15 09:00:00

aero[1:10 , c(1, 3, 4)] # 1, 3, 4번째 열만 출력

##                   Date Wind direction Wind speed
## 1  2020-12-15 00:00:00           21.5        0.2
## 2  2020-12-15 01:00:00          118.1        0.1
## 3  2020-12-15 02:00:00          125.1        0.3
## 4  2020-12-15 03:00:00           30.6        0.6
## 5  2020-12-15 04:00:00           95.3        0.5
## 6  2020-12-15 05:00:00          323.0        0.7
## 7  2020-12-15 06:00:00           77.4        0.2
## 8  2020-12-15 07:00:00          104.3        0.5
## 9  2020-12-15 08:00:00           13.9        0.0
## 10 2020-12-15 09:00:00          126.7        0.2

aero[1:10 , 2:5] # 2부터 다섯번째 열까지 출력

##    Temperature Wind direction Wind speed PM2.5
## 1        -12.3           21.5        0.2    10
## 2        -13.1          118.1        0.1    12
## 3        -13.6          125.1        0.3    13
## 4        -12.4           30.6        0.6    14
## 5        -13.7           95.3        0.5    13
## 6        -12.2          323.0        0.7    11
## 7        -12.8           77.4        0.2     9
## 8        -14.1          104.3        0.5     8
## 9        -14.9           13.9        0.0     8
## 10       -13.4          126.7        0.2     9

aero[1:10 , seq(from = 1, to = ncol(aero), by = 2)] #첫번째 열부터 마지막 열까지 홀수만 출력

##                   Date Wind direction PM2.5        EC    NO3-     Na+      K+
## 1  2020-12-15 00:00:00           21.5    10 0.6255070 1.39549 0.03549 0.03115
## 2  2020-12-15 01:00:00          118.1    12 0.6197119 1.49809 0.03399 0.02605
## 3  2020-12-15 02:00:00          125.1    13 0.5682204 1.54819 0.03159 0.03265
## 4  2020-12-15 03:00:00           30.6    14 0.5837869 1.53099 0.03029 0.03495
## 5  2020-12-15 04:00:00           95.3    13 0.5339275 1.34259 0.02999 0.03495
## 6  2020-12-15 05:00:00          323.0    11 0.4581114 1.17579 0.02659 0.03425
## 7  2020-12-15 06:00:00           77.4     9 0.3653754 0.94399 0.01959 0.03385
## 8  2020-12-15 07:00:00          104.3     8 0.4000881 0.79839 0.01839 0.02505
## 9  2020-12-15 08:00:00           13.9     8 0.4107523 0.83159 0.01929 0.02015
## 10 2020-12-15 09:00:00          126.7     9 0.5356505 0.89959 0.01999 0.01815
##       Ca2+
## 1  0.01698
## 2  0.01718
## 3  0.01488
## 4  0.01028
## 5  0.01728
## 6  0.01998
## 7  0.02048
## 8  0.02018
## 9  0.02758
## 10 0.02538

aero[1:10 , c("PM2.5", "OC", "EC")] # 다음과 같은 열 이름만 출력력

##    PM2.5        OC        EC
## 1     10 1.9441572 0.6255070
## 2     12 2.0489682 0.6197119
## 3     13 1.8051302 0.5682204
## 4     14 2.0836662 0.5837869
## 5     13 1.7662002 0.5339275
## 6     11 1.4007012 0.4581114
## 7      9 0.8442722 0.3653754
## 8      8 0.9623322 0.4000881
## 9      8 0.9659122 0.4107523
## 10     9 1.3125702 0.5356505

• grep: 변수명 중에서 특정 문자로 시작하거나 끝나거나 포함하고 있는 것을 추출할 때

argument	설명
pattern	특정한 문자로 시작하거나 끝나거나 포함하고 있는지를 문자형으로 표현한다. a로 시작하는 것은 “^a”, a로 끝나는 것은 “a$”, a를 포함하는 것은 “a”로 지정한다.x 문자형 벡터를 지정한다.
value	논리형으로 TRUE으로 pattern를 만족하는 벡터의 값을 반환하고,FALSE를 지정하면 pattern를 만족하는 벡터의 위치(index)를 반환한다.

grep("O", colnames(aero), value = TRUE) # 헤더이름 중 O가 들어가는 이름

## [1] "OC"    "SO42-" "NO3-"

grep("O", colnames(aero), value = FALSE) # 헤더이름 중 O가 들어가는 순서 

## [1] 6 8 9

aero[1:10, grep("O", colnames(aero), value = TRUE)] # 헤더이름 중 O가 들어가는 이름 중 10열까지 출력

##           OC   SO42-    NO3-
## 1  1.9441572 1.40232 1.39549
## 2  2.0489682 1.45922 1.49809
## 3  1.8051302 1.48532 1.54819
## 4  2.0836662 1.54722 1.53099
## 5  1.7662002 1.45522 1.34259
## 6  1.4007012 1.39652 1.17579
## 7  0.8442722 1.27212 0.94399
## 8  0.9623322 1.19452 0.79839
## 9  0.9659122 1.18652 0.83159
## 10 1.3125702 1.13932 0.89959

• substr: 문자 중에서 일부를 추출. 엑셀의 mid 함수와 동일

argument	설명
x	문자형 벡터를 지정한다.
start	추출할 문자의 첫 번째 위치를 정수형으로 지정.
stop	추출한 문자의 마지막 위치를 정수형으로 지정.

substr(colnames(aero), start = 1, stop = 2) # 변수명 중에서 처음 두 개의 문자를 추출

##  [1] "Da" "Te" "Wi" "Wi" "PM" "OC" "EC" "SO" "NO" "Cl" "Na" "NH" "K+" "Mg" "Ca"

substr(colnames(aero), 1, 2) # 변수명 중에서 처음 두 개의 문자를 추출

##  [1] "Da" "Te" "Wi" "Wi" "PM" "OC" "EC" "SO" "NO" "Cl" "Na" "NH" "K+" "Mg" "Ca"

aero[1:10 , substr(colnames(aero), start = 1, stop = 2) == "SO"] # 변수명 중에서 처음 두 개의 문자를 추출한 것이 “SO”와같은 열을 추출 (1-10행)

##  [1] 1.40232 1.45922 1.48532 1.54722 1.45522 1.39652 1.27212 1.19452 1.18652
## [10] 1.13932

원하는 케이스의 행 추출

엑셀의 필터링 기능

• aero 데이터에서 풍향 자료를 기준으로 북풍/남풍으로 구분해 볼 수 있습니다. 즉 풍향 (0-360도)이 270보다 크거나 90보다 작으면 북풍, 그렇지 않으면 남풍일 것입니다.
• 우선 모든 자료를 남풍 (South)로 할당해준 후 북풍 사례만을 골라서 북풍 (North)로 바꿔줍니다.

aero$main.wd = "North"
aero[is.na(aero$`Wind direction`) == F & aero$`Wind direction` >=  90 & aero$`Wind direction` <  270, "main.wd"] = "South"

south = aero[aero$main.wd == "South" , ] # 풍상풍하열 값 중 'u'만 추출하여 upwind에 할당
north = aero[aero$main.wd == "North" , ]

• factor를 이용하면 그룹에 따른 boxplot을 그릴 수 있다.

boxplot(aero$PM2.5 ~ aero$main.wd) #factor를 이용한 boxplot

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

aero[aero$PM2.5 >= 5, ] # Organic이 5 ug m-3 이상인 자료
aero[(aero$main.wd == 'North') & (aero$PM2.5 >= 35), ] # Organic이 5 ug m-3 이상이고 풍상지역인 자료
aero[(aero$main.wd == 'South') & (aero$PM2.5 >= 35), ] # Organic이 5 ug m-3 이상이거나 풍상지역인 자료
aero[(aero$main.wd == 'South') | (aero$PM2.5 >= 35), 2:3] # Organic이 5 ug m-3 이상이거나 풍상지역인 자료 중 2번째, 3번째 컬럼
aero[(aero$main.wd == 'South') | (aero$PM2.5 >= 355), c("OC", "EC")] # Organic이 5 ug m-3 이상이거나 풍상지역인 자료 중 "ORG", "SO4" 컬럼

시계열 자료 처리

• aero의 첫 번째 행 (Date)의 형식은 모두 YYYY-mm-dd HH:MI:SS의 형태임
• 엑셀로 부를땐, 시간 파일로 잘 읽히는 경우가 많음…

aero$Date[1]

• 위 형태는 R에서 시간 형태가 아닌 텍스트 형태로 인식했기 때문에, ‘lubridate’ 라이브러리의 ‘as.POSIXct()’ 함수를 이용하여 시간 형태로 변경시켜줘야 함.

argument	설명
x	변경할 변수
tz	사용자가 원하는 time zone. 미지정 시 “GMT”. 한국은 “Asia/Seoul”.

aero$Time = as.POSIXct(aero$Date, tz = "Asia/Seoul")
aero$Time[1]

## [1] "2020-12-15 09:00:00 KST"

• 뒤에 KST가 붙은 것을 알 수 있다. 이제 R은 해당 텍스트를 시간형식으로 인식한다.
• (엑셀로 불러들인 경우) 여전히 UTC이다. 이럴 경우 as.POSIXlt()를 사용하여 KST로 바꿔준다.
• 주의할 점은 UTC -> KST로 가면서 +9의 시간이 더해지기 때문에 -9 시간을 고려해주어야한다.

aero$Time = as.POSIXlt(aero$Date, tz = "Asia/Seoul")
aero$Time[1]

## [1] "2020-12-15 09:00:00 KST"

aero$Time = aero$Time - hours(9)
aero$Time[1]

## [1] "2020-12-15 KST"

• meteo 자료 부르기

meteo = read_excel(path =  "sample_data.xlsx", sheet = "Sheet2", col_names = TRUE)#, fileEncoding = "WINDOWS-1252")

• meteo 자료 형태 보기

str(meteo)

## tibble [768 × 16] (S3: tbl_df/tbl/data.frame)
##  $ Date   : POSIXct[1:768], format: "2020-12-15" "2020-12-15" ...
##  $ time   : POSIXct[1:768], format: "1899-12-31 00:00:00" "1899-12-31 01:00:00" ...
##  $ prec   : chr [1:768] "○" "○" "○" "○" ...
##  $ prec15 : chr [1:768] "0" "0" "0" "0" ...
##  $ prec60 : chr [1:768] "0" "0" "0" "0" ...
##  $ prec3h : chr [1:768] "0" "0" "0" "0" ...
##  $ prec6h : chr [1:768] "0" "0" "0" "0" ...
##  $ prec12h: chr [1:768] "0" "0" "0" "0" ...
##  $ precd  : chr [1:768] "0" "0" "0" "0" ...
##  $ temp   : chr [1:768] "-12.3" "-13.1" "-13.6" "-12.4" ...
##  $ wd     : num [1:768] 54.1 102.2 148.8 143.1 334 ...
##  $ wdd    : chr [1:768] "NE" "-" "SSE" "SE" ...
##  $ ws     : num [1:768] 0.7 0.1 0.5 0.6 0.3 0.8 0 0.1 0 0.3 ...
##  $ wd10   : num [1:768] 21.5 118.1 125.1 30.6 95.3 ...
##  $ wdd10  : chr [1:768] "NNE" "-" "SE" "NNE" ...
##  $ ws10   : num [1:768] 0.2 0.1 0.3 0.6 0.5 0.7 0.2 0.5 0 0.2 ...

• temp, wd, wd10 (10 m 풍향)등이 숫자가 아닌 문자로 표시되어 있다. 바꿔줘야함.

meteo$temp = as.numeric(meteo$temp)

## Warning: 강제형변환에 의해 생성된 NA 입니다

meteo$prec60 = as.numeric(meteo$prec60)

## Warning: 강제형변환에 의해 생성된 NA 입니다

• meteo 중 필요한 컬럼만 취하기!

meteo.slct = meteo[,c(1,2,5,10,11,13,14,16)]

• 시간자료로 병합하기

meteo.slct$Date[1:10]

##  [1] "2020-12-15 UTC" "2020-12-15 UTC" "2020-12-15 UTC" "2020-12-15 UTC"
##  [5] "2020-12-15 UTC" "2020-12-15 UTC" "2020-12-15 UTC" "2020-12-15 UTC"
##  [9] "2020-12-15 UTC" "2020-12-15 UTC"

meteo.slct$time[1:10]

##  [1] "1899-12-31 00:00:00 UTC" "1899-12-31 01:00:00 UTC"
##  [3] "1899-12-31 02:00:00 UTC" "1899-12-31 03:00:00 UTC"
##  [5] "1899-12-31 04:00:00 UTC" "1899-12-31 05:00:00 UTC"
##  [7] "1899-12-31 06:00:00 UTC" "1899-12-31 07:00:00 UTC"
##  [9] "1899-12-31 08:00:00 UTC" "1899-12-31 09:00:00 UTC"

meteo.slct$Time = meteo.slct$Date + hours(hour(meteo.slct$time))

meteo.slct$Time = as.POSIXct(meteo.slct$Time, tz = "Asia/Seoul")
meteo.slct$Time = meteo.slct$Time - hours(9)
meteo.slct$Time[1]

## [1] "2020-12-15 KST"

• aero와 meteo 자료 합치기 (merge 함수 사용; 엑셀의 vlookup): 기준이 되는 열을 지정후 두 자료를 병합. 엑셀과 달리 기준이 되는 열을 여러개 지정가능함.

argument	설명
x, y	합칠 데이터 셋
by	기준이 될 열 이름
all	겹치지 않는 행이 존재할 시 살려둘지, 삭제할지를 정함
all.x	x 데이터 셋의 겹치지 않는 행을 살려둠 (TRUE일 경우)
all.y	y 데이터 셋의 겹치지 않는 행을 살려둠 (TRUE일 경우)

aero.meteo = merge(aero, meteo.slct, by = "Time", all = T) #시간단위 평균의 gas와 aero 자료를 병합

• 일평균 계산 (엑셀의 부분합): aggregate 함수를 사용하여 일평균의 부분합을 계산함.

argument	설명
x	계산할 항목
by	그룹화할 항목
FUN	계산할 함수 (mean, median, max, min, …)
na.action	자료 중에 NA가 포함되었을 시 처리 방법

d.aero = aggregate(aero.meteo[,3:25], by = list(format(aero.meteo$Time, "%Y-%m-%d")), FUN = mean, na.rm = T) # 부분합 결과를 할당

## Warning in mean.default(X[[i]], ...): 인자가 수치형 또는 논리형이 아니므로 NA를
## 반환합니다
## Warning in mean.default(X[[i]], ...): 인자가 수치형 또는 논리형이 아니므로 NA를
## 반환합니다
## Warning in mean.default(X[[i]], ...): 인자가 수치형 또는 논리형이 아니므로 NA를
## 반환합니다
## Warning in mean.default(X[[i]], ...): 인자가 수치형 또는 논리형이 아니므로 NA를
## 반환합니다
## Warning in mean.default(X[[i]], ...): 인자가 수치형 또는 논리형이 아니므로 NA를
## 반환합니다
## Warning in mean.default(X[[i]], ...): 인자가 수치형 또는 논리형이 아니므로 NA를
## 반환합니다
## Warning in mean.default(X[[i]], ...): 인자가 수치형 또는 논리형이 아니므로 NA를
## 반환합니다
## Warning in mean.default(X[[i]], ...): 인자가 수치형 또는 논리형이 아니므로 NA를
## 반환합니다
## Warning in mean.default(X[[i]], ...): 인자가 수치형 또는 논리형이 아니므로 NA를
## 반환합니다
## Warning in mean.default(X[[i]], ...): 인자가 수치형 또는 논리형이 아니므로 NA를
## 반환합니다
## Warning in mean.default(X[[i]], ...): 인자가 수치형 또는 논리형이 아니므로 NA를
## 반환합니다
## Warning in mean.default(X[[i]], ...): 인자가 수치형 또는 논리형이 아니므로 NA를
## 반환합니다
## Warning in mean.default(X[[i]], ...): 인자가 수치형 또는 논리형이 아니므로 NA를
## 반환합니다
## Warning in mean.default(X[[i]], ...): 인자가 수치형 또는 논리형이 아니므로 NA를
## 반환합니다
## Warning in mean.default(X[[i]], ...): 인자가 수치형 또는 논리형이 아니므로 NA를
## 반환합니다
## Warning in mean.default(X[[i]], ...): 인자가 수치형 또는 논리형이 아니므로 NA를
## 반환합니다
## Warning in mean.default(X[[i]], ...): 인자가 수치형 또는 논리형이 아니므로 NA를
## 반환합니다
## Warning in mean.default(X[[i]], ...): 인자가 수치형 또는 논리형이 아니므로 NA를
## 반환합니다
## Warning in mean.default(X[[i]], ...): 인자가 수치형 또는 논리형이 아니므로 NA를
## 반환합니다
## Warning in mean.default(X[[i]], ...): 인자가 수치형 또는 논리형이 아니므로 NA를
## 반환합니다
## Warning in mean.default(X[[i]], ...): 인자가 수치형 또는 논리형이 아니므로 NA를
## 반환합니다
## Warning in mean.default(X[[i]], ...): 인자가 수치형 또는 논리형이 아니므로 NA를
## 반환합니다
## Warning in mean.default(X[[i]], ...): 인자가 수치형 또는 논리형이 아니므로 NA를
## 반환합니다
## Warning in mean.default(X[[i]], ...): 인자가 수치형 또는 논리형이 아니므로 NA를
## 반환합니다
## Warning in mean.default(X[[i]], ...): 인자가 수치형 또는 논리형이 아니므로 NA를
## 반환합니다
## Warning in mean.default(X[[i]], ...): 인자가 수치형 또는 논리형이 아니므로 NA를
## 반환합니다
## Warning in mean.default(X[[i]], ...): 인자가 수치형 또는 논리형이 아니므로 NA를
## 반환합니다
## Warning in mean.default(X[[i]], ...): 인자가 수치형 또는 논리형이 아니므로 NA를
## 반환합니다
## Warning in mean.default(X[[i]], ...): 인자가 수치형 또는 논리형이 아니므로 NA를
## 반환합니다
## Warning in mean.default(X[[i]], ...): 인자가 수치형 또는 논리형이 아니므로 NA를
## 반환합니다
## Warning in mean.default(X[[i]], ...): 인자가 수치형 또는 논리형이 아니므로 NA를
## 반환합니다
## Warning in mean.default(X[[i]], ...): 인자가 수치형 또는 논리형이 아니므로 NA를
## 반환합니다

d.aero$Group.1[1]

## [1] "2020-12-15"

d.aero$Group.1 = as.POSIXct(d.aero$Group.1, tz= "Asia/Seoul") # mi.aero$Group.1을 시간형식으로 변환 as.POSIXlt 함수 사용용
d.aero$Group.1[1]

## [1] "2020-12-15 KST"

내부 데이터 저장

• 텍스트 파일은 주로 txt 또는 csv로 저장함 (write.csv 함수 사용)

argument	설명
x	저장할 데이터 프레임
file	저장할 파일명
append	기존 파일 뒤에 병합하여 저장하는 옵션, 논리형 (T 또는 F)
sep	각 열의 구분자 (seperator)
row.names	행 이름 (주로 NA로 지정)
col.names	열 이름 (주로 지정하지 않음)
fileEncoding	파일 encoding 형식 (주로 지정하지 않음)

write.csv(file = "sample_data_dd.csv", d.aero, row.names = F) # 행번호를 없애기 위해 F를 씁니다.
write.csv(file = "sample_data.csv", aero.meteo, row.names = F) # 행번호를 없애기 위해 F를 씁니다.

실전편

기본 설정

setwd("~/Rclass/")

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

getwd()

## [1] "D:/OneDrive - hufs.ac.kr/Workspace_hufs/Rclass/지구환경과"

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

다음으로는 분석에 필요한 라이브러리를 로드해줍니다.

library("lubridate")
library("plyr")

강원 고성 측정자료를 불러오기 (CSV 파일) 불러오기

DAT = read.csv(file = “파일명”, header = TRUE)

cfc11 = read.csv(file = "CFC11_Hour_Background_GS.csv", header = TRUE)
cfc12 = read.csv(file = "CFC12_Hour_Background_GS.csv", header = TRUE)
cfc113 = read.csv(file = "CFC113_Hour_Background_GS.csv", header = TRUE)
co2 = read.csv(file = "CO2_Hour_Background_GS.csv", header = TRUE)
ch4 = read.csv(file = "CH4_Hour_Background_GS.csv", header = TRUE)
n2o = read.csv(file = "N2O_Hour_Background_GS.csv", header = TRUE)

불러온 자료를 확인해보기

head(cfc11)

##                  date cfc11
## 1 2012-01-01 00:00:00   240
## 2 2012-01-01 01:00:00   240
## 3 2012-01-01 02:00:00   241
## 4 2012-01-01 03:00:00   240
## 5 2012-01-01 04:00:00   241
## 6 2012-01-01 05:00:00   242

head(cfc12)

##                  date cfc12
## 1 2012-01-01 00:00:00   525
## 2 2012-01-01 01:00:00   520
## 3 2012-01-01 02:00:00   525
## 4 2012-01-01 03:00:00   524
## 5 2012-01-01 04:00:00   521
## 6 2012-01-01 05:00:00   522

head(cfc113)

##                  date cfc113
## 1 2012-01-01 00:00:00   82.4
## 2 2012-01-01 01:00:00   89.7
## 3 2012-01-01 02:00:00   88.7
## 4 2012-01-01 03:00:00   85.7
## 5 2012-01-01 04:00:00   83.9
## 6 2012-01-01 05:00:00   86.6

head(co2)

##                  date   co2
## 1 2012-01-16 00:00:00 402.8
## 2 2012-01-16 01:00:00 403.5
## 3 2012-01-16 02:00:00 403.3
## 4 2012-01-16 03:00:00 403.3
## 5 2012-01-16 04:00:00 403.3
## 6 2012-01-16 05:00:00 403.7

head(ch4)

##                  date    ch4
## 1 2012-01-16 01:00:00 1921.7
## 2 2012-01-16 02:00:00 1919.7
## 3 2012-01-16 03:00:00 1918.7
## 4 2012-01-16 04:00:00 1911.5
## 5 2012-01-16 05:00:00 1913.5
## 6 2012-01-16 06:00:00 1912.5

head(n2o)

##                  date   n2o
## 1 2012-01-01 01:00:00 315.5
## 2 2012-01-01 02:00:00 316.4
## 3 2012-01-01 03:00:00    NA
## 4 2012-01-01 04:00:00    NA
## 5 2012-01-01 05:00:00 318.1
## 6 2012-01-01 06:00:00 318.4

시간 자료를 인식시켜주리

cfc11$date = as.POSIXct(cfc11$date, tz = "Asia/Seoul")
cfc12$date = as.POSIXct(cfc12$date, tz = "Asia/Seoul")
cfc113$date = as.POSIXct(cfc113$date, tz = "Asia/Seoul")
co2$date = as.POSIXct(co2$date, tz = "Asia/Seoul")
ch4$date = as.POSIXct(ch4$date, tz = "Asia/Seoul")
n2o$date = as.POSIXct(n2o$date, tz = "Asia/Seoul")
n2o$date[1]

## [1] "2012-01-01 01:00:00 KST"

자료를 하나로 합치기

ghgs = merge(co2, ch4, by = "date", all = T)
ghgs = merge(ghgs, n2o, by = "date", all = T)
ghgs = merge(ghgs, cfc11, by = "date", all = T)
ghgs = merge(ghgs, cfc12, by = "date", all = T)
ghgs = merge(ghgs, cfc113, by = "date", all = T)
head(ghgs)

##                  date co2 ch4   n2o cfc11 cfc12 cfc113
## 1 2012-01-01 00:00:00  NA  NA    NA   240   525   82.4
## 2 2012-01-01 01:00:00  NA  NA 315.5   240   520   89.7
## 3 2012-01-01 02:00:00  NA  NA 316.4   241   525   88.7
## 4 2012-01-01 03:00:00  NA  NA    NA   240   524   85.7
## 5 2012-01-01 04:00:00  NA  NA    NA   241   521   83.9
## 6 2012-01-01 05:00:00  NA  NA 318.1   242   522   86.6

그래프 그리기

• 기본적으로 그래프는 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(ghgs$co2 ~ ghgs$date) # (1): 일반적인 점 그래프

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

plot(ghgs$co2 ~ ghgs$date, type = "l", ylim = c(360,450)) # (3): (2) + y축 범위 변경

plot(ghgs$co2 ~ ghgs$date, type = "l", ylim = c(360,450), col = "darkorange") # (4): (3) + 색깔 변경

plot(ghgs$co2 ~ ghgs$date, type = "l", ylim = c(360,450), col = "darkorange", xlab = "Time", ylab = "CO2") # (5): (4) + 이름 변경

plot(ghgs$co2 ~ ghgs$date, type = "l", ylim = c(360,450), col = "darkorange", xlab = "Time", ylab = expression("CO"[2]*" (ppm)")) # (6): (5) + 특수문자 및 아래첨자 
{plot(ghgs$co2 ~ ghgs$date, type = "l", ylim = c(360,450), col = "darkorange", xlab = "Time", ylab = expression("CO"[2]*" (ppm)")) 
points(ghgs$co2 ~ ghgs$date, col = "darkorange", pch = 22, cex = 0.1)
legend("topright", legend = expression("CO"[2]), col = "darkorange", pch = 22, cex = 0.7)} # (7): (6) + 심볼 추가 + 레전드 추가 

{plot(ghgs$cfc11 ~ ghgs$date, type = "l", ylim = c(0,300), col = "darkorange", xlab = "Time", ylab = "CFC (ppt)") 
points(ghgs$cfc11 ~ ghgs$date, col = "darkorange", pch = 22, cex = 0.7)

points(ghgs$cfc113 ~ ghgs$date, col = "blue", pch = 21, cex = 0.7) # (8): (7) + SO4 자료 추가
lines(ghgs$cfc113 ~ ghgs$date, col = "blue")
legend("topright", legend = c("CFC-11", "CFC-113"), col = c("darkorange", "blue"), pch = c(22, 21), cex = 0.7)} 

• y축라벨과 숫자의 간격을 줄이는 방법

{plot(ghgs$cfc11 ~ ghgs$date, type = "l", ylim = c(100,300), col = "darkorange", xlab = "Time", ylab = "")
  mtext(text = "CFC-11 (ppt)", 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_cfc11.png", width = 5, height = 3, units = "in", res = 200) # 저장할 스케치북 열기
  plot(ghgs$cfc11 ~ ghgs$date, type = "l", ylim = c(0,100), col = "darkorange", xlab = "Time", ylab = "")
  mtext(text = "CFC-11 (ppt)", line = 2, side = 2)
  dev.off() # 스케치북 닫기
} # (7): (6) 그래프 저장

## png 
##   2

산점도 및 추세선 추가

• lm (linear model): 선형 회귀직선
• abline: 선 추가 명령어

sl = lm(ghgs$ch4 ~ ghgs$co2) # y ~ x의 방정식 (y = bx + a)
(sl2 = summary(sl)) # sl 변수의 종합적 통계

## 
## Call:
## lm(formula = ghgs$ch4 ~ ghgs$co2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -163.642  -19.741   -0.869   17.401  199.879 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 1.183e+03  7.991e+00   148.0   <2e-16 ***
## ghgs$co2    1.873e+00  1.949e-02    96.1   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 33.86 on 26252 degrees of freedom
##   (결측으로 인하여 70178개의 관측치가 삭제되었습니다.)
## Multiple R-squared:  0.2602, Adjusted R-squared:  0.2602 
## F-statistic:  9235 on 1 and 26252 DF,  p-value: < 2.2e-16

{plot(ghgs$ch4 ~ ghgs$co2, type = "p", col = "darkorange", xlab = "", ylab = "") #ORG와 SO4 산점도. 축 라벨은 없음
mtext(text = expression("CH"[4]*" (ppb)"), line = 2, side = 2) # y축 라밸 추가
mtext(text = expression("CO"[2]*" (ppm)"), 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) # 추세식 표기
} 

Openair

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

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

install.packages("openair")

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

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

summary plot

• 시계열 자료의 전체적인 경향을 파악하기 좋은 그래프. 자료 및 결측치의 분포, 평균, 중간값 등에 대한 정보를 제공함함.
• 중요한 것은 시간은 date, 풍속은 ws, 풍향은 wd로 인식함

summaryPlot(ghgs) # 강원 고성 summary plot

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

windrose

• 풍향, 풍속에 대한 정보를 제공

windRose(ghgs) # 기본 바랑장미
windRose(ghgs, type = "month", layout = c(2, 1)) # 월별 바랑장미

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

pollution rose

• 오염장미

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

Polar plot

• 풍향/풍속에 따른 오염물질의 평균 농도를 보여줌

polarPlot(ghgs, pollutant = "co2") 
polarPlot(ghgs, pollutant = "co2", statistic = "cpf", percentile = 75) # Conditional Probability Function
polarPlot(ghgs, pollutant = "co2", statistic = "cpf", percentile = c(0,10)) # Conditional Probability Function
polarPlot(ghgs, pollutant = "co2", statistic = "cpf", percentile = c(50,60)) # Conditional Probability Function
polarPlot(ghgs, pollutant = "co2", type = "daylight") # ”year”, ”hour”, ”month”, ”season”, ”weekday”, ”weekend”, ”monthyear”, ”daylight”

Calendar plot

• 달력에 원하는 오염물질의 농도등을 표시하는 그림. 보다 직관적이기에 발표자료에 많이 사용됨.

calendarPlot(ghgs[year(ghgs$date) == 2020,], pollutant = "co2") # 고성 CO2 농도 기준, 표시

Trend 분석 (smoothTrend, TheilSen)

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

TheilSen(ghgs, pollutant = "co2", deseason = TRUE,  ylab = "CO2", main = "Theil-Sen Trend") # Theil-sen 기울기 계산산

## Taking bootstrap samples. Please wait.

TheilSen(ghgs, pollutant = "co2", deseason = TRUE,  ylab = "CO2", main = "Theil-Sen Trend", type = "season") # 계절절을 기준으로 분류

## Taking bootstrap samples. Please wait.

### timeVariation 분석

• 월별, 시간별, 요일별 분석을 수행 가능

# 기본 그림
timeVariation(ghgs, pollutant = "co2")

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

그 외 유용한 함수

scatterPlot(ghgs, x = "co2", y = "ch4", method = "density", col = "jet", type = "season")

## (loaded the KernSmooth namespace)

aero.meteo = rollingMean(ghgs, pollutant = "co2", hours = 8, new.name = "ghgs.8", data.thresh = 75) # 8시간 이동 평균

trendLevel(ghgs, pollutant = "co2", cols = "jet", y = "hour")