지난 포스팅에서는 R data.table 패키지에서 선형회귀 모델의 오른쪽 부분의 변수 조합을 .SD, .SDcols, lapply(), sapply()를 사용해서 간단하게 단 몇줄의 코드로 처리하는 방법(rfriend.tistory.com/609)을 소개하였습니다.

 

이번 포스팅에서는 R data.table 패키지에서 '조건이 있는 상태에서 Key를 기준으로 데이터셋을 Left Join 하는 방법 (Conditional Joins)'을 소개하겠습니다.  그리고 base R이나 dplyr 대비 data.table의 조건이 있는 경우의 데이터셋끼리 병합이 얼마나 간단한지 비교를 해보겠습니다.  이번 포스팅은 R data.table의 Vignettes 을 참조하였습니다.

 

(1) data.table을 이용한 조건이 있는 경우의 Left Join

(2) dplyr을 이용한 조건이 있는 경우의 Left Join 비교

 

 

 

먼저, data.table 패키지를 불러오고, 예제로 사용할 데이터로 Lahman 패키지에 들어있는 야구 투구 통계 데이터인 'Pitching' 데이터셋과 야구 팀 성적 통계 데이터인 'Teams' 데이터셋을 참조해서 data.table과 data.frame으로 만들어보겠습니다.

 

library(data.table)

## Lahman database on baseball
#install.packages("Lahman")
library(Lahman)

data("Pitching")

## data.frame
Pitching_df <- data.frame(Pitching)


## data.table
setDT(Pitching)

str(Pitching)
# Classes 'data.table' and 'data.frame':	47628 obs. of  30 variables:
#   $ playerID: chr  "bechtge01" "brainas01" "fergubo01" "fishech01" ...
# $ yearID  : int  1871 1871 1871 1871 1871 1871 1871 1871 1871 1871 ...
# $ stint   : int  1 1 1 1 1 1 1 1 1 1 ...
# $ teamID  : Factor w/ 149 levels "ALT","ANA","ARI",..: 97 142 90 111 90 136 111 56 97 136 ...
# $ lgID    : Factor w/ 7 levels "AA","AL","FL",..: 4 4 4 4 4 4 4 4 4 4 ...
# $ W       : int  1 12 0 4 0 0 0 6 18 12 ...
# $ L       : int  2 15 0 16 1 0 1 11 5 15 ...
# $ G       : int  3 30 1 24 1 1 3 19 25 29 ...
# $ GS      : int  3 30 0 24 1 0 1 19 25 29 ...
# $ CG      : int  2 30 0 22 1 0 1 19 25 28 ...
# $ SHO     : int  0 0 0 1 0 0 0 1 0 0 ...
# $ SV      : int  0 0 0 0 0 0 0 0 0 0 ...
# $ IPouts  : int  78 792 3 639 27 3 39 507 666 747 ...
# $ H       : int  43 361 8 295 20 1 20 261 285 430 ...
# $ ER      : int  23 132 3 103 10 0 5 97 113 153 ...
# $ HR      : int  0 4 0 3 0 0 0 5 3 4 ...
# $ BB      : int  11 37 0 31 3 0 3 21 40 75 ...
# $ SO      : int  1 13 0 15 0 0 1 17 15 12 ...
# $ BAOpp   : num  NA NA NA NA NA NA NA NA NA NA ...
# $ ERA     : num  7.96 4.5 27 4.35 10 0 3.46 5.17 4.58 5.53 ...
# $ IBB     : int  NA NA NA NA NA NA NA NA NA NA ...
# $ WP      : int  7 7 2 20 0 0 1 15 3 44 ...
# $ HBP     : int  NA NA NA NA NA NA NA NA NA NA ...
# $ BK      : int  0 0 0 0 0 0 0 2 0 0 ...
# $ BFP     : int  146 1291 14 1080 57 3 70 876 1059 1334 ...
# $ GF      : int  0 0 0 1 0 1 1 0 0 0 ...
# $ R       : int  42 292 9 257 21 0 30 243 223 362 ...
# $ SH      : int  NA NA NA NA NA NA NA NA NA NA ...
# $ SF      : int  NA NA NA NA NA NA NA NA NA NA ...
# $ GIDP    : int  NA NA NA NA NA NA NA NA NA NA ...
# - attr(*, ".internal.selfref")=<externalptr> 


data("Teams")
setDT(Teams)
str(Teams)
# Classes 'data.table' and 'data.frame':	2925 obs. of  48 variables:
#   $ yearID        : int  1871 1871 1871 1871 1871 1871 1871 1871 1871 1872 ...
# $ lgID          : Factor w/ 7 levels "AA","AL","FL",..: 4 4 4 4 4 4 4 4 4 4 ...
# $ teamID        : Factor w/ 149 levels "ALT","ANA","ARI",..: 24 31 39 56 90 97 111 136 142 8 ...
# $ franchID      : Factor w/ 120 levels "ALT","ANA","ARI",..: 13 36 25 56 70 85 91 109 77 9 ...
# $ divID         : chr  NA NA NA NA ...
# $ Rank          : int  3 2 8 7 5 1 9 6 4 2 ...
# $ G             : int  31 28 29 19 33 28 25 29 32 58 ...
# $ Ghome         : int  NA NA NA NA NA NA NA NA NA NA ...
# $ W             : int  20 19 10 7 16 21 4 13 15 35 ...
# $ L             : int  10 9 19 12 17 7 21 15 15 19 ...
# $ DivWin        : chr  NA NA NA NA ...
# $ WCWin         : chr  NA NA NA NA ...
# $ LgWin         : chr  "N" "N" "N" "N" ...
# $ WSWin         : chr  NA NA NA NA ...
# $ R             : int  401 302 249 137 302 376 231 351 310 617 ...
# $ AB            : int  1372 1196 1186 746 1404 1281 1036 1248 1353 2571 ...
# $ H             : int  426 323 328 178 403 410 274 384 375 753 ...
# $ X2B           : int  70 52 35 19 43 66 44 51 54 106 ...
# $ X3B           : int  37 21 40 8 21 27 25 34 26 31 ...
# $ HR            : int  3 10 7 2 1 9 3 6 6 14 ...
# $ BB            : int  60 60 26 33 33 46 38 49 48 29 ...
# $ SO            : int  19 22 25 9 15 23 30 19 13 28 ...
# $ SB            : int  73 69 18 16 46 56 53 62 48 53 ...
# $ CS            : int  16 21 8 4 15 12 10 24 13 18 ...
# $ HBP           : int  NA NA NA NA NA NA NA NA NA NA ...
# $ SF            : int  NA NA NA NA NA NA NA NA NA NA ...
# $ RA            : int  303 241 341 243 313 266 287 362 303 434 ...
# $ ER            : int  109 77 116 97 121 137 108 153 137 166 ...
# $ ERA           : num  3.55 2.76 4.11 5.17 3.72 4.95 4.3 5.51 4.37 2.9 ...
# $ CG            : int  22 25 23 19 32 27 23 28 32 48 ...
# $ SHO           : int  1 0 0 1 1 0 1 0 0 1 ...
# $ SV            : int  3 1 0 0 0 0 0 0 0 1 ...
# $ IPouts        : int  828 753 762 507 879 747 678 750 846 1548 ...
# $ HA            : int  367 308 346 261 373 329 315 431 371 573 ...
# $ HRA           : int  2 6 13 5 7 3 3 4 4 3 ...
# $ BBA           : int  42 28 53 21 42 53 34 75 45 63 ...
# $ SOA           : int  23 22 34 17 22 16 16 12 13 77 ...
# $ E             : int  243 229 234 163 235 194 220 198 218 432 ...
# $ DP            : int  24 16 15 8 14 13 14 22 20 22 ...
# $ FP            : num  0.834 0.829 0.818 0.803 0.84 0.845 0.821 0.845 0.85 0.83 ...
# $ name          : chr  "Boston Red Stockings" "Chicago White Stockings" "Cleveland Forest Citys" "Fort Wayne Kekiongas" ...
# $ park          : chr  "South End Grounds I" "Union Base-Ball Grounds" "National Association Grounds" "Hamilton Field" ...
# $ attendance    : int  NA NA NA NA NA NA NA NA NA NA ...
# $ BPF           : int  103 104 96 101 90 102 97 101 94 106 ...
# $ PPF           : int  98 102 100 107 88 98 99 100 98 102 ...
# $ teamIDBR      : chr  "BOS" "CHI" "CLE" "KEK" ...
# $ teamIDlahman45: chr  "BS1" "CH1" "CL1" "FW1" ...
# $ teamIDretro   : chr  "BS1" "CH1" "CL1" "FW1" ...
# - attr(*, ".internal.selfref")=<externalptr>

 

 

(1) data.table을 이용한 조건이 있는 경우의 Left Join

 

(1-1) 참여 경기 수 조건 하에 팀별 연도별 평균자책점 순위 (rank_in_team)

 

(조건) G > 5 : 6개 이상의 경기(G, game)에 참여하여 공을 던진 투수에 한하여,

(연산) ERA (평균 자책점, Earned Run Average)의 순위를 구해서 rank_in_team 변수를 만들되,

(그룹 기준) 'teamID' & 'yearID' 그룹 별로 ERA 순위(frank(ERA))를 구하여라.

 

조건, 연산, 그룹 기준이 모두 Pitching[조건, 연산, 그룹 기준] 의 구문으로 해서 단 한줄로 모두 처리가 가능합니다. 간결함의 끝판왕이라고나 할까요!

 

## -- rank in team using frank() function by teamID & yearID groups
## to exclude pitchers with exceptional performance in a few games,
##   subset first; then define rank of pitchers within their team each year
##   (in general, we should put more care into the 'ties.method' of frank)
Pitching[G > 5, rank_in_team := frank(ERA), by = .(teamID, yearID)]

head(Pitching)
#    playerID yearID stint teamID lgID  W  L  G GS CG SHO SV IPouts   H  ER HR BB SO BAOpp   ERA
# 1: bechtge01   1871     1    PH1   NA  1  2  3  3  2   0  0     78  43  23  0 11  1    NA  7.96
# 2: brainas01   1871     1    WS3   NA 12 15 30 30 30   0  0    792 361 132  4 37 13    NA  4.50
# 3: fergubo01   1871     1    NY2   NA  0  0  1  0  0   0  0      3   8   3  0  0  0    NA 27.00
# 4: fishech01   1871     1    RC1   NA  4 16 24 24 22   1  0    639 295 103  3 31 15    NA  4.35
# 5: fleetfr01   1871     1    NY2   NA  0  1  1  1  1   0  0     27  20  10  0  3  0    NA 10.00
# 6: flowedi01   1871     1    TRO   NA  0  0  1  0  0   0  0      3   1   0  0  0  0    NA  0.00
#    IBB WP HBP BK  BFP GF   R SH SF GIDP rank_in_team
# 1:  NA  7  NA  0  146  0  42 NA NA   NA           NA
# 2:  NA  7  NA  0 1291  0 292 NA NA   NA            1
# 3:  NA  2  NA  0   14  0   9 NA NA   NA           NA
# 4:  NA 20  NA  0 1080  1 257 NA NA   NA            1
# 5:  NA  0  NA  0   57  0  21 NA NA   NA           NA
# 6:  NA  0  NA  0    3  1   0 NA NA   NA           NA

 

 

(1-2) 팀내 순위 조건하에 다른 데이터셋에 앴는 팀 성적 (team_performance) Left join

 

(조건) Pitching 데이터셋에서 팀/연도별로 순위가 1등인 투수에 한해 (rank_in_team == 1)

(연산) Teams 데이터셋의 순위(Rank)를 가져다가 Pitching 데이터셋에 team_performance 이름의 변수로 만들되,

(병합 기준) Pitching 과 Teams 데이터셋의 'teamID'와 'yearID'를 기준으로 매칭하시오.

 

 

Pitching[rank_in_team == 1, # condition for rows
         team_performance := Teams[.SD, Rank, 
                                   on = c('teamID', 'yearID')] # left join by keys
         ]


head(Pitching)
#    playerID yearID stint teamID lgID  W  L  G GS CG SHO SV IPouts   H  ER HR BB SO BAOpp   ERA
# 1: bechtge01   1871     1    PH1   NA  1  2  3  3  2   0  0     78  43  23  0 11  1    NA  7.96
# 2: brainas01   1871     1    WS3   NA 12 15 30 30 30   0  0    792 361 132  4 37 13    NA  4.50
# 3: fergubo01   1871     1    NY2   NA  0  0  1  0  0   0  0      3   8   3  0  0  0    NA 27.00
# 4: fishech01   1871     1    RC1   NA  4 16 24 24 22   1  0    639 295 103  3 31 15    NA  4.35
# 5: fleetfr01   1871     1    NY2   NA  0  1  1  1  1   0  0     27  20  10  0  3  0    NA 10.00
# 6: flowedi01   1871     1    TRO   NA  0  0  1  0  0   0  0      3   1   0  0  0  0    NA  0.00
#    IBB WP HBP BK  BFP GF   R SH SF GIDP rank_in_team team_performance
# 1:  NA  7  NA  0  146  0  42 NA NA   NA           NA               NA
# 2:  NA  7  NA  0 1291  0 292 NA NA   NA            1                4
# 3:  NA  2  NA  0   14  0   9 NA NA   NA           NA               NA
# 4:  NA 20  NA  0 1080  1 257 NA NA   NA            1                9
# 5:  NA  0  NA  0   57  0  21 NA NA   NA           NA               NA
# 6:  NA  0  NA  0    3  1   0 NA NA   NA           NA               NA

 

 

 

(2) dplyr을 이용한 조건이 있는 경우의 Left Join 비교

 

(2-1) 참여 경기 수 조건하에 팀별 연도별 평균자책점 순위 (rank_in_team)

이제 위의 (1-1)에서 data.table로 수행했던 것과 동일한 과업을 dplyr 패키지로 수행해보겠습니다.

 

(조건) G > 5 : 6개 이상의 경기(G, game)에 참여하여 공을 던진 투수에 한하여,

(연산) ERA (평균 자책점, Earned Run Average)의 순위를 구해서 rank_in_team 변수를 만들되,

(그룹 기준) 'teamID' & 'yearID' 그룹 별로 ERA 순위(frank(ERA))를 구하여라.

 

1단계에서 G > 5 라는 조건으로 filter() 를 하여 그룹별로 rank_in_team 을 구해서 별도의 'rank_in_team_df' data.frame을 만든 후에, --> 2단계에서 이를 원본 'Pitching_df'에 left_join() 을 해서 left join 병합을 해주었습니다. G > 5 라는 조건(condition)이 들어감으로써 2단계로 나누어서 진행이 되다보니 코드가 길어졌습니다.

 

 

## -- doing the same operation using 'dplyr'
library(dplyr)

## -- rank in team by teamID & yearID group using dense_rank() window function
rank_in_team_df <- Pitching_df %>% 
  filter(G > 5) %>% 
  group_by(teamID, yearID) %>% 
  mutate(rank_in_team = dense_rank(ERA)) %>% 
  select(playerID, lgID, teamID, yearID, stint, rank_in_team)


## left outer join {dplyr}
Pitching_df <- left_join(Pitching_df, 
                         rank_in_team_df, 
                         by = c('playerID', 'lgID', 'teamID', 'yearID', 'stint'))

# ## left outer join {base}
# Pitching_df <- merge(x = Pitching_df, 
#                      y = rank_in_team_df, 
#                      by = c('playerID', 'lgID', 'teamID', 'yearID', 'stint'), 
#                      all.x = TRUE, all.y = FALSE)


head(Pitching_df)
# playerID yearID stint teamID lgID  W  L  G GS CG SHO SV IPouts   H  ER HR BB SO BAOpp   ERA IBB
# 1 bechtge01   1871     1    PH1   NA  1  2  3  3  2   0  0     78  43  23  0 11  1    NA  7.96  NA
# 2 brainas01   1871     1    WS3   NA 12 15 30 30 30   0  0    792 361 132  4 37 13    NA  4.50  NA
# 3 fergubo01   1871     1    NY2   NA  0  0  1  0  0   0  0      3   8   3  0  0  0    NA 27.00  NA
# 4 fishech01   1871     1    RC1   NA  4 16 24 24 22   1  0    639 295 103  3 31 15    NA  4.35  NA
# 5 fleetfr01   1871     1    NY2   NA  0  1  1  1  1   0  0     27  20  10  0  3  0    NA 10.00  NA
# 6 flowedi01   1871     1    TRO   NA  0  0  1  0  0   0  0      3   1   0  0  0  0    NA  0.00  NA
# WP HBP BK  BFP GF   R SH SF GIDP rank_in_team
# 1  7  NA  0  146  0  42 NA NA   NA           NA
# 2  7  NA  0 1291  0 292 NA NA   NA            1
# 3  2  NA  0   14  0   9 NA NA   NA           NA
# 4 20  NA  0 1080  1 257 NA NA   NA            1
# 5  0  NA  0   57  0  21 NA NA   NA           NA
# 6  0  NA  0    3  1   0 NA NA   NA           NA

 

 

 

(2-2) 팀내 순위 조건하에 다른 데이터셋에 앴는 팀 성적 (team_performance) Left join

 

이제 위의 (1-2)에서 data.table로 했던 것과 똑같은 과업을 dplyr로 해보겠습니다.

 

(조건) Pitching 데이터셋에서 팀/연도별로 순위가 1등인 투수에 한해 (rank_in_team == 1) (연산) Teams 데이터셋의 순위(Rank)를 가져다가 Pitching 데이터셋에 team_performance 이름의 변수로 만들되,
(병합 기준) Pitching 과 Teams 데이터셋의 'teamID'와 'yearID'를 기준으로 매칭하시오.

 

아래처럼 dplyr로 하게 되면 '팀/연도별로 순위가 1등인 투수에 한해(rank_in_team==1)' 라는 조건을 충족시키기 위해서 조건절을 포함한 ifelse() 절을 한번 더 수행해줘야 합니다.

 

## -- merging team performance
Pitching_df <- left_join(Pitching_df, 
                         Teams[, c('teamID', 'yearID', 'Rank')], 
                         by = c('teamID', 'yearID'))

## condition: rank_in_team == 1
Pitching_df$team_performance <- ifelse(Pitching_df$rank_in_team == 1, 
                                       Pitching_df$Rank, # if TRUE
                                       NA) # if FALSE

head(Pitching_df)
# bechtge01   1871     1    PH1   NA  1  2  3  3  2   0  0     78  43  23  0 11  1    NA  7.96  NA
# 2 brainas01   1871     1    WS3   NA 12 15 30 30 30   0  0    792 361 132  4 37 13    NA  4.50  NA
# 3 fergubo01   1871     1    NY2   NA  0  0  1  0  0   0  0      3   8   3  0  0  0    NA 27.00  NA
# 4 fishech01   1871     1    RC1   NA  4 16 24 24 22   1  0    639 295 103  3 31 15    NA  4.35  NA
# 5 fleetfr01   1871     1    NY2   NA  0  1  1  1  1   0  0     27  20  10  0  3  0    NA 10.00  NA
# 6 flowedi01   1871     1    TRO   NA  0  0  1  0  0   0  0      3   1   0  0  0  0    NA  0.00  NA
# WP HBP BK  BFP GF   R SH SF GIDP rank_in_team Rank team_performance
# 1  7  NA  0  146  0  42 NA NA   NA           NA    1               NA
# 2  7  NA  0 1291  0 292 NA NA   NA            1    4                4
# 3  2  NA  0   14  0   9 NA NA   NA           NA    5               NA
# 4 20  NA  0 1080  1 257 NA NA   NA            1    9                9
# 5  0  NA  0   57  0  21 NA NA   NA           NA    5               NA
# 6  0  NA  0    3  1   0 NA NA   NA           NA    6               NA

 

 

 

위에서 소개했던 똑같은 과업을 수행하는 data.table과 dplyr의 조건있는 데이터셋 병합(conditional left join) 예제 코드를 나란히 제시해서 비교해보면 아래와 같습니다. data.table이 dplyr (혹은 base R) 대비 조건있는 데이터셋 병합의 경우 비교할 수 없을 정도로 훨~씬 코드가 간결합니다!

 

구분

data.table

dplyr

(1)

conditional

left join

by key

w/ a list

library(data.table)

Pitching[G > 5, # condition

         rank_in_team := frank(ERA), 

         by = .(teamID, yearID)]

library(dplyr)

rank_in_team_df <- Pitching_df %>% 

  filter(G > 5) %>% 

  group_by(teamID, yearID) %>% 

  mutate(rank_in_team = dense_rank(ERA)) %>% 

  select(playerID, lgID, teamID, yearID, stint, rank_in_team)

## left outer join {dplyr}

Pitching_df <- left_join(

  Pitching_df, 

  rank_in_team_df, 

  by = c('playerID', 'lgID', 'teamID', 'yearID', 'stint'))

(2)

conditional

left join

by key

w/ a data.table

Pitching[rank_in_team == 1,          team_performance := Teams[           .SD, Rank, 
         on = c('teamID', 'yearID')]
         ]

## -- merging team performance

Pitching_df <- left_join(

  Pitching_df, 

  Teams[, c('teamID', 'yearID', 'Rank')], 

  by = c('teamID', 'yearID'))

## condition: rank_in_team == 1

Pitching_df$team_performance <- ifelse(

  Pitching_df$rank_in_team == 1, 

  Pitching_df$Rank, # if TRUE

  NA) # if FALSE

 

 

[ Reference ]

* R data.table vignettes 'Using .SD for Data Analysis'
  : cran.r-project.org/web/packages/data.table/vignettes/datatable-sd-usage.html

 

 

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

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

 

728x90
반응형
Posted by Rfriend
,

Rfriend님의
글이 좋았다면 응원을 보내주세요!

지난 포스팅에서는 R data.table 패키지에서 .SD 가 무엇이고, .SD와 .SDcols 를 사용해서 특정 칼럼만을 가져오는 방법(을 소개하였습니다. 그리고 lapply()도 같이 사용해서 여러개의 칼럼의 데이터 유형을 변환하는 방법도 같이 설명하였습니다. (rfriend.tistory.com/608)

 

이번 포스팅에서는 지난 포스팅에 이어서, R data.table의 .SD, .SDcols와 lapply(), combn() 함수를 같이 사용하여 '선형회귀 모델의 오른쪽 부분에 여러 변수의 모든 가능한 조합을 적용하여 다루는 방법'을 소개하겠습니다.

 

(1) 설명변수 간 모든 가능한 조합 만들기: combn(x, m, simplify = TRUE, ...)

(2) 설명변수 간 모든 가능한 조합으로 선형회귀모델 적합하여 회귀계수 가져오기

(3) 설명변수 조합별 회귀계수로 막대그래프 그리기: barplot()

 

만약 모델에서 사용하는 칼럼의 개수가 몇 개 안된다면 수작업으로 일일이 변수 이름을 나열해서 적어주면 됩니다. 하지만, 모델에서 사용하는 칼럼 개수가 수십개인 경우(** 몇 개의 칼럼 간의 조합을 사용하는 경우 기하급수적으로 칼럼 개수가 늘어남)에는 키보드를 쳐가면서 일일이 변수 이름(조합)을 나열하기가 시간도 오래 걸리고, 오타 에러를 유발할 위험도 높아집니다. 이럴 경우에 이번에 소개하는 내용이 매우 유용하고 간편하게 사용될 수 있습니다.

 

먼저, data.table 패키지를 불러오고, 예제로 사용할 데이터로 Lahman 패키지에 들어있는 야규 투구 통계 데이터인 'Pitching' 데이터셋을 Data.Table로 만들어보겠습니다.

 

library(data.table)

## Lahman database on baseball
#install.packages("Lahman")
library(Lahman)
data("Pitching")
setDT(Pitching)
str(Pitching)
# Classes 'data.table' and 'data.frame':	47628 obs. of  30 variables:
#   $ playerID: chr  "bechtge01" "brainas01" "fergubo01" "fishech01" ...
# $ yearID  : int  1871 1871 1871 1871 1871 1871 1871 1871 1871 1871 ...
# $ stint   : int  1 1 1 1 1 1 1 1 1 1 ...
# $ teamID  : Factor w/ 149 levels "ALT","ANA","ARI",..: 97 142 90 111 90 136 111 56 97 136 ...
# $ lgID    : Factor w/ 7 levels "AA","AL","FL",..: 4 4 4 4 4 4 4 4 4 4 ...
# $ W       : int  1 12 0 4 0 0 0 6 18 12 ...
# $ L       : int  2 15 0 16 1 0 1 11 5 15 ...
# $ G       : int  3 30 1 24 1 1 3 19 25 29 ...
# $ GS      : int  3 30 0 24 1 0 1 19 25 29 ...
# $ CG      : int  2 30 0 22 1 0 1 19 25 28 ...
# $ SHO     : int  0 0 0 1 0 0 0 1 0 0 ...
# $ SV      : int  0 0 0 0 0 0 0 0 0 0 ...
# $ IPouts  : int  78 792 3 639 27 3 39 507 666 747 ...
# $ H       : int  43 361 8 295 20 1 20 261 285 430 ...
# $ ER      : int  23 132 3 103 10 0 5 97 113 153 ...
# $ HR      : int  0 4 0 3 0 0 0 5 3 4 ...
# $ BB      : int  11 37 0 31 3 0 3 21 40 75 ...
# $ SO      : int  1 13 0 15 0 0 1 17 15 12 ...
# $ BAOpp   : num  NA NA NA NA NA NA NA NA NA NA ...
# $ ERA     : num  7.96 4.5 27 4.35 10 0 3.46 5.17 4.58 5.53 ...
# $ IBB     : int  NA NA NA NA NA NA NA NA NA NA ...
# $ WP      : int  7 7 2 20 0 0 1 15 3 44 ...
# $ HBP     : int  NA NA NA NA NA NA NA NA NA NA ...
# $ BK      : int  0 0 0 0 0 0 0 2 0 0 ...
# $ BFP     : int  146 1291 14 1080 57 3 70 876 1059 1334 ...
# $ GF      : int  0 0 0 1 0 1 1 0 0 0 ...
# $ R       : int  42 292 9 257 21 0 30 243 223 362 ...
# $ SH      : int  NA NA NA NA NA NA NA NA NA NA ...
# $ SF      : int  NA NA NA NA NA NA NA NA NA NA ...
# $ GIDP    : int  NA NA NA NA NA NA NA NA NA NA ...
# - attr(*, ".internal.selfref")=<externalptr>

 

(1) 설명변수 간 모든 가능한 조합 만들기: combn(x, m, simplify = TRUE, ...)

 

먼저, 코드 이해를 돕기 위해서 utils 패키지의 combn(x, m, simplify = TRUE, ...) 함수에 대해서 설명하겠습니다.

 

combn() 함수는 투입되는 객체 x 내 원소들 간의 m 개로 이루어진 모든 가능한 조합을 만들어줍니다. 이때 simplify = TRUE (default 설정) 이면 바로 아래의 예처럼 행렬(matrix), 배열(array) 형태로 간소화해서 조합의 결과를 나타내주며, simplify = FALSE 로 설정해주면 두번째 예처럼 원소 간 조합의 결과를 리스트(list) 형태로 해서 위에서 아래로 길게 나타내줍니다.

 

## combn {utils}
## : combn(x, m, FUN = NULL, simplify = TRUE, ...)
## : Generate all combinations of the elements of 'x' taken 'm' at a time
combn(x = letters[1:4], m = 2)
# [,1] [,2] [,3] [,4] [,5] [,6]
# [1,] "a"  "a"  "a"  "b"  "b"  "c" 
# [2,] "b"  "c"  "d"  "c"  "d"  "d"


## simplify = FALSE : returns the combinations in a list format
combn(x = letters[0:4], m = 2, simplify = FALSE)
# [[1]]
# [1] "a" "b"
# 
# [[2]]
# [1] "a" "c"
# 
# [[3]]
# [1] "a" "d"
# 
# [[4]]
# [1] "b" "c"
# 
# [[5]]
# [1] "b" "d"
# 
# [[6]]
# [1] "c" "d"

 

 

앞에서 combn() 함수로 인풋 x 의 모든 원소 간 조합을 만들 수 있다는 것을 알았으므로, 이제 c('yearID', 'teamID', 'G', 'L') 의 4개 설명변수 간 모든 조합을 만들어보겠습니다. 이때 m=0~4 (설명변수 개수) 까지 바꾸어가면서 lapply() 로 combn() 함수에 x = c('yearID', 'teamID', 'G', 'L') 의 4개 원소를 가진 인풋을 적용해서 가능한 모든 조합을 만들어보겠습니다. 그리고 조합의 결과를 unlist() 해서 아래의 models 결과처럼 각 조합이 문자형 벡터(character vector)로 해서 차곡 차곡 리스트에 묶여있도록 하겠습니다.

 

## this generates a list of the 2^k possible extra variables
## for models of the form ERA ~ G + (...)
extra_var <- c('yearID', 'teamID', 'G', 'L')
models <- unlist(
  lapply(0L:length(extra_var), combn, x = extra_var, simplify = FALSE), 
  recursive = FALSE # the function will not recurse beyond the first level items in x
)

models
# [[1]]
# character(0)
# 
# [[2]]
# [1] "yearID"
# 
# [[3]]
# [1] "teamID"
# 
# [[4]]
# [1] "G"
# 
# [[5]]
# [1] "L"
# 
# [[6]]
# [1] "yearID" "teamID"
# 
# [[7]]
# [1] "yearID" "G"     
# 
# [[8]]
# [1] "yearID" "L"     
# 
# [[9]]
# [1] "teamID" "G"     
# 
# [[10]]
# [1] "teamID" "L"     
# 
# [[11]]
# [1] "G" "L"
# 
# [[12]]
# [1] "yearID" "teamID" "G"     
# 
# [[13]]
# [1] "yearID" "teamID" "L"     
# 
# [[14]]
# [1] "yearID" "G"      "L"     
# 
# [[15]]
# [1] "teamID" "G"      "L"     
# 
# [[16]]
# [1] "yearID" "teamID" "G"      "L"

 

 

위의 lapply() 와 combn() 함수를 같이 쓴 코드가 잘 이해가 안가신다면, 바로 아래처럼 m=0~4 까지 일일이 반복적으로 길게 쓴 코드를 lapply() 사용해서 좀더 간결하게 쓴 것이 위의 코드라고 생각하면 되겠습니다.

 

## or, equivalently by manual
unlist(
  list(
    combn(x=extra_var, m=0, simplify=F), 
    combn(x=extra_var, m=1, simplify=F), 
    combn(x=extra_var, m=2, simplify=F), 
    combn(x=extra_var, m=3, simplify=F), 
    combn(x=extra_var, m=4, simplify=F)
  ), 
  recursive = FALSE
)

 

 

 

(2) 설명변수 간 모든 가능한 조합으로 선형회귀모델 적합하여 회귀계수 가져오기

 

R에서 선형회귀모형은 lm(y ~ x, data) 형태의 구문을 사용합니다. 그리고 coef(lm(y ~ x, data)) 로 적합된 회귀모형의 회귀계수에 접근할 수 있습니다.

아래의 Pitching[ , coef(lm(ERA ~ ., data = .SD))['W'], .SDcols = c('W', rhs)] 코드가 무척 어려워보일 수 있는데요, 하나씩 설명해보자면요, data = .SD, .SDcols = c('W', rhs) 는 Pitching 데이터셋을 참조해서 .SDcols 에 있는 칼럼(즉, 'W' 변수와 (1)번에서 만든 models 설명변수 조합에 해당하는 rhs 변수들)만 선별해서 가져오라는 뜻입니다. 이들 설명변수 조합을 사용해 lm(ERA ~ .) 로 회귀모형을 적합하며, coef(lm(ERA ~ ., data=.SD))['W'] 는 변수 'W'의 적합된 추정 회귀계수를 가져오라는 뜻입니다.

## -- coefficient of 'W' from linear regression model
## using ERA ~ . and data = .SD, then varying which
##   columns are included in .SD allows us to perform this
##   iteration over 16 models succinctly.
##   coef(.)['W'] extracts the W coefficient from each model fit
lm_coef = sapply(models, function(rhs) {
  Pitching[ , coef(lm(ERA ~ ., data = .SD))['W'], .SDcols = c('W', rhs)]
})



## coefficients of 'W' from the linear regression models above
lm_coef
# W           W           W           W           W           W           W           W 
# -0.20643825 -0.20877606 -0.20763604 -0.11362170 -0.18675857 -0.20689003 -0.09819310 -0.18763789 
# W           W           W           W           W           W           W           W 
# -0.10765791 -0.18214290 -0.12478923 -0.09595359 -0.18197728 -0.11646014 -0.11793832 -0.11143365

 

 

위의 sapply() 함수와 .SD, .SDcols 를 사용하지 않고 각 설명변수 16개 조합별로 일일이 수작업으로 회귀모형을 나열해서 적으려면 아래처럼 단순 작업을 여러번 해야 합니다. 시간도 오래걸리고, 또 오타 에러를 낼 위험이 높습니다.

 

## all combinations of linear regression models by manually
Pitching[ , coef(lm(ERA ~ ., data = .SD))['W'], .SDcols = c('W')]
Pitching[ , coef(lm(ERA ~ ., data = .SD))['W'], .SDcols = c('W', 'yearID')]
Pitching[ , coef(lm(ERA ~ ., data = .SD))['W'], .SDcols = c('W', 'teamID')]
Pitching[ , coef(lm(ERA ~ ., data = .SD))['W'], .SDcols = c('W', 'G')]
Pitching[ , coef(lm(ERA ~ ., data = .SD))['W'], .SDcols = c('W', 'L')]
Pitching[ , coef(lm(ERA ~ ., data = .SD))['W'], .SDcols = c('W', 'yearID', 'teamID')]
Pitching[ , coef(lm(ERA ~ ., data = .SD))['W'], .SDcols = c('W', 'yearID', 'G')]
Pitching[ , coef(lm(ERA ~ ., data = .SD))['W'], .SDcols = c('W', 'yearID', 'L')]
Pitching[ , coef(lm(ERA ~ ., data = .SD))['W'], .SDcols = c('W', 'teamID', 'G')]
Pitching[ , coef(lm(ERA ~ ., data = .SD))['W'], .SDcols = c('W', 'teamID', 'L')]
Pitching[ , coef(lm(ERA ~ ., data = .SD))['W'], .SDcols = c('W', 'G', 'L')]
Pitching[ , coef(lm(ERA ~ ., data = .SD))['W'], .SDcols = c('W', 'yearID', 'teamID', 'G')]
Pitching[ , coef(lm(ERA ~ ., data = .SD))['W'], .SDcols = c('W', 'yearID', 'teamID', 'L')]
Pitching[ , coef(lm(ERA ~ ., data = .SD))['W'], .SDcols = c('W', 'yearID', 'G', 'L')]
Pitching[ , coef(lm(ERA ~ ., data = .SD))['W'], .SDcols = c('W', 'teamID', 'G', 'L')]
Pitching[ , coef(lm(ERA ~ ., data = .SD))['W'], .SDcols = c('W', 'yearID', 'teamID', 'G', 'L')]

 

 

(3) 설명변수 조합별 회귀계수로 막대그래프 그리기: barplot()

 

위의 (2)번에서 구한 각 설명변수 16개 조합별 변수 'W'의 회귀계수를 가지고 barplot() 함수를 사용해서 막대그래프(bar plot)를 그려보겠습니다. 

여기서도 'W' 회귀계수가 적합되었을 때 각 모델별로 사용되었던 설명변수 조합을 names.arg = sapply(models, paste, collapse = '/') 로 해서 넣어주었습니다.

 

## -- barplot
# here are 16 visually distinct colors, taken from the list of 20 here:
#   https://sashat.me/2017/01/11/list-of-20-simple-distinct-colors/
col16 = c('#e6194b', '#3cb44b', '#ffe119', '#0082c8',
          '#f58231', '#911eb4', '#46f0f0', '#f032e6',
          '#d2f53c', '#fabebe', '#008080', '#e6beff',
          '#aa6e28', '#fffac8', '#800000', '#aaffc3')

par(oma = c(2, 0, 0, 0))

barplot(lm_coef, names.arg = sapply(models, paste, collapse = '/'), 
        main = "Wins coefficient\n With various covariates", 
        col = col16, las = 2L, cex.names = .8)
        

 

 

참고로, paste() 함수로 'models' 리스트에 들어있는 변수 이름들의 조합을 하나의 합칠 때, collpase = '/' 를 설정해주면 아래처럼 각 리스트 원소 내 여러개의 변수 이름을 하나로 합치면서 사이에 '/'를 넣어줍니다. (paste 함수에서 sep와 collapse 의 차이점은 여기를 참고하세요.)

sapply(models, paste, collapse = '/')
# [1] ""                  "yearID"            "teamID"            "G"                 "L"                
# [6] "yearID/teamID"     "yearID/G"          "yearID/L"          "teamID/G"          "teamID/L"         
# [11] "G/L"               "yearID/teamID/G"   "yearID/teamID/L"   "yearID/G/L"        "teamID/G/L"       
# [16] "yearID/teamID/G/L"

 

 

[ Reference ]

* R data.table vignettes 'Using .SD for Data Analysis'

  : cran.r-project.org/web/packages/data.table/vignettes/datatable-sd-usage.html

 

 

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

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

 

 

 

728x90
반응형
Posted by Rfriend
,

Rfriend님의
글이 좋았다면 응원을 보내주세요!

R data.table 패키지를 사용하다 보면 base R이나 tydeverse 에서는 사용하지 않는, data.table에서만 사용하는 .SD, .SDcols, .N 등의 매개변수를 볼 수 있는데요, 이게 data.table 패키지를 굉장히 이질적이고 코드가 이해가 되지 않는 어려운 프로그래밍 언어라는 첫인상을 주는 것 같습니다.  사실 이 첫번째 관문만 무사히 넘으면 data.table 의 강력함과 간결함에 매료될만도 한데 말이지요.

 

이번 포스팅에서는 data.table의 vignettes을 참조하여 R의 data.table 패키지에서

 

(1) .SD는 무엇인가? (What is .SD in data.table?)

(2) .SDcols 로 일부 칼럼만 가져오기 (Column Subsetting using .SDcols)

(3) lapply()와 .SDcols로 칼럼 유형 변환하기 (Column Type Conversion)

(4) 패턴이 일치하는 특정 칼럼만 가져오기 (Column subsetting using pattern-based matching)

 

에 대해서 소개하겠습니다.

 

 

(1) .SD는 무엇인가? (What is .SD in data.table package?)

 

data.table 패키지에서 .SD 는 "데이터의 부분집합, 데이터 자기 자신, 데이터의 자기 참조" ( 'Subset, Selfsame, or Self-reference of the Data')를 의미한다고 이해할 수 있습니다.  다시 말하자면, .SD는 data.table 자기 자신에 대한 재귀적인 참조 (a reflexive reference to the data.table), data.table 그 자체 (data.table itself) 를 의미합니다.

 

 

R의 Lahman database에 있는 야구 팀과 투수의 통계 데이터를 사용해서 예를 들어보겠습니다.  setDT(Pitching) 은 Lists와 DataFrame을 참조해서 Data.Table로 만들어줍니다.

 

먼저, 아래의 예는 .SD 를 사용해서 'Pitching' data.table 자체를 재귀적으로 참조해오는 예입니다. (복사가 아니라 참조임. not copy, but referece to the Pitching data.table)

 

library(data.table)

## .SD stands for Subset, Selfsame, or Self-reference of the Data.
## .SD is a reflexive reference to the data.table itself. 
## .SD is helpful for chaining together "queries" (extractions/subsets/etc using [). 


## Lahman database on baseball
install.packages("Lahman")
library(Lahman)

data("Pitching")
setDT(Pitching)
str(Pitching)
# Classes 'data.table' and 'data.frame':	47628 obs. of  30 variables:
#   $ playerID: chr  "bechtge01" "brainas01" "fergubo01" "fishech01" ...
# $ yearID  : int  1871 1871 1871 1871 1871 1871 1871 1871 1871 1871 ...
# $ stint   : int  1 1 1 1 1 1 1 1 1 1 ...
# $ teamID  : Factor w/ 149 levels "ALT","ANA","ARI",..: 97 142 90 111 90 136 111 56 97 136 ...
# $ lgID    : Factor w/ 7 levels "AA","AL","FL",..: 4 4 4 4 4 4 4 4 4 4 ...
# $ W       : int  1 12 0 4 0 0 0 6 18 12 ...
# $ L       : int  2 15 0 16 1 0 1 11 5 15 ...
# $ G       : int  3 30 1 24 1 1 3 19 25 29 ...
# $ GS      : int  3 30 0 24 1 0 1 19 25 29 ...
# $ CG      : int  2 30 0 22 1 0 1 19 25 28 ...
# $ SHO     : int  0 0 0 1 0 0 0 1 0 0 ...
# $ SV      : int  0 0 0 0 0 0 0 0 0 0 ...
# $ IPouts  : int  78 792 3 639 27 3 39 507 666 747 ...
# $ H       : int  43 361 8 295 20 1 20 261 285 430 ...
# $ ER      : int  23 132 3 103 10 0 5 97 113 153 ...
# $ HR      : int  0 4 0 3 0 0 0 5 3 4 ...
# $ BB      : int  11 37 0 31 3 0 3 21 40 75 ...
# $ SO      : int  1 13 0 15 0 0 1 17 15 12 ...
# $ BAOpp   : num  NA NA NA NA NA NA NA NA NA NA ...
# $ ERA     : num  7.96 4.5 27 4.35 10 0 3.46 5.17 4.58 5.53 ...
# $ IBB     : int  NA NA NA NA NA NA NA NA NA NA ...
# $ WP      : int  7 7 2 20 0 0 1 15 3 44 ...
# $ HBP     : int  NA NA NA NA NA NA NA NA NA NA ...
# $ BK      : int  0 0 0 0 0 0 0 2 0 0 ...
# $ BFP     : int  146 1291 14 1080 57 3 70 876 1059 1334 ...
# $ GF      : int  0 0 0 1 0 1 1 0 0 0 ...
# $ R       : int  42 292 9 257 21 0 30 243 223 362 ...
# $ SH      : int  NA NA NA NA NA NA NA NA NA NA ...
# $ SF      : int  NA NA NA NA NA NA NA NA NA NA ...
# $ GIDP    : int  NA NA NA NA NA NA NA NA NA NA ...
# - attr(*, ".internal.selfref")=<externalptr>



## .SD on Ungrouped Data
## In terms of subsetting, .SD is a subset of the data, the set itself. 

Pitching[, .SD]
# playerID yearID stint teamID lgID  W  L  G GS CG SHO SV IPouts   H  ER HR BB SO BAOpp
# 1: bechtge01   1871     1    PH1   NA  1  2  3  3  2   0  0     78  43  23  0 11  1    NA
# 2: brainas01   1871     1    WS3   NA 12 15 30 30 30   0  0    792 361 132  4 37 13    NA
# 3: fergubo01   1871     1    NY2   NA  0  0  1  0  0   0  0      3   8   3  0  0  0    NA
# 4: fishech01   1871     1    RC1   NA  4 16 24 24 22   1  0    639 295 103  3 31 15    NA
# 5: fleetfr01   1871     1    NY2   NA  0  1  1  1  1   0  0     27  20  10  0  3  0    NA
# ---                                                                                       
#   47624: zamorda01   2019     1    NYN   NL  0  1 17  0  0   0  0     26  10   5  1  5  8 0.294
# 47625: zeuchtj01   2019     1    TOR   AL  1  2  5  3  0   0  0     68  22  12  2 11 20 0.250
# 47626: zimmejo02   2019     1    DET   AL  1 13 23 23  0   0  0    336 145  86 19 25 82 0.311
# 47627: zimmeky01   2019     1    KCA   AL  0  1 15  0  0   0  0     55  28  22  2 19 18 0.337
# 47628: zobribe01   2019     1    CHN   NL  0  0  1  0  0   0  0      3   0   0  0  2  1 0.000
# ERA IBB WP HBP BK  BFP GF   R SH SF GIDP
# 1:  7.96  NA  7  NA  0  146  0  42 NA NA   NA
# 2:  4.50  NA  7  NA  0 1291  0 292 NA NA   NA
# 3: 27.00  NA  2  NA  0   14  0   9 NA NA   NA
# 4:  4.35  NA 20  NA  0 1080  1 257 NA NA   NA
# 5: 10.00  NA  0  NA  0   57  0  21 NA NA   NA
# ---                                           
#   47624:  5.19   1  0   1  0   41  3   5  0  1    0
# 47625:  4.76   0  2   0  0   99  0  13  0  0    1
# 47626:  6.91   2  3   6  0  504  0  89  3  4    5
# 47627: 10.80   0  2   0  0  102  3  22  0  0    1
# 47628:  0.00   0  0   0  0    5  1   0  0  0    0

 

 

Pitching[, .SD] 는 단순히 Pitching data.table 자체를 그대로 반환하는데요, identical() 로 두 data.table이 동일한지 여부를 확인해보면 TRUE 입니다.

 

identical(Pitching, Pitching[ , .SD])
# [1] TRUE

 

 

(2) .SDcols 로 일부 칼럼만 가져오기 (Column Subsetting using .SDcols)

 

.SD 를 사용하는 가장 일반적인 예는 .SDcols와 함께 일부 칼럼의 부분집합을 선택해서 가져오는 것입니다. 

 

아래 예에서는 Pitching[, .SD] 로 Pitching data.table 자체를 재귀적으로 참조한 후에 --> Pitching[, .SD, .SDcols = c("playerID", "W", "L", "G")] 처럼 .SDcols 에 칼럼 이름을 넣어줘서 일부 칼럼만 선택적으로 가져와 보겠습니다.

 

library(data.table)

library(Lahman)
data("Pitching")

## Coerce lists and Data.Frames to Data.Table by Reference
setDT(Pitching)


## ** Column Subsetting: .SDcols **
Pitching[, .SD, .SDcols = c("playerID", "W", "L", "G")]

# playerID  W  L  G
# 1: bechtge01  1  2  3
# 2: brainas01 12 15 30
# 3: fergubo01  0  0  1
# 4: fishech01  4 16 24
# 5: fleetfr01  0  1  1
# ---                   
#   47624: zamorda01  0  1 17
# 47625: zeuchtj01  1  2  5
# 47626: zimmejo02  1 13 23
# 47627: zimmeky01  0  1 15
# 47628: zobribe01  0  0  1



 

(3) 칼럼 유형 변환하기 (Column Type Conversion)

 

이번에는 lapply()와 .SDcols를 사용해서 data.table에서 여러개 칼럼의 유형을 한꺼번에 변환해보겠습니다.

 

먼저, setDT(Teams) 로 'Teams' DataFrame을 참조하여 Data.Table로 만들고, c('teamIDBR', 'teamIDlahman45', 'teamIDretro') 의 3개 칼럼이 '문자형 인지 여부(is.character)'를 확인해보겠습니다. 3개 칼럼 모두 is.character()가 TRUE 이므로 문자형 맞군요.

 

library(data.table)

## Lahman database on baseball
install.packages("Lahman")
library(Lahman)
data(Teams)

## Coerce lists and Data.Frames to Data.Table by Reference
setDT(Teams)

str(Teams)
# Classes 'data.table' and 'data.frame':	2925 obs. of  48 variables:
#   $ yearID        : int  1871 1871 1871 1871 1871 1871 1871 1871 1871 1872 ...
# $ lgID          : Factor w/ 7 levels "AA","AL","FL",..: 4 4 4 4 4 4 4 4 4 4 ...
# $ teamID        : Factor w/ 149 levels "ALT","ANA","ARI",..: 24 31 39 56 90 97 111 136 142 8 ...
# $ franchID      : Factor w/ 120 levels "ALT","ANA","ARI",..: 13 36 25 56 70 85 91 109 77 9 ...
# $ divID         : chr  NA NA NA NA ...
# $ Rank          : int  3 2 8 7 5 1 9 6 4 2 ...
# $ G             : int  31 28 29 19 33 28 25 29 32 58 ...
# $ Ghome         : int  NA NA NA NA NA NA NA NA NA NA ...
# $ W             : int  20 19 10 7 16 21 4 13 15 35 ...
# $ L             : int  10 9 19 12 17 7 21 15 15 19 ...
# $ DivWin        : chr  NA NA NA NA ...
# $ WCWin         : chr  NA NA NA NA ...
# $ LgWin         : chr  "N" "N" "N" "N" ...
# $ WSWin         : chr  NA NA NA NA ...
# $ R             : int  401 302 249 137 302 376 231 351 310 617 ...
# $ AB            : int  1372 1196 1186 746 1404 1281 1036 1248 1353 2571 ...
# $ H             : int  426 323 328 178 403 410 274 384 375 753 ...
# $ X2B           : int  70 52 35 19 43 66 44 51 54 106 ...
# $ X3B           : int  37 21 40 8 21 27 25 34 26 31 ...
# $ HR            : int  3 10 7 2 1 9 3 6 6 14 ...
# $ BB            : int  60 60 26 33 33 46 38 49 48 29 ...
# $ SO            : int  19 22 25 9 15 23 30 19 13 28 ...
# $ SB            : int  73 69 18 16 46 56 53 62 48 53 ...
# $ CS            : int  16 21 8 4 15 12 10 24 13 18 ...
# $ HBP           : int  NA NA NA NA NA NA NA NA NA NA ...
# $ SF            : int  NA NA NA NA NA NA NA NA NA NA ...
# $ RA            : int  303 241 341 243 313 266 287 362 303 434 ...
# $ ER            : int  109 77 116 97 121 137 108 153 137 166 ...
# $ ERA           : num  3.55 2.76 4.11 5.17 3.72 4.95 4.3 5.51 4.37 2.9 ...
# $ CG            : int  22 25 23 19 32 27 23 28 32 48 ...
# $ SHO           : int  1 0 0 1 1 0 1 0 0 1 ...
# $ SV            : int  3 1 0 0 0 0 0 0 0 1 ...
# $ IPouts        : int  828 753 762 507 879 747 678 750 846 1548 ...
# $ HA            : int  367 308 346 261 373 329 315 431 371 573 ...
# $ HRA           : int  2 6 13 5 7 3 3 4 4 3 ...
# $ BBA           : int  42 28 53 21 42 53 34 75 45 63 ...
# $ SOA           : int  23 22 34 17 22 16 16 12 13 77 ...
# $ E             : int  243 229 234 163 235 194 220 198 218 432 ...
# $ DP            : int  24 16 15 8 14 13 14 22 20 22 ...
# $ FP            : num  0.834 0.829 0.818 0.803 0.84 0.845 0.821 0.845 0.85 0.83 ...
# $ name          : chr  "Boston Red Stockings" "Chicago White Stockings" "Cleveland Forest Citys" "Fort Wayne Kekiongas" ...
# $ park          : chr  "South End Grounds I" "Union Base-Ball Grounds" "National Association Grounds" "Hamilton Field" ...
# $ attendance    : int  NA NA NA NA NA NA NA NA NA NA ...
# $ BPF           : int  103 104 96 101 90 102 97 101 94 106 ...
# $ PPF           : int  98 102 100 107 88 98 99 100 98 102 ...
# $ teamIDBR      : chr  "BOS" "CHI" "CLE" "KEK" ...
# $ teamIDlahman45: chr  "BS1" "CH1" "CL1" "FW1" ...
# $ teamIDretro   : chr  "BS1" "CH1" "CL1" "FW1" ...
# - attr(*, ".internal.selfref")=<externalptr> 


## check whether columns are character type or not
col_lists = c('teamIDBR', 'teamIDlahman45', 'teamIDretro')
Teams[, sapply(.SD, is.character), .SDcols = col_lists]
# teamIDBR teamIDlahman45    teamIDretro 
# TRUE           TRUE           TRUE


head(unique(Teams[[col_lists[1L]]]))
# [1] "BOS" "CHI" "CLE" "KEK" "NYU" "ATH"

 

 

이제 위의 3개 문자형 칼럼들을 요인형(factor type)으로 lapply()와 .SDcols 를 사용해서 한꺼번에 데이터 유형을 변환해보겠습니다.

 

이때 주의를 해야 할 것이 있는데요, 아래에 색깔을 칠해 놓은 것처럼, (col_lists) 처럼 괄호 () 로 싸주어야 col_lists 안의 3개 칼럼 이름을 인식해서 할당을 해줍니다. (괄호를 안쳐주면 'col_lists' 라는 이름으로 할당해버립니다.)

 

Teams[ , (col_lists) := lapply(.SD, factor), .SDcols = col_lists]

 

## Converting columns to factor by adding the := assignment operator.
## we must wrap fkt in parentheses () to force data.table to interprete this as column names.
Teams[ , (col_lists) := lapply(.SD, factor), .SDcols = col_lists]

 

col_lists의 3개 칼럼의 데이터 유형이 무엇인지 확인해보면, 요인형(factor)로 잘 변환이 되었네요.

sapply(Teams[, .SD, .SDcols = col_lists], class)
# teamIDBR teamIDlahman45    teamIDretro 
# "factor"       "factor"       "factor"



head(unique(Teams[[col_lists[1L]]]))
# [1] BOS CHI CLE KEK NYU ATH
# 101 Levels: ALT ANA ARI ATH ATL BAL BLA BLN BLU BOS BRA BRG BRO BSN BTT BUF BWW CAL CEN CHC ... WSN

 

 

(4) 패턴이 일치하는 특정 칼럼만 가져오기 (Column subsetting using pattern-based matching)

 

data.table의 .SDcols는 패턴 매칭을 지원합니다. 아래의 예에서는 'Teams' Data.Table에서 'team' 이라는 단어가 들어가 있는 칼럼 이름(.SDcols = patterns('team'))을 선별해서 가져와 보겠습니다.

 

names(Teams)

# [1] "yearID"         "lgID"           "teamID"         "franchID"       "divID"         
# [6] "Rank"           "G"              "Ghome"          "W"              "L"             
# [11] "DivWin"         "WCWin"          "LgWin"          "WSWin"          "R"             
# [16] "AB"             "H"              "X2B"            "X3B"            "HR"            
# [21] "BB"             "SO"             "SB"             "CS"             "HBP"           
# [26] "SF"             "RA"             "ER"             "ERA"            "CG"            
# [31] "SHO"            "SV"             "IPouts"         "HA"             "HRA"           
# [36] "BBA"            "SOA"            "E"              "DP"             "FP"            
# [41] "name"           "park"           "attendance"     "BPF"            "PPF"           
# [46] "teamIDBR"       "teamIDlahman45" "teamIDretro" 



## pattern-based matching of columns in .SDcols to select all columns 
## which contain team back to factor. 
Teams[ , .SD, .SDcols = patterns('team')]

# teamID teamIDBR teamIDlahman45 teamIDretro
# 1:    BS1      BOS            BS1         BS1
# 2:    CH1      CHI            CH1         CH1
# 3:    CL1      CLE            CL1         CL1
# 4:    FW1      KEK            FW1         FW1
# 5:    NY2      NYU            NY2         NY2
# ---                                           
#   2921:    SLN      STL            SLN         SLN
# 2922:    TBA      TBR            TBA         TBA
# 2923:    TEX      TEX            TEX         TEX
# 2924:    TOR      TOR            TOR         TOR
# 2925:    WAS      WSN            MON         WAS

 

아래의 예는 team_idx = grep('team', names(Teams), value=TRUE)로, 먼저 names(Teams)를 통해 얻은 전체 칼럼 이름들 중에서 'team'이 들어가 있는 칼럼 이름을 찾아서 team_idx 에 저장을 해주었습니다.  다음으로, 위의 (3)번에서 했던 lapply()와 .SD, .SDcols를 이용하여 칼럼 이름에 'team'이 들어간 모든 칼럼의 데이터 유형을 요인형(factor type)으로 일괄 변환해주었습니다.

 

team_idx = grep('team', names(Teams), value = TRUE)
team_idx

# [1] "teamID"         "teamIDBR"       "teamIDlahman45" "teamIDretro"



Teams[ , (team_idx) := lapply(.SD, factor), .SDcols = team_idx]
sapply(Teams[, .SD, .SDcols = team_idx], class)

# teamID       teamIDBR teamIDlahman45    teamIDretro 
# "factor"       "factor"       "factor"       "factor"

 

[ Reference ]

* R data.table vignettes 'Using .SD for Data Analysis'

  : cran.r-project.org/web/packages/data.table/vignettes/datatable-sd-usage.html

 

 

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

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

 

728x90
반응형
Posted by Rfriend
,

Rfriend님의
글이 좋았다면 응원을 보내주세요!

지난번 포스팅에서는 R dplyr 패키지의 group_by(), mutate(), filter(), select() 함수를 사용해서 단위 그룹별 관측치 개수별로 DataFrame을 구분해서 생성하는 방법을 소개(https://rfriend.tistory.com/606)하였습니다.


이번 포스팅에서는 R data.table 패키지를 사용해서 동일하게 단위 그룹별 관측치 개수별로 DataTable을 구분해서 생성하는 방법을 소개하겠습니다.


(1) R data.table을 사용해서 그룹 별 관측치 개수 세기

(2) R data.table을 사용해서 그룹 별 관측치 개수 별로 data.table 구분해서 생성하기





만약 데이터셋이 대용량이라면 R dplyr 패키지를 사용하는 것보다 R data.table 패키지를 사용하는 것이 속도 면이나 메모리 효율성 면에서 유리합니다.



먼저 data.table 패키지를 불러오고, 예제로 사용할 'id' 변수를 관측치 개수를 세는 단위 기준으로 삼는 간단한 data.table을 만들어보겠습니다.



library(data.table)


id <- c("A", "B", "C", "C", "D", "E", "E", "E")
x1 <- c(3, 2, 1, 1, 3, 3, 2, 4)
x2 <- c(60, 20, 30, 10, 70, 10, 20, 30)

DT <- data.table(id, x1, x2)
DT
# id x1 x2
# 1:  A  3 60
# 2:  B  2 20
# 3:  C  1 30
# 4:  C  1 10
# 5:  D  3 70
# 6:  E  3 10
# 7:  E  2 20
# 8:  E  4 30





  (1) R data.table을 사용해서 그룹 별 관측치 개수 세기


data.table 패키지에서 관측치 개수를 세는 것은 '.N' 매개변수를 사용하며, 그룹별 연산은 'by' 매개변수를 사용합니다. 그리고 새로운 변수를 만들어서 이름을 부여할 때는  'new_col_name := operation_method' 처럼 ':=' 로 할당을 해주면 됩니다.  


'id' 기준 그룹별로 관측치 개수를 세어서 'n' 이라는 이름으로 새로운 칼럼을 생성해보았습니다. id가 'A', 'B',  D' 는 관측치가 각 1개씩이며, 'C'는 관측치가 2개, 'E'는 관측치가 3개임을 알 수 있습니다.



## number of rows by 'id'
DT[, n := .N, by = id]
DT
# id x1 x2 n
# 1:  A  3 60 1
# 2:  B  2 20 1
# 3:  C  1 30 2
# 4:  C  1 10 2
# 5:  D  3 70 1
# 6:  E  3 10 3
# 7:  E  2 20 3
# 8:  E  4 30 3





  (2) R data.table을 사용해서 그룹 별 관측치 개수 별로 data.table 구분해서 생성하기


data.table 패키지도 dplyr 처럼 코드를 연속해서 이어서 쓸 수 있습니다. 위의 (1)번에서 구한 그룹별 관측치 개수를 조건으로 해서, 관측치 개수가 1개인 데이터셋([n == 1])과, 그룹별 관측치 개수가 2개 이상인 데이터셋(n >=2])을 조건절을 이어써줘서 선별해보겠습니다.



DT[, n:=.N, by = id][n == 1,]
# id x1 x2 n
# 1:  A  3 60 1
# 2:  B  2 20 1
# 3:  D  3 70 1

DT[, n:=.N, by = id][n >= 2,]
# id x1 x2 n
# 1:  C  1 30 2
# 2:  C  1 10 2
# 3:  E  3 10 3
# 4:  E  2 20 3
# 5:  E  4 30 3




그룹별 관측치 개수를 조건절로 사용해서 원하는 그룹별 관측치별 개수 1개, 2개 이상 관측치 데이터셋을 구분하고 나면, 이제 원래 데이터셋에 있었던 칼럼인 c("id", "x1", "x2") 만 칼럼을 선택해서 가져오고 나머지는 버리도록 하겠습니다.



DT_1 = DT[, n := .N, by = id][n == 1,][, c("id", "x1", "x2")]
DT_1
# id x1 x2
# 1:  A  3 60
# 2:  B  2 20
# 3:  D  3 70

DT_2 = DT[, n := .N, by = id][n >= 2,][,c("id", "x1", "x2")]
DT_2
# id x1 x2
# 1:  C  1 30
# 2:  C  1 10
# 3:  E  3 10
# 4:  E  2 20
# 5:  E  4 30




이전의 dplyr 보다 data.table 패키지의 코드가 처음에는 낮설어 보일 수 있기는 합니다만, 코드의 길이나 표현이 좀더 간결하고 속도/성능도 우수하므로 사용해보시길 권해드립니다.



[Reference]

* R data.table 의 DT[i, j, by]에서 행 subset, 열 select 하고 계산하기: https://rfriend.tistory.com/566

* R data.table의 by 구문으로 그룹별 집계하기: https://rfriend.tistory.com/567

* R data.table의 Reference Semantics: https://rfriend.tistory.com/573



많은 도움이 되었기를 바랍니다.
행복한 데이터 과학자 되세요!  :-)



728x90
반응형
Posted by Rfriend
,

Rfriend님의
글이 좋았다면 응원을 보내주세요!

이번 포스팅에서는 데이터 전처리 과정 중에서


(1) 기준이 되는 단위나 그룹 별로 관측치의 개수를 구해서

(2) 그룹 별 관측치 개수가 1개인 그룹의 DataFrame 과 2개 이상인 그룹의 DataFrame을 구분


해서 생성하는 방법을 소개하겠습니다.


dplyr 패키지를 이용해서 chain operator (%>%) 로 한꺼번에 코드를 짜면 편리합니다. dplyr 패키지의 group_by(), mutate(), filter(), select() 등의 여러개의 함수를 사용하였습니다.


(그룹 별 관측치 개수를 구해서 먼저 DataFrame으로 만들어놓고, 그 다음에 merge 함수를 이용해서 그룹별로 관측치 개수를 원본 DataFrame에 합쳐준 후에, 조건문으로 그룹별 관측치 개수별로 DataFrame을 구분해서 생성해주는 방식으로 step by step 진행해도 됩니다.)





먼저, dplyr 패키지를 불러오고, 간단한 예제 DataFrame을 만들어보겠습니다. 이때 칼럼 'id'가 관측치 개수를 세는 기준이 되는 그룹이 되겠습니다.



library(dplyr)

id <- c("A", "B", "C", "C", "D", "E", "E", "E")
x1 <- c(3, 2, 1, 1, 3, 3, 2, 4)
x2 <- c(60, 20, 30, 10, 70, 10, 20, 30)

df <- data.frame(id, x1 x2)
df
# id x1 x2
# 1  A  3 60
# 2  B  2 20
# 3  C  1 30
# 4  C  1 10
# 5  D  3 70
# 6  E  3 10
# 7  E  2 20
# 8  E  4 30

 




  (1) 단위 그룹별로 관측치 개수 구하기


(a) dplyr 패키지의 group_by() 함수에 관측치 개수를 세는 단위 기준이 되는 그룹변수 'id' 를 넣어주고,


(b) mutate() 함수의 n() 매개변수로 관측치 개수(number of rows)를 세어서 'n'이라는 이름의 새로운 변수를 만들어줍니다.


(c) 그리고 filter() 함수를 사용해서 위의 (b)에서 새로 만든 그룹별 관측치 개수 'n'에 대해 조건 '관측치 개수가 1개인 조건 (n == 1)' 을 만족하는 행만 걸러내는 원리입니다.



df %>%
  group_by(id) %>%
  mutate(n = n()) %>%
  filter(n == 1)


# # A tibble: 3 x 4
# # Groups:   id [3]
# id       x1    x2     n
# <chr> <dbl> <dbl> <int>
# 1 A         3    60     1
# 2 B         2    20     1
# 3 D         3    70     1






  (2) 단위 그룹별로 관측치 개수가 1개 vs. 2개 이상인 그룹 구분해서 DataFrame 만들기


우리가 원하는 최종 산출물은 단위 그룹별 관측치 개수는 필요하지 않으므로 (d) select(id, x1, x2) 함수를 사용해서 원본 데이터에 있는 변수들만 선택해서 가져오겠습니다.


그리고 그룹별 관측치 개수가 1개인 관측치만으로 이루어진 'df_1' 라는 이름의 DataFrame과, 그룹별 관측치 개수가 2개 이상인 관측치들로 이루어진 'df_2' 라는 이름의 DataFrame을 각각 구분해서 생성해보겠습니다.


dplyr 패키지의 체인 연산(%>%) 으로 한꺼번에 코드를 짤 수 있으므로 코드가 단순하고 가독성이 좋습니다.



## data.frame which has only 1 observation by x1 category
df_1 <- df %>%
  group_by(id) %>%
  mutate(n = n()) %>%
  filter(n == 1) %>%
  select(id, x1, x2)

df_1
# id       x1    x2
# <chr> <dbl> <dbl>
# 1 A         3    60
# 2 B         2    20
# 3 D         3    70


# data.frame which has more than 2 observations by x1 categories
df_2 <- df %>%
  group_by(id) %>%
  mutate(n = n()) %>%
  filter(n >= 2) %>%
  select(id, x1, x2)

df_2
# id       x1    x2
# <chr> <dbl> <dbl>
# 1 C         1    30
# 2 C         1    10
# 3 E         3    10
# 4 E         2    20
# 5 E         4    30



[ Reference ]

* R dplyr 패키지 기본함수 소개: https://rfriend.tistory.com/234

* R dplyr 패키지의 새로운 변수만들기 소개: https://rfriend.tistory.com/235

* R dplyr 패키지의 chain operator (or pipe operator) 소개: https://rfriend.tistory.com/236

* R dplyr 패키지의 그룹별 행의 개수 세기 소개: https://rfriend.tistory.com/240



다음 포스팅에서는 이번 포스팅과 동일한 결과를 얻기위해 R data.table 패키지를 사용해서 그룹별 관측치 개수별로 data.table을 구분해서 생성하는 방법(https://rfriend.tistory.com/607)을 소개하겠습니다.


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

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



728x90
반응형
Posted by Rfriend
,

Rfriend님의
글이 좋았다면 응원을 보내주세요!

지리적 레스터 데이터 (geographic raster data) 는 일반적으로 (a) 레스터 헤더(raster header)(b) 행렬 (matrix)로 구성이 됩니다.

- (a) 레스터 헤더 (raster header): 좌표 참조 시스템(CRS, Coordinate Reference System), 시작점(the origin)과 범위 (the extent)를 정의함
- (b) 행렬 (matrix): 동일한 크기의 픽셀 또는 셀(pixel, or cell)을 표현. 픽셀 ID(pixel IDs)와 픽셀 값(pixel values).


보통 행렬의 시작점(the origin, or starting point)은 행렬의 좌측 하단 구석에 위치한 좌표를 의미하는데요, R의 raster 패키지의 레스터 헤더에서는 시작점의 기본값으로 좌측 상단 구석에 위치한 좌표를 시작점으로 사용하므로 주의가 필요합니다.

레스터 헤더에서 범위 (the extent) sms 행의 개수, 열의 개수, 셀의 크기 해상도로 정의합니다. 각각의 셀에 접근하거나 수정하려면 시작점(the origin)으로 부터의 셀 ID를 사용하거나 또는 명시적으로 행과 열을 지정하면 됩니다.

레스터 데이터의 행렬 표현법은 각 셀의 네 개 구석의 좌표를 명시적으로 저장하지 않으며, 대신 시작점(the origin)만 저장하고 나머지는 시작점으로부터의 행과 열의 ID를 가지고 각 셀에 접근하는 방식이므로 벡터 데이터 표현과 비교해서 상대적으로 효율적이고 속도가 빠릅니다. (예: 벡터 사각형 폴리곤의 경우 각 셀별로 5개 점의 좌표를 저장해야 함.) 하지만 벡터 데이터의 각 도형별로 여러개의 값을 가질 수 있는 반면에, 레스터 데이터의 경우 각 셀별로 하나의 값만을 가질 수 있습니다.

 

 

[레스터 데이터 유형 (1) 셀 ID (Cell IDs), (2) 셀 값 (Cell Values) ]

 

 

지리공간 벡터 데이터 처리 및 분석에 sf 패키지를 사용했었는데요, 레스터 데이터 처리 및 분석은 R raster 패키지를 사용합니다. 

spDataLarge 패키지에 내장되어 있는 미국 유타지역의 Zion 국립공원 지역의 레스터 샘플 데이터를 raster() 함수를 사용해서 불러오고, 레스터 데이터의 속성정보들을 살펴보겠습니다. 혹시 raster, spDataLarge, rgdal 패키지를 사전에 설치하지 않았다면 install.packages() 로 패키지를 먼저 설치해주세요.

raster() 함수로 불러온 레스터 데이터셋 이름 'srtm_raster' 을 입력해주면 클래스(class) 차원(dimentions), 해상도(resolution), 범위(extent), 좌표 참조 시스템 (Coordinates Reference System), 출처(Source), 이름(names), 최소/최대 값(min, max values) 속성 정보를 볼 수 있습니다.

 

# =================================
# R GeoSpatial Data Analysis 
# Raster Data using raster package
# =================================

library(raster)
library(spDataLarge)

install.packages("rgdal")
library(rgdal)

## raster dataset from the spDataLarge package, 
## a few raster objects and one vector object covering an area of the Zion National Park (Utah, USA).
raster_filepath = system.file("raster/srtm.tif", package = "spDataLarge")
srtm_raster = raster(raster_filepath)
                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                          
srtm_raster
# class      : RasterLayer 
# dimensions : 457, 465, 212505  (nrow, ncol, ncell)
# resolution : 0.0008333333, 0.0008333333  (x, y)
# extent     : -113.2396, -112.8521, 37.13208, 37.51292  (xmin, xmax, ymin, ymax)
# crs        : +proj=longlat +datum=WGS84 +no_defs 
# source     : /Library/Frameworks/R.framework/Versions/4.0/Resources/library/spDataLarge/raster/srtm.tif 
# names      : srtm 
# values     : 1024, 2892  (min, max)

 

 

이들 레스터 데이테의 속성 정보에 함수를 사용해서 일일이 접근이 가능합니다. 먼저, 좌표 참조 시스템(Coordinate Reference System, CRS)는 crs() 함수를 통해 확인할 수 있습니다.

 

## crs() : coordinate reference system
crs(srtm_raster)
# CRS arguments: +proj=longlat +datum=WGS84 +no_defs

 

 

dim() 함수는 레스터 데이터의 차원인 행의 개수(number of rows), 열의 개수 (number of columns), 층의 개수 (number of layers)를 한꺼번에 볼 수 있게 해줍니다. 

ncol() 은 열의 개수, nrow() 는 행의 개수, nlayers() 는 층의 개수를 개별적으로 접근할 수 있는 함수입니다.

 

## dim() : number of rows, columns and layers
dim(srtm_raster)
# [1] 457 465   1


nrow(srtm_raster)
# [1] 457

ncol(srtm_raster)
# [1] 465

nlayers(srtm_raster)
# [1] 1

 

 

ncell() 함수는 레스터 데이터의 전체 셀의 개수 (number of cells, or pixels) 를 알려주며, 이는 (셀의 개수 = 행의 개수 X 열의 개수) 를 의미합니다.

res() 함수는 레스터 데이터의 공간 해상도 (the raster's spatial resolution)를 나타냅니다.

 

## ncell() : number of cells (pixels) = number of rows * number of columns
ncell(srtm_raster)
# [1] 212505


## res() : the raster's spatial resolution
res(srtm_raster)
# [1] 0.0008333333 0.0008333333

 

 

extent() 함수는 공간의 x와 y 좌표의 최소, 최대값의 범위를 나타냅니다. 

반면에 summary() 함수는 각 셀의 특성 값의 분위수값과 평균의 요약 통계량을 나타냅니다.

 

## extent() : spatial extent
extent(srtm_raster)
# class      : Extent 
# xmin       : -113.2396 
# xmax       : -112.8521 
# ymin       : 37.13208 
# ymax       : 37.51292


## summary() : Summary of the values of a Raster* object (quartiles and mean)
summary(srtm_raster)
# srtm
# Min.    1024
# 1st Qu. 1535
# Median  1837
# 3rd Qu. 2115
# Max.    2892
# NA's       0
# Warning message:
# In .local(object, ...) :
# summary is an estimate based on a sample of 1e+05 cells (47.06% of all cells)

 

 

inMemory() 함수는 레스터 데이터가 메모리에 저장되어 있는지 (기본 설정), 아니면  디스크에 저장되어 있는지를 확인할 때 사용합니다. (블리언 반환)

 

## inmemory() : whether the raster data is stored in memory (the default) or on disk
inMemory(srtm_raster)
# [1] FALSE

 

 

 

레스터 데이터를 시각화하는데는 여러개의 패키지와 함수가 있는데요, 가장 기본적으로 쉽고 빠르게 시각화하는 방법은 plot() 함수를 사용하는 것입니다.

 

## plotting
plot(srtm_raster, main = 'basic raster plot')

 

[ Reference ]

* Raster data: https://geocompr.robinlovelace.net/spatial-class.html#raster-data

 

 

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

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

 

 

728x90
반응형
Posted by Rfriend
,

Rfriend님의
글이 좋았다면 응원을 보내주세요!

[R] 반복문 프로그램 진행 경과 막대로 나타내기 (progress bar)



R에서 for loop 반복문을 실행하거나 데이터를 다운로드 하다보면 전체 수행 회수 중에서 현재 어디까지 진행이 된 것인지 중간 중간 확인해보고 싶을 때가 있습니다. 특히 연산이 오래걸리고 for loop 반복회수가 많거나, 대용량 데이터를 다운로드 해야할 경우라면 얼마나 진행이 되었고, 얼마나 더 기다려야 하는건지 중간에 확인할 수 없다면 무척 답답할 것입니다.


이번 포스팅에서는 R의 반복문이나 다운로드의 진행 경과 (progress status) 를 출력해주는 방법을 소개하겠습니다.


'progress' 패키지를 활용하여

(1) 순환문 진행상태를 막대(bar), 비율(percent), 추정 완료시간 출력하기

(2) 순환문 진행상태를 현재까지 수행 개수(current), 총 수행 (예정) 개수(total), 현재까지 상세 소요 시간(elapsedfull) 출력하기

(3) 다운로드의 진행 경과 출력하기 (download progress)


'randomForest' 패키지를 사용하여

(4) Random Forest 분석의 진행상태 출력하기



[ R 'progress' 패키지를 사용한 진행 경과 막대 출력하기 ]







  (1) 순환문 진행상태를 막대(bar), 비율(percent), 추정 완료시간 출력하기


R의 'progress' 패키지는 진행 경과를 막대 형태로 출력해주는데 있어 다양한 매개변수를 사용하여 원하는 형식으로 진행 상태를 출력(Configurable progress bars)할 수 있게 해줍니다.

아래의 (1) ~ (3) 까지는 R의 'progress' 패키지를 사용합니다. 처음 사용하는 분이라면 'progress' 패키지 설치부터 해주세요.



## 'progress' Package: Configurable progress bars

install.packages("progress")
library(progress)

 



'progress' 패키지의 진행 상태 막대 (Progress bar) 는 R6 객체이며, progress_bar$new() 를 사용해서 진행상태 막대 R6 객체를 생성할 수 있습니다.


progress_bar$new() 안의 total 은 진행상태를 확인하는 총 회수 (tick 의 사전적 의미는 시계가 '똑딱 똑딱 움직이는 소리'를 나타냄) 를 지정해주는 매개변수입니다. 이 total 값이 알려져있지 않을 때는 NA 를 사용하면 되며, 기본값은 total=100 입니다.


그리고, 아래의 10,000회를 반복하는 for loop 예제문에서 Sys.sleep() 는 괄호 안에 지정한 시간만큼 R 실행을 잠깐 멈추라는 뜻입니다. 아래 예의 for loop 반복문 안에는 특별히 연산을 수행하는 것이 없으므로 Sys.sleep(1 / 1000) * 10,000 회 수행하는 만큼의 시간이 걸리겠네요.


기본 설정값(default)만 사용한 결과, 아래처럼 진행 경과 막대(progress bar)와 비율(percent) 이 매우 간결한 형태로 출력이 되었습니다.


그리고 100% 모두 진행이 되면 콘솔 창에서 진행 경과 막대 출력 결과가 사라집니다.(clear=TRUE 가 기본 설정이므로)



pb <- progress_bar$new(total = 10000)
for (i in 1:10000) {
  pb$tick()
  Sys.sleep(1 / 1000) # Suspend execution of R expressions for a specified time interval
}





앞서 진행 경과 막대(progress bar)가 설정가능하다("Configurable")고 말씀드렸는데요, 이는progress_bar$new() 로 R6 객체를 만들 때 "format 매개변수"에 다양한 "설정 값 (Token)" 들 중에서 원하는 설정 값을 선택해서 넣어주면 됩니다.

아래의 예에서는 format 매개변수에 Token 값으로
  • :bar   >> 진행 경과 막대 출력
  • :percent  >> 진행 경과 비율 출력
  • :eta   >> 진행 완료 추정 시간 출력
을 사용하였습니다. 그러면 아래의 화면 캡쳐처럼 for loop 문의 진행경과가 막대 형태와 퍼센트로 출력이 되고, 추정 완료 시간도 같이 출력됩니다.

clear = FALSE 로 매개변수를 설정해주면 진행 경과가 100%가 된 이후에도 '진행 경과 막대' 출력 결과가 사라지지 않고 그대로 남아있게 됩니다.


## format: The format of the progress bar with :bar, :percent, :eta tokens
pb <- progress_bar$new(
  format = " Progress: [:bar] :percent, Estimated completion time: :eta",
  total = 10000, # totalnumber of ticks to complete (default 100)
  clear = FALSE, # whether to clear the progress bar on completion (default TRUE)
  width= 80 # width of the progress bar
  )


for (i in 1:10000) {
  pb$tick() # increases the number of ticks by one (or another specified value)
  Sys.sleep(1 / 1000)
}





위의 progress_bar$new() 에서 width 매개변수는 '진행 경과 막대'의 폭을 설정해줄 때 사용합니다. R의 기본 설정 폭의 값은 options('width') 로 확인해볼 수 있는데요, 115 이군요. 위의 예에서는 width=80 으로서 기본 설정값보다는 좀더 폭이 좁게 조정해 본 것입니다.



options('width')
# $width
# [1] 115

 




  (2) 순환문 진행상태를 현재까지 tick 개수(current), 총 tick 개수(total),

      현재까지 상세 소요 시간(elapsedfull) 출력하기


이번에는 진행 상태를 현재 tick 의 회수 (current), 총 tick 의 회수 (total), 그리고 현재까지 실제 소요 시간(elapsedfull) 을 출력해보겠습니다.

  • :current       >> 현재까지 tick 개수
  • :total          >> 총 tick 개수
  • :elapsedfull  >> 현재까지 소요된 상세 시간 (hh:mm:ss format)


## format: The format of the progress bar with :current, :total, :elapsedfull tokens
pb <- progress_bar$new(
  format = "Current tick number :current / Total tick number :total, Elapsed time :elapsedfull",
  total = 10000, # totalnumber of ticks to complete (default 100)
  clear = FALSE, # whether to clear the progress bar on completion (default TRUE)
  width= 80 # width of the progress bar
)


for (i in 1:10000) {
  pb$tick() # increases the number of ticks by one (or another specified value)
  Sys.sleep(1 / 1000)
}






  (3) 다운로드의 진행 경과 출력하기 (download progress)


progress_bar$new() 의 R6 객체에 format 설정을 통해서 다운로드 할 때 파일 크기는 얼마이고, 그중에서 몇 바이트를 다운로드 진행했는지도 경과 막대로 표시할 수 있습니다.


  • :rate       >> 다운로드 비율, 초당 Bytes
  • :elapsed  >> 소요 시간 (단위: 초)



## Download (or other) rates
pb <- progress_bar$new(
  format = " downloading foobar at :rate, got :bytes in :elapsed",
  clear = FALSE, total = NA, width = 60)

f <- function() {
  for (i in 1:100) {
    pb$tick(sample(1:100 * 1000, 1))
    Sys.sleep(2/100)
  }
  pb$tick(1e7)
  invisible()
}
f()






  (4) randomForest' 패키지를 사용하여 Random Forest 분석의 진행상태 출력하기


Decision Tree를 여러개 수행해서 결과를 평균내거나 다수결로 취하는 ensemble 기법인 Random Forest 의 경우, R의 random forest 패키지의 do.trace 옵션을 사용하면 진행 경과를 출력할 수 있습니다.

do.trace=TRUE 또는 do.trance=integer (로그를 콘솔에 남기기 원하는 간격) 의 형식으로 입력해주시면 됩니다.

아래 예시 코드는 ntree=10000 으로 해놓고, 1000회 별로 콘솔에 로그를 남기게 됩니다.



library("randomForest")
set.seed(1)

rf = randomForest(Species~., data=iris,
    ntree=10000,
    proximity=T,
    do.trace=1000) # <--- 추가


ntree OOB 1 2 3
1000: 4.67% 0.00% 6.00% 8.00%
2000: 4.00% 0.00% 6.00% 6.00%
3000: 4.00% 0.00% 6.00% 6.00%
4000: 4.00% 0.00% 6.00% 6.00%
5000: 4.00% 0.00% 6.00% 6.00%
6000: 4.00% 0.00% 6.00% 6.00%
7000: 4.00% 0.00% 6.00% 6.00%
8000: 4.00% 0.00% 6.00% 6.00%
9000: 4.67% 0.00% 6.00% 8.00%
10000: 4.67% 0.00% 6.00% 8.00%

 



[ Reference ]

* Pakcage 'progress' : https://cran.r-project.org/web/packages/progress/progress.pdf


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

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



728x90
반응형
Posted by Rfriend
,

Rfriend님의
글이 좋았다면 응원을 보내주세요!

지난번 포스팅에서는 R의 sf 패키지를 사용하여 "sf 클래스 (The sf class)"를 만드는 방법을 소개하였습니다.


이번 포스팅에서는 지리공간 벡터 데이터(GeoSpatial Vector data)를 분석할 때 사용하는 "sf 패키지"와 "sp 패키지"의 관계, "sf 클래스"와 "sp 클래스" 간 변환에 대해서 소개하겠습니다.


(1) R의 sf 패키지와 sp 패키지의 관계/ 역사

(2) sf 클래스를 sp 클래스로 전환하기: as(sf_class, Class="Spatial")

(3) sp 클래스를 sf 클래스로 전환하기: st_as_sf(sp_class)



[ R 지리공간 데이터 분석 패키지: sf package와 sp package 클래스 간 전환 ]



  (1) R의 sf 패키지와 sp 패키지의 관계/ 역사


sf 패키지는 {sp 패키지 + R 인터페이스 with (GDAL, GEOS, PROJ) + s2 패키지} 의 패키지 기능들을 승계한 종합 패키지로서, 지리공간 벡터 데이터를 단순 지리특성 (Simple features) 으로 인코딩할 수 있게 해줍니다. [1]



>> sp 패키지 : Classes and Methods for Spatial Data [2], [3]


(sf 패키지의 부모격이 되는) sp 패키지는 2005년 Pebesma 와 Bivand 가 개발해서 처음으로 공개하였습니다. 2005년 전에는 지리적인 좌표가 일반적으로 여느 숫자와 다름없이 처리되었다면, sp 패키지가 나오고 부터는 지리공간 점(Point), 선(LineString), 면/다각형(Polygon), 격자(Grid) 와 속성(Attributes)을 지원하는 클래스(Classes)와 메소드(Methods) 로 처리되는 방식으로 바뀌게 됩니다.


sp 패키지는 테두리 상자(bounding box), 좌표 참조 시스템(CRS, Coordinate Reference System), 속성(Attributes) 등의 정보를 S4 클래스 시스템을 사용해서 "Spatial" 객체안의 슬롯에 저장을 합니다. 이를 통해 지리공간 데이터에 대한 연산 작업을 할 수 있게 해줍니다. 또한, sp 패키지는 지리공간 데이터에 대해 summary()나 plot() 함수와 같은 R 에 내장된 함수도 사용할 수 있게 해줍니다.



>> sf 패키지 : Simple Features for R


sf 패키지는 Edzer Pebesma, Roger Bivand 등이 2016년 10월에 최초로 오픈소스로 공개하였으며, R로 단순 지리특성 기하 (Simple Feature Geometry) 형태로 지리 벡터 데이터를 인코딩하는 표준화된 방법을 지원합니다. sf 패키지는 sp 패키지의 기능을 승계하였으며, 이에 더해 지리공간 데이터를 읽고 쓰는 'GDAL', 지리적 연산을 할 때 사용하는 'GEOS', 지도의 투영 변환(projection conversions)과 데이터 변환(datum transformations)을 위한 'PROJ' 와 R과의 인터페이스를 제공합니다. 그리고 선택적으로 지리적 좌표에 대한 구면 기하 연산 (spherical geometry operations) 을 위해 's2' 패키지를 사용합니다.


단순 지리특성 모델 (Simple Features Model)을 지원하는 "sf 패키지"를 사용하면 좋은 점들로는[4],

  • 지리공간 벡터 데이터를 빠르게 읽고 쓸 수 있음
  • 지리공간 벡터 데이터 시각화 성능의 고도화
    (tmap, leaflet, mapview 지리공간 데이터 시각화 패키지가 sf 클래스 지원)
  • 대부분의 연산에서 sf 객체는 DataFrame 처럼 처리가 가능함
  • sf 함수들은 '%>%' 연산자 (chain operator) 와 함께 사용할 수 있고, R의 tidyverse 패키지들과도 잘 작동함
    (sp 패키지도 spdplyr 패키지를 설치하면 dplyr의  %>% 체인 연산자와 기능을 사용할 수 있음)
  • sf 함수이름은 'st_' 로 시작하여 상대적으로 일관성이 있고 직관적임

등을 들 수 있습니다.


sf 패키지의 장점이 이렇게 많으므로 지리공간 벡터 데이터를 처리하고 분석하고 시각화하는데 있어 sf 패키지를 사용하지 않을 이유가 없습니다!




  (2) sf 클래스를 sp 클래스로 전환하기: as(sf_class, Class="Spatial")


예제로 spData 패키지에 들어있는 'world' 라는 이름의 "sf", "data.frame" 클래스의 데이터셋을 사용하겠습니다.



library(sp) # Classes and methods for spatial data
library(sf) # Support for simple features, a standardized way to encode spatial vector data.
library(spData) # load geographic data

names(world)
# [1] "iso_a2"    "name_long" "continent" "region_un" "subregion" "type"      "area_km2"  "pop"      
# [9] "lifeExp"   "gdpPercap" "geom"

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




"sf" 클래스를 "sp" 클래스로 변환(converting the sf class into the sp class) 할 때는 as(sf_class, Class = "Spatial") 함수와 매개변수를 사용합니다.(sp 패키지는 S4의 Spatial 클래스를 사용)



## Converting sf class into sp's Spatial objects using as(sf_class, Class="Spatial))
world_sp_class = as(world, Class = "Spatial")

class(world_sp_class)
# [1] "SpatialPolygonsDataFrame"
# attr(,"package")
# [1] "sp"

 




  (3) sp 클래스를 sf 클래스로 전환하기: st_as_sf(sp_class)


반대로, "sp" 클래스를 "sf" 클래스로 변환(converting the sp class into sf class) 할 때는 sf 패키지의 st_as_sf() 함수를 사용합니다.



## Converting Spatial objects into sf class using st_as_sf()
world_sf_class = st_as_sf(world_sp_class)


class(world_sf_class)
# [1] "sf"         "data.frame"




[ Reference ]

[1] sf package: https://cran.r-project.org/web/packages/sf/index.html

[2] sp package: https://cran.r-project.org/web/packages/sp/index.html

[3] The history of R-spatial: https://geocompr.robinlovelace.net/intro.html

[4] Why simple features?: https://geocompr.robinlovelace.net/spatial-class.html



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

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




728x90
반응형
Posted by Rfriend
,

Rfriend님의
글이 좋았다면 응원을 보내주세요!

지난번 포스팅에서는 R의 sf 패키지를 사용해서


- 벡터 데이터의 기하 유형 (Geometry types of vector data)을 단순 지리특성 기하로 정의하기

- 벡터 데이터의 여러개 단순 지리특성 기하(sfg)를 하나의 단순 지리특성 칼럼(sfc) 객체로 합치기

- st_sfc() 함수로 좌표 참조 시스템(CRS, Coordinate Reference Systems) 을 설정하고 확인하기


에 대해서 소개하였습니다.


이번 포스팅에서는 앞서의 포스팅들에서 소개했던 내용을 모두 포괄하는 sf패키지의 "sf 클래스 (The sf class)" 를 만드는 방법을 소개하겠습니다.


(1) sf 패키지의 "sf 클래스(the sf class)"는 무엇인가?

(2) sf 클래스 만들기: st_sf()

(3) sf 클래스의 클래스, 기하유형, CRS 확인하기: class(), st_geometry_type(), st_crs)






  (1) sf 패키지의 "sf 클래스(the sf class)"는 무엇인가?


R의 sf 패키지의 "sf 클래스(the sf class)"는 위의 그림에서 설명하는 바와 같이,


>> 지리공간 벡터 데이터 (GeoSpatial Vector Data)의 단순 지리특성 기하 (sfg, Simple feature geometry) 객체들을 리스트로 하나의 객체로 합친 단순 지리특성 칼럼 (sfc, Simple feature columns) 과,

>> 지리공간 벡터의 단순 지리특성 기하별로 속성 값 (예: 이름, 용도, 면적, 가격, 인구수, 온도, 습도 등)들을 모아놓은 data.frame 데이터 구조를


하나의 "sf 클래스" ( = sfc + data.frame) 로 합쳐 놓은 것입니다.


실제 지리공간 벡터 데이터를 분석한다고 했을 때는 sfc 객체의 기하 (Geometry) 정보와 data.frame 의 속성 (Non-geometric Attributes) 정보가 모두 있어야지, 기하 유형(점, 선, 면/다각형 등) 별로 속성 값을 분석할 수 있겠지요.


R의 sf 패키지로 지리공간 벡터 데이터 분석을 한다고 했을 때 새로운 용어가 자꾸 나오다 보니 위계 체계 (hierarchy), 자료 유형과 구조가 혼란스러울 수 있는데요, 위의 그림의 위계 체계, Input/Output 관계를 참고하시기 바랍니다.




  (2) sf 클래스 만들기: st_sf()



(2-1) 기하 객체 (Geometry object)


아래의 예에서는 먼저, library(sf) 로 sf 패키지를 로딩하고, 지리공간 벡터 데이터(Vector data)의 점 기하 유형 (Point geometry types) 들로 이루어진 단순 지리특성 기하 (sfg, Simple feature geometry) 객체를 만들어 보겠습니다.


다음으로 st_sfc() 함수를 사용하여 앞서 만든 3개의 점(points) 단순 지리특성 기하(sfg) 객체를 하나의 단순 지리특성 칼럼 (sfc, Simple feature columns) 객체로 합치고, 이때 좌표 참조 시스템(CRS, Coordinate Reference Systems) 은 EPSG 코드 정의를 사용해서 WGS84 (World Geodetic System 1984) 인 crs = 4326 을 설정해주겠습니다.



library(sf)

## (1) sfg(Simple feature geometry) object
## points: the coordinate of (longitude, latitude)
seoul_point_sfg_1 = st_point(c(127.059, 37.511))
seoul_point_sfg_2 = st_point(c(127.063, 37.512))
seoul_point_sfg_3 = st_point(c(127.073, 37.516))


## (2) sfc(Simple feature columns) object
seoul_points_sfc = st_sfc(seoul_point_sfg_1, seoul_point_sfg_2, seoul_point_sfg_3,
                          crs = 4326) # EPSG code, WGS84

seoul_points_sfc
# Geometry set for 3 features
# geometry type:  POINT
# dimension:      XY
# bbox:           xmin: 127.059 ymin: 37.511 xmax: 127.073 ymax: 37.516
# geographic CRS: WGS 84
# POINT (127.059 37.511)
# POINT (127.063 37.512)
# POINT (127.073 37.516)

 




(2-2) 비기하 data.frame 속성 (data.frame Non-geometry Attributes)


이번에는 서울의 COEX, GBC(현대자동차 Global Business Center), JansilStadium(잠실 종합 운동장) 의 3개 POI (Point of Interest)별 면적(area)과 가격(price) 속성(Attributes) 정보를 가지는 data.frame을 만들어보겠습니다. (아래 속성 데이터 값은 모두 그냥 가짜로 입력한 예시 값입니다.)



## (3) data.frame object with non-geometric attributes
seoul_df_attrib = data.frame(
  poi = c("COEX", "GBC", "JamsilStadium"), # Point of Interest
  area = c(20, 25, 15),
  price = c(80, 100, 50)
)

seoul_df_attrib
#             poi area price
# 1          COEX   20    80
# 2           GBC   25   100
# 3 JamsilStadium   15    50





(2-3) sf 클래스 = sfc 기하 객체(sfc Geometry) + data.frame 속성 (Non-geometry Attributes)


마지막으로, st_sf() 함수를 사용하여 위의 (2-1)에서 정의한 벡터 데이터의 기하 정보와, (2-2)에서 정의한 비기하 속성값 data.frame (Non-geometry Attributes) 을 합쳐서 sf 클래스 (the sf class) 를 만들어보겠습니다.



## (4) sf object using st_sf() function
## : sfc object (geometry) + data.frame object (non-geometric attributes)
seoul_sf = st_sf(seoul_df_attrib, # data.frame object (non-geometric attributes)
                 geometry = seoul_points_sfc) # sfc object (geometry)

seoul_sf
# Simple feature collection with 3 features and 3 fields
# geometry type:  POINT
# dimension:      XY
# bbox:           xmin: 127.059 ymin: 37.511 xmax: 127.073 ymax: 37.516
# geographic CRS: WGS 84
# poi area price               geometry
# 1          COEX   20    80 POINT (127.059 37.511)
# 2           GBC   25   100 POINT (127.063 37.512)
# 3 JamsilStadium   15    50 POINT (127.073 37.516)





  (3) sf 클래스의 클래스, 기하 유형, CRS 확인하기: class(), st_geometry_type(), st_crs)


sf 클래스의 클래스를 class() 함수를 사용해서 확인해보면, 기하(Geometry) 정보를 가지고 있는 "sf" 클래스 객체와 속성 정보 값을 가지고 있는 "data.frame" 클래스의 두개 클래스 객체로 구성되어 있음을 알 수 있습니다.


## check the class
class(seoul_sf)
# [1] "sf"         "data.frame"




st_geometry_type() 함수를 사용해서 기하 유형을 확인해보면, 이 예제에서는 점(Point) 3개로 구성되어 있음을 알 수 있습니다.



## check the geometry types
st_geometry_type(seoul_sf)
# [1] POINT POINT POINT
# 18 Levels: GEOMETRY POINT LINESTRING POLYGON MULTIPOINT MULTILINESTRING MULTIPOLYGON ... TRIANGLE




좌표 참조 시스템(CRS, Coordinate Reference Systems)은 st_crs() 함수로 확인할 수 있는데요, 위의 (1)번에서 단순 지리특성 칼럼(sfc, Simple feature columns)을 생성할 때 crs = 4326  으로 설정을 해주었기 때문에 User input: EPSG:4326 으로 등록되어 있음을 확인할 수 있습니다.



## check the CRS
st_crs(seoul_sf)
# Coordinate Reference System:
#   User input: EPSG:4326
# wkt:
#   GEOGCRS["WGS 84",
#           DATUM["World Geodetic System 1984",
#                 ELLIPSOID["WGS 84",6378137,298.257223563,
#                           LENGTHUNIT["metre",1]]],
#           PRIMEM["Greenwich",0,
#                  ANGLEUNIT["degree",0.0174532925199433]],
#           CS[ellipsoidal,2],
#           AXIS["geodetic latitude (Lat)",north,
#                ORDER[1],
#                ANGLEUNIT["degree",0.0174532925199433]],
#           AXIS["geodetic longitude (Lon)",east,
#                ORDER[2],
#                ANGLEUNIT["degree",0.0174532925199433]],
#           USAGE[
#             SCOPE["unknown"],
#             AREA["World"],
#             BBOX[-90,-180,90,180]],
#           ID["EPSG",4326]]




이렇게 만든 sf 클래스를 대상으로 leaflet 패키지를 사용하여 웹 기반으로 상호작용하는 동적 지도를 만드는 방법은 https://rfriend.tistory.com/593 를 참조하세요.

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

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


728x90
반응형
Posted by Rfriend
,

Rfriend님의
글이 좋았다면 응원을 보내주세요!

지난번 포스팅에서는 R의 sf패키지의 st_sfc() 함수를 사용해서 여러개의 기본 지리특성 기하(sfg, Simple feature geometry)들을 하나의 단순 지리특성 칼럼(sfc, Simple feature columns) 객체로 합치는 방법을 소개하였습니다.

 

지리공간 데이터 분석을 할 때 제일 먼저 해야 하고 또 확인해야 하는 것이 좌표계, 좌표 참조 시스템(Coordinate Reference Systems, CRS)을 설정하고 확인하는 것입니다.

 

좌표 참조 시스템(CRS) (또는 좌표계, 공간 참조 시스템, Spatial Reference System, SRS)좌표 기반으로 국소, 지역, 또는 전지구적으로 지리상의 객체를 위치시킬 때 사용되는 참조 시스템입니다. 좌표 참조 시스템은 특정 지도 투영법 (map projection)을 정의하고, 서로 다른 좌표 참조 시스템 간의 변환을 할 수 있도록 해줍니다.

 

위도, 경도 좌표의 기준이 어떻게 되는지, 누가 언제 정의한 좌표 체계인지를 정의해 놓은 것입니다.

 

오픈소스 PostGIS 에서는 Spatial Reference System Table에 SRID 별로 'auth_name', 'auth_srid', 'srtext', 'proj4text' 등의 정보가 저장되어 있는데요, R에서도 좌표계, 좌표 참조 시스템(CRS)이 있어서 좌표계 간에 변환을 할 수 있게 해줍니다.

 

 

 

참고로, 전세계를 한번에 나타내야 할 때 많이 사용하는 전지구 좌표계로는 WGS84(EPSG:4326), Bessel 1841(EPSG:4004, EPSG:4162), GRS80(EPSG)4019, EPGS:4737), Google Mercator(EPSG:3857) 등이 있습니다.

 

 

[전지구 좌표계]
*WGS84 경위도: GPS가 사용하는 좌표계EPSG:4326, EPSG:4166 (Korean 1995)+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs 
*Bessel 1841 경위도: 한국과 일본에 잘 맞는 지역타원체를 사용한 좌표계EPSG:4004, EPSG:4162 (Korean 1985)+proj=longlat +ellps=bessel +no_defs +towgs84=-115.80,474.99,674.11,1.16,-2.31,-1.63,6.43
*GRS80 경위도: WGS84와 거의 유사EPSG:4019, EPSG:4737 (Korean 2000)+proj=longlat +ellps=GRS80 +no_defs
*Google Mercator: 구글지도/빙지도/야후지도/OSM 등 에서 사용중인 좌표계EPSG:3857(공식), 

EPSG:900913(통칭)

+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs

출처: https://www.osgeo.kr/17 [OSGeo(Open Source GeoSpatial) 한국어 지부 - OSGeo Korean Chapter]

 

한국 주요 좌표계 EPSG 코드 및 proj4 인자 정리는 https://www.osgeo.kr/17 를 참고하세요.

 

 

이번 포스팅에서는 R의 sf패키지를 사용해서 단순 지리특성 칼럼(sfc, Simple feature columns)에서

 

(1) 특정 좌표 참조 시스템(CRS, Coordinate Reference Systems) 지정해주기

(2) 좌표 참조 시스템(CRS, Coordinate Reference Systems) 확인: st_crs()

(3) raw proj4string 을 사용해서 좌표 참조 시스템 지정해주기

 

를 소개하겠습니다.

 

 

  (1) 특정 좌표 참조 시스템(CRS, Coordinate Reference Systems) 지정해주기

 

기본 지리특성 칼럼(sfc, Simple feature columns)의 좌표 참조 시스템의 기본 값은 NA (Not Available) 입니다.

 

 

library(sf)

## Coordinate Reference Systems (CRS)
## default CRS is NA (Not Available)
point_sfg_1 = st_point(c(1, 2))
point_sfg_2 = st_point(c(3, 4))
points_sfc = st_sfc(point_sfg_1, point_sfg_2)

st_crs(points_sfc)
# Coordinate Reference System: NA

 

 

 

st_sfc() 함수의 crs 매개변수를 사용해서 특정 좌표 참조 시스템(CRS)을 설정해줄 수 있습니다.

가령, 아래 예에서는 EPSG 코드 4326 번으로 좌표 참조 시스템을 설정하여 보았습니다.

 

* EPSG 코드의 EPSG는 European Petroleum Survey Group 의 앞글자를 딴 줄임말입니다. EPSG 그룹은 좌표 시스템 정보 데이터베이스와 지도 투영법과 데이터에 대한 표준 문서 자료를 발행하는 곳입니다.

 

* WGS84 는 World Geodetic System 1984 의 앞글자를 딴 줄임말입니다. WGS84 는 전지구 좌표계로서, 우리가 매일 사용하는 GPS(Global Positioning System), 지도 제작, 측지학 등에서 사용합니다.

 

 

 

## EPSG: European Petroleum Survey Group
## : They publish a database of coordinate system information
##   plus some very good related documents on map projections and datums

## Specify CRS using 'crs' argument
point_sfc_wgs = st_sfc(point_sfg_1, point_sfg_2, crs = 4326) # EPSG code

 

point_sfc_wgs
# Geometry set for 2 features
# geometry type:  POINT
# dimension:      XY
# bbox:           xmin: 1 ymin: 2 xmax: 3 ymax: 4
# geographic CRS: WGS 84
# POINT (1 2)
# POINT (3 4)

 

 

 

 

  (2) 좌표 참조 시스템(CRS, Coordinate Reference Systems) 확인: st_crs()

 

좌표 참조 시스템 확인은 st_crs() 함수를 사용합니다.

 

 

st_crs(point_sfc_wgs)
# Coordinate Reference System:
#   User input: EPSG:4326
# wkt:
#   GEOGCRS["WGS 84",
#           DATUM["World Geodetic System 1984",
#                 ELLIPSOID["WGS 84",6378137,298.257223563,
#                           LENGTHUNIT["metre",1]]],
#           PRIMEM["Greenwich",0,
#                  ANGLEUNIT["degree",0.0174532925199433]],
#           CS[ellipsoidal,2],
#           AXIS["geodetic latitude (Lat)",north,
#                ORDER[1],
#                ANGLEUNIT["degree",0.0174532925199433]],
#           AXIS["geodetic longitude (Lon)",east,
#                ORDER[2],
#                ANGLEUNIT["degree",0.0174532925199433]],
#           USAGE[
#             SCOPE["unknown"],
#             AREA["World"],
#             BBOX[-90,-180,90,180]],
#           ID["EPSG",4326]]

 

 

 

 

  (3) raw proj4string 을 사용해서 좌표 참조 시스템 지정해주기

 

raw proj4string 을 직접 입력해서 좌표 참조 시스템을 지정해주려면 st_sfc() 함수의 crs 매개변수에 직접 따옴표 안에 proj, datum 을 입력해주면 됩니다. 아래 예는 WGS84를 raw proj4string로 해서 좌표 참조 시스템을 등록한 것입니다.

 

 

## to specify CRS using a raw proj4string
point_sfc_wgs_2 = st_sfc(point_sfg_1, point_sfg_2,
                         crs = "+proj=longlat +datum=WGS84 +no_defs")

st_crs(point_sfc_wgs_2)
# Coordinate Reference System:
#   User input: +proj=longlat +datum=WGS84 +no_defs
# wkt:
#   GEOGCRS["unknown",
#           DATUM["World Geodetic System 1984",
#                 ELLIPSOID["WGS 84",6378137,298.257223563,
#                           LENGTHUNIT["metre",1]],
#                 ID["EPSG",6326]],
#           PRIMEM["Greenwich",0,
#                  ANGLEUNIT["degree",0.0174532925199433],
#                  ID["EPSG",8901]],
#           CS[ellipsoidal,2],
#           AXIS["longitude",east,
#                ORDER[1],
#                ANGLEUNIT["degree",0.0174532925199433,
#                          ID["EPSG",9122]]],
#           AXIS["latitude",north,
#                ORDER[2],
#                ANGLEUNIT["degree",0.0174532925199433,
#                          ID["EPSG",9122]]]]

 

 

 

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

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

 

 

728x90
반응형
Posted by Rfriend
,

Rfriend님의
글이 좋았다면 응원을 보내주세요!