RMD와 rpubs에서 plotly 사용여부 연습
yoonani
2016년 8월 8일
R Markdown을 한동안 안 쓰니 많이 잊어먹었습니다. 다시 한번 R Markdown을 시작하려고 합니다. 연습삼아 내용은 OECD.Stat에서 온실가스 배출 자료를 받아서 저장한 파일로부터 우리나라의 온실가스 배출량을 얻어오고 같은 기간 서울의 년평균 자료를 기상자료개방포털의 기온분석으로부터 다운로드 받은 자료를 사용하였습니다.
자료 준비하기
- OECD로부터 자료를 다운로드 받은 파일을 읽어온다. (관련링크)
ghg <- read.csv("../data/oecd_GHG.csv", header=T)
str( ghg )## 'data.frame': 997 obs. of 17 variables:
## $ COU : Factor w/ 46 levels "AUS","AUT","BEL",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ Country : Factor w/ 46 levels "Australia","Austria",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ POL : Factor w/ 1 level "GHG": 1 1 1 1 1 1 1 1 1 1 ...
## $ Pollutant : Factor w/ 1 level "Greenhouse gases": 1 1 1 1 1 1 1 1 1 1 ...
## $ VAR : Factor w/ 1 level "TOTAL": 1 1 1 1 1 1 1 1 1 1 ...
## $ Variable : Factor w/ 1 level "Total emissions excluding LULUCF": 1 1 1 1 1 1 1 1 1 1 ...
## $ YEA : int 1990 1991 1992 1993 1994 1995 1996 1997 1998 1999 ...
## $ Year : int 1990 1991 1992 1993 1994 1995 1996 1997 1998 1999 ...
## $ Unit.Code : Factor w/ 1 level "T_CO2_EQVT": 1 1 1 1 1 1 1 1 1 1 ...
## $ Unit : Factor w/ 1 level "Tonnes of CO2 equivalent": 1 1 1 1 1 1 1 1 1 1 ...
## $ PowerCode.Code : int 3 3 3 3 3 3 3 3 3 3 ...
## $ PowerCode : Factor w/ 1 level "Thousands": 1 1 1 1 1 1 1 1 1 1 ...
## $ Reference.Period.Code: logi NA NA NA NA NA NA ...
## $ Reference.Period : logi NA NA NA NA NA NA ...
## $ Value : num 418623 418674 423080 423765 424093 ...
## $ Flag.Codes : Factor w/ 3 levels "","E; I","I": 1 1 1 1 1 1 1 1 1 1 ...
## $ Flags : Factor w/ 3 levels "","Estimated value; Incomplete data",..: 1 1 1 1 1 1 1 1 1 1 ...names( ghg )## [1] "COU" "Country"
## [3] "POL" "Pollutant"
## [5] "VAR" "Variable"
## [7] "YEA" "Year"
## [9] "Unit.Code" "Unit"
## [11] "PowerCode.Code" "PowerCode"
## [13] "Reference.Period.Code" "Reference.Period"
## [15] "Value" "Flag.Codes"
## [17] "Flags"- 분석에 필요한 변수들만으로 재구성하고 적절한 변수명(ghg.s)으로 저장한다.
ghg.s <- ghg[c("COU", "Year", "Value", "Flag.Codes")]
str(ghg.s)## 'data.frame': 997 obs. of 4 variables:
## $ COU : Factor w/ 46 levels "AUS","AUT","BEL",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ Year : int 1990 1991 1992 1993 1994 1995 1996 1997 1998 1999 ...
## $ Value : num 418623 418674 423080 423765 424093 ...
## $ Flag.Codes: Factor w/ 3 levels "","E; I","I": 1 1 1 1 1 1 1 1 1 1 ...- 우리나라의 자료를 불러와 새로운 데이터 프레임으로 만든다. OECD 자료에서 변수명 “COU”는 국가코드로 우리나라의 경우 “KOR”이고, 이를 이용하여 subset() 함수로 우리나라 자료만 선택해서 kor.ghg에 저장한다.
kor.ghg <- subset(ghg.s, COU=="KOR")
head(kor.ghg)## COU Year Value Flag.Codes
## 376 KOR 1990 292259
## 377 KOR 1991 315369
## 378 KOR 1992 342117
## 379 KOR 1993 377444
## 380 KOR 1994 402352
## 381 KOR 1995 433906- 서울의 평균 온도 자료를 가져온다. 기상자료개방포털로부터 다운로드 받은 csv 파일은 처음 몇 줄에 걸쳐 데이터에 대한 설명을 하고 있으므로 건너뛰고 변수명 또한 따로 생성하기 위해 skip=7을 이용해 파일의 7줄을 건너뛰고 8번째 줄부터 바로 데이터로 읽는다(header=F)
se <- read.csv("../data/se_1990_2013.csv", skip=7, header=F)- 변수명은 Year(년도), obs.p(관측지점 코드), ave.T(평균온도), min.T(최저온도), max.T(최고온도)가 되도록 한다.
names( se ) <- c("Year", "obs.p", "ave.T", "min.T", "max.T")
str(se)## 'data.frame': 24 obs. of 5 variables:
## $ Year : int 1990 1991 1992 1993 1994 1995 1996 1997 1998 1999 ...
## $ obs.p: int 108 108 108 108 108 108 108 108 108 108 ...
## $ ave.T: num 12.8 12.3 12.5 12 13.5 12.2 12.2 12.9 13.8 13.2 ...
## $ min.T: num -17.1 -15 -9.7 -11.2 -12.6 -11.1 -13.8 -13.7 -15.4 -12.3 ...
## $ max.T: num 35.5 34.7 32.8 31.9 38.4 33.7 35.4 36.1 32.8 35.4 ...- OECD 자료와 평균기온 자료를 하나의 데이터 프레임으로 만듭니다. 두 자료는 Year 변수를 갖고 있어 이 변수를 기준으로 merge() 함수를 이용하여 합칩니다.
kor <- merge(kor.ghg, se, by="Year")
str(kor)## 'data.frame': 24 obs. of 8 variables:
## $ Year : int 1990 1991 1992 1993 1994 1995 1996 1997 1998 1999 ...
## $ COU : Factor w/ 46 levels "AUS","AUT","BEL",..: 28 28 28 28 28 28 28 28 28 28 ...
## $ Value : num 292259 315369 342117 377444 402352 ...
## $ Flag.Codes: Factor w/ 3 levels "","E; I","I": 1 1 1 1 1 1 1 1 1 1 ...
## $ obs.p : int 108 108 108 108 108 108 108 108 108 108 ...
## $ ave.T : num 12.8 12.3 12.5 12 13.5 12.2 12.2 12.9 13.8 13.2 ...
## $ min.T : num -17.1 -15 -9.7 -11.2 -12.6 -11.1 -13.8 -13.7 -15.4 -12.3 ...
## $ max.T : num 35.5 34.7 32.8 31.9 38.4 33.7 35.4 36.1 32.8 35.4 ...head(kor)## Year COU Value Flag.Codes obs.p ave.T min.T max.T
## 1 1990 KOR 292259 108 12.8 -17.1 35.5
## 2 1991 KOR 315369 108 12.3 -15.0 34.7
## 3 1992 KOR 342117 108 12.5 -9.7 32.8
## 4 1993 KOR 377444 108 12.0 -11.2 31.9
## 5 1994 KOR 402352 108 13.5 -12.6 38.4
## 6 1995 KOR 433906 108 12.2 -11.1 33.7년도별 증감비율 계산
온실가스 배출량과 평균기온의 전년도 대비 증감비율을 계산해 봅시다. 이를 위해 벡터 변수내의 값의 증감을 계산하는 방법에 대해 먼저 알아봅시다.
다음의 코드에서 변수 x는 1, 3, 2, 7, 5, 10의 6개의 원소로 구성된 벡터로 변화량은 2, -1, 5, -2, 5 임을 알 수 있습니다. 이것을 R로 계산하기 위해 x.s는 처음의 다섯개의 원소를 가진 벡터, x.l은 첫번째 원소를 제외한 벡터로 x.l에서 x.s를 빼면 차이를 구할 수 있습니다.
( x <- c(1, 3, 2, 7, 5, 10) )## [1] 1 3 2 7 5 10( x.s <- x[-6] )## [1] 1 3 2 7 5( x.l <- x[-1] )## [1] 3 2 7 5 10( x.l - x.s )## [1] 2 -1 5 -2 5은근히 복잡합니다. R을 개발하는 분들이 가만히둘리 없습니다. 다음의 diff() 함수가 이 역할을 합니다.
( x.d <- diff(x) )## [1] 2 -1 5 -2 5여기서 한가지 주의할 점은 diff() 함수를 구하면 결과 벡터의 원소의 갯수가 기존 벡터보다 하나가 적다는 것입니다. 이는 2번째에서 1번째를 뺀 값, 3번째에서 2번째를 뺀 값, … , n번째에서 n-1번째를 뺀 값을 구하므로 길이가 1 줄어든 다는 점 잘 기억해주시기 바랍니다.
length(x)## [1] 6length( diff(x) )## [1] 5diff()를 이용하여 증감분을 구하였으니, 이제 전년도 대비 증감율을 구해 봅시다. 먼저 다음의 코드를 보겠습니다.
diff(kor$Value) / kor$Value[-length(kor$Value)]## [1] 0.079073698 0.084814931 0.103259996 0.065991246 0.078423868
## [6] 0.078565404 0.064365935 -0.140329921 0.088195732 0.070437890
## [11] 0.025699267 0.041377801 0.019895007 0.016575358 0.006813854
## [16] 0.009799197 0.028008362 0.023384510 0.005656268 0.099082618
## [21] 0.042228568 0.005369922 0.014956934우리나라의 온실가스 배출량은 kor$Value에 저장되어 있고 diff() 함수를 이용하여 증감분을 구했습니다. 그 다음 증감율을 구하기 위해 직전 년도 대비 값을 가져오는 코드로 kor$Value[-length(kor$Value)] 를 사용하였습니다. 여기서 생성된 값은 온실가스 배출량 자료중 마지막 자료는 제외하였습니다. 이렇게 마지막 제외한 것은 diff()를 통해 직전 값과의 차이를 가져오고 직전 값과의 차이는 직전 값 대비 증감율을 구하기 위함입니다. 이렇게 해야 두 벡터의 길이가 동일해 집니다.
온실가스 배출량의 증감율과 서울 평균기온의 변화의 증감율을 다음과 같이 구합니다.
ghg.diff <- diff(kor$Value) / kor$Value[-length(kor$Value)]
ave.T.diff <- diff(kor$ave.T) / kor$ave.T[-length(kor$ave.T)]그래프 작성
먼저 R의 기본 plot을 그림을 그려봤습니다. 여기서 y축 구간을 구하기 위해 온실가스 배출량의 증감율의 최소값과 평균온도 증감율의 최소값 중 더 작은 값으로 y축의 하한으로 하였고, 상한은 온실가스 배출량의 증감율의 최대값과 평균온도 증감율의 최대값 중 더 큰 값으로 하였습니다.
plot(1991:2013, ave.T.diff, type="l", col="blue", xlab="Year", ylab="increasing rate", ylim=c( min(min(ave.T.diff), min(ghg.diff)), max(max(ave.T.diff), max(ghg.diff))))
lines(1991:2013, ghg.diff, type="l", col="red")
abline(h=0, lty=2)
너무도 많이 사용하는 ggplot으로도 그려보았습니다. 이를 위해 먼저 Hadley-wickham이 제안한(?) tidy data로 만들기 위해 reshape 패키지의 melt를 이용하여 자료를 재구성하고 ggplot으로 그림을 그렸습니다.
library(ggplot2)
library(reshape2)
gt.kor <- data.frame(Year=1990:2012, ghg=ghg.diff, temp=ave.T.diff)
gt.m <- melt(gt.kor, id="Year")
p <- ggplot(gt.m, aes(Year, value)) + geom_line(aes(colour=variable)) + geom_abline(intercept=0, slope=0)
p
또한 plotly를 이용하여 동적인 그래프를 그려보았습니다.(오늘 테스트 글의 목적으로 ploty로 작성한 그래프가 rpubs.com 에서 잘 보이는지 확인해 보고자 합니다.)
library(plotly)##
## Attaching package: 'plotly'## The following object is masked from 'package:ggplot2':
##
## last_plot## The following object is masked from 'package:graphics':
##
## layoutgg <- ggplotly(p)
gg별 의미 없이 연습삼아 작성해 본 것입니다. 기후에 대해 이야기하기에 우리나라 데이터만 쓰는 것도 무리라고 생각하니 그래프의 의미를 고민하지는 말아주세요 :)
감사합니다.
자료 출처
- 우리나라 온실가스 배출량 : [OECD.stat], http://stats.oecd.org/
- 서울시 년간 평균온도 : [기상자료개방포털], https://data.kma.go.kr/stcs/grnd/grndTaList.do?pgmNo=70