이번 포스팅에서는 R에서 칼럼의 순서를 역순으로 재정렬하는 3가지 방법을 소개하겠습니다. 

 

(1) 위치 인덱싱 (position indexing)

(2) 수작업으로 칼럼 이름 인덱싱 (column name indexing manually)

(3) rev() 함수를 써서 역순 재정렬하여 칼럼 이름 인덱싱 (column name indexing using rev() function)

 

 

먼저 예제로 사용할 X1~X10의 10개의 칼럼으로 이루어진 간단한 DataFrame을 만들어보겠습니다. 

 

## sample DataFrame with 10 columns
df <- data.frame(matrix(1:30, nrow=3))

print(df)
# X1 X2 X3 X4 X5 X6 X7 X8 X9 X10
# 1  1  4  7 10 13 16 19 22 25  28
# 2  2  5  8 11 14 17 20 23 26  29
# 3  3  6  9 12 15 18 21 24 27  30

 

 

이제 X1~X5 까지는 순서를 역순으로 재정열하고, X6~X10은 원래의 순서를 그대로 유지하게끔 재정렬 해보겠습니다. 

(before: X1, X2, X3, X4, X5, X6, X7, X8, X9, X10

 after :    X5, X4, X3, X2, X1, X6, X7, X8, X9, X10) 

 

 

(1) 위치 인덱싱 (position indexing)

 

## X~X5를 역순으로 바꾸고, X6~X10은 그대로 두기
## (방법1) 위치(position) indexing
df2 <- df[c(5:1, # reverse
            6:10)]
print(df2)
# X5 X4 X3 X2 X1 X6 X7 X8 X9 X10
# 1 13 10  7  4  1 16 19 22 25  28
# 2 14 11  8  5  2 17 20 23 26  29
# 3 15 12  9  6  3 18 21 24 27  30

 

 

 

(2) 수작업으로 칼럼 이름 인덱싱 (column name indexing manually)

 

## (방법2) 칼럼 이름으로 indexing
df3 <- df[c(paste0(rep("X", 5), 5:1), # reverse, manually
            paste0(rep("X", 5), 6:10))]
print(df3)
# X5 X4 X3 X2 X1 X6 X7 X8 X9 X10
# 1 13 10  7  4  1 16 19 22 25  28
# 2 14 11  8  5  2 17 20 23 26  29
# 3 15 12  9  6  3 18 21 24 27  30

 

 

 

(3) rev() 함수를 써서 역순 재정렬하여 칼럼 이름 인덱싱 (column name indexing using rev() function)

 

## (방법3) 칼럼 이름으로 indexing, vec() 함수로 칼럼 이름 역순으로 재정렬
## column names
col_vec <- names(df)
print(col_vec)
# [1] "X1"  "X2"  "X3"  "X4"  "X5"  "X6"  "X7"  "X8"  "X9"  "X10"

## reverse column names
rev(col_vec)
# [1] "X10" "X9"  "X8"  "X7"  "X6"  "X5"  "X4"  "X3"  "X2"  "X1"


df4 <- df[c(rev(col_vec[1:5]), # reverse, using rev() function
            col_vec[6:10])]

print(df4)
# X5 X4 X3 X2 X1 X6 X7 X8 X9 X10
# 1 13 10  7  4  1 16 19 22 25  28
# 2 14 11  8  5  2 17 20 23 26  29
# 3 15 12  9  6  3 18 21 24 27  30

 

 

이번 포스팅이 많은 도움이 되었기를 바랍니다. 

행복한 데이터 과학자 되세요. 

반응형
Posted by Rfriend

댓글을 달아 주세요

이번 포스팅에서는 국가통계포털 사이트에서 받은 2020년, 2021년도 실업률과 취업자 수 통계 데이터를 가지고 R의 dplyr과 ggplot2 패키지를 사용해서 아래의 데이터 전처리 및 시각화하는 방법을 소개하겠습니다. 

 

1. 취업자 수 증가율(%) 변수 계산 (전년 동월 대비)

2. 실업률과 취업자 수 증가율 변수의 평균, 분산, 표준편차, 중앙값, 최대값, 최소값 계산

3. 실업률과 취업자 수 증가율 변수의 시계열 그래프 그리기

4. 실업률과 취업자 수 증가율 변수의  히스토그램 그리기 (히스토그램의 구간은 10개)

 

 

먼저, 국가통계포털 사이트에서 받은 2020년, 2021년도 실업률과 취업자 수 통계 데이터를 입력해서 DataFrame을 만들어보겠습니다. 데이터 자료 구조를 어떻게 해서 만드는지 유심히 봐주세요. 

 

## making a dataframe
df <- data.frame(
  month=c("01", "02", "03", "04", "05", "06", "07", "08", "09", "10", "11", "12"), 
  unemploy_rate_2020=c(4.1, 4.1, 4.2, 4.2, 4.5, 4.3, 4.0, 3.1, 3.6, 3.7, 3.4, 4.1), 
  unemploy_rate_2021=c(5.7, 4.9, 4.3, 4.0, 4.0, 3.8, 3.2, 2.6, 2.7, 2.8, 2.6, 3.5), 
  employed_num_2020=c(26800, 26838, 26609, 26562, 26930, 27055, 27106, 27085, 27012, 27088, 27241, 26526), 
  employed_num_2021=c(25818, 26365, 26923, 27214, 27550, 27637, 27648, 27603, 27683, 27741, 27795, 27298)
)

print(df)
#    month unemploy_rate_2020 unemploy_rate_2021 employed_num_2020 employed_num_2021
# 1     01                4.1                5.7             26800             25818
# 2     02                4.1                4.9             26838             26365
# 3     03                4.2                4.3             26609             26923
# 4     04                4.2                4.0             26562             27214
# 5     05                4.5                4.0             26930             27550
# 6     06                4.3                3.8             27055             27637
# 7     07                4.0                3.2             27106             27648
# 8     08                3.1                2.6             27085             27603
# 9     09                3.6                2.7             27012             27683
# 10    10                3.7                2.8             27088             27741
# 11    11                3.4                2.6             27241             27795
# 12    12                4.1                3.5             26526             27298

 

 

 

1. 취업자 수 증가율(%) 변수 계산 (전년 동월 대비)

 

dplyr 패키지로 새로운 변수를 생성하는 방법은 https://rfriend.tistory.com/235 를 참고하세요. 

dplyr 패키지의 chain operation, pipe operator %>% 사용 방법은 https://rfriend.tistory.com/236 를 참고하세요. 

 

## 1. 취업자 수 증가율(%) 변수 계산 (전년 동월 대비)
library(dplyr)
df2 <- df %>% 
  transform(
    employed_inc_rate = 100*(employed_num_2021 - employed_num_2020)/employed_num_2020) # percentage


print(df2)
#    month unemploy_rate_2020 unemploy_rate_2021 employed_num_2020 employed_num_2021 employed_inc_rate
# 1     01                4.1                5.7             26800             25818         -3.664179
# 2     02                4.1                4.9             26838             26365         -1.762426
# 3     03                4.2                4.3             26609             26923          1.180052
# 4     04                4.2                4.0             26562             27214          2.454634
# 5     05                4.5                4.0             26930             27550          2.302265
# 6     06                4.3                3.8             27055             27637          2.151174
# 7     07                4.0                3.2             27106             27648          1.999557
# 8     08                3.1                2.6             27085             27603          1.912498
# 9     09                3.6                2.7             27012             27683          2.484081
# 10    10                3.7                2.8             27088             27741          2.410662
# 11    11                3.4                2.6             27241             27795          2.033699
# 12    12                4.1                3.5             26526             27298          2.910352

 

 

 

2. 실업률과 취업자 수 증가율 변수의 평균, 분산, 표준편차, 중앙값, 최대값, 최소값 계산

 

dplyr 패키지로 데이터의 요약통계량을 계산하는 방법은 https://rfriend.tistory.com/235 를 참고하세요. 

여러개의 패키지별로 그룹별 요약통계량을 계산하는 방법은 https://rfriend.tistory.com/125 를 참고하세요. 

 

## 2. 실업률과 취업자 수 증가율 변수의 평균, 분산, 표준편차, 중앙값, 최대값, 최소값 계산
df2 %>% 
  summarise(
    unemploy_rate_2021_mean = mean(unemploy_rate_2021), 
    unemploy_rate_2021_var = var(unemploy_rate_2021), 
    unemploy_rate_2021_sd = sd(unemploy_rate_2021), 
    unemploy_rate_2021_median = median(unemploy_rate_2021), 
    unemploy_rate_2021_max = max(unemploy_rate_2021), 
    unemploy_rate_2021_min = min(unemploy_rate_2021)
  )

# unemploy_rate_2021_mean unemploy_rate_2021_var unemploy_rate_2021_sd 
#                   3.675              0.9547727             0.9771247
# 
# unemploy_rate_2021_median unemploy_rate_2021_max unemploy_rate_2021_min
#                      3.65                    5.7                    2.6


df2 %>% 
  summarise(
    employed_inc_rate_mean = mean(employed_inc_rate), 
    employed_inc_rate_var = var(employed_inc_rate), 
    employed_inc_rate_sd = sd(employed_inc_rate), 
    employed_inc_rate_median = median(employed_inc_rate), 
    employed_inc_rate_max = max(employed_inc_rate), 
    employed_inc_rate_min = min(employed_inc_rate)
  )

# employed_inc_rate_mean employed_inc_rate_var employed_inc_rate_sd 
#               1.367697              3.970439             1.992596
# 
# employed_inc_rate_median employed_inc_rate_max employed_inc_rate_min
#                 2.092436              2.910352             -3.664179

 

 

3. 실업률과 취업자 수 증가율 변수의 시계열 그래프 그리기

 

ggplot2 로 시계열 그래프 그리기는 https://rfriend.tistory.com/73 를 참고하세요. 

 

## 3. 실업률과 취업자 수 증가율 변수의 시계열 그래프 그리기 
library(ggplot2)

ggplot(df2, aes(x=month, y=unemploy_rate_2021, group=1)) +
  geom_line() +
  ylim(0, max(df2$unemploy_rate_2021)) +
  ggtitle("Time Series Plot of Unemployment Rate, Year 2021")

Time Series Plot of Unemployment Rate

 

 

ggplot(df2, aes(x=month, y=employed_inc_rate, group=1)) +
  geom_line() +
  ylim(min(df2$employed_inc_rate), max(df2$employed_inc_rate)) +
  ggtitle("Time Series Plot of Employment Increase Rate, Year 2021")

Time Series Plot of Employment Increase Rate

 

 

 

4. 실업률과 취업자 수 증가율 변수의  히스토그램 그리기 (히스토그램의 구간은 10개)

 

ggplot2 패키지로 히스토그램 그리기는 https://rfriend.tistory.com/67 를 참고하세요. 

 

ggplot(df2, aes(x=employed_inc_rate)) + 
  geom_histogram(bins=10) + 
  ggtitle("Histogram of Unemployment Rate, Year 2021")

Histogram of Unemployment Rate

 

 

ggplot(df2, aes(x=employed_inc_rate)) + 
  geom_histogram(bins=10) + 
  ggtitle("Histogram of Employment Incease Rate, Year 2021")

Histogram of Employment Increase Rate

 

 

이번 포스팅이 많은 도움이 되었기를 바랍니다. 

행복한 데이터 과학자 되세요!  :-)

 

반응형
Posted by Rfriend

댓글을 달아 주세요

이번 포스팅에서는 R data.frame에서 여러개의 칼럼 이름을 '변경 전 칼럼 이름 : 변경 후 칼럼 이름'의 매핑 테이블 (old_column_name : new_column_name mapping table) 을 이용해서 한꺼번에 변경하는 방법을 소개하겠습니다. data.frame에 칼럼 개수가 엄청 많고, 특정 칼럼에 대해서 선별적으로 칼럼 이름을 변경하고 싶을 때 전:후 칼럼 이름 매핑 테이블을 사용하는 이번 포스팅의 방법을 사용하면 편리합니다. 

 

renaming column names using mapping table in R data.frame

 

 

(1) 모든 칼럼을 순서대로 칼럼 이름을 변경하고 싶은 경우

 

참고로, R 에서 names(), rename() 등의 함수를 이용해서 칼럼 이름을 변경하는 방법은 https://rfriend.tistory.com/41 를 참고하세요. 

 

 

먼저, "X1" ~ "X10" 까지의 10개 칼럼을 가지는 예제 data.frame 을 만들어보겠습니다. 

 

## -- creating a sample data.frame with 10 columns
df <- data.frame(matrix(1:30, nrow=3))

print(df)
# X1 X2 X3 X4 X5 X6 X7 X8 X9 X10
# 1  1  4  7 10 13 16 19 22 25  28
# 2  2  5  8 11 14 17 20 23 26  29
# 3  3  6  9 12 15 18 21 24 27  30

 

 

다음으로, '변경 전 칼럼 이름 : 변경 후 칼럼 이름' 매핑 테이블을 만들어보겠습니다. 아래 예제에서는 변경 전 칼럼 이름 "X1"~"X10" 을 --> 변경 후 칼럼 이름 "var1"~"var10" 의 매핑 테이블 data.frame을 만들었습니다. (특정 칼럼만 선별적으로 변경하고 싶으면 해당 칼럼의 "변경 전 : 변경 후 매핑 테이블"을 만들면 됩니다.)

 

## -- creating a key(old column name):value(new column name) mapping table
old_col_nm <- names(df)
print(old_col_nm)
# [1] "X1"  "X2"  "X3"  "X4"  "X5"  "X6"  "X7"  "X8"  "X9"  "X10"

col_cnt <- ncol(df) # 10
new_col_nm <- paste0(c(rep("var", col_cnt)), 1:col_cnt)
print(new_col_nm)
# [1] "var1"  "var2"  "var3"  "var4"  "var5"  "var6"  "var7"  "var8"  "var9"  "var10"

df_col_dict <- data.frame("old_col_nm" = old_col_nm, "new_col_nm" = new_col_nm)
print(df_col_dict)
# old_col_nm new_col_nm
# 1          X1       var1
# 2          X2       var2
# 3          X3       var3
# 4          X4       var4
# 5          X5       var5
# 6          X6       var6
# 7          X7       var7
# 8          X8       var8
# 9          X9       var9
# 10        X10      var10

 

 

 

마지막으로, dplyr 패키지의 rename_at() 함수를 사용해서 "변경 전 칼럼 이름(old_col_nm)"을 "변경 후 칼럼 이름(new_col_nm)" 으로 변경해 보겠습니다. 

 

## -- changing data.frame's column names using key(old_col):value(new_col) mapping table
library(dplyr)
df_new <- df %>% 
  rename_at(vars(as.character(df_col_dict$old_col_nm)), 
            ~ as.character(df_col_dict$new_col_nm))

print(df_new)
# var1 var2 var3 var4 var5 var6 var7 var8 var9 var10
# 1    1    4    7   10   13   16   19   22   25    28
# 2    2    5    8   11   14   17   20   23   26    29
# 3    3    6    9   12   15   18   21   24   27    30

 

 

 

(2) 특정 칼럼만 선별적으로 이름을 바꾸고 싶은 경우

 

아래의 'col_dict' 테이블을 칼럼 이름을 변경하고자 하는 특정 칼럼의 old_col_nm : new_col_nm 으로 만들어서 적용하면 됩니다.

가령, 기존의 c1~c5'까지의 칼럼들 중에서 'c2', 'c4' 의 2개 칼럼만 선별적으로 변경하고 싶으면 아래처럼 'col_dict' 테이블을 만들어서 적용하면 돼요.

 

old_col_nm = c("c2", "c4")
new_col_nm = c("v2", "v4")

col_dict <- data.frame("old" = old_col_nm, "new" = new_col_nm)
print(col_dict)
# old new
# 2 c2 v2
# 4 c4 v4


library(dplyr)
c_df_new <- c_df %>%
rename_at(vars(as.character(col_dict$old)), ~ as.character(col_dict$new))

print(c_df_new)
# c1 v2 c3 v4 c5
# 1 1 4 7 10 13
# 2 2 5 8 11 14
# 3 3 6 9 12 15

 

이번 포스팅이 많은 도움이 되었기를 바랍니다. 

행복한 데이터 과학자 되세요!  :-)

 

반응형
Posted by Rfriend

댓글을 달아 주세요

이번 포스팅에서는 R의 data.table 패키지의 dcast() 함수를 사용해서 문자열(string)을 대상으로 데이터를 재구조화할 때 집계 함수 (aggregation function) 로서 

 

  (1) 문자열 원소의 개수 (length)

  (2) 문자열을 콤마로 구분해서 붙여쓰기

  (3) 첫번째 문자열만 가져오기

 

하는 방법을 소개하겠습니다. 

 

 

 

 

먼저 간단한 예제 data.table을 만들어보겠습니다. 

 

##---------------------------------
## R data.table dcast() for string
##---------------------------------

#install.packages("data.table")
library(data.table)

x1 <- c('g1', 'g1', 'g1', 'g1', 'g2', 'g2', 'g2', 'g2')
x2 <- c('v1', 'v2', 'v3', 'v3', 'v1', 'v2', 'v2', 'v3')
x3 <- c('a', 'b', 'c', 'd', 'e', 'f', 'g', 'h')

dt <- data.table(x1, x2, x3)
print(dt)
#    x1  x2  x3
# 1: g1  v1   a
# 2: g1  v2   b
# 3: g1  v3   c
# 4: g1  v3   d
# 5: g2  v1   e
# 6: g2  v2   f
# 7: g2  v2   g
# 8: g2  v3   h

 

 

R 의 data.table 패키지의 dcast() 함수를 사용해 데이터를 재구조화하면서 문자열을 대상으로 집계(value.var)를 할 때 집계 함수 (aggregation function) 을 명시적으로 적어주지 않으면 아래와 같은 경고 메시지가 발생합니다. 

 

Warning message: Aggregate function missing, defaulting to 'length'

 

이것은 문자열을 대상으로는 합계(sum), 최소값(min), 최대값(max), 평균(mean) 등의 숫자형을 대상으로 하는 요약통계량을 사용할 수 없기 때문입니다. 

 

##-- warning message
##: Aggregate function missing, defaulting to 'length'
dcast(dt, x1 ~ x2, 
      value.var = "x3")

# Aggregate function missing, defaulting to 'length'
#    x1  v1  v2  v3
# 1: g1   1   1   2
# 2: g2   1   2   1

 

 

따라서 dcast() 함수로 데이터를 재고조화시 문자열을 대상으로 집계를 한다면

  (1) 문자열 원소의 개수 (length)

  (2) 문자열을 콤마로 구분해서 붙여쓰기

  (3) 첫번째 문자열만 가져오기

와 같이 문자열에 맞는 집계함수를 지정해주어야 합니다. 

 

 

 

(1) dcast() 함수로 데이터셋 재구조화 시 문자열을 원소의 개수 (length) 로 집계

 

문자열 대상 집계일 때는 default 설정이 원소의 개수 (length) 이므로 위와 결과는 동일합니다만, 이번에는 경고 메시지가 안떴습니다. 

 

##-- (1) counting the number of values as an aggregation function for string values
dcast(dt, x1 ~ x2, 
      fun.aggregate = length, 
      value.var = "x3")

#    x1  1  2  3
# 1: g1  1  1  2
# 2: g2  1  2  1

 

 

 

(2) dcast() 함수로 데이터셋 재구조화 시 문자열을 콤마로 구분해서 붙여쓰기

 

dcast() 로 재구조화 시 하나의 셀 안에 여러개의 원소가 존재하게 될 경우, 이들 문자열 원소들을 콤마로 구분해서 옆으로 나란히 붙여서 집계하는 사용자 정의 함수를 fun.aggregate 매개변수란에 써주었습니다. 

 

##-- (2) concatenation as an aggregation function for string values
dcast(dt, x1 ~ x2, 
      fun.aggregate = function(x) if (length(x)==1L) x else paste(x, collapse=","), 
      value.var = "x3")

#    x1  1     2     3
# 1: g1  a     b   c,d
# 2: g2  e   f,g     h

 

 

 

(3) dcast() 함수로 데이터셋 재구조화 시 첫번째 문자열만 가져오기

 

dcast() 로 재구조화 시 하나의 셀 안에 여러개의 원소가 존재하게 될 경우, 이들 복수개의 원소들 중에서 첫번째 원소만 가져오는 사용자정의함수를 fun.aggregate 매개변수란에 작성해주었습니다. 

 

##-- (3) keeping the first value as an aggregation function for string values
dcast(dt, x1 ~ x2, 
      fun.aggregate = function(x) if (length(x)==1L) x else x[1], 
      value.var = "x3")
      
#    x1  1  2  3
# 1: g1  a  b  c
# 2: g2  e  f  h

 

 

이번 포스팅이 많은 도움이 되었기를 바랍니다. 

행복한 데이터 과학자 되세요!  :-)

 

반응형
Posted by Rfriend

댓글을 달아 주세요

주식을 하는 분들은 아마도 대표적인 시계열 데이터인 주가의 이동평균, 누적평균 그래프에 이미 익숙할 것입니다. 

 

이번 포스팅에서는 R의 zoo 패키지의 rollapply() 라는 window function

 

(1) Rolling Windows 를 사용해서 시계열 데이터의 이동 평균 구하기

     (average of time series using rolling windows)

(2) Expanding Windows 를 사용해서 시계열 데이터의 누적 평균 구하기

     (average of time series using expanding windows)

 

방법을 소개하겠습니다. 

 

 

[ 이동 평균 (average using Rolling Windows) vs. 누적 평균 (average using Expanding Windows) ]

moving average (rolling windows) vs. cumulative average (expanding windows) using R zoo rollapply() function

 

시계열 데이터를 전처리하고 분석할 때 Window Function 을 자주 사용하는데요, 

 - Rolling Windows : 특정 window width (예: 10분, 1시간, 1일 등) 를 유지한채 측정 단위시간별로 이동하면서 분석 

 - Expanding Windows : 처음 시작 시점은 고정한 채, 시간이 흐름에 따라 신규로 포함되는 데이터까지 누적해서 분석

하는 차이가 있습니다. 바로 위에 Rolling Windows 와 Expanding Windows 를 도식화 해놓은 자료를 보면 금방 이해가 될거예요. 

 

만약 시계열 데이터에 추세(trend) 나 계절성 (seasonality) 이 있다면 Rolling Windows 가 적당하며, 시계열 데이터에 추세나 계절성이 없이 안정적(stable) 이다면 Expanding Windows 를 사용해서 더 많은 데이터를 이용해서 요약 통계량을 계산하는게 유리할 수 있겠습니다. 

 

시계열 예측 모델링할 때는 Rolling Windows 를 사용해서 모델 성능을 검증합니다. 

 

 

R 의 zoo 패키지의 rollapply() 함수를 사용할 것이므로, zoo 패키지를 먼저 설치하고 임포팅합니다. 

그리고 예제로 사용할 간단한 시계열 데이터를 만들어보겠습니다. 추세와 노이즈가 있는 시계열 데이터 입니다. 

 

## ------------
## Wimdow functions in Time Series
## (1) Rolling window
## (2) Expanding window
## R zoo's rollapply(): https://www.rdocumentation.org/packages/zoo/versions/1.8-9/topics/rollapply
## ------------

install.packages("zoo")
library(zoo)

## generating a time series with trend and noise
set.seed(1) # for reproducibility
x <- rnorm(n=100, mean=0, sd=10) + 1:100

plot(x, type='l', 
     main="time series plot with trend and noise")

time series with trend and noise

 

 

 

(1) Rolling Windows 를 사용해서 시계열 데이터의 이동 평균 구하기

     (average of time series using rolling windows)

 

zoo 패키지의 rollapply() 함수에서

- width 매개변수는 'window width' 를 설정할 때 사용합니다.

- FUN 매개변수에는 원하는 함수를 지정해줄 수 있으므로 매우 강력하고 유연하게 사용할 수 있습니다. 아래 예에서는 평균(mean)과 최대값(max) 을 계산하는 함수를 사용해보았습니다.

- align 은 데이터의 기준을 정렬할 때 왼쪽("left"), 중앙("centered", default 설정), 오른쪽("right") 중에서 지정할 수 있습니다. 이때 align="left"로 설정해주면 자칫 잘못하면 미래의 데이터를 가져다가 요약 통계량을 만드는 실수 (lookahead) 를 할 수도 있으므로, 만약 예측 모델링이 목적이라면 lookahead 를 하는건 아닌지 유의해야 합니다. 

- partial=TRUE 로 설정하면 양쪽 끝부분에 window width 의 개수에 데이터 포인트 개수가 모자라더라도 있는 데이터만 가지고 부분적으로라도 함수의 통계량을 계산해줍니다. 

 

## (1) Rolling Windows

## (1-1) moving average
f_avg_rolling_win <- rollapply(
  data=zoo(x), 
  width=10, # window width
  FUN=function(w) mean(w), 
  # 'align' specifies whether the index of the result should be left-aligned 
  # or right-aligned or centered (default) 
  # compared to the rolling window of observations. 
  align="right", 
  # If 'partial=TRUE', then the subset of indexes 
  # that are in range are passed to FUN.
  partial=TRUE)

## (1-2) moving max
f_max_rolling_win <- rollapply(
  zoo(x), 
  10, 
  function(w) max(w), 
  align="right", 
  partial=TRUE)

plot(x, col="gray", lwd=1, type="l", main="Average and Max using Rolling Window")
lines(f_avg_rolling_win, col="blue", lwd=2, lty="dotted")
lines(f_max_rolling_win, col="red", lwd=2, lty="dashed")
legend("topleft", 
       c("Average with Rolling Windows", "Max with Rolling Windows"), 
       col = c("blue", "red"), 
       lty = c("dotted", "dashed"))

moving average and max using the rolling windows

 

 

 

 

(2) Expanding Windows 를 사용해서 시계열 데이터의 누적 평균 구하기

     (average of time series using expanding windows)

 

R 에서 zoo 패키지의 rollapply() 함수로 Expanding Windwos 를 사용하려면 width = seq_along(x) 를 지정해주면 누적으로 함수를 계산해줍니다. 

 

아래 예에서는 누적으로 평균과 최대값을 계산해서 시각화 한건데요, 우상향 하는 추세가 있는 시계열이다보니 누적으로 평균을 구하면 시계열 초반의 낮은 값들까지 모두 포함이 되어서 누적평균 값이 최근 값들을 제대로 따라가지 못하고 있습니다.

반면, 누적으로 최대값을 계산한 값은 중간에 소폭 값이 줄어들더라도 계산 시점까지 누적으로 최대값을 계산하므로, 항상 우상향하는 누적 최대값을 보여주고 있습니다.

(위의 (1)번의 이동평균, 이동최대값과 (2) 누적평균, 누적최대값을 비교해서 보세요.)

 

# (2) Expanding Windows

## (2-1) cumulative average
f_avg_expanding_win <- rollapply(
  data=zoo(x), 
  width=seq_along(x), # expanding windows
  FUN=function(w) mean(w), # average
  align="right", 
  partial=TRUE)

## (2-2) cumulative max
f_max_expanding_win <- rollapply(
  zoo(x), 
  seq_along(x), # expanding windows
  function(w) max(w), # max
  align="right", 
  partial=TRUE)

## plotting
plot(x, col="gray", lwd=1, type="l", main="Average and Max using Expanding Window")
lines(f_avg_expanding_win, col="blue", lwd=2, lty="dotted")
lines(f_max_expanding_win, col="red", lwd=2, lty="dashed")
legend("topleft", 
       c("Average with Expanding Windows", "Max with Expanding Windows"), 
       col = c("blue", "red"), 
       lty = c("dotted", "dashed"))

 

 

[ Reference ]

- R zoo's rollapply(): https://www.rdocumentation.org/packages/zoo/versions/1.8-9/topics/rollapply

 

 

이번 포스팅이 많은 도움이 되었기를 바랍니다. 

행복한 데이터 과학자 되세요!  :-)

 

반응형
Posted by Rfriend

댓글을 달아 주세요

지난번 포스팅에서는 지리공간 레스터 객체 데이터로 부터 일부를 가져오기 (Raster subsetting) 하는 방법(rfriend.tistory.com/629)을 소개하였습니다.

 

이번 포스팅에서는 지리공간 레스터 객체 데이터로 부터 요약 통계량(summary statistics), 기술 통계량(descriptive statistics)을 계산하고, 시각화(visualization) 하는 방법을 알아보겠습니다.

 

(1) 레스터 객체 요약 통계량 구하기 (summarizing raster objects)

(2) 레스터 객체 시각화 하기 (visualizing raster objects)

 

 

 

 

(1) 레스터 객체 요약 통계량 구하기 (summarizing raster objects)

 

먼저, raster 패키지를 불러오고, 예제로 사용하기 위해 6행 6열, 총 36개의 픽셀로 이루어진 간단한 레스터 객체를 만들어보겠습니다.

 

레스터 객체를 출력해보면 아래의 예처럼 차원(dimensions), 해상도(resolution), 그리고 X와 Y축의 최소값과 최대값의 범위(extent)에 대한 요약 통계량을 볼 수 있습니다.

 

## ================================
## GeoSpatial data analysis using R
## : Summarizing raster objects
## ================================

library(raster)

## sample raster object 
elev <- raster(nrows = 6, ncols = 6, 
               xmn = -1.5, xmx = 1.5, ymn = -1.5, ymx = 1.5, 
               vals = 1:36)

elev 
# class      : RasterLayer 
# dimensions : 6, 6, 36  (nrow, ncol, ncell)
# resolution : 0.5, 0.5  (x, y)
# extent     : -1.5, 1.5, -1.5, 1.5  (xmin, xmax, ymin, ymax)
# crs        : +proj=longlat +datum=WGS84 +no_defs 
# source     : memory
# names      : layer 
# values     : 1, 36  (min, max)
]

 

 

레스터 객체에 대해 summary(raster_obj) 함수를 사용하면 레스터 객체의 픽셀 속성값에 대한 최소값(Minimum value), 1사분위수(1st Quantile), 중위값(Median), 3사분위수(3rd Quantile), 최대값(Maximum value), 결측값 개수(NA's) 에 대한 기술통계량(descriptive statistics)을 확인할 수 있습니다.

 

## summary() function for raster objects
summary(elev)

# layer
# Min.     1.00
# 1st Qu.  9.75
# Median  18.50
# 3rd Qu. 27.25
# Max.    36.00
# NA's     0.00

 

 

레스터 객체의 픽셀 속성값에 대해 cellStats(raster_obj, summary_function) 함수를 사용하여 평균(mean), 분산(variance), 표준편차(standard deviation), 사분위수(quantile), 합계(summation) 를 구할 수 있습니다.

 

## -- cessStats() function for raster objects

## mean
cellStats(elev, mean)
# [1] 18.5


## variance
cellStats(elev, var)
# [1] 111


## standard deviation
cellStats(elev, sd)
# [1] 10.53565


## quantile
cellStats(elev, quantile)
# 0%   25%   50%   75%  100% 
# 1.00  9.75 18.50 27.25 36.00


## sum
cellStats(elev, sum)
# [1] 666

 

 

RasterBlack Calss, RasterStack Class와 같이 여러개의 층을 가지는 레스터 클래스 객체 (multi-layered raster classes objects) 에 대해서 앞서 소개한 요약 통계량을 구하는 함수를 적용하면 각 층별로 따로 따로 요약 통계량이 계산됩니다.

 

 

 

(2) 레스터 객체 시각화 하기 (visualizing raster objects)

 

## -- visualization of raster objects

## plot()
plot(elev, main = "raster objects")

 

 

레스터 객체의 픽셀 속성값에 대해 히스토그램(histogram for raster objects)을 그리려면 hist() 함수를 사용하면 됩니다. 

 

## histogram
hist(elev, main = "histogram for raster objects")

 

레스터 객체의 픽셀 속성값에 대해서 raster 패키지의 raster::density() 함수를 사용하여 추정 밀도 곡선(smoothed density estimates curve) 을 그릴 수 있습니다.

 

## raster::density() : density plot (smoothed density estimates)
density(elev, main = "density plot for raster objects")

 

 

레스터 객체의 픽셀 속성값에 대해 boxplot() 함수를 사용하면 상자그림(box plot for raster objects)을 그릴 수 있습니다.

 

## boxplot
boxplot(elev, main = "box plot for raster objects")

 

 

만약 레스터 객체에 대해 시각화를 하려고 하는데 안된다면 values() 또는 getValues() 함수를 사용해서 레스터 객체로 부터 픽셀의 속성값을 반환받은 결과에 대해서 시각화를 하면 됩니다.

 

 

 

[Reference]

[1] "Geocomputation with R - Attritube data operations"
    : geocompr.robinlovelace.net/attr.html

 

Geocomputation with R

Geocomputation with R is for people who want to analyze, visualize and model geographic data with open source software. It is based on R, a statistical programming language that has powerful data processing, visualization, and geospatial capabilities. The

geocompr.robinlovelace.net

 

이번 포스팅이 많은 도움이 되었기를 바랍니다.

행복한 데이터 과학자 되세요!  :-)

 

반응형
Posted by Rfriend

댓글을 달아 주세요

지난번 포스팅에서는 지리공간 레스터 객체 데이터를 시작하면서, 문자형을 값으로 가지는 범주형 변수를 요인형으로 변환하여 --> 요인형 변수를 속성정보로 하여 레스터 객체를 만드는 방법을 소개하였습니다. (rfriend.tistory.com/628

 

이번 포스팅에서는 레스터 객체 데이터의 일부를 가져오는 방법(Raster subsetting)을 소개하겠습니다.

 

(1) 레스터 객체의 각 픽셀(pixcels, cells) 중 일부 행(rows), 열(column) 가져오기

    (subsetting rows, columns from laster objects) 

(2) 다층 레스터 객체(multi-layered laster objects)로 부터 일부 층(layers) 가져오기

    (subsetting layers from multi-layered laster objects, a raster brick or stack)

(3) 레스터 객체의 일부 값을 수정하기 (modifying values in raster objects)

 

 

 

 

레스터 객체의 일부분 가져오기(Raster subsetting)는 Base R 의 '[' 연산자와 raster 패키지의 raster::subset() 함수를 사용해서 할 수 있습니다. 먼저 Base R의 '[' 연산자를 사용한 일부분 가져오기를 살펴보겠습니다. 

 

 

(1) 레스터 객체의 각 픽셀(pixcels, cells) 중 일부 행(rows), 열(column) 가져오기

    (subsetting rows, columns from laster objects)

 

예제로 사용하기 위해, 6개의 행과 6개의 열의 총 36개 픽셀로 구성되어 있고 1~36의 정수를 속성 값으로 가지는 레스터 객체를 만들어보겠습니다. 

 

## ========================================================
## R Raster subsetting
## - reference: https://geocompr.robinlovelace.net/attr.html
## ========================================================

library(raster)

## sample raster object
elev <- raster(nrows = 6, ncols = 6, 
              xmn = -1.5, xmx = 1.5, ymn = -1.5, ymx = 1.5, 
              vals = 1:36)

elev
# class      : RasterLayer 
# dimensions : 6, 6, 36  (nrow, ncol, ncell)
# resolution : 0.5, 0.5  (x, y)
# extent     : -1.5, 1.5, -1.5, 1.5  (xmin, xmax, ymin, ymax)
# crs        : +proj=longlat +datum=WGS84 +no_defs 
# source     : memory
# names      : layer 
# values     : 1, 36  (min, max)

 

 

레스터 객체의 행과 열의 일부분을 가져오는 것은 Base R의 '[' 연산자('[' operator) 를 사용합니다. raster_object[i, j] 또는 raster_object[Cell_IDs] 구문으로 i행과 j 열을 가져올 수 있습니다.  레스터 객체로 부터 일부분 가져오기는 '(a) 지리정보가 아닌 속성 정보를 가져오기(non-spatial subsetting)'와, '(b) 지리정보 가져오기 (spatial subsetting)' 의 두가지 유형으로 나눌 수 있는데요, 이번 포스팅은 '(a) 레스터 객체로부터 지리정보가 아닌 속성 정보를 가져오기'에 대해서만 다루겠습니다. 

 

elev[1, 1] 은 evel 레스터 객체로 부터 1행과 1열에 위치한 픽셀(pixel, cell) 에 위치한 속성 값을 가져온 것입니다. 

elev[1] 은 Cell ID 1번에 위치한 속성 값을 가져온 것입니다. (R 레스터는 좌측 상단에서 부터 1행 1열, Cell ID 1번이 시작하므로 elev[1, 1] 과 elev[1] 은 동일한 위치의 픽셀 속성 값을 가져옴)

 

## -- (1) Raster subsetting is done with the base R operator '['
## , which accepts a variety of inputs:
## - non-spatial subsetting: row-column indexing, cell IDs
## - spatial subsetting: coordinates, another spatial object

## subsetting the value of the top left pixel in the raster object 'elev'
## (1-1) row-column indexing
elev[1, 1]
# [1] 1

## (1-2) cell IDs
elev[1]
# [1] 1

 

 

레스터 객체의 전체 속성 값을 Cell ID의 순서대로 모두 가져오고 싶으면 values(), getValues() 함수 또는 [] 연산자를 사용하면 됩니다. 

 

## -- extracting all values or compelte rows
## (a) values()
values(elev)
# [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
# [26] 26 27 28 29 30 31 32 33 34 35 36


## or, eqivalently
# (b) getValues()
getValues(elev) 

# (c) []
elev[]

 

 

 

(2) 다층 레스터 객체(multi-layered laster objects)로 부터 일부 층(layers) 가져오기

    (subsetting layers from multi-layered laster objects, a raster brick or stack)

 

다층 레스터 객체(multi-layered laster objects)로는 RasterBrick 클래스와 RasterStack 클래스가 있습니다. (자세한 설명은 rfriend.tistory.com/617 를 참고하세요.)  이들 다층 레스터 객체로 부터 특정 층을 가져오려면 raster 패키지의 subset() 함수를 사용합니다. 

 

먼저, 예제로 사용할 수 있도록 grain 이라는 레스터 객체를 하나 더 만들고, 이를 (1)에서 만들었던 elev 레스터 객체에 stack() 함수로 쌓아서 다층 레스터 객체를 만들어보겠습니다

 

## -- for multi-layered raster objects 'stack' or 'brick', 
## values(), getValues() returns the cell value(s) for each layer. 
## 3 raster classes: https://rfriend.tistory.com/617

grain_order <- c("clay", "silt", "sand")
grain_char <- sample(grain_order, 36, replace = TRUE)
grain_fact <- factor(grain_char, levels = grain_order)
grain <- raster(nrows = 6, ncols = 6, res = 0.5, 
                xmn = -1.5, xmx = 1.5, ymn = -1.5, ymx = 1.5,
                vals = grain_fact)

grain
# class      : RasterLayer 
# dimensions : 6, 6, 36  (nrow, ncol, ncell)
# resolution : 0.5, 0.5  (x, y)
# extent     : -1.5, 1.5, -1.5, 1.5  (xmin, xmax, ymin, ymax)
# crs        : +proj=longlat +datum=WGS84 +no_defs 
# source     : memory
# names      : layer 
# values     : 1, 3  (min, max)
# attributes :
#   ID VALUE
# 1  clay
# 2  silt
# 3  sand


## making a RaserStack class object using stack() function
r_stack <- stack(elev, grain)
names(r_stack) <- c("elev", "grain")

r_stack
# class      : RasterStack 
# dimensions : 6, 6, 36, 2  (nrow, ncol, ncell, nlayers)
# resolution : 0.5, 0.5  (x, y)
# extent     : -1.5, 1.5, -1.5, 1.5  (xmin, xmax, ymin, ymax)
# crs        : +proj=longlat +datum=WGS84 +no_defs 
# names      : elev, grain 
# min values :    1,     1 
# max values :   36,     3

 

 

values() 로 r_stack 안의 각 층별 속성값을 조회해보면 아래처럼 "elev" 층과 "grain" 층의 각 Cell 별로 들어있는 속성 값이 들어있음을 알 수 있습니다. (values(r_stack), getValues(r_stack), r_stack[] 모두 동일하게 각 층의 모든 속성값을 반환합니다.) 

 

values(r_stack)
#getValues(r_stack)  # or equivalently
#r_stack[]           # or equivalently

# elev grain
# [1,]    1     3
# [2,]    2     3
# [3,]    3     1
# [4,]    4     3
# [5,]    5     3
# [6,]    6     1
# [7,]    7     1
# [8,]    8     3
# [9,]    9     1
# [10,]   10     1
# [11,]   11     2
# [12,]   12     3
# [13,]   13     2
# [14,]   14     2
# [15,]   15     1
# [16,]   16     3
# [17,]   17     1
# [18,]   18     1
# [19,]   19     2
# [20,]   20     3
# [21,]   21     1
# [22,]   22     3
# [23,]   23     1
# [24,]   24     2
# [25,]   25     1
# [26,]   26     2
# [27,]   27     3
# [28,]   28     2
# [29,]   29     3
# [30,]   30     1
# [31,]   31     1
# [32,]   32     3
# [33,]   33     1
# [34,]   34     2
# [35,]   35     1
# [36,]   36     1

 

 

다층 레스터 객체에서 특정 층을 가져오려면 raster::subset(raster_object, "layer_name") 구문을 사용할 수 있습니다.  아래 예에서는 r_stack 다층 레스터 객체에서 "elev" 층(layer)을 가져온 것입니다. 

 

## (2-1) raster::subset()
raster::subset(r_stack, "elev")
# class      : RasterLayer 
# dimensions : 6, 6, 36  (nrow, ncol, ncell)
# resolution : 0.5, 0.5  (x, y)
# extent     : -1.5, 1.5, -1.5, 1.5  (xmin, xmax, ymin, ymax)
# crs        : +proj=longlat +datum=WGS84 +no_defs 
# source     : memory
# names      : elev 
# values     : 1, 36  (min, max)


## (2-2) '[[' operator
r_stack[["elev"]]


## (2-3) '$' operator
r_stack$elev

 

 

 

(3) 레스터 객체의 일부 값을 수정하기 (modifying values in raster objects)

 

만약 레스터 객체의 일부 값을 새로운 값으로 덮어쓰기, 업데이트, 수정을 하고 싶으면 앞에서 소개한 subsetting 기능을 활용해서 원하는 위치의 행과 열, 픽셀 ID(Pixel ID, Cell ID) 또는 층(layer) 의 subsetting 한 부분에 새로운 속성 값을 할당(<-, =) 하면 됩니다. 

 

아래의 첫번째 예는 elev[1, 1]  <- 0 은 elev 레스터 객체의 1행 1열 위치의 셀에 새로운 값 '0'을 할당한 것입니다.  

아래의 두번째 예는 evel[1, 1:3] <- 0 은 elev 레스터 객체의 1행의 1~3열 위치 셀에 새로운 값 '0'을 할당하여 속성값을 수정한 것입니다. 

 

## -- (3) modifying the existing values by overwriting 
##        with a subsetting operation
elev[1, 1] <- 0

elev[]
# [1]  0  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
# [26] 26 27 28 29 30 31 32 33 34 35 36


## multiple cells can be modified with a subsetting operation. 
elev[1, 1:3] <- 0
values(elev)
# [1]  0  0  0  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
# [26] 26 27 28 29 30 31 32 33 34 35 36

 

 

[Reference]

* Geo-computation with R

  - 'Attribute data operations': geocompr.robinlovelace.net/attr.html

 

 

이번 포스팅이 많은 도움이 되었기를 바랍니다.

행복한 데이터 과학자 되세요!  :-)

 

반응형
Posted by Rfriend

댓글을 달아 주세요

지난번 포스팅까지는 R 지리공간 벡터 데이터 (geometry vector objects) 에 대한 속성정보 가져오기 (rfriend.tistory.com/622), 그룹별 집계하기 (rfriend.tistory.com/624), 두 개의 테이블 간 Join 하기 (rfriend.tistory.com/625), Key 값에 대한 문자열 매핑을 테이블 Join 하기 (rfriend.tistory.com/626), 새로운 속성정보 만들기와 지리정보 제거하기 (rfriend.tistory.com/627)에 대해서 소개하였습니다.

 

이번 포스팅부터는 레스터 객체 데이터셋을 조작하는 방법(manipulating raster objects dataset)을 몇 개의 포스팅으로 나누어서 소개하겠습니다.

 

그 중에서도 이번 포스팅에서는 R 문자형을 요인형으로 변환하여 레스터 객체 데이터의 속성으로 만드는 방법(making raster objects attributes by converting character to factor type)을 알아보겠습니다.

 

 

R raster objects' attritubes with numeric, integer, logical and factor types. No support for character.

 

 

 

 

R의 레스터 객체(raster objects)는 데이터 속성으로 숫자형(numeric), 정수형(integer), 논리형(logical), 요인형(factor) 데이터 유형을 지원하며, 문자형(character)은 지원하지 않습니다. 따라서 문자형으로 이루어진 범주형 변수 값(categorical variables' values)을 가지고 레스터 객체의 속성(attritubes)을 만들고 싶으면 (1) 먼저 문자형을 요인형으로 변환 (또는 논리형으로 변환)하고, --> (2) 요인형 값을 속성 값으로 해서 레스터 객체를 만들어야 합니다.

 

 

먼저, 레스터 객체 데이터에 대한 이해를 돕기위해, 정수형(integer) 값을 속성 값으로 가지는 레스터 객체를 raster 패키지의 raster() 함수를 사용해서 만들어보겠습니다. 아래의 예에서는 6 x 6 = 36개의 픽셀(pixels, cells)에,  x와 y의 좌표값의 최소~최대값 범위의 좌표에, 1~36까지의 정수를 속성 값으로 가지는 레스터 객체 데이터 입니다.

 

R 레스터 객체 데이터(Raster object in R)에 대한 보다 자세한 소개는 rfriend.tistory.com/589rfriend.tistory.com/605 를 참고하세요.

 

## ========================================================
## R GeoSpatial data analysis
## : Manipulating Raster Objects
## [reference] https://geocompr.robinlovelace.net/attr.html
## ========================================================

## raster data represent continuous surfaces. 
## Because of their unique structure, 
## subsetting and other operations on raster datasets work in a different way. 

library(raster)

## -- creating an example raster dataset
elev <- raster(nrows = 6, # integer > 0. Number of rows
               ncols = 6, # integer > 0. Number of columns
               #res = 0.5, # numeric vector of length 1 or 2 to set the resolution
               xmn = -1.5, # minimum x coordinate (left border)
               xmx = 1.5,  # maximum x coordinate (right border)
               ymn = -1.5, # minimum y coordinate (bottom border)
               ymx = 1.5,  # maximum y coordinate (top border)
               vals = 1:36) # values for the new RasterLayer

elev
# class      : RasterLayer 
# dimensions : 6, 6, 36  (nrow, ncol, ncell)
# resolution : 0.5, 0.5  (x, y)
# extent     : -1.5, 1.5, -1.5, 1.5  (xmin, xmax, ymin, ymax)
# crs        : +proj=longlat +datum=WGS84 +no_defs 
# source     : memory
# names      : layer 
# values     : 1, 36  (min, max)


plot(elev, main = 'raster datasets with numeric valeus')

 

 

이제 본론으로 넘어가서, 문자형으로 이루어진 범주형 변수를 요인형으로 변환한 후에, 이를 다시 레스터 객체의 속성으로 하여 레스터 객체를 만드는 방법을 소개하겠습니다.

 

 

(1)

 

문자형(character)을 요인형(factor)으로 변환할 때는 R base 패키지의 factor() 를 사용합니다. 이때 요인의 수준(levels)은 범주형 변수 내 유일한(unique) 문자열을 오름차순으로 정렬(sorting in an increasing order)하여 부여가 되는데요, 만약 요인의 수준(levels)을 특정 순서에 맞게 분석가가 수작업으로 설정을 하고 싶다면 levels 매개변수를 사용해서 직접 요인 수준을 입력을 해주면 됩니다.

 

아래 예에서는 모래의 굵기의 수준에 따라서 "점토(clay)", "미세모래(silt)", "모래(sand)"의 순서로 levels 매개변수를 이용해 요인의 수준을 설정해서, factor() 함수로 문자형으로 구성된 범주형 변수('grain_order')를 요인형 변수('grain_fact')으로 변환해준 것입니다. 

<--> 만약 factor() 함수로 요인형으로 변환할 때 levels 매개변수로 요인 수준의 순서를 지정해주지 않는다면, default 인 유일한 문자형 값들의 오름차순 정렬에 따라서 ["clay", "sand", "silt"] 의 순서로 설정이 되었을 것입니다.

 

## -- raster objects can contain values of class numeric, integer, logical or factor, 
## but not character.
## Raster objects can also contain categorical values of 
## class logical or factor variables in R.


grain_order <- c("clay", "silt", "sand")

## random sampling with replacement
set.seed(1004)
grain_char <- sample(grain_order, 36, replace = TRUE)
grain_char
# [1] "sand" "silt" "clay" "clay" "clay" "sand" "silt" "silt" "silt" "clay" "sand" "silt" "sand" "sand"
# [15] "sand" "clay" "sand" "clay" "sand" "clay" "sand" "sand" "silt" "clay" "sand" "silt" "sand" "clay"
# [29] "clay" "clay" "sand" "sand" "clay" "sand" "sand" "clay"


## converting character into factor
grain_fact <- factor(grain_char, 
                     # ordered factor: clay < silt < sand in terms of grain size.
                     levels = grain_order) 

grain_fact
# [1] sand silt clay clay clay sand silt silt silt clay sand silt sand sand sand clay sand clay sand clay
# [21] sand sand silt clay sand silt sand clay clay clay sand sand clay sand sand clay
# Levels: clay silt sand

 

 

 

(2) 요인형 값을 속성으로 한 레스터 객체 데이터 만들기 (creating raster objects with factor values)

 

위의 (1)번에서 문자열 값을 요인형으로 변환한 값을 가지고 raster 패키지의 raster() 함수를 사용해서 레스터 객체를 만들어보겠습니다.  행의 수(nrows) 6개, 열의 수(ncols) 6개의 총 36개 픽셀(pixcels, cells)을 가지는 레스터 객체에 raster(nrows = 6, ncols = 6, res = 0.5, vals = grain_fact)('grain_fact' 는 요인형 값을 가지는 벡터임) 을 입력하였습니다.

 

레스터 객체 데이터에 대해서 plot() 함수로 간단하게 시각화를 할 수 있습니다. 이때 요인형의 수준 값(levels)인 clay = 1, silt = 2, sand = 3 의 정수값을 가지고 시각화를 하게 됩니다.

 

## -- creating a raster dataset with factor variable in R
grain_raster <- raster(nrows = 6, ncols = 6, 
                       #res = 0.5, 
                       xmn = -1.5, xmx = 1.5, ymn = -1.5, ymx = 1.5, 
                       vals = grain_fact)

grain_raster
# class      : RasterLayer 
# dimensions : 6, 6, 36  (nrow, ncol, ncell)
# resolution : 0.5, 0.5  (x, y)
# extent     : -1.5, 1.5, -1.5, 1.5  (xmin, xmax, ymin, ymax)
# crs        : +proj=longlat +datum=WGS84 +no_defs 
# source     : memory
# names      : layer 
# values     : 1, 3  (min, max)
# attributes :
#   ID VALUE
# 1  clay
# 2  silt
# 3  sand


plot(grain_raster, main = 'raster dataset with factor(categorical) values')

 

 

이번에는 levels() 와 cbind() 함수를 사용해서 기존의 요인의 수준에다가 wetness = c("wet", "moist", "dry") 라는 새로운 요인 수준을 추가(adding new factor levels to the attritube table) 하여 보겠습니다. 

 

## -- Use the function levels() for retrieving and 
##    adding new factor levels to the attribute table:
levels(grain_raster)
# [[1]]
# ID VALUE
# 1  1  clay
# 2  2  silt
# 3  3  sand

## adding new factor levels
levels(grain_raster)[[1]] = cbind(levels(grain_raster)[[1]], 
                                  wetness = c("wet", "moist", "dry"))


## retrieving factor levels
levels(grain_raster)
# [[1]]
# ID VALUE wetness
# 1  1  clay     wet
# 2  2  silt   moist
# 3  3  sand     dry

 

 

레스터 객체 데이터의 "attritubes"라는 이름의 속성 테이블(attritube table)에 요인형 수준 값이 저장되어 있으며, 여기에 새로운 요인형 수준 값이 추가되는 것입니다.  ratify() 함수를 사용하면 이 속성 테이블(attritube table)에 저장되어 있는 값을 조회할 수 있습니다.

 

## -- The raster object stores the corresponding look-up table 
## or “Raster Attribute Table” (RAT) as a data frame in a new slot named attributes, 
## which can be viewed with ratify(grain)

ratify(grain_raster)
# class      : RasterLayer 
# dimensions : 6, 6, 36  (nrow, ncol, ncell)
# resolution : 0.5, 0.5  (x, y)
# extent     : -1.5, 1.5, -1.5, 1.5  (xmin, xmax, ymin, ymax)
# crs        : +proj=longlat +datum=WGS84 +no_defs 
# source     : memory
# names      : layer 
# values     : 1, 3  (min, max)
# attributes :
#   ID
# 1
# 2
# 3

 

 

레스터 객체의 각 픽셀(pixcel, cell)은 단 하나의 값만을 가질 수 있습니다. 이 예제에서는 각 픽셀이 하나의 요인 수준 값을 가지고 있는데요, 아래 예는 레스터 객체의 1번, 10번, 26번 픽셀의 요인 수준 값이 3, 1, 2 이고, 이들 요인 수준 값에 대응하는 속성 정보를 factorValues() 함수를 사용하여 속성 테이블(attritube table) 로 부터 매핑하여 보여주는 것입니다.

 

## -- looking up the attributes in the corresponding attribute table
grain_raster[c(1, 10, 26)]
# [1] 3 1 2

factorValues(grain_raster, grain_raster[c(1, 10, 26)])
# VALUE wetness
# 1  sand     dry
# 2  clay     wet
# 3  silt   moist

 

[Reference]

* Geo-Computation with R
  - 'Attribute data operations': geocompr.robinlovelace.net/attr.html

 

 

다음번 포스팅에서는 지리공간 레스터 데이터 가져오는 방법(Raster subsetting, rfriend.tistory.com/629) 에 대하여 소개하겠습니다.

 

이번 포스팅이 많은 도움이 되었기를 바랍니다.

행복한 데이터 과학자 되세요!

 

 

반응형
Posted by Rfriend

댓글을 달아 주세요

지난번 포스팅에서는 지리 벡터 데이터의 두 테이블의 Join Key 를 정규표현식을 사용한 패턴 매칭을 통해 일치를 시킨 후에 두 데이터 테이블을 Join 하는 방법(rfriend.tistory.com/626)을 소개하였습니다.  

 

이번 포스팅에서는 sf 클래스의 지리 벡터 데이터에서 Base R, dplyr, tidyr 의 R 패키지를 이용하여 기존 속성을 가지고 새로운 속성을 만드는 방법과, 지리 정보(spatial information)을 제거하는 방법을 소개하겠습니다. 

 

(1) Base R 로 지리 벡터 데이터에 새로운 속성 만들기

(2) dplyr 로 지리 벡터 데이터에 새로운 속성 만들기 : mutate(), transmute()

(3) tidyr 로 지리 벡터 데이터의 기존 속성을 합치거나 분리하기 : unite(), separate()

(4) dplyr 로 지리 벡터 데이터의 속성 이름 바꾸기 : rename(), setNames()

(5) 지리 벡터 데이터에서 지리 정보 제거하기 : st_drop_geometry()

 

 

 

(1) Base R 로 지리 벡터 데이터에 새로운 속성 만들기

 

먼저 예제로 사용할 sf 클래스 객체 데이터셋으로는, spData 패키지에 내장된 "world" 데이터셋을 사용하겠습니다. "world" 데이터셋에는 177개 국가의 지리기하 geometry 정보를 포함하고 있으며, sf 클래스를 사용하여 지리공간 데이터 처리 및 분석을 할 수 있습니다. 원래의 데이터를 덮어쓰지 않기 위해 "world_new" 라는 복사 데이터셋을 추가로 하나 더 만들었습니다.  그리고 데이터 전처리를 위해 dplyr, tidyr 패키지를 불러오겠습니다. 

 

## =============================================================
## GeoSpatial data analysis using R
## : Creating vector attributes and removing spatial information
## : reference: https://geocompr.robinlovelace.net/attr.html
## =============================================================

library(sf)
library(spData)
library(dplyr) # mutate(), transmute()
library(tidyr) # unite(), separate()

str(world)
# tibble [177 x 11] (S3: sf/tbl_df/tbl/data.frame)
# $ iso_a2   : chr [1:177] "FJ" "TZ" "EH" "CA" ...
# $ name_long: chr [1:177] "Fiji" "Tanzania" "Western Sahara" "Canada" ...
# $ continent: chr [1:177] "Oceania" "Africa" "Africa" "North America" ...
# $ region_un: chr [1:177] "Oceania" "Africa" "Africa" "Americas" ...
# $ subregion: chr [1:177] "Melanesia" "Eastern Africa" "Northern Africa" "Northern America" ...
# $ type     : chr [1:177] "Sovereign country" "Sovereign country" "Indeterminate" "Sovereign country" ...
# $ area_km2 : num [1:177] 19290 932746 96271 10036043 9510744 ...
# $ pop      : num [1:177] 8.86e+05 5.22e+07 NA 3.55e+07 3.19e+08 ...
# $ lifeExp  : num [1:177] 70 64.2 NA 82 78.8 ...
# $ gdpPercap: num [1:177] 8222 2402 NA 43079 51922 ...
# $ geom     :sfc_MULTIPOLYGON of length 177; first list element: List of 3
# ..$ :List of 1
# .. ..$ : num [1:8, 1:2] 180 180 179 179 179 ...
# ..$ :List of 1
# .. ..$ : num [1:9, 1:2] 178 178 179 179 178 ...
# ..$ :List of 1
# .. ..$ : num [1:5, 1:2] -180 -180 -180 -180 -180 ...
# ..- attr(*, "class")= chr [1:3] "XY" "MULTIPOLYGON" "sfg"
# - attr(*, "sf_column")= chr "geom"
# - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA NA NA NA NA
# ..- attr(*, "names")= chr [1:10] "iso_a2" "name_long" "continent" "region_un" ...


## do not overwrite the original data
world_new = world

 

 

sf 클래스 객체의 속성(Attributes)에 대해서 Base R 의 기본 구문인 '$' 와 '=', '<-' 을 사용해서 새로운 속성을 만들 수 있습니다.

가령, 국가별 인구("pop")를 면적("area_km2")로 나누어서 '(km2 면적 당) 인구밀도("pop_dens")' 라는 새로운 속성을 만들어 보겠습니다.  새로운 속성 칼럼인 "pop_dens"가 생성이 된 후에도 "world_new" 데이터셋은 sf 클래스 객체의 특성을 계속 유지하고 있습니다. 따라서 plot(world_new["pop_dens"] 로 시각화를 하면, 아래와 같이 다면 (multi-polygons) 기하로 표현되는 세계 국가별 지도 위에, 기존 칼럼을 사용해 계산해서 새로 만든 속성 칼럼인 '인구밀도("pop_dens")' 를 시각화할 수 있습니다. 

 

## creating new column using Base R
world_new$pop_dens = world_new$pop / world_new$area_km2


## attribute data operations preserve the geometry of the simple features
## plot() will visualize map using geometry
plot(world_new['pop_dens'])

 

 

(2) dplyr 로 지리 벡터 데이터에 새로운 속성 만들기 : mutate(), transmute()


위 (1)번의 Base R 대비 dplyr 를 사용하면, 한꺼번에 여러개의 신규 속성 칼럼을 생성할 수 있고, 체인('%>%')을 이용해서 파이프 연산자로 이전 작업 결과를 다음 작업으로 넘겨주는 방식으로 작업 흐름을 코딩할 수 있어 가독성이 뛰어난 코드를 짤 수 있는 장점이 있습니다. 

mutate() 함수는 기존 데이터셋에 새로 만든 데이터셋을 추가해주는 반면에, transmute() 함수는 기존의 칼럼은 모두 제거를 하고 새로 만든 칼럼과 지리기하 geometry 칼럼만을 반환합니다. (이때 sf 클래스 객체의 경우 '지리기하 geometry 칼럼'이 껌딱지처럼 달라붙어서 계속 따라 다닌다는 점이 중요합니다!)

 

 

## creating new columns using dplyr
world %>% 
  mutate(pop_dens = pop / area_km2)

# Simple feature collection with 177 features and 11 fields
# geometry type:  MULTIPOLYGON
# dimension:      XY
# bbox:           xmin: -180 ymin: -90 xmax: 180 ymax: 83.64513
# geographic CRS: WGS 84
# # A tibble: 177 x 12
# iso_a2 name_long  continent  region_un subregion type  area_km2     pop lifeExp gdpPercap                        geom pop_dens
# * <chr>  <chr>      <chr>      <chr>     <chr>     <chr>    <dbl>   <dbl>   <dbl>     <dbl> <MULTIPOLYGON [arc_degree]>    <dbl>
# 1 FJ     Fiji       Oceania    Oceania   Melanesia Sove~   1.93e4  8.86e5    70.0     8222. (((180 -16.06713, 180 -16.~    45.9 
# 2 TZ     Tanzania   Africa     Africa    Eastern ~ Sove~   9.33e5  5.22e7    64.2     2402. (((33.90371 -0.95, 34.0726~    56.0 
# 3 EH     Western S~ Africa     Africa    Northern~ Inde~   9.63e4 NA         NA         NA  (((-8.66559 27.65643, -8.6~    NA   
# 4 CA     Canada     North Ame~ Americas  Northern~ Sove~   1.00e7  3.55e7    82.0    43079. (((-122.84 49, -122.9742 4~     3.54
# 5 US     United St~ North Ame~ Americas  Northern~ Coun~   9.51e6  3.19e8    78.8    51922. (((-122.84 49, -120 49, -1~    33.5 
# 6 KZ     Kazakhstan Asia       Asia      Central ~ Sove~   2.73e6  1.73e7    71.6    23587. (((87.35997 49.21498, 86.5~     6.33
# 7 UZ     Uzbekistan Asia       Asia      Central ~ Sove~   4.61e5  3.08e7    71.0     5371. (((55.96819 41.30864, 55.9~    66.7 
# 8 PG     Papua New~ Oceania    Oceania   Melanesia Sove~   4.65e5  7.76e6    65.2     3709. (((141.0002 -2.600151, 142~    16.7 
# 9 ID     Indonesia  Asia       Asia      South-Ea~ Sove~   1.82e6  2.55e8    68.9    10003. (((141.0002 -2.600151, 141~   140.  
# 10 AR     Argentina  South Ame~ Americas  South Am~ Sove~   2.78e6  4.30e7    76.3    18798. (((-68.63401 -52.63637, -6~    15.4 
# # ... with 167 more rows

 

 

transmute() 함수를 사용해서 '(km2 면적당) 인구밀도 ("pop_dens")' 의 신규 속성을 만들면, 기존의 속성 칼럼들은 모두 제거 되고 새로운 '인구밀도' 속성 칼럼과 '지리기하 geom' 칼럼만 남게 됩니다. 

 

## -- transmute() drops all other existing columns except for the sticky geometry column
world %>% 
  transmute(pop_dens = pop / area_km2)

# Simple feature collection with 177 features and 1 field
# geometry type:  MULTIPOLYGON
# dimension:      XY
# bbox:           xmin: -180 ymin: -90 xmax: 180 ymax: 83.64513
# geographic CRS: WGS 84
# # A tibble: 177 x 2
# pop_dens                                                                                    geom
# *    <dbl>                                                             <MULTIPOLYGON [arc_degree]>
# 1    45.9  (((180 -16.06713, 180 -16.55522, 179.3641 -16.80135, 178.7251 -17.01204, 178.5968 -1...
# 2    56.0  (((33.90371 -0.95, 34.07262 -1.05982, 37.69869 -3.09699, 37.7669 -3.67712, 39.20222 ...
# 3    NA    (((-8.66559 27.65643, -8.665124 27.58948, -8.6844 27.39574, -8.687294 25.88106, -11....
# 4     3.54 (((-122.84 49, -122.9742 49.00254, -124.9102 49.98456, -125.6246 50.41656, -127.4356...
# 5    33.5  (((-122.84 49, -120 49, -117.0312 49, -116.0482 49, -113 49, -110.05 49, -107.05 49,...
# 6     6.33 (((87.35997 49.21498, 86.59878 48.54918, 85.76823 48.45575, 85.72048 47.45297, 85.16...
# 7    66.7  (((55.96819 41.30864, 55.92892 44.99586, 58.50313 45.5868, 58.68999 45.50001, 60.239...
# 8    16.7  (((141.0002 -2.600151, 142.7352 -3.289153, 144.584 -3.861418, 145.2732 -4.373738, 14...
# 9   140.   (((141.0002 -2.600151, 141.0171 -5.859022, 141.0339 -9.117893, 140.1434 -8.297168, 1...
# 10    15.4  (((-68.63401 -52.63637, -68.25 -53.1, -67.75 -53.85, -66.45 -54.45, -65.05 -54.7, -6...
# # ... with 167 more rows

 

 

 

(3) tidyr 로 지리 벡터 데이터의 기존 속성을 합치거나 분리하기 : unite(), separate()

 

tidyr 패키지에 있는 unite(data, col, lll, sep = "_", remove = TRUE) 함수를 사용하면 기존의 속성 칼럼을 합쳐서 새로운 속성 칼럼을 만들 수 있습니다. 이때 remove = TRUE 매개변수를 설정해주면 기존의 합치려는 두 개의 칼럼은 제거되고, 새로 만들어진 칼럼만 남게 됩니다. 

 

## -- unite(): pasting together existing columns.
library(tidyr)

world_unite = world %>% 
  unite("con_reg", continent:region_un, 
        sep = ":", 
        remove = TRUE)

world_unite
# Simple feature collection with 177 features and 9 fields
# geometry type:  MULTIPOLYGON
# dimension:      XY
# bbox:           xmin: -180 ymin: -90 xmax: 180 ymax: 83.64513
# geographic CRS: WGS 84
# # A tibble: 177 x 10
# iso_a2 name_long   con_reg     subregion    type    area_km2     pop lifeExp gdpPercap                                    geom
# <chr>  <chr>       <chr>       <chr>        <chr>      <dbl>   <dbl>   <dbl>     <dbl>             <MULTIPOLYGON [arc_degree]>
# 1 FJ     Fiji        Oceania:Oc~ Melanesia    Sovere~   1.93e4  8.86e5    70.0     8222. (((180 -16.06713, 180 -16.55522, 179.3~
# 2 TZ     Tanzania    Africa:Afr~ Eastern Afr~ Sovere~   9.33e5  5.22e7    64.2     2402. (((33.90371 -0.95, 34.07262 -1.05982, ~
# 3 EH     Western Sa~ Africa:Afr~ Northern Af~ Indete~   9.63e4 NA         NA         NA  (((-8.66559 27.65643, -8.665124 27.589~
# 4 CA     Canada      North Amer~ Northern Am~ Sovere~   1.00e7  3.55e7    82.0    43079. (((-122.84 49, -122.9742 49.00254, -12~
# 5 US     United Sta~ North Amer~ Northern Am~ Country   9.51e6  3.19e8    78.8    51922. (((-122.84 49, -120 49, -117.0312 49, ~
# 6 KZ     Kazakhstan  Asia:Asia   Central Asia Sovere~   2.73e6  1.73e7    71.6    23587. (((87.35997 49.21498, 86.59878 48.5491~
# 7 UZ     Uzbekistan  Asia:Asia   Central Asia Sovere~   4.61e5  3.08e7    71.0     5371. (((55.96819 41.30864, 55.92892 44.9958~
# 8 PG     Papua New ~ Oceania:Oc~ Melanesia    Sovere~   4.65e5  7.76e6    65.2     3709. (((141.0002 -2.600151, 142.7352 -3.289~
# 9 ID     Indonesia   Asia:Asia   South-Easte~ Sovere~   1.82e6  2.55e8    68.9    10003. (((141.0002 -2.600151, 141.0171 -5.859~
# 10 AR     Argentina   South Amer~ South Ameri~ Sovere~   2.78e6  4.30e7    76.3    18798. (((-68.63401 -52.63637, -68.25 -53.1, ~
# # ... with 167 more rows

 

 

tidyr 패키지의 separate(data, col, sep = "[^:alnum:]]+", remove = TRUE, ...) 함수는 기존에 존재하는 칼럼을 구분자(sep)를 기준으로 두 개의 칼럼으로 분리(separation, splitting)를 해줍니다. 이때 remove = TRUE 매개변수를 설정해주면 기존의 분리하려는 원래 속성 칼럼은 제거가 됩니다.

역시, unite() 또는 separate() 함수 적용 후의 결과 데이터셋은 sf 클래스 객체로서, 제일 뒤에 지리기하 geometry 칼럼이 계속 붙어 있습니다. 

 

## -- separate(): splitting one column into multiple columns
## : using either a regular expression or character positions. 
world_separate <- world_unite %>% 
  separate(con_reg, 
           c("continent", "region_un"), 
           sep = ":")


world_separate
# Simple feature collection with 177 features and 10 fields
# geometry type:  MULTIPOLYGON
# dimension:      XY
# bbox:           xmin: -180 ymin: -90 xmax: 180 ymax: 83.64513
# geographic CRS: WGS 84
# # A tibble: 177 x 11
#    iso_a2 name_long  continent  region_un subregion   type   area_km2     pop lifeExp gdpPercap                              geom
#    <chr>  <chr>      <chr>      <chr>     <chr>       <chr>     <dbl>   <dbl>   <dbl>     <dbl>       <MULTIPOLYGON [arc_degree]>
#  1 FJ     Fiji       Oceania    Oceania   Melanesia   Sover~   1.93e4  8.86e5    70.0     8222. (((180 -16.06713, 180 -16.55522,~
#  2 TZ     Tanzania   Africa     Africa    Eastern Af~ Sover~   9.33e5  5.22e7    64.2     2402. (((33.90371 -0.95, 34.07262 -1.0~
#  3 EH     Western S~ Africa     Africa    Northern A~ Indet~   9.63e4 NA         NA         NA  (((-8.66559 27.65643, -8.665124 ~
#  4 CA     Canada     North Ame~ Americas  Northern A~ Sover~   1.00e7  3.55e7    82.0    43079. (((-122.84 49, -122.9742 49.0025~
#  5 US     United St~ North Ame~ Americas  Northern A~ Count~   9.51e6  3.19e8    78.8    51922. (((-122.84 49, -120 49, -117.031~
#  6 KZ     Kazakhstan Asia       Asia      Central As~ Sover~   2.73e6  1.73e7    71.6    23587. (((87.35997 49.21498, 86.59878 4~
#  7 UZ     Uzbekistan Asia       Asia      Central As~ Sover~   4.61e5  3.08e7    71.0     5371. (((55.96819 41.30864, 55.92892 4~
#  8 PG     Papua New~ Oceania    Oceania   Melanesia   Sover~   4.65e5  7.76e6    65.2     3709. (((141.0002 -2.600151, 142.7352 ~
#  9 ID     Indonesia  Asia       Asia      South-East~ Sover~   1.82e6  2.55e8    68.9    10003. (((141.0002 -2.600151, 141.0171 ~
# 10 AR     Argentina  South Ame~ Americas  South Amer~ Sover~   2.78e6  4.30e7    76.3    18798. (((-68.63401 -52.63637, -68.25 -~
# # ... with 167 more rows

 

 

 

(4) dplyr 로 지리 벡터 데이터의 속성 이름 바꾸기 : rename(), setNames()

 

dplyr 패키지로 sf 객체이자 data frame 인 데이터셋의 특정 속성 변수 이름을 바꾸려면 rename(data, new_name = old_name) 함수를 사용합니다. 

 

## -- dplyr::rename()
names(world)
# [1] "iso_a2"    "name_long" "continent" "region_un" "subregion" "type"      "area_km2"  "pop"       "lifeExp"   "gdpPercap"
# [11] "geom"


world_nm <- world %>% 
  rename(name = name_long)

names(world_nm)
# [1] "iso_a2"    "name"      "continent" "region_un" "subregion" "type"      "area_km2"  "pop"       "lifeExp"   "gdpPercap"
# [11] "geom"

 

 

만약 여러개의 속성 칼럼 이름을 한꺼번에 변경, 혹은 부여하고 싶다면 stats 패키지의 setNames(object = nm, nm) 함수를 dplyr의 체인 파이프 라인에 같이 사용하면 편리합니다. 

 

##-- setNames() changes all column names at once
new_names <- c("i", "n", "c", "r", "s", "t", "a", "p", "l", "gP", "geom")

world_setnm <- world %>% 
  setNames(new_names)

names(world_setnm)
# [1] "i"    "n"    "c"    "r"    "s"    "t"    "a"    "p"    "l"    "gP"   "geom"

 

 

 

(5) 지리 벡터 데이터에서 지리 정보 제거하기 : st_drop_geometry()

 

위의 (1)~(4)번까지의 새로운 속성 칼럼 생성/ 이름 변경 등의 작업을 하고 난 결과 데이터셋에는 항상 '지리기하 geometry' 칼럼이 착 달라붙어서 따라 다녔습니다. 이게 지리공간 데이터 분석을 할 때 '속성(attributes)' 정보와 '지리기하 geometry' 정보를 같이 연계하여 분석해야 하는 경우에는 매우 유용합니다. 하지만, 경우에 따라서는 '속성'정보만 필요하고 '지리기하 geometry' 정보는 필요하지 않을 수도 있을 텐데요, 이럴때 많은 저장 용량과 처리 부하를 차지하는 '지리기하 geometry' 칼럼은 제거하고 싶을 것입니다. 이때 sf 패키지의 st_drop_geometry() 함수를 사용해서 sf 클래스 객체에 들어있는 '지리기하 geometry' 정보를 제거할 수 있습니다. 

 

## -- st_drop_geometry(): remove the geometry

## "sf" class
class(world)
# [1] "sf"         "tbl_df"     "tbl"        "data.frame"

## "data.frame"
world_data <- world %>% st_drop_geometry()
class(world_data)
# [1] "tbl_df"     "tbl"        "data.frame"

 

 

[Reference]

[1] Geo-computation with R - Attribute data operations
     : geocompr.robinlovelace.net/attr.html

 

 

이번 포스팅이 많은 도움이 되었기를 바랍니다. 

행복한 데이터 과학자 되세요! :-)

 

 

반응형
Posted by Rfriend

댓글을 달아 주세요

지난번 포스팅에서는 두 개의 sf 클래스 객체의 지리 벡터 데이터 테이블을 R dplyr 패키지의 함수를 사용하여 Mutating Joins, Filtering Joins, Nesting Joins 하는 방법을 소개하였습니다(rfriend.tistory.com/625). 

 

이번 포스팅에서는 여기서 특수한 경우로 조금 더 깊이 들어가서, 두 테이블을 Join 하는 기준이 되는 Key 칼럼이 문자열로 되어 있고, 데이터 표준화가 미흡한 문제로 인해 정확하게 매칭이 안되어서 Join 이 안되는 경우에, R의 stringr 패키지를 사용해 정규 표현식의 문자열 매칭(a string matching using regular expression)으로 Key 값을 변환하여 두 테이블을 Join 하는 방법을 소개하겠습니다. 

 

 

 

먼저, 전세계 국가별 지리기하와 속성 정보를 모아놓은 sf 클래스 객체의 지리 벡터 데이터셋인 "world" 와, 2016년과 2017년 국가별 커피 생산량을 집계한 data frame 인 "coffee_data" 의 두 개 데이터셋을 spData 로 부터 가져오겠습니다. 

그리고 두 개 테이블 Join 을 위해 dplyr 패키지를 불러오고, 정규 표현식을 이용한 문자열 매칭을 위해 stringr 패키지를 불러오겠습니다. 

 

"world" 데이터셋은 177개의 행(국가)과 11개의 열(속성(attritubes)과 지리기하 칼럼(gemgraphy column)) 으로 이루어져 있습니다. "coffee_data"는 47개의 행과 3개의 열로 구성되어 있습니다. 

 

## =========================================================
## inner join using a string matching
## - reference: https://geocompr.robinlovelace.net/attr.html
## =========================================================

library(sf)
library(spData)
library(dplyr)
library(stringr) # for a string matching

## -- two geography vector dataset tables : world, coffee_data
## -- (a) world: World country pologons in spData
names(world)
# [1] "iso_a2"  "name_long" "continent" "region_un" "subregion" "type"  "area_km2"  "pop"  "lifeExp"   "gdpPercap"
# [11] "geom"

dim(world)
# [1] 177  11


## -- (b) coffee_data: World coffee productiond data in spData
## : estimated values for coffee production in units of 60-kg bags in each year
names(coffee_data)
# [1] "name_long"      "coffee_production_2016" "coffee_production_2017"

dim(coffee_data)
# [1] 47  3

 

 

(1) 두 테이블 inner join 하기: inner_join(x, y, by)

 

"world"와 "coffee_data"의 두개 데이터 테이블을 inner join 해보면 45개의 행(즉, 국가)과 13개의 열(= "world"로 부터 11개의 칼럼 + "coffee_data"로 부터 2개의 칼럼) 으로 이루어진 Join 결과를 반환합니다. 

위에서 "coffee_data" 데이터셋이 47개의 행으로 이루어졌다고 했는데요, inner join 한 결과는 행이 45개로서 2개가 서로 차이가 나는군요. 

 

## -- inner join
world_coffee_inner = inner_join(x = world, 
                                y = coffee_data, 
                                by = "name_long")

## or shortly
world_coffee_inner = inner_join(world, coffee_data)
# Joining, by = "name_long"


dim(world_coffee_inner)
# [1] 45 13


nrow(world_coffee_inner)
# [1] 45

 

 

(2) 두 문자열의 원소 차이 알아보고 문자열 매칭으로 찾아보기: setdiff(), str_subset()

 

Join 전과 후에 어느 국가에서 차이가 나는지 확인해 보기 위해 setdiff() 함수를 사용해서 Join의 Key로 사용이 되었던 'name_long' (긴 국가 이름)에 대해 "coffee_data" 와 "world" 데이터의 원소 간 차이를 구해보았습니다. 그랬더니 ["Congo, Dem. Rep. of", "Others"] 의 2개 'name_long' 에서 차이가 있네요. 

 

다음으로, "world" 의 'name_long' 칼럼의 원소 중에서  "Dem"으로 시작하고 "Congo"를 포함하고 있는 문자열을 stringr 패키지의 str_subset(string, pattern) 함수를 사용해 정규 표현식의 문자열 매칭으로 찾아보겠습니다. "world" 데이터셋의 'name_long' 칼럼에는 "Democratic Republic of the Congo" 라는 이름으로 데이터가 들어가 있네요. ("coffee_data"  데이터셋에는 "Confo, Dem. Rep. of" 라고 들어가 있다보니, 서로 같은 국가임에도 left_join() 을 했을 때 서로 정확하게 매칭이 안되어 Join 이 안되었습니다.)

 

참고로, str_subset() 은 x[str_detect(x, pattern)] 의 wrapper 입니다. 그리고 grep(pattern, x, value = TRUE) 와 동일한 역할을 수행합니다. 

 

## setdiff(): calculates the set difference of subsets of two data frames
setdiff(coffee_data$name_long, world$name_long)
# [1] "Congo, Dem. Rep. of" "Others"


## string matching (regex) function from the stringr package
str_subset(world$name_long, "Dem*.+Congo")
# [1] "Democratic Republic of the Congo"

 

 

 

(3) 문자열 매칭으로 Key 값 업데이트 하고, 다시 두 테이블 inner join 하기

 

이제 Join Key로 사용하는 'name_long' 칼럼에서 "Congo" 국가에 대한 표기가 "world" 와 "coffee_data" 의 두 개 데이터셋이 서로 조금 다르다는 이유로 Join 이 안된다는 문제를 해결해 보겠습니다. 

grepl(pattern, x) 함수로 "coffee_data" 데이터셋의 'name_long' 칼럼에서 "Congo" 가 들어있는 행을 찾아서, 그 행의 값의 str_subset() 의 정규표현식 문자열 매칭으로 찾은 (str_subset(world$name_long, "Dem*.+Congo") 이름인 "Demogratic Republic of the Congo" 라는 이름으로 대체를 해보겠습니다. 이렇게 하면 "world"와 "coffee_data"에 있는 "Congo" 국가의 긴 이름이 동일하게 "Demogratic Republic of Congo"로 되어 Join 이 제대로 될 것입니다. 

 

## updating 'name_long' values using a string matching
coffee_data$name_long[grepl("Congo", coffee_data$name_long)] = 
  str_subset(world$name_long, "Dem*.+Congo")

## inner join again using an updated key
world_coffee_match = inner_join(world, coffee_data)
#> Joining, by = "name_long"

nrow(world_coffee_match)
#> [1] 46

 

 

 참고로, R에서 문자열 패턴 매칭을 할 때 grepl(pattern, x) 은 패턴 매칭되는 여부에 대해 TRUE, FALSE 로 블러언 값을 반환하는 반면에, grep(pattern, x) 은 패턴 매칭이 되는(TRUE) 위치 인덱스(Position Index)를 반환합니다. 

 

## -- grepl: pattern matching and returns boolean
grepl("Congo", coffee_data$name_long)
# [1] FALSE FALSE FALSE FALSE FALSE FALSE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
# [21] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
# [41] FALSE FALSE FALSE FALSE FALSE FALSE FALSE


## -- grep: pattern matching and returns position
grep("Congo", coffee_data$name_long)
# [1] 7

 

 

 

(4) Join 할 때 테이블 쓰는 순서의 중요성

 

dplyr 패키지로 두 테이블을 Join 할 때 왼쪽(x, LHS, Left Hand Side)에 써주는 테이블의 데이터 구조로 Join 한 결과를 반환합니다. 즉, Join 할 테이블을 써주는 순서가 중요합니다. 

가령, 아래의 예에서는 "world" 가 'sf' 클래스의 지리 벡터 객체이고, 'coffee_data'는 tydiverse의 tibble, data.frame 객체입니다. left_join(world, coffee_data) 로 'world' 의 'sf' 지리 벡터 객체를 Join 할 때 왼쪽(LHS, x)에 먼저 써주면 Join 한 결과도 'sf' 클래스의 지리 벡터 객체가 됩니다.(R이 지리공간 벡터 데이터임을 알고 'sf' 클래스를 적용한 지리공간 데이터 처리 및 분석이 가능함). 

반면에, left_join(coffee_data, world) 로 'coffee_data'의 'data.frame'을 Join 할 때 왼쪽(LHS, x)에 먼저 써주면 Join 한 결과도 'data.frame' 객체가 반환됩니다. (지리공간 'sf' 클래스가 더이상 아님) 

 

## starting with a non-spatial dataset and 
## adding variables from a simple features object.
## the result is not another simple feature object, 
## but a data frame in the form of a tidyverse tibble: 
## the output of a join tends to match its first argument.

## -- (a) 'sf' object first, then returns 'sf' object.
world_coffee = left_join(world, coffee_data)
#> Joining, by = "name_long"

class(world_coffee)
# [1] "sf"         "tbl_df"     "tbl"        "data.frame"


## -- (b) 'data.frame' object first, then returns 'data.frame' object.
coffee_world = left_join(coffee_data, world)
#> Joining, by = "name_long"

class(coffee_world)
#> [1] "tbl_df"     "tbl"        "data.frame"

 

 

(5) data.frame을 'sf' 클래스 객체로 변환하기

 

'sf' 패키지의 st_as_df() 함수를 사용하면 data.frame 을 'sf' 클래스 객체로 변환할 수 있습니다. 

## -- converting data.frame to 'sf' class object

st_as_sf(coffee_world)

# imple feature collection with 47 features and 12 fields (with 2 geometries empty)
# geometry type:  MULTIPOLYGON
# dimension:      XY
# bbox:           xmin: -117.1278 ymin: -33.76838 xmax: 156.02 ymax: 35.49401
# geographic CRS: WGS 84
# # A tibble: 47 x 13
# name_long coffee_producti~ coffee_producti~ iso_a2 continent region_un subregion type  area_km2     pop lifeExp gdpPercap
# <chr>                <int>            <int> <chr>  <chr>     <chr>     <chr>     <chr>    <dbl>   <dbl>   <dbl>     <dbl>
#   1 "Angola"                NA               NA AO     Africa    Africa    Middle A~ Sove~ 1245464.  2.69e7    60.9     6257.
# 2 "Bolivia"                3                4 BO     South Am~ Americas  South Am~ Sove~ 1085270.  1.06e7    68.4     6325.
# 3 "Brazil"              3277             2786 BR     South Am~ Americas  South Am~ Sove~ 8508557.  2.04e8    75.0    15374.
# 4 "Burundi"               37               38 BI     Africa    Africa    Eastern ~ Sove~   26239.  9.89e6    56.7      803.
# 5 "Cameroo~                8                6 CM     Africa    Africa    Middle A~ Sove~  460325.  2.22e7    57.1     3196.
# 6 "Central~               NA               NA CF     Africa    Africa    Middle A~ Sove~  621860.  4.52e6    50.6      597.
# 7 "Congo, ~                4               12 NA     NA        NA        NA        NA         NA  NA         NA         NA 
# 8 "Colombi~             1330             1169 CO     South Am~ Americas  South Am~ Sove~ 1151883.  4.78e7    74.0    12716.
# 9 "Costa R~               28               32 CR     North Am~ Americas  Central ~ Sove~   53832.  4.76e6    79.4    14372.
# 10 "C\u00f4~              114              130 CI     Africa    Africa    Western ~ Sove~  329826.  2.25e7    52.5     3055.
# # ... with 37 more rows, and 1 more variable: geom <MULTIPOLYGON [arc_degree]>

 

다음번 포스팅에서는  '지리공간 벡터 데이터에서 새로운 속성을 만들고 지리공간 정보를 제거하는 방법'에 대해서 알아보겠습니다. 

 

[Reference]

- Geocomputation with R, 'Attribute data operations': geocompr.robinlovelace.net/attr.html

 

 

이번 포스팅이 많은 도움이 되었기를 바랍니다. 

행복한 데이터 과학자 되세요!  :-)

 

반응형
Posted by Rfriend

댓글을 달아 주세요