R, Python 분석과 프로그래밍의 친구 (by R Friend)

지난번 포스팅에서는 지리공간 벡터 데이터의속성 정보의 행과 열을 가져오기 (rfriend.tistory.com/622) 에 대해서 알아보았습니다. 그리고 지리공간 벡터 데이터는 지리 정보(geometry column) + 속성 정보(attributes)로 이루어지며, 속성 정보의 일부를 subset 하여 가져오면 지리정보(geometry column)가 같이 따라 온다는 것도 소개하였습니다.

이번 포스팅에서는 지리공간 벡터 데이터의 속성 정보의 요약 정보를 그룹별로 집계하는 방법(aggregating geographic vector attritubes by group)을 R 패키지별로 나누어서 소개하겠습니다.

(1) Base R의 aggregate() 함수를 사용해서 지리 벡터 데이터의 속성 정보를 그룹별로 집계하기

(2) dplyr 의 summarise() 함수를 사용해서 지리 벡터 데이터의 속성 정보를 그룹별로 집계하기

(3) data.table 로 지리 벡터 데이터의 속성 정보를 그룹별로 집계하기

(1) Base R의 aggregate() 함수를 사용해서 지리 벡터 데이터의 속성 정보를 그룹별로 집계하기

그리고, Base R, dplyr, data.table 의 3개 패키지 내 함수를 사용해서 아래 문제의 그룹별 속성 정보 집계를 해보겠습니다.

대륙(continent) 그룹 별로 인구(pop) 속성 데이터의 합(summation)을 구하고, 국가의 수 (number of countries) 를 집계한 후, --> 인구 합이 가장 큰 상위 3개(top 3) 대륙을 선별하시오.

예제로 spData 패키지에 들어있는 'world' 라는 이름의, 10개 속성(attributes)과 1개 지리기하 칼럼(geometry column)을 가지고 있는 세계 국가에 대한 지리 벡터 데이터를 사용하겠습니다. 그리고 속성 정보 집계를 하는데 사용할 dplyr, data.table 패키지도 불러오겠습니다.

## =================================
## GeoSpatial Data Analysis using R
## : Vector attribute aggregation
## =================================

library(sf)
library(spData)     # for dataset
library(dplyr)      # for data aggregation
library(data.table) # for data aggregation

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

Base R 패키지에 기본으로 내장되어 있는 aggregate(x ~ group, FUN, data, ...) 함수를 사용해서 '대륙(continent)' 그룹별로 속성 정보인 '인구(pop)' 의 합(sum) 을 집계해보겠습니다. 이렇게 합계를 집계하면 data.frame 을 반환하였으며, 집계된 결과에 지리 기하(geometry) 정보는 없습니다.

## aggregation operations summarize dataset by a 'grouping variable'.

## (1) aggregation using aggregate() function from Base R
## : returns a non-spatial data.frame
world_agg_pop = aggregate(pop ~ continent, 
                          FUN = sum, 
                          data = world, 
                          na.rm = TRUE)

str(world_agg_pop)
# 'data.frame':	6 obs. of  2 variables:
#   $ continent: chr  "Africa" "Asia" "Europe" "North America" ...
# $ pop      : num  1.15e+09 4.31e+09 6.69e+08 5.65e+08 3.78e+07 ...

class(world_agg_pop)
# [1] "data.frame"


world_agg_pop
# continent        pop
# 1        Africa 1154946633
# 2          Asia 4311408059
# 3        Europe  669036256
# 4 North America  565028684
# 5       Oceania   37757833
# 6 South America  412060811

'sf' 객체에 대해 aggregate() 함수를 적용하면 'sf' 클래스에 있는 aggregate.sf() 함수가 자동으로 활성화되어 적용이 됩니다. world['pop'] 의 'sf' 클래스를 x 로 하여 aggregate(x = sf_object, by = list(), FUN, ...) 구문으로 아래와 같이 대륙 그룹별로 인구의 합계를 집계하면 "data.frame"을 확장한 "sf" 클래스의 집계된 객체를 반환하며, "data.frame"의 제일 끝에 "geometry" 칼럼이 리스트로 같이 집계(aggregation)가 됩니다. (즉, 각 대륙별로 국가들의 다각형 면(polygon)을 모아놓은 리스트가 "geometry"에 생김).

## aggregate.sf() is activated automatically when 'x' is an sf object and a 'by' argument is provided.
## : returns a spatial 'sf' class and data.frame
world_agg_pop2 = aggregate(
  x = world['pop'], 
  by = list(world$continent), # a lost of grouping elements
  FUN = sum,  
  na.rm = TRUE)

str(world_agg_pop2)
# Classes 'sf' and 'data.frame':	8 obs. of  3 variables:
#   $ Group.1 : chr  "Africa" "Antarctica" "Asia" "Europe" ...
# $ pop     : num  1.15e+09 0.00 4.31e+09 6.69e+08 5.65e+08 ...
# $ geometry:sfc_GEOMETRY of length 8; first list element: List of 2
# ..$ :List of 1
# .. ..$ : num [1:356, 1:2] 32.8 32.6 32.5 32.2 31.5 ...
# ..$ :List of 1
# .. ..$ : num [1:49, 1:2] 49.5 49.8 50.1 50.2 50.5 ...
# ..- attr(*, "class")= chr [1:3] "XY" "MULTIPOLYGON" "sfg"
# - attr(*, "sf_column")= chr "geometry"
# - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: 3 2
# ..- attr(*, "names")= chr [1:2] "Group.1" "pop"


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

names(world_agg_pop2)
# [1] "Group.1"  "pop"      "geometry"

world_agg_pop2
# Simple feature collection with 8 features and 2 fields
# Attribute-geometry relationship: 0 constant, 1 aggregate, 1 identity
# geometry type:  GEOMETRY
# dimension:      XY
# bbox:           xmin: -180 ymin: -90 xmax: 180 ymax: 83.64513
# geographic CRS: WGS 84
# Group.1        pop                       geometry
# 1                  Africa 1154946633 MULTIPOLYGON (((32.83012 -2...
# 2              Antarctica          0 MULTIPOLYGON (((-163.7129 -...
# 3                    Asia 4311408059 MULTIPOLYGON (((120.295 -10...
# 4                  Europe  669036256 MULTIPOLYGON (((-51.6578 4....
# 5           North America  565028684 MULTIPOLYGON (((-61.68 10.7...
# 6                 Oceania   37757833 MULTIPOLYGON (((169.6678 -4...
# 7 Seven seas (open ocean)          0 POLYGON ((68.935 -48.625, 6...
# 8           South America  412060811 MULTIPOLYGON (((-66.95992 -...

참고로, world['pop'] 는 "sf" 객체이기 때문에  aggregate() 함수 적용 시 자동으로  ''sf" 클래스의 aggregate.sf() 함수가 자동으로 활성화되어 적용되고 집계 결과가 "sf" 객체로 반환되는 반면에, world$pop 는 숫자형 벡터(numeric vector) 이므로 aggregate() 함수를 적용하면 "data.frame" 으로 집계 결과를 반환하게 됩니다.

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

world['pop']
# 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                                                                                    geom
# <dbl>                                                             <MULTIPOLYGON [arc_degree]>
# 1    885806 (((180 -16.06713, 180 -16.55522, 179.3641 -16.80135, 178.7251 -17.01204, 178.5968 -1...
# 2  52234869 (((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  35535348 (((-122.84 49, -122.9742 49.00254, -124.9102 49.98456, -125.6246 50.41656, -127.4356...
# 5 318622525 (((-122.84 49, -120 49, -117.0312 49, -116.0482 49, -113 49, -110.05 49, -107.05 49,...
# 6  17288285 (((87.35997 49.21498, 86.59878 48.54918, 85.76823 48.45575, 85.72048 47.45297, 85.16...
# 7  30757700 (((55.96819 41.30864, 55.92892 44.99586, 58.50313 45.5868, 58.68999 45.50001, 60.239...
# 8   7755785 (((141.0002 -2.600151, 142.7352 -3.289153, 144.584 -3.861418, 145.2732 -4.373738, 14...
# 9 255131116 (((141.0002 -2.600151, 141.0171 -5.859022, 141.0339 -9.117893, 140.1434 -8.297168, 1...
# 10  42981515 (((-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


## -- vs. numeric
class(world$pop)
# [1] "numeric"

world$pop[1:20]
# [1]    885806  52234869        NA  35535348 318622525  17288285  30757700   7755785 255131116  42981515  17613798  73722860
# [13]  13513125  46024250  37737913  13569438  10572466  10405844 143819666    382169

이제 우리가 풀고자 하는 요구사항에 답해보겠습니다. 요구사항은 "대륙 그룹별로 속성정보인 인구(pop)의 합과 국가의 수를 집계하고, 인구 합을 기준으로 상위 3개의 대륙을 찾으라"는 것이었습니다. aggregate(x ~ group, data, FUN) 구문으로 그룹별 집계를 하되, 집계해야 할 것이 (a) 인구의 "합(sum)"과, & (b) 국가의 "개수(number of rows)" 의 2개 이므로 FUN = function(x) c(pop = sum(x), n = length(x)) 의 함수를 사용해서 합과 개수를 구하도록 했습니다. 그리고 do.call(data.frame, ...) 을 사용해서 리스트로 반환된 집계 결과를 data.frame으로 변환을 해주었습니다.

그리고 "인구 합을 기준으로 상위 3개 대륙"을 찾기 위해서 order(-world_pop$pop.pop) 로 내림차순 정렬('-' 부호가 내림차순 정렬을 의미함. order 는 정렬 후의 위치 인덱스를 반환함.)을 한 인덱스를 사용해서 정렬을 한 후에 [1:3, ] 으로 상위 3개 행만을 가져왔습니다.

## group by aggregation: sum of pop, number of rows
world_pop = do.call(data.frame, aggregate(pop ~ continent, 
                      data = world, 
                      FUN = function(x) c(pop = sum(x), n = length(x))))




## sorting by pop.pop in descending order and select top 3 
world_pop[order(-world_pop$pop.pop),][1:3,]

#   continent    pop.pop pop.n
# 2      Asia 4311408059    45
# 1    Africa 1154946633    48
# 3    Europe  669036256    37

(2) dplyr 의 summarise() 함수를 사용해서 지리 벡터 데이터의 속성 정보를 그룹별로 집계하기

위의 (1)번에서는 Base R의 aggregate() 함수를 사용하였는데요, dplyr 패키지의 group_by() 로 기준이 되는 그룹을 지정해주고 summarise() (또는 summarize()) 함수를 사용해서 다양한 집계 함수를 사용하여 Base R 보다 가독성이 좋고 속도가 빠른 코드를 짤 수 있습니다.

아래는 dplyr 패키지의 체인(%>%) 으로 파이프 연산자(pipe operator)를 사용해서 대륙 그룹별로 인구의 합계를 구해본 것입니다. 이때, 'sf' 객체에 대해 summarise() 함수를 적용하면 집계 결과는 자동으로 "data.frame"을 확장한 "sf" 클래스의 객체로 반환을 하고, 제일 마지막 칼럼에 "geom" 지리기하 리스트를 포함하게 됩니다.

## (2) aggregation using summarise() function from dplyr
library(dplyr)

world_agg_pop3 = world %>% 
  group_by(continent) %>% 
  summarise(sum(pop, na.rm = TRUE))

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

names(world_agg_pop3)
# [1] "continent" "pop"   "geom"

dplyr 의 summarise() 함수를 쓰면 (a) 집계한 결과의 칼럼 이름을 새로 부여할 수 있고, (b) 여러개의 집계 함수(예: 합, 개수, 평균, 표준편차, 최대, 최소 등)를 동시에 사용할 수 있다는 편리한 장점이 있습니다.

아래 예에서는 대륙별로 인구에 대한 합(sum(pop, na.rm=TRUE))과 국가의 개수(n()) 를 구하고, 집계된 후 결과의 칼럼 이름도 "pop_sum", "n_countries" 로 새로 부여해 본 것입니다.

## summarise(sum(), n()) {dplyr}
world %>% 
  group_by(continent) %>% 
  summarise(pop_sum = sum(pop, na.rm = TRUE), 
            n_countries = n())

# `summarise()` ungrouping output (override with `.groups` argument)
# Simple feature collection with 8 features and 3 fields
# geometry type:  MULTIPOLYGON
# dimension:      XY
# bbox:           xmin: -180 ymin: -90 xmax: 180 ymax: 83.64513
# geographic CRS: WGS 84
# # A tibble: 8 x 4
#   continent                pop_sum n_countries                                                                          geom
# <chr>                      <dbl>       <int>                                                   <MULTIPOLYGON [arc_degree]>
# 1 Africa                1154946633          51 (((32.83012 -26.74219, 32.58026 -27.47016, 32.46213 -28.30101, 32.20339 -28.~
# 2 Antarctica                     0           1 (((-48.66062 -78.04702, -48.1514 -78.04707, -46.66286 -77.83148, -45.15476 -~
# 3 Asia                  4311408059          47 (((120.295 -10.25865, 118.9678 -9.557969, 119.9003 -9.36134, 120.4258 -9.665~
# 4 Europe                 669036256          39 (((-51.6578 4.156232, -52.24934 3.241094, -52.55642 2.504705, -52.93966 2.12~
# 5 North America          565028684          18 (((-61.68 10.76, -61.105 10.89, -60.895 10.855, -60.935 10.11, -61.77 10, -6~
# 6 Oceania                 37757833           7 (((169.6678 -43.55533, 170.5249 -43.03169, 171.1251 -42.51275, 171.5697 -41.~
# 7 Seven seas (open oce~          0           1 (((68.935 -48.625, 69.58 -48.94, 70.525 -49.065, 70.56 -49.255, 70.28 -49.71~
# 8 South America          412060811          13 (((-66.95992 -54.89681, -67.29103 -55.30124, -68.14863 -55.61183, -68.63999 ~

이제 위에서 주어진 요구사항인 "대륙 그룹별로 인구의 합계와 국가 수를 집계하고, 인구합을 기준으로 상위 3개 대륙을 찾으라"는 문제를 dplyr 의 체인 파이프 연산자를 써서 풀어보겠습니다. top_n(n = 3, wt = pop) 함수로 인구합을 기준으로 상위 3개를 구하였고, arrange(desc(pop)) 로 인구합을 기준으로 내림차순 정렬하였습니다. 이때 지리기하 칼럼은 생략하고 속성(attritubes)에 대한 집계 결과만 얻고자 하므로 제일 마지막에 st_drop_geometry() 함수로 "geometry" 칼럼을 제거하였습니다.

## finding the world's 3 most populous continents by chaining functions using dplyr
world %>% 
  dplyr::select(pop, continent) %>% 
  group_by(continent) %>% 
  summarise(pop = sum(pop, na.rm = TRUE), n_countries = n()) %>% 
  top_n(n = 3, wt = pop) %>% # select top 3 rows by 'pop' varialbe
  arrange(desc(pop)) %>%     # sorting in a descending order
  st_drop_geometry()         # drop geometry column
  
# `summarise()` ungrouping output (override with `.groups` argument)
# # A tibble: 3 x 3
# continent        pop n_countries
# * <chr>          <dbl>       <int>
#   1 Asia      4311408059          47
# 2 Africa    1154946633          51
# 3 Europe     669036256          39  

(3) data.table 로 지리 벡터 데이터의 속성 정보를 그룹별로 집계하기

마지막으로 data.table 패키지를 사용해서 "대륙 그룹별로 인구 합계와 국가 수를 집계하고, 인구 합계 기준 상위 3개 대륙을 찾으라"는 문제를 풀어보겠습니다.

먼저, world_dt = data.table(world) 로 'world' sf 객체를 data.table 객체로 변환하였습니다. 그리고 data.table 의 체인 파이프를 사용해서 필요한 작업 흐름을 이어서 작성해보았는데요, 역시 data.table 이 코드가 가장 간결하고, 또 속도도 가장 빠릅니다. 아래 코드에서 .(pop = sum(pop, na.rm = T), n = .N) 은 인구 합계(sum(pop, na.rm=T)와 행의 개수(.N) 를 대륙 그룹별로(by = list(continent)) 집계한 후에, [order(-pop)] 로 인구 합계를 기준으로 내림차순 정렬('-' 부호는 내림차순을 의미함)을 하고, [1:3]으로 상위 3개 대륙의 집계 결과만 가져온 것입니다. 코드가 '간결함' 그 자체입니다!

단, 이때 최종 집계 결과는 data.table 로 반환이 되며, 지리기하 정보는 없습니다.

## (3) aggregating vector attributes using data.table package
library(data.table)

world_dt = data.table(world)

world_dt[, .(pop = sum(pop, na.rm = TRUE), n = .N), 
         by = list(continent)][order(-pop)][1:3]

## or equivalently above, using chaining
world_dt[, 
         .(pop = sum(pop, na.rm = TRUE), n = .N), # aggregate
         by = list(continent)][    # by groups
           order(-pop)][           # sorting in descending order
             1:3]                  # subset rows (1:3)

# continent        pop  n
# 1:      Asia 4311408059 47
# 2:    Africa 1154946633 51
# 3:    Europe  669036256 39

만약, 그룹별 지리기하 벡터 데이터를 집계한 결과에 지리기하 칼럼(geometry column)을 같이 리스트로 묶어서 속성정보와 함께 집계하고 싶으면 data.table 객체에서 칼럼을 가져올 때 world_dt[, .(pop = sum(pop, na.rm = T), n = .N, geom = geom) 처럼 "geom" 칼럼도 같이 가져오고, 이후의 집계 및 정렬, 인덱싱은 동일합니다.

## with geometry column (lists)
world_dt[, .(pop = sum(pop, na.rm = TRUE), n = .N, geom = geom), 
         by = list(continent)][order(-pop)][1:3]
# continent        pop  n     geom
# 1:      Asia 4311408059 47  <XY[1]>
# 2:      Asia 4311408059 47  <XY[1]>
# 3:      Asia 4311408059 47 <XY[13]>

[Reference]

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

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

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

Tistory 의 코드 블록에 'R' 언어도 사용할 수 있게 해달라고 며칠전에 제안 문의를 드렸었는데요, 드디어 'R' 언어가 추가되었네요. 잘 사용하겠습니다. Tistory, 감사합니다. 꾸벅~!  

지리공간 벡터 데이터는 (a) 점(point), 선(linestring), 면(polygon) 등의 리스트로 이루어진 지리공간 데이터(geographic data)와, 지리공간 데이터를 제외한 (b) 속성 데이터(attribute data)로 구성이 됩니다. 속성 데이터(attritubes)를 사용해서 각 지리공간 별 특성 (예: 지역 이름, 면적, 인구, 예상수명 등)을 표현하게 됩니다. 

이번 포스팅에서는 벡터 데이터에서 속성 정보를 가져오는 여러가지 방법을 소개하겠습니다.  'sf' 객체는 R data.frame 에서 사용하는 클래스를 지원(!!!)하므로 Base R 이나 혹은 dplyr 패키지를 사용하여 데이터 전처리하는 작업에 이미 능숙한 분이라면 이번 포스팅은 특별한 것 없이 매우 쉽게, 복습하는 기분으로 보실 수 있을 거예요. 

(1) 'sf' 객체에서 속성 정보만 가져오기: st_drop_geometry()

(2) Base R 구문으로 벡터 데이터 속성 정보의 행과 열 가져오기
    (subsetting attributes from geographic vector data using Base R)

(3) dplyr 로 벡터 데이터 속성 정보의 행과 열 가져오기
    (subsetting attributes from geographic vector data using dplyr)

(4) 한개 칼럼만 가져온 결과를 벡터로 반환하기 (extracting a single vector)

먼저, 지리공간 데이터 객체를 다루기 위한 sf 패키지와 데이터 전처리에 사용할 dplyr 패키지를 불러오고, 예제로 사용할 세계 나라 데이터 ('world') 를 spData 패키지로 부터 가져오겠습니다. 

'world' 벡터 데이터셋은 10개의 속성 데이터 (10 attributes) 와 1개의 지리기하 칼럼 (1 geometry column) 으로 구성이 되어 있으며, 총 177개의 행 (즉, 국가) 에 대한 정보가 들어 있습니다. 

## ==============================
## R GeoSpatial data analysis
## - Vector Attribute data operations
## - reference: https://geocompr.robinlovelace.net/attr.html
## ==============================

library(sf)
library(dplyr)
library(spData)

## vector dataset 'world': 10 attritubes and 1 geometry column
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" ...

dim(world)
# [1] 177  11

(1) 'sf' 객체에서 속성 정보만 가져오기: st_drop_geometry()

지리공간 'sf' 객체에는 항상 점(points)/선(lines)/면(polygons) 등의 지리기하 데이터를 리스트로 가지고 있는 geometry 칼럼이 항상 따라 다닙니다. 만약 'sf' 객체로 부터 이 geometry 칼럼을 제거하고 나머지 속성 정보만으로 DataFrame 을 만들고 싶다면 sf 패키지의 st_drop_geometry() 메소드를 사용합니다. geometry 칼럼의 경우 지리기하 점/선/면 등의 리스트 정보를 가지고 있기 때문에 메모리 점유가 크기 때문에, 사용할 필요가 없다면 geometry 칼럼은 제거하고 속성 정보만 dataframe 으로 만들어서 분석을 진행하는게 좋겠습니다. 

(아래의 (2)번 부터 소개하는 부분 가져오기(subset) 방법을 사용하면 여전히 geometry 칼럼이 그림자차럼 같이 따라와서 여전히 달라붙어 있을 것입니다.) 

## Extracting the attribute data of an 'sf' object
world_df <- st_drop_geometry(world)
class(world_df)
# [1] "tbl_df"     "tbl"        "data.frame"

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

(2) Base R 구문으로 벡터 데이터 속성 정보의 행과 열 가져오기
    (subsetting attributes from geographic vector data using Base R)

R data.frame 에서 i 행(row)과 j 열(column)을 가져올 때는 df[i, j] , subset(), $ 구문을 사용합니다.  이때  (a) i 행과  j 열에는 정수로 위치(subset rows and columns by position) 을 사용하거나, (b) j 행의 이름 (subset columns by name)을 사용할 수 있으며, (c) 논리 벡터 (logical vector)를 사용해서 i 행의 부분집합(subset rows by a logical vector) 을 가져올 수 있습니다. 

여기서 다시 한번 짚고 넘어갈 점은, 특정 행과 열의 부분집합을 가져오면 제일 끝에 geometry 칼럼이 꼭 껌딱지처럼 달라붙어서 같이 따라온다는 점입니다.

(2-a-1) 위치를 지정해서 지리공간 벡터 데이터로 부터 부분 행 가져오기 (subset rows by position)

## -- Vecter arrtibute subsetting
## subsetting columns returns results with geometry column as well.

## -- subset rows by position
world[1:3, ]

# Simple feature collection with 6 features and 10 fields
# geometry type:  MULTIPOLYGON
# dimension:      XY
# bbox:           xmin: -180 ymin: -18.28799 xmax: 180 ymax: 83.23324
# geographic CRS: WGS 84
# # A tibble: 6 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.555~
#   2 TZ     Tanzania   Africa     Africa    Eastern ~ Sover~   9.33e5  5.22e7    64.2     2402. (((33.90371 -0.95, 34.07262 -~
#   3 EH     Western S~ Africa     Africa    Northern~ Indet~   9.63e4 NA         NA         NA  (((-8.66559 27.65643, -8.6651~

(2-a-2) 위치를 지정해서 지리공간 벡터 데이터로 부터 부분 열 가져오기 (subset columns by position)

## -- subset columns by position
world[, 1:3]

# Simple feature collection with 177 features and 3 fields
# geometry type:  MULTIPOLYGON
# dimension:      XY
# bbox:           xmin: -180 ymin: -90 xmax: 180 ymax: 83.64513
# geographic CRS: WGS 84
# # A tibble: 177 x 4
# iso_a2 name_long        continent                                                                                   geom
# <chr>  <chr>            <chr>                                                                <MULTIPOLYGON [arc_degree]>
# 1 FJ     Fiji             Oceania       (((180 -16.06713, 180 -16.55522, 179.3641 -16.80135, 178.7251 -17.01204, 178.5968~
# 2 TZ     Tanzania         Africa        (((33.90371 -0.95, 34.07262 -1.05982, 37.69869 -3.09699, 37.7669 -3.67712, 39.202~
# 3 EH     Western Sahara   Africa        (((-8.66559 27.65643, -8.665124 27.58948, -8.6844 27.39574, -8.687294 25.88106, -~
# 4 CA     Canada           North America (((-122.84 49, -122.9742 49.00254, -124.9102 49.98456, -125.6246 50.41656, -127.4~
# 5 US     United States    North America (((-122.84 49, -120 49, -117.0312 49, -116.0482 49, -113 49, -110.05 49, -107.05 ~
# 6 KZ     Kazakhstan       Asia          (((87.35997 49.21498, 86.59878 48.54918, 85.76823 48.45575, 85.72048 47.45297, 85~
# 7 UZ     Uzbekistan       Asia          (((55.96819 41.30864, 55.92892 44.99586, 58.50313 45.5868, 58.68999 45.50001, 60.~
# 8 PG     Papua New Guinea Oceania       (((141.0002 -2.600151, 142.7352 -3.289153, 144.584 -3.861418, 145.2732 -4.373738,~
# 9 ID     Indonesia        Asia          (((141.0002 -2.600151, 141.0171 -5.859022, 141.0339 -9.117893, 140.1434 -8.297168~
# 10 AR     Argentina        South America (((-68.63401 -52.63637, -68.25 -53.1, -67.75 -53.85, -66.45 -54.45, -65.05 -54.7,~
# # ... with 167 more rows

(2-b) 칼럼 이름을 사용해서 지리공간 벡터 데이터로 부터 부분 열 가져오기 (subset columns by name)

## -- subset columns by name
world[, c("name_long", "area_km2")]

# Simple feature collection with 177 features and 2 fields
# geometry type:  MULTIPOLYGON
# dimension:      XY
# bbox:           xmin: -180 ymin: -90 xmax: 180 ymax: 83.64513
# geographic CRS: WGS 84
# # A tibble: 177 x 3
# name_long         area_km2                                                                                    geom
# <chr>                <dbl>                                                             <MULTIPOLYGON [arc_degree]>
# 1 Fiji                19290. (((180 -16.06713, 180 -16.55522, 179.3641 -16.80135, 178.7251 -17.01204, 178.5968 -1...
# 2 Tanzania           932746. (((33.90371 -0.95, 34.07262 -1.05982, 37.69869 -3.09699, 37.7669 -3.67712, 39.20222 ...
# 3 Western Sahara      96271. (((-8.66559 27.65643, -8.665124 27.58948, -8.6844 27.39574, -8.687294 25.88106, -11....
# 4 Canada           10036043. (((-122.84 49, -122.9742 49.00254, -124.9102 49.98456, -125.6246 50.41656, -127.4356...
# 5 United States     9510744. (((-122.84 49, -120 49, -117.0312 49, -116.0482 49, -113 49, -110.05 49, -107.05 49,...
# 6 Kazakhstan        2729811. (((87.35997 49.21498, 86.59878 48.54918, 85.76823 48.45575, 85.72048 47.45297, 85.16...
# 7 Uzbekistan         461410. (((55.96819 41.30864, 55.92892 44.99586, 58.50313 45.5868, 58.68999 45.50001, 60.239...
# 8 Papua New Guinea   464520. (((141.0002 -2.600151, 142.7352 -3.289153, 144.584 -3.861418, 145.2732 -4.373738, 14...
# 9 Indonesia         1819251. (((141.0002 -2.600151, 141.0171 -5.859022, 141.0339 -9.117893, 140.1434 -8.297168, 1...
# 10 Argentina         2784469. (((-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

(2-c) 논리 벡터를 사용해서 지리공간 벡터 데이터로 부터 부분 행 가져오기 (subset rows by a logical vectors)

## -- subset using 'logical vectors'
small_area_bool <- world$area_km2 < 10000
summary(small_area_bool)
#    Mode   FALSE    TRUE 
# logical     170       7

small_countries <- world[small_area_bool, ]
nrow(small_countries)
# [1]  7

## or equivalently above
small_countries <- world[world$area_km2 < 10000, ]

## -- using base R subset() function
small_countries <- subset(world, area_km2 < 10000)

(3) dplyr 로 벡터 데이터 속성 정보의 행과 열 가져오기
    (subsetting attributes from geographic vector data using dplyr)

R의 dplyr 패키지를 사용해서 지리공간 벡터 데이터의 속성 정보를 처리하면 Base R 대비 코드의 가독성이 좋고, 속도가 빠르다는 장점이 있습니다.  dplyr 패키지에서 체인('%>%')으로 파이프 연산자 (pipe operator) 를 사용하면 작업 흐름 (work flow) 를 자연스럽게 따라가면서 코딩을 할 수 있어서 코드의 가독성이 상당히 좋습니다.  그리고 dplyr 은 소스코드가 C++로 되어 있어 속도도 상당히 빠른 편입니다.  

dplyr 의 select(), slice(), filter(), pull() 함수를 사용해서 지리공간 벡터 데이터의 속성 정보의 부분 행과 열을 가져와 보겠습니다. 

(3-1) dplyr::select(sf, name) 함수를 사용하여 이름으로 특정 열 선택하기 (selects columns by name)

이때, select() 함수는 dplyr 패키지 뿐만 아니라 raster 패키지에도 존재하므로 dplyr::select() 처럼 select() 함수 이름 앞에 패키지 이름을 명시적으로 같이 표기해 줌으로써 불명확성을 없앨 수 있습니다. dplyr::select(data.frame 이름, 칼럼1, 칼럼2, ...) 의 구문으로 원하는 열의 속성 정보만 가져올 수 있습니다. (지리공간 벡터 데이터라고 해서 select() 구문에는 특별한 것이 없으며, subset 결과의 끝에 geometry 칼럼이 달라 붙어서 따라와 있다는 점이 다릅니다.)

## -- using 'dplyr' package
## : select(), slice(), filter(), pull()

## select(): selects columns by name
## : "geom" column remains. 
world1 = dplyr::select(world, name_long, area_km2)

names(world1)
# [1] "name_long" "area_km2"  "geom"

(3-2) dplyr::select(sf, position) 함수를 사용하여 위치로 특정 열 선택하기 (select columns by position)

## select() columns by position
world1 = dplyr::select(world, 2, 7)

names(world1)
# [1] "name_long" "area_km2"  "geom"

(3-3) dplyr::select(sf, name1:name2) 함수를 사용하여 범위 내의 모든 열 선택하기 (a range of columns by ':')

':' 연산자를 사용하여 두개의 위치(position1:position2)나 칼럼 이름(name1:name2)들의 범위 사이에 들어있는 모든 열을 한꺼번에 가져올 수 있습니다. 선택해서 가져오고 싶은 열이 많고 범위로 표현할 수 있을 경우 유용합니다. 

## select() also allows subsetting of a range of columns with the help of the : operator
world2 = dplyr::select(world, name_long:area_km2)

names(world2)
# [1] "name_long" "continent" "region_un" "subregion" "type"   "area_km2"  "geom"

(3-4) dplyr::select(sf, -name) 으로 특정 열을 제거하기 (omit specific columns with the - operator)

## Omit specific columns with the - operator
# all columns except subregion and area_km2 (inclusive)
world3 = dplyr::select(world, -region_un, -subregion, -area_km2, -gdpPercap)

names(world3)
# [1] "iso_a2"    "name_long" "continent" "type"      "pop"       "lifeExp"   "geom" 

(3-5) dplyr::select(sf, name_new = name_old) 로 선택해서 가져온 열 이름을 바꾸기 (subset & rename columns)

## select() lets you subset and rename columns at the same time
world4 = dplyr::select(world, name_long, population = pop)

names(world4)
# [1] "name_long"  "population" "geom"

(3-6) dplyr::select(sf, contains(string)) 로 특정 문자열을 포함한 칼럼을 선택하기

dplyr 의 select() 함수는 contains(), starts_with(), ends_with(), num_range() 와 같은 도우미 함수 (helper functions) 를 사용해서 매우 강력하고 편리하게 특정 칼럼을 선택해서 가져올 수 있습니다. 

## select() works with ‘helper functions’ for advanced subsetting operations, 
## including contains(), starts_with() and num_range()
## contains(): Contains a literal string
world5 = dplyr::select(world, contains('Ex'))

names(world5)
# [1] "lifeExp" "geom"

(3-7) dplyr::select(sf, starts_with(string)) 로 특정 문자열로 시작하는 칼럼을 선택하기

## starts_with(): Starts with a prefix
world6 = dplyr::select(world, starts_with('life'))

names(world6)
# [1] "lifeExp" "geom"

(3-8) dplyr::select(sf, ends_with(string)) 로 특정 문자열로 끝나는 칼럼을 선택하기

## ends_with(): Ends with a suffix
world7 = dplyr::select(world, ends_with('Exp'))

names(world7)
# [1] "lifeExp" "geom"

(3-9) dplyr::select(df, matches(regular_expression)) 으로 정규 표현식을 만족하는 특정 칼럼을 선택하기

아래 예는 정규 표현식(regular expression) "[p]" 를 matches() 도움미 함수 안에 사용해서, 칼럼 이름에 대/소문자 'p' 문자열이 들어있는 칼럼을 선택해서 가져온 것입니다.  

## matches(): Matches a regular expression.
world8 = dplyr::select(world, matches("[p]"))

names(world8)
# [1] "type"      "pop"       "lifeExp"   "gdpPercap" "geom"

(3-10) dplyr::select(df, num_range("x", num1:num2)) 으로 숫자 범위 안의 모든 칼럼 선택하기

만약 칼럼 이름이 "x1", "x2", "x3", "x4", "x5" 처럼 접두사(prefix)가 동일하고 뒤에 순차적인 정수가 붙어있는 형태라면 특정 정수 범위의 칼럼들을 선택해서 가져올 때 num_range() 도우미 함수를 select() 함수 안에 사용할 수 있습니다.  "world" sf 객체의 데이터셋은 칼럼 이름이 특정 접두사로 동일하게 시작하지 않으므로, 아래에는 간단한 data.frame 을 만들어서 예를 들어보였습니다. 

## -- num_range(): Matches a numerical range like x01, x02, x03. 

x1 = c(1:5)
x2 = c(6:10)
x3 = c(11:15)
x4 = c(16:20)
x5 = c(21:25)

df = data.frame(x1, x2, x3, x4, x5)
df
# x1 x2 x3 x4 x5
# 1  1  6 11 16 21
# 2  2  7 12 17 22
# 3  3  8 13 18 23
# 4  4  9 14 19 24
# 5  5 10 15 20 25

## subset using num_range()
df2 = dplyr::select(df, num_range("x", 1:3))
df2
# x1 x2 x3
# 1  1  6 11
# 2  2  7 12
# 3  3  8 13
# 4  4  9 14
# 5  5 10 15

(3-11) dplyr::slice(sf, position1:position2) 으로 특정 행 잘라서 가져오기 (slice rows by position)

## -- slice(): is the row-equivalent of select()
dplyr::slice(world, 3:5)

# Simple feature collection with 3 features and 10 fields
# geometry type:  MULTIPOLYGON
# dimension:      XY
# bbox:           xmin: -171.7911 ymin: 18.91619 xmax: -8.665124 ymax: 83.23324
# geographic CRS: WGS 84
# # A tibble: 3 x 11
# iso_a2 name_long continent region_un subregion type  area_km2     pop lifeExp gdpPercap
# <chr>  <chr>     <chr>     <chr>     <chr>     <chr>    <dbl>   <dbl>   <dbl>     <dbl>
#   1 EH     Western ~ Africa    Africa    Northern~ Inde~   9.63e4 NA         NA         NA 
# 2 CA     Canada    North Am~ Americas  Northern~ Sove~   1.00e7  3.55e7    82.0    43079.
# 3 US     United S~ North Am~ Americas  Northern~ Coun~   9.51e6  3.19e8    78.8    51922.
# # ... with 1 more variable: geom <MULTIPOLYGON [arc_degree]>

(3-12) dplyr::filter(sf, logical_vector) 로 조건을 만족하는 특정 행 걸러내기 (filter rows with conditions)

## -- filter()
# Countries with a population longer than 1 B.
world9 = filter(world, pop > 1000000000)

# world9
# Simple feature collection with 2 features and 10 fields
# geometry type:  MULTIPOLYGON
# dimension:      XY
# bbox:           xmin: 68.17665 ymin: 7.965535 xmax: 135.0263 ymax: 53.4588
# geographic CRS: WGS 84
# # A tibble: 2 x 11
# iso_a2 name_long continent region_un subregion type  area_km2    pop lifeExp gdpPercap
# * <chr>  <chr>     <chr>     <chr>     <chr>     <chr>    <dbl>  <dbl>   <dbl>     <dbl>
#   1 IN     India     Asia      Asia      Southern~ Sove~ 3142892. 1.29e9    68.0     5385.
# 2 CN     China     Asia      Asia      Eastern ~ Coun~ 9409830. 1.36e9    75.9    12759.
# # ... with 1 more variable: geom <MULTIPOLYGON [arc_degree]>

(3-13)  dplyr 의 체인('%>%')을 사용한 파이프 연산자로 filter(), select(), slice() 의 작업 흐름 만들기 (workflow by pipe operator)

아래에는 'world' 의 sf 객체 데이터셋으로 부터 먼저 대륙이 아시아인 국가를 걸러내고 (world %>% filter(continent == "Asia"), 이 결과를 받아서 긴 이름과 인구 칼럼을 선택하고 (%>% dplyr::select(name_long, pop)), 이 결과를 받아서 1행에서 3행까지 잘라서 가져오기(%>% slice(1:3)) 를 체인(chain %>%)으로 연결해서 파이프 연산자(pipe operator)로 작업 흐름에 따라 코딩을 한 것입니다. 

두번째 코드는 위와 똑같은 결과를 얻기 위해 함수 안의 함수인 중첩 함수 (nested functions) 로 표현해 본 것입니다. 앞서의 dplyr 의 chain(%>%)을 사용한 파이프 연산자(pipe operator)로 쓴 코드에 비해서 중첩 함수로 쓴 코드는 상대적으로 가독성이 떨어집니다. 

##-- pipe operator %>%: It enables expressive code: the output of a previous function 
## becomes the first argument of the next function, enabling chaining.
world10 = world %>%
  filter(continent == "Asia") %>%
  dplyr::select(name_long, pop) %>%
  slice(1:3)

world10
# Simple feature collection with 3 features and 2 fields
# geometry type:  MULTIPOLYGON
# dimension:      XY
# bbox:           xmin: 46.46645 ymin: -10.35999 xmax: 141.0339 ymax: 55.38525
# geographic CRS: WGS 84
# # A tibble: 3 x 3
# name_long        pop                                                                         geom
# <chr>          <dbl>                                                  <MULTIPOLYGON [arc_degree]>
# 1 Kazakhstan  17288285 (((87.35997 49.21498, 86.59878 48.54918, 85.76823 48.45575, 85.72048 47.452~
# 2 Uzbekistan  30757700 (((55.96819 41.30864, 55.92892 44.99586, 58.50313 45.5868, 58.68999 45.5000~
# 3 Indonesia  255131116 (((141.0002 -2.600151, 141.0171 -5.859022, 141.0339 -9.117893, 140.1434 -8.~


## -- The alternative to %>% is nested function calls, which is harder to read:
world11 = slice(
  dplyr::select(
    filter(world, continent == "Asia"),
    name_long, pop),
  1:3)

(4) dplyr::pull(sf, column_name) 으로 한개 칼럼만 가져온 결과를 벡터로 반환하기 (extracting a single vector)

일반적으로 dplyr 에서 특정 칼럼을 선택(select)해서 가져오면 data.frame 으로 결과를 반환합니다. 만약, 한개의 특정 칼럼만을 벡터(a single vector)로 가져오고 싶다면 명시적으로 pull() 함수를 사용해줘야 합니다. (혹은 df[ , "column_name"] 도 벡터 반환)

## Most dplyr verbs return a data frame. 
## To extract a single vector, one has to explicitly use the pull() command
# create throw-away data frame
d = data.frame(pop = 1:10, area = 1:10)

# (a) -- return data frame object when selecting a single column
d[, "pop", drop = FALSE] # equivalent to d["pop"]
select(d, pop)

# pop
# 1    1
# 2    2
# 3    3
# 4    4
# 5    5
# 6    6
# 7    7
# 8    8
# 9    9
# 10  10


# (b) -- return a vector when selecting a single column
d[, "pop"]
pull(d, pop)

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

아래의 예는 pull() 함수를 명시적으로 사용해서 인구("pop") 한개 열의 1행~5행의 관측치를 가져와서 벡터(vetor objects)로 반환한 것입니다. 

## -- vector objects

world$pop[1:5]
# [1]    885806  52234869        NA  35535348 318622525

pull(world, pop)[1:5]
# [1]    885806  52234869        NA  35535348 318622525

다음번 포스팅에서는 지리기하 벡터 데이터에서 속성 정보에 대해 그룹별로 집계(aggregating geographic vector attributes by group)하는 방법을 소개하겠습니다. (rfriend.tistory.com/624)

[Reference]

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

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

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

지난번 포스팅에서는 좌표계 (CRS, Coordinate Reference System) 로서 지리 좌표계(Geographic CRS)와 투영 좌표계(Projected CRS) 그리고 R에서 좌표계를 조회하고 설정하는 방법에 대해서 알아보았습니다. (rfriend.tistory.com/618)

이번 포스팅에서는 좌표계 (CRS) 정보 안에 들어있는 공간의 단위 (Spatial Units)에 대해서 알아보겠습니다. 지도를 제작하거나 볼 때 측정 단위 (measurement units)가 미터(meters) 인지 혹은 피트(feets) 인지 명시적으로 표현하고 정확하게 확인할 필요가 있습니다. 그래야지 벡터의 지리적 데이터나 레스터의 픽셀에서 측정되는 단위라는 맥락(context)를 알 수 있고, 실제 지표면과 지도 표현 간의 관계, 거리를 알 수 있고, 또 거리나 면적 등을 계산할 수 있기 때문입니다.

(1) 지리공간 벡터 데이터에서의 측정 단위 (Units in Vector data)

(2) 지리공간 레스터 데이터에서의 측정 단위 (Units in Raster data)

(1) 지리공간 벡터 데이터에서의 측정 단위 (Units in Vector data)

sf 객체의 지리공간 벡터 데이터는 단위에 대해서 native support 이여서, 다른 외부 모듈이나 확장 프로그램을 설치하지 않아도 sf 객체 내에 단위가 들어가 있습니다. 그래서 sf 객체 벡터 데이터에 대해서 연산을 하게 되면 units 패키지에 의해 정의된 "단위 속성(units arrtibute)"도 같이 반환해주어서 단위로 인한 혼란을 미연에 방지해줍니다. (대부분의 좌표계는 미터(meters)를 사용하지만, 일부는 피트(feets)를 사용하기 때문에 단위가 혼란스러울 수 있음. raster 패키지는 단위가 native support가 아님.)

아래 예제는 R의 spData 패키지에 들어있는 "world" 데이터셋에서 '대한민국 (South Korea, Republic of Korea)' 의 벡터 데이터를 가져와서, sf 패키지의 st_area() 함수로 면적을 계산해본 것입니다.

## ========================
## Units in Geospatial data
## ========================

## -- getting 'world' data from spData package
library(spData)
names(world)
# [1] "iso_a2"    "name_long" "continent" "region_un" "subregion" "type"      "area_km2"  "pop"       "lifeExp"  
# [10] "gdpPercap" "geom"

south_korea = world[world$name_long == "Republic of Korea", ]

south_korea
# Simple feature collection with 1 feature and 10 fields
# geometry type:  MULTIPOLYGON
# dimension:      XY
# bbox:           xmin: 126.1174 ymin: 34.39005 xmax: 129.4683 ymax: 38.61224
# geographic CRS: WGS 84
# # A tibble: 1 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_degre>
# 1 KR   Republic ~ Asia   Asia   Eastern ~ Sove~ 99044. 5.07e7    81.7  33426. (((126.1748 37.74969, 12~

plot(south_korea[1])

sf 패키지의 st_area() 함수로 벡터 데이터의 면적으로 계산하면, 결과값의 뒤에 [m^2] 이라고 해서 2차원 공간 상의 "제곱미터" 단위가 같이 반환된 것을 알 수 있습니다.

## -- calculating the area of South Korea
library(sf)
install.packages("lwgeom")
library(lwgeom)

st_area(south_korea)
# 99044366830 [m^2]

attributes() 함수로 벡터 데이터에 대한 면적 계산 결과에 대한 속성 정보를 확인해 볼 수 있습니다. "units" 클래스에 "m^2" (제곱미터) 단위라고 나오는 군요.

## -- attributes of units
attributes(st_area(south_korea))
# $units
# $numerator
# [1] "m" "m"
# 
# $denominator
# character(0)
# 
# attr(,"class")
# [1] "symbolic_units"
# 
# $class
# [1] "units"

위의 대한민국 면적 계산의 단위가 "제곱미터 [m^2]" 이다보니 결과값의 자리수가 너무 길게 나와서 보기에 어렵습니다.  그래서 보기에 편하도록 면적 계산의 단위를 "제곱킬로미터 [km^2]"로 변경하려면 units 패키지의 set_units(sf_object, units) 함수로 단위를 설정해주면 변경이 됩니다.

만약 set_units() 함수로 단위를 재설정하는게 아니라, 기존의 면적 단위인 '제곱미터(m^2)' 로 계산된 결과값을 1,000,000 으로 나누게 되면 결과값은 맞더라도 단위가 '제곱미터(m^2)' 로 그대로여서, 우리가 원하던 단위인 '제곱킬로미터(km^2)' 가 아니게 되므로 주의가 필요합니다. 즉, 단위를 변경해서 계산하고 결과값을 얻고 싶다면 units 패키지의 st_units() 함수로 단위를 명확하게 설정해주는게 필요합니다.

## -- cautious; the units is [m^2], not [km^2]
st_area(south_korea) / 1000000
# 99044.37 [m^2]


##-- set_units() : setting the units
library(units)
units::set_units(st_area(south_korea), km^2)
# 99044.37 [km^2]

(2) 지리공간 레스터 데이터에서의 측정 단위 (Units in Raster data)

앞에서 소개한 벡터 데이터를 다루는 sf 패키지는 단위가 native support 여서 조회나 계산 결과를 반환할 때 단위(units)를 속성으로 반환한다고 하였습니다. 하지만 레스터 데이터를 다루는 raster 패키지는 단위에 대해서 native support 가 아직까지는 아니므로, 단위에 대해서 혼란스러울 수 있으므로 조심해야 합니다.

아래의 예는 spDataLarge 패키지에 들어있는 strm.tif 파일을 raster() 함수로 읽어왔습니다. 이 데이터는 st_crs() 함수로 좌표계를 확인해보면 "WGS84 투영"을 사용하므로, 십진수 각도 (decimal degrees as units) 를 단위로 사용합니다. res() 함수로 해상도를 확인해보면, 단지 숫자형 벡터 (numeric vector) 만 반환할 뿐, 단위에 대한 속성 정보는 없습니다! (no units attributes).

## -- units in raster data
library(raster)
library(spDataLarge)

raster_filepath = system.file("raster/srtm.tif", package = "spDataLarge") 
new_raster = raster(raster_filepath)


## -- getting CRS
st_crs(new_raster)
# 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]]]]


## -- WGS84 projection with decimal degrees as units
## -- cautious; simply returns numeric vector without any unit.
res(new_raster)
# [1] 0.0008333333 0.0008333333

만약 우리가 아래의 예처럼 UTM 투영을 사용한다면, 이에 따라서 단위가 바뀔 텐데요, res()로 해상도를 살펴보면 역시 단지 숫자형 벡터만 반환할 뿐, 단위 속성 정보는 없습니다.  따라서 R에서 raster 패키지로 레스터 데이터를 다룰 때는 단위에 대해서 특별히 조심할 필요가 있겠습니다.

## -- if we used the UTM projection, the units would change.
repr = projectRaster(new_raster, crs = "+init=epsg:26912")

## -- no units attributes, just only returns numeric vector
res(repr)
# [1] 73.8 92.5

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

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

지난번 포스팅에서는 3가지 종류의 레스터 클래스 객체로서 RasterLayer 클래스, RasterBrick 클래스, RasterStack 클래스가 있고, 각 특성은 무엇인지, 그리고 R raster 패키지를 사용해서 각 레스터 클래스 객체를 어떻게 만들 수 있는지 살펴보았습니다. (rfriend.tistory.com/617)

또 이전 포스팅 rfriend.tistory.com/596 에서 좌표계가 무엇이고, 벡터 데이터를 다룰 때 사용하는 sf 패키지를 사용해서 좌표계를 지정하고 확인하는 방법을 소개한 적이 있습니다.

이번 포스팅에서는 지리 공간 데이터 분석에서 가장 기본이 되고 또 처음에 확인을 해보아야 하는 좌표계, 좌표 참조 시스템(CRS, Coordinate Reference Systems)에 대해서 좀더 자세하게 소개하겠습니다. (이하, CRS는 "좌표계"로 부르겠습니다.)

(1) 좌표계, 좌표 참조 시스템 (CRS, Coordinate Reference Systems)은 무엇인가?

(2) 지리 좌표계 (Geographic Coordinate Reference Systems)

(3) 투영(투사) 좌표계 (Projected Coordinate Reference Systems)

(4) R에서의 좌표계 (Coordinate Reference Systems in R)

(1) 좌표계 (CRS, Coordinate Reference Systems)는 무엇인가?

좌표계 (Coordinate Reference Systems) 는 지리공간 데이터를 저장/처리/분석하는데 있어서 지구의 표면을 어떻게 데이터의 공간 원소들로 표현할 지를 정의합니다.  좌표계는 좌표를 기반으로 국소, 지역, 또는 전지구적으로 지리상의 객체를 위치시킬 때 사용하는 참조 시스템입니다. 좌표계는 특정 지도 투영법(map projection)을 정의하고, 서로 다른 좌표계 간의 변환 (CRS transformtations)을 할 수 있도록 해줍니다.

좌표계는 지리 좌표계 (Geographic CRSs) 이거나 또는 투영 좌표계 (Projected CRSs) 인데요, 아래에 각각에 대해서 하나씩 설명하겠습니다.

(2) 지리 좌표계 (Geographic Coordinate Reference Systems)

지리 좌표계는 지구의 표면을 위도(Latitude)와 경도(Longitude) 의 두 개의 값을 사용해서 위치를 정의합니다. 위도(Latitude)는 적도면(Equator plane)으로 부터 남과 북 방향으로의 거리를 각도로 나타낸 것입니다. 경도(Longitude)는 본초 자오선 면(the Prime Meridian plane)으로 부터 동-서 방향으로의 거리를 각도로 나타낸 것입니다. 즉, 지리 좌표계에서의 거리는 미터(meters)가 아니라 각도(anglular distance)를 사용해서 측정합니다.

[  위도와 경도 (Latitude and Longitude) ]

지리 좌표계에서의 지구 표면은 구면의 표면(Spherical surface) 또는 타원의 표면(Ellipsoidal surface)으로 표현됩니다. 

구면 모델 (Spherical models) 은 지구가 주어진 반경(radius)으로 부터 완벽한 구면이라고 가정합니다. 구면 모델은 단순하고 계산이 빠르다는 장점이 있지만, 부정확하기 때문에 많이 이용되지는 않습니다. (실제의 지구 표면은 아래의 그림처럼 울퉁불퉁하며, 완벽한 구면이 아니기 때문입니다.)

타원 모델(Ellipsoidal models)은 적도 반지름(the equatorial radius)과 극 반지름(the polar radius)의 두 개의 모수를 사용해서 지구 표면을 정의합니다. 실제 지구는 아래의 그림처럼 위아래가 눌려서 옆으로 찌그러져 있으므로 타원 모델로 지구 표면을 표현하는 것이 구면 모델보다 더 적합합니다. (적도 반지름이 극 반지름보다 약 11.5 km 더 길이가 깁니다.)

[ 지구의 표면, 지오이드(GEOID), 타원(Ellipse), and 구(Sphere) ]

타원면(Ellipsoids)은 좌표계(CRSs)의 광의의 구성요소 중의 하나인 'the datum' 의 한 부분입니다. 'the datum'은 PROJ CRS 라이브러리 안의 'ellps' 모수로 무슨 타원면을 사용할지에 대한 정보와, 데카르트 좌표(Cartesian coordinates)와 지구 표면 위치 간의 정확한 관계에 대한 정보를 포함하고 있습니다.  이러한 추가적인 세부사항들은 'proj4string' 표기의 'towgs84' 매개변수에 저장이 되어 있습니다.

이러한 타원면의 the datum 정보들은 지구 표면의 지역적인 변이를 허용합니다. (가령, 큰 산맥이 있다면 지역 좌표계에 반영할 수 있음.)

the datum 에는 지역적 데이터(local datum)와 지심에서 바라 본 데이터(geocentric datum)의 두 종류가 있습니다.  'NAD83'과 같은 지역적 데이터(local datum)에서는 타원 표면이 특정 위치의 표면과 정렬이 되어 이동이 되어 있습니다. 반면에, 'WGS84'와 같은 지심에서 바라 본 데이터 (geocentric datum) 에서는 지구 중력의 중심이 중심 위치가 되며, 투사의 정확도가 특정 위치에 최적화되어 있지는 않습니다.

(3) 투영(투사) 좌표계 (Projected Coordinate Reference Systems)

투영 좌표계(Projected Coordinate Reference Systems)는 암묵적으로 평평한 표면 위의 데카르트 좌표를 기반으로 합니다. 투영 좌표계는 원점(origin), x축과 y축, 그리고 미터(meters)와 같은 선형 측정 단위를 가지고 있습니다. 모든 투영 좌표계는 지리 좌표계(geographic CRSs)를 기반으로 하며, 3차원의 지구 표면을 투영 좌표계의 x축(동-서)와 y축(남-북)의 값으로 변환하는 지도 투사에 의존합니다.

평평한 평면을 가정하고 투영 변환을 하기 때문에 면적, 방향, 거리, 모양 등에 대한 왜곡이 발생하게 됩니다. 투영된 좌표계는 이들 중 단지 하나 혹은 두개의 특성 만을 보존할 수 있습니다.  투영 후 좌표계가 보존하는 특성이 무엇이냐에 따라서 보통 이름이 지어지는데요, 가령, 등면적 보존 면적 (equal-area preserve area), 방위각 보존 방향 (azimuthal preserve direction), 등거리 보존 거리 (equidistant preserve distance), 등각 보존 국소 형태 (conformal preserve local shape) 등과 같이 부릅니다.

투영의 유형과 원리(projection types and principles)에는 크게 평면(planar), 원뿔(conic), 실린더(cylindrical) 의 3가지 그룹이 있습니다. (아래 그림 참조)

(a) 평면 투영 (planar projection) 은 데이터를 지구와 접촉하는 하나의 점이나 접촉하는 선과 함께 평평한 표면에 투사를 합니다. 평면 투영법은 일반적으로 극 지방(polar regions)을 매핑할 때 사용됩니다. 투영의 기준(시작)이 되는 표준 점(standard point)으로 부터 멀어질 수록 데이터 왜곡이 심해집니다.

(b) 원뿔 투영 (conic projection)에서는 지구의 표면이 원뿔과 접촉하는 하나의 선이나 두개의 선들로 원뿔에 투영이 됩니다. 원뿔과 만나는 선인 표준 선(standard line)에서 왜곡이 가장 작고, 표준 선으로부터 멀어질 수록 왜곡이 심해집니다. 따라서 중간 위도 지역(mid-latitude areas)의 투영법으로 가장 적합합니다.

(c) 원통 투영 (cylindrical projection) 은 지구의 표면을 원통과 접촉하는 하나의 선이나 두개의 선들을 원통에 투사합니다. 원통 투영은 보통 지구 전체(entire world)를 매핑할 때 많이 사용합니다.

[ 지도 투영의 3가지 유형(Projection types)과 각 투영 유형별 왜곡(Distortion) ]

(4) R에서의 좌표계 (CRSs in R)

R에서 좌표계를 표현할 때는 (a) epsg 코드 (epsg code)나 또는 (b) proj4string 정의 (proj4string definition)를 사용합니다. 이들 두 가지의 좌표계 표현 방법별로 장단점이 있습니다.

(a) epsg 코드는 더 짧고, 따라서 기억하기 쉽습니다. 그리고 epsg 코드는 단 하나의 잘 정의된 좌표계를 지칭합니다. (EPSG: European Petroleum Survey Group, 지도 투영과 datums 에 대한 좌표계 정보 데이터베이스를 제공함)

반면에, (b) proj4string 정의는 투사 유형이나 datum, 타원체 등의 다른 모수들을 구체화할 수 있는 유연성이 있다는 장점이 있는 대신에, 사용자가 구체화를 해야 하므로 길고 복잡하며 기억하기 어렵습니다.

지리공간 R 패키지들은 PROJ 라이브러리를 사용하여 광범위한 좌표계(CRSs)를 지원합니다. R의 rgdal::make_EPSG() 함수를 사용하면 R이 지원하는 투영 좌표계의 리스트를 data.frame 형태로 반환받아 확인할 수 있습니다. 아래에 조회해본 결과, 현재 6,609개의 좌표계를 지원하고 있으며, EPSG 코드 (code), 설명 (note), proj4string 정의 (prj4), 투영 방법 (prj_method) 의 4개 칼럼으로 좌표계 정보를 제공합니다.  가령 'EPSG코드가 4326' 인 행을 조회해보면 note에 'WGS 84' (World Geodetic System 1984) 라고 설명이 나오는군요.

## ========================================
## Coordinate Reference Systems (CRS) in R
## ========================================

crs_rgdal = rgdal::make_EPSG()

str(crs_rgdal)
# 'data.frame':	6609 obs. of  4 variables:
#   $ code      : int  3819 3821 3822 3823 3824 3887 3888 3889 3906 4000 ...
# $ note      : chr  "HD1909" "TWD67" "TWD97" "TWD97" ...
# $ prj4      : chr  "+proj=longlat +ellps=bessel +no_defs +type=crs" "+proj=longlat +ellps=aust_SA +no_defs +type=crs" "+proj=geocent +ellps=GRS80 +units=m +no_defs +type=crs" "+proj=longlat +ellps=GRS80 +no_defs +type=crs" ...
# $ prj_method: chr  "(null)" "(null)" "(null)" "(null)" ...


head(crs_rgdal, 10)
# code      note                                                   prj4 prj_method
# 1  3819    HD1909         +proj=longlat +ellps=bessel +no_defs +type=crs     (null)
# 2  3821     TWD67        +proj=longlat +ellps=aust_SA +no_defs +type=crs     (null)
# 3  3822     TWD97 +proj=geocent +ellps=GRS80 +units=m +no_defs +type=crs     (null)
# 4  3823     TWD97          +proj=longlat +ellps=GRS80 +no_defs +type=crs     (null)
# 5  3824     TWD97          +proj=longlat +ellps=GRS80 +no_defs +type=crs     (null)
# 6  3887      IGRS +proj=geocent +ellps=GRS80 +units=m +no_defs +type=crs     (null)
# 7  3888      IGRS          +proj=longlat +ellps=GRS80 +no_defs +type=crs     (null)
# 8  3889      IGRS          +proj=longlat +ellps=GRS80 +no_defs +type=crs     (null)
# 9  3906  MGI 1901         +proj=longlat +ellps=bessel +no_defs +type=crs     (null)
# 10 4000 MOLDREF99 +proj=geocent +ellps=GRS80 +units=m +no_defs +type=crs     (null)


## EPSG 4326 : WGS 84
crs_rgdal[crs_rgdal$code == 4326,]
# code   note                                          prj4 prj_method
# 261 4326 WGS 84 +proj=longlat +datum=WGS84 +no_defs +type=crs     (null)

• (4-1) 벡터 데이터의 좌표계 확인: sf 패키지의 st_crs() 함수
• (4-2) 레스터 지리 데이터의 좌표계 확인: raster 패키지의 projection() 함수

(4-1) 벡터 데이터의 좌표계 (CRSs in Vector geographic data)

벡터 지리 데이터에 대해서는 sf 패키지의 st_crs() 함수를 사용해서 좌표계를 확인할 수 있습니다.

아래의 예에서는 spDataLarge 패키지에 들어있는 Zion 국립 공원의 경계를 다각형면(Polygon)으로 나타내는 zion.gpkg 벡터 데이터를 st_read() 함수로 불러와서, st_crs() 함수로 좌표계를 조회해 본 것입니다. 좌표계가 'User input: UTM Zone 12, Northern Hemisphere' 라고 나오는군요.

## st_read() : read vector dataset in R sf package
library(sf)

vector_filepath = system.file("vector/zion.gpkg", package = "spDataLarge")
vector_zion = st_read(vector_filepath)
# Reading layer `zion' from data source `/Library/Frameworks/R.framework/Versions/4.0/Resources/library/spDataLarge/vector/zion.gpkg' using driver `GPKG'
# Simple feature collection with 1 feature and 11 fields
# geometry type:  POLYGON
# dimension:      XY
# bbox:           xmin: 302903.1 ymin: 4112244 xmax: 334735.5 ymax: 4153087
# projected CRS:  UTM Zone 12, Northern Hemisphere


## -- st_crs() : get CRS in vector dataset
st_crs(vector_zion)
# Coordinate Reference System:
#   User input: UTM Zone 12, Northern Hemisphere 
# wkt:
#   BOUNDCRS[
#     SOURCECRS[
#       PROJCRS["UTM Zone 12, Northern Hemisphere",
#               BASEGEOGCRS["GRS 1980(IUGG, 1980)",
#                           DATUM["unknown",
#                                 ELLIPSOID["GRS80",6378137,298.257222101,
#                                           LENGTHUNIT["metre",1,
#                                                      ID["EPSG",9001]]]],
#                           PRIMEM["Greenwich",0,
#                                  ANGLEUNIT["degree",0.0174532925199433]]],
#               CONVERSION["UTM zone 12N",
#                          METHOD["Transverse Mercator",
#                                 ID["EPSG",9807]],
#                          PARAMETER["Latitude of natural origin",0,
#                                    ANGLEUNIT["degree",0.0174532925199433],
#                                    ID["EPSG",8801]],
#                          PARAMETER["Longitude of natural origin",-111,
#                                    ANGLEUNIT["degree",0.0174532925199433],
#                                    ID["EPSG",8802]],
#                          PARAMETER["Scale factor at natural origin",0.9996,
#                                    SCALEUNIT["unity",1],
#                                    ID["EPSG",8805]],
#                          PARAMETER["False easting",500000,
#                                    LENGTHUNIT["Meter",1],
#                                    ID["EPSG",8806]],
#                          PARAMETER["False northing",0,
#                                    LENGTHUNIT["Meter",1],
#                                    ID["EPSG",8807]],
#                          ID["EPSG",16012]],
#               CS[Cartesian,2],
#               AXIS["(E)",east,
#                    ORDER[1],
#                    LENGTHUNIT["Meter",1]],
#               AXIS["(N)",north,
#                    ORDER[2],
#                    LENGTHUNIT["Meter",1]]]],
#     TARGETCRS[
#       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["latitude",north,
#                    ORDER[1],
#                    ANGLEUNIT["degree",0.0174532925199433]],
#               AXIS["longitude",east,
#                    ORDER[2],
#                    ANGLEUNIT["degree",0.0174532925199433]],
#               ID["EPSG",4326]]],
#     ABRIDGEDTRANSFORMATION["Transformation from GRS 1980(IUGG, 1980) to WGS84",
#                            METHOD["Position Vector transformation (geog2D domain)",
#                                   ID["EPSG",9606]],
#                            PARAMETER["X-axis translation",0,
#                                      ID["EPSG",8605]],
#                            PARAMETER["Y-axis translation",0,
#                                      ID["EPSG",8606]],
#                            PARAMETER["Z-axis translation",0,
#                                      ID["EPSG",8607]],
#                            PARAMETER["X-axis rotation",0,
#                                      ID["EPSG",8608]],
#                            PARAMETER["Y-axis rotation",0,
#                                      ID["EPSG",8609]],
#                            PARAMETER["Z-axis rotation",0,
#                                      ID["EPSG",8610]],
#                            PARAMETER["Scale difference",1,
#                                      ID["EPSG",8611]]]]

만약 좌표계가 비어있거나 잘못 입력되어 있는 경우 st_set_crs(vector_object, EPSG code) 구문으로 좌표계를 설정할 수 있습니다. 단, 아래의 경고 문구처럼 "st_set_crs() 함수는 좌표계를 변경하는 것이 투영 데이터를 변환하는 것은 아니며, 투영 데이터 변환을 하려면 st_transform() 함수를 이용 (st_crs<- : replacing crs does not reproject data; use st_transform for that)"해야 합니다.

## -- st_set_crs() : setting a CRS (coordinate reference system)

vector_zion_2 = st_set_crs(vector_zion, 4326) # set CRS with EPSG 4326 code
# Warning message:
#   st_crs<- : replacing crs does not reproject data; use st_transform for that

[ 벡터 데이터: 지리 좌표계 (WGS 84; 왼쪽) vs. 투영 좌표계 (NAD83/ UTM zone 12N; 오른쪽) ]

(4-2) 레스터 데이터에서 좌표계 (CRSs in Raster data)

레스터 모델의 객체에 대해서는 raster 패키지의 projection() 함수를 사용해서 좌표계를 확인하거나 설정할 수 있습니다.

아래 예는 R spDataLarge 패키지에 포함되어 있는 strm.tif 파일을 raster() 함수를 사용해 레스터 데이터로 읽어와서, projection() 함수로 좌표계를 확인해본 것입니다.

## -- raster::projection() : get or set CRS in raster* objects
library(raster)

raster_filepath = system.file("raster/srtm.tif", package = "spDataLarge")
new_raster = raster(raster_filepath)

projection(new_raster) # get CRS in raster objects
# [1] "+proj=longlat +datum=WGS84 +no_defs"

레스터 데이터에 대해서 좌표계를 새로 설정할 때도 역시 projection() 함수를 사용합니다. 벡터 데이터의 경우 좌표계를 설정할 때 'EPSG 코드'나 'Proj4string 정의' 모두 사용 가능한 반면에, 레스터 데이터는 'Proj4string 정의'만 사용할 수 있습니다.

## -- set CRS for raster objects
projection(new_raster) = "+proj=utm +zone=12 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 
                            +units=m +no_defs" 
# Warning message:
#   In showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj = prefer_proj) :
#   Discarded datum Unknown based on GRS80 ellipsoid in Proj4 definition

projection(new_raster)
# [1] "+proj=utm +zone=12 +ellps=GRS80 +units=m +no_defs"

[ Reference ]

Geographic data in R - 'CRSs' : geocompr.robinlovelace.net/spatial-class.html

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

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

지난번 포스팅에서는 레스터 데이터의 유형으로 (a) 셀 ID (Cell IDs), (b) 셀 값 (Cell Values) 에 대해서 알아보고, RasterLayer 예제 데이터를 R raster 패키지의 raster() 함수로 불러와서 여러 속성 정보를 살펴보았습니다. (rfriend.tistory.com/605)

이번 포스팅에서는 3가지의 레스터 클래스 (Raster Classes) 의 특장점과 언제 사용하면 좋을지에 대해서 소개하겠습니다.

(1) RasterLayer class
(2) RasterBrick class
(3) RasterStack class

(4) 언제 어떤 레스터 클래스를 사용하는 것이 좋은가?

[ 3가지의 레스터 클래스 (Raster Classes) ]

(1) RasterLayer Class

RasterLayer Class는 레스터 객체 중에서 가장 간단한 형태의 클래스로서, 단 한 개의 층으로 구성되어 있습니다. 

지난번 포스팅에서 소개했던 것처럼, RasterLayer Class 객체를 만드는 가장 쉬운 방법은 기존의 RasterLayer Class 객체 파일을 읽어오는 것입니다.  아래 예에서는 raster 패키지의 raster() 함수를 사용해서 spDataLarge 패키지에 내장되어 있는 srtm.tif 레스터 층 클래스 객체를 읽어와서 raster_layer 라는 이름의 단 한개의 층만을 가진 RasterLayer Class 객체를 만들어보겠습내다 .  nlayers() 함수로 층의 개수를 살펴보면 '1'개 인 것을 확인할 수 있습니다.

## ==========================================
## R Raster Classes
## : RasterLayer, RasterBrick and RasterStack
## ==========================================

library(raster)
library(spDataLarge)

## -- (1) RasterLayer
## : The RasterLayer class represents the simplest form of a raster object, 
##   and consists of only one layer.

## read-in a raster file from disk or from a server
raster_filepath = system.file("raster/srtm.tif", package = "spDataLarge")
raster_filepath
# [1] "/Library/Frameworks/R.framework/Versions/4.0/Resources/library/spDataLarge/raster/srtm.tif"

raster_layer = raster(raster_filepath)


raster_layer
# 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)


## number of layers
nlayers(raster_layer)
# [1] 1

raster 패키지는 rgdal 패키지의 도움을 받아 수많은 드라이버를 지원합니다. 

아래는 raster::writeFormats() 로 확인해 본, raster 패키지가 지원하는 "쓰기 포맷 (Write Formats)" 들 입니다.

## write Formats
raster::writeFormats()
#      name        long_name                                
# [1,] "raster"    "R-raster"                               
# [2,] "SAGA"      "SAGA GIS"                               
# [3,] "IDRISI"    "IDRISI"                                 
# [4,] "IDRISIold" "IDRISI (img/doc)"                       
# [5,] "BIL"       "Band by Line"                           
# [6,] "BSQ"       "Band Sequential"                        
# [7,] "BIP"       "Band by Pixel"                          
# [8,] "ascii"     "Arc ASCII"                              
# [9,] "CDF"       "NetCDF"                                 
# [10,] "ADRG"      "ARC Digitized Raster Graphics"          
# [11,] "BMP"       "MS Windows Device Independent Bitmap"   
# [12,] "BT"        "VTP .bt (Binary Terrain) 1.3 Format"    
# [13,] "BYN"       "Natural Resources Canada's Geoid"       
# [14,] "CTable2"   "CTable2 Datum Grid Shift"               
# [15,] "EHdr"      "ESRI .hdr Labelled"                     
# [16,] "ELAS"      "ELAS"                                   
# [17,] "ENVI"      "ENVI .hdr Labelled"                     
# [18,] "ERS"       "ERMapper .ers Labelled"                 
# [19,] "GPKG"      "GeoPackage"                             
# [20,] "GS7BG"     "Golden Software 7 Binary Grid (.grd)"   
# [21,] "GSBG"      "Golden Software Binary Grid (.grd)"     
# [22,] "GTX"       "NOAA Vertical Datum .GTX"               
# [23,] "GTiff"     "GeoTIFF"                                
# [24,] "HDF4Image" "HDF4 Dataset"                           
# [25,] "HFA"       "Erdas Imagine Images (.img)"            
# [26,] "IDA"       "Image Data and Analysis"                
# [27,] "ILWIS"     "ILWIS Raster Map"                       
# [28,] "INGR"      "Intergraph Raster"                      
# [29,] "ISCE"      "ISCE raster"                            
# [30,] "ISIS2"     "USGS Astrogeology ISIS cube (Version 2)"
# [31,] "ISIS3"     "USGS Astrogeology ISIS cube (Version 3)"
# [32,] "KRO"       "KOLOR Raw"                              
# [33,] "LAN"       "Erdas .LAN/.GIS"                        
# [34,] "Leveller"  "Leveller heightfield"                   
# [35,] "MBTiles"   "MBTiles"                                
# [36,] "MRF"       "Meta Raster Format"                     
# [37,] "NGW"       "NextGIS Web"                            
# [38,] "NITF"      "National Imagery Transmission Format"   
# [39,] "NTv2"      "NTv2 Datum Grid Shift"                  
# [40,] "NWT_GRD"   "Northwood Numeric Grid Format .grd/.tab"
# [41,] "PAux"      "PCI .aux Labelled"                      
# [42,] "PCIDSK"    "PCIDSK Database File"                   
# [43,] "PCRaster"  "PCRaster Raster File"                   
# [44,] "PDF"       "Geospatial PDF"                         
# [45,] "PDS4"      "NASA Planetary Data System 4"           
# [46,] "PNM"       "Portable Pixmap Format (netpbm)"        
# [47,] "RMF"       "Raster Matrix Format"                   
# [48,] "ROI_PAC"   "ROI_PAC raster"                         
# [49,] "RRASTER"   "R Raster"                               
# [50,] "RST"       "Idrisi Raster A.1"                      
# [51,] "SAGA"      "SAGA GIS Binary Grid (.sdat, .sg-grd-z)"
# [52,] "SGI"       "SGI Image File Format 1.0"              
# [53,] "Terragen"  "Terragen heightfield"                   
# [54,] "VICAR"     "MIPL VICAR file"                        
# [55,] "netCDF"    "Network Common Data Format" 

그리고 아래는 rgdal::gdalDrivers() 로 확인해 본, rgdal 패키지에서 지원하는 140 여개의 드라이버들입니다. 드라이버별로 생성(create), 복사(copy) 가능 여부는 아래에서 확인이 가능합니다.

## gdal Drivers
rgdal::gdalDrivers()
#                    name                                                long_name create  copy isRaster
# 1               AAIGrid                                      Arc/Info ASCII Grid  FALSE  TRUE     TRUE
# 2                  ACE2                                                     ACE2  FALSE FALSE     TRUE
# 3                  ADRG                            ARC Digitized Raster Graphics   TRUE FALSE     TRUE
# 4                   AIG                                     Arc/Info Binary Grid  FALSE FALSE     TRUE
# 5                   ARG                                Azavea Raster Grid format  FALSE  TRUE     TRUE
# 6                AirSAR                                AirSAR Polarimetric Image  FALSE FALSE     TRUE
# 7                   BAG                               Bathymetry Attributed Grid  FALSE  TRUE     TRUE
# 8                BIGGIF                       Graphics Interchange Format (.gif)  FALSE FALSE     TRUE
# 9                   BLX                                     Magellan topo (.blx)  FALSE  TRUE     TRUE
# 10                  BMP                     MS Windows Device Independent Bitmap   TRUE FALSE     TRUE
# 11                  BSB                              Maptech BSB Nautical Charts  FALSE FALSE     TRUE
# 12                   BT                      VTP .bt (Binary Terrain) 1.3 Format   TRUE FALSE     TRUE
# 13                  BYN                         Natural Resources Canada's Geoid   TRUE FALSE     TRUE
# 14                  CAD                                           AutoCAD Driver  FALSE FALSE     TRUE
# 15                 CALS                                            CALS (Type 1)  FALSE  TRUE     TRUE
# 16                 CEOS                                               CEOS Image  FALSE FALSE     TRUE
# 17                COASP                          DRDC COASP SAR Processor Raster  FALSE FALSE     TRUE
# 18                  COG                        Cloud optimized GeoTIFF generator  FALSE  TRUE     TRUE
# 19                COSAR               COSAR Annotated Binary Matrix (TerraSAR-X)  FALSE FALSE     TRUE
# 20                  CPG                                          Convair PolGASP  FALSE FALSE     TRUE
# 21                  CTG                           USGS LULC Composite Theme Grid  FALSE FALSE     TRUE
# 22              CTable2                                 CTable2 Datum Grid Shift   TRUE FALSE     TRUE
# 23                 DAAS          Airbus DS Intelligence Data As A Service driver  FALSE FALSE     TRUE
# 24              DERIVED               Derived datasets using VRT pixel functions  FALSE FALSE     TRUE
# 25                DIMAP                                               SPOT DIMAP  FALSE FALSE     TRUE
# 26                DIPEx                                                    DIPEx  FALSE FALSE     TRUE
# 27                 DOQ1                                     USGS DOQ (Old Style)  FALSE FALSE     TRUE
# 28                 DOQ2                                     USGS DOQ (New Style)  FALSE FALSE     TRUE
# 29                 DTED                                    DTED Elevation Raster  FALSE  TRUE     TRUE
# 30              E00GRID                                 Arc/Info Export E00 GRID  FALSE FALSE     TRUE
# 31              ECRGTOC                                          ECRG TOC format  FALSE FALSE     TRUE
# 32                EEDAI                              Earth Engine Data API Image  FALSE FALSE     TRUE
# 33                 EHdr                                       ESRI .hdr Labelled   TRUE  TRUE     TRUE
# 34                  EIR                                        Erdas Imagine Raw  FALSE FALSE     TRUE
# 35                 ELAS                                                     ELAS   TRUE FALSE     TRUE
# 36                 ENVI                                       ENVI .hdr Labelled   TRUE FALSE     TRUE
# 37                  ERS                                   ERMapper .ers Labelled   TRUE FALSE     TRUE
# 38                 ESAT                                     Envisat Image Format  FALSE FALSE     TRUE
# 39                 FAST                                        EOSAT FAST Format  FALSE FALSE     TRUE
# 40                  FIT                                                FIT Image  FALSE  TRUE     TRUE
# 41              FujiBAS                                   Fuji BAS Scanner Image  FALSE FALSE     TRUE
# 42                  GFF Ground-based SAR Applications Testbed File Format (.gff)  FALSE FALSE     TRUE
# 43                  GIF                       Graphics Interchange Format (.gif)  FALSE  TRUE     TRUE
# 44                  GMT                                   GMT NetCDF Grid Format  FALSE  TRUE     TRUE
# 45                 GPKG                                               GeoPackage   TRUE  TRUE     TRUE
# 46       GRASSASCIIGrid                                         GRASS ASCII Grid  FALSE FALSE     TRUE
# 47                 GRIB                             GRIdded Binary (.grb, .grb2)  FALSE  TRUE     TRUE
# 48                GS7BG                     Golden Software 7 Binary Grid (.grd)   TRUE  TRUE     TRUE
# 49                 GSAG                        Golden Software ASCII Grid (.grd)  FALSE  TRUE     TRUE
# 50                 GSBG                       Golden Software Binary Grid (.grd)   TRUE  TRUE     TRUE
# 51                  GSC                                              GSC Geogrid  FALSE FALSE     TRUE
# 52                  GTX                                 NOAA Vertical Datum .GTX   TRUE FALSE     TRUE
# 53                GTiff                                                  GeoTIFF   TRUE  TRUE     TRUE
# 54                  GXF                             GeoSoft Grid Exchange Format  FALSE FALSE     TRUE
# 55               GenBin                           Generic Binary (.hdr Labelled)  FALSE FALSE     TRUE
# 56                 HDF4                       Hierarchical Data Format Release 4  FALSE FALSE     TRUE
# 57            HDF4Image                                             HDF4 Dataset   TRUE FALSE     TRUE
# 58                 HDF5                       Hierarchical Data Format Release 5  FALSE FALSE     TRUE
# 59            HDF5Image                                             HDF5 Dataset  FALSE FALSE     TRUE
# 60                  HF2                               HF2/HFZ heightfield raster  FALSE  TRUE     TRUE
# 61                  HFA                              Erdas Imagine Images (.img)   TRUE  TRUE     TRUE
# 62                 HTTP                                    HTTP Fetching Wrapper  FALSE FALSE     TRUE
# 63                  IDA                                  Image Data and Analysis   TRUE FALSE     TRUE
# 64  IGNFHeightASCIIGrid                  IGN France height correction ASCII Grid  FALSE FALSE     TRUE
# 65                ILWIS                                         ILWIS Raster Map   TRUE  TRUE     TRUE
# 66                 INGR                                        Intergraph Raster   TRUE  TRUE     TRUE
# 67                 IRIS                             IRIS data (.PPI, .CAPPi etc)  FALSE FALSE     TRUE
# 68                 ISCE                                              ISCE raster   TRUE FALSE     TRUE
# 69                  ISG                      International Service for the Geoid  FALSE FALSE     TRUE
# 70                ISIS2                  USGS Astrogeology ISIS cube (Version 2)   TRUE FALSE     TRUE
# 71                ISIS3                  USGS Astrogeology ISIS cube (Version 3)   TRUE  TRUE     TRUE
# 72           JAXAPALSAR               JAXA PALSAR Product Reader (Level 1.1/1.5)  FALSE FALSE     TRUE
# 73                 JDEM                                      Japanese DEM (.mem)  FALSE FALSE     TRUE
# 74          JP2OpenJPEG               JPEG-2000 driver based on OpenJPEG library  FALSE  TRUE     TRUE
# 75                 JPEG                                                JPEG JFIF  FALSE  TRUE     TRUE
# 76      KMLSUPEROVERLAY                                        Kml Super Overlay  FALSE  TRUE     TRUE
# 77                  KRO                                                KOLOR Raw   TRUE FALSE     TRUE
# 78                  L1B                     NOAA Polar Orbiter Level 1b Data Set  FALSE FALSE     TRUE
# 79                  LAN                                          Erdas .LAN/.GIS   TRUE FALSE     TRUE
# 80                  LCP                        FARSITE v.4 Landscape File (.lcp)  FALSE  TRUE     TRUE
# 81               LOSLAS                        NADCON .los/.las Datum Grid Shift  FALSE FALSE     TRUE
# 82             Leveller                                     Leveller heightfield   TRUE FALSE     TRUE
# 83                  MAP                                         OziExplorer .MAP  FALSE FALSE     TRUE
# 84              MBTiles                                                  MBTiles   TRUE  TRUE     TRUE
# 85                  MEM                                         In Memory Raster   TRUE FALSE     TRUE
# 86                  MFF                                        Vexcel MFF Raster   TRUE  TRUE     TRUE
# 87                 MFF2                                 Vexcel MFF2 (HKV) Raster   TRUE  TRUE     TRUE
# 88                  MRF                                       Meta Raster Format   TRUE  TRUE     TRUE
# 89                 MSGN                           EUMETSAT Archive native (.nat)  FALSE FALSE     TRUE
# 90                  NDF                                        NLAPS Data Format  FALSE FALSE     TRUE
# 91             NGSGEOID                              NOAA NGS Geoid Height Grids  FALSE FALSE     TRUE
# 92                  NGW                                              NextGIS Web   TRUE  TRUE     TRUE
# 93                 NITF                     National Imagery Transmission Format   TRUE  TRUE     TRUE
# 94                 NTv1                                    NTv1 Datum Grid Shift  FALSE FALSE     TRUE
# 95                 NTv2                                    NTv2 Datum Grid Shift   TRUE FALSE     TRUE
# 96              NWT_GRC               Northwood Classified Grid Format .grc/.tab  FALSE FALSE     TRUE
# 97              NWT_GRD                  Northwood Numeric Grid Format .grd/.tab   TRUE  TRUE     TRUE
# 98                  OZI                                   OziExplorer Image File  FALSE FALSE     TRUE
# 99                 PAux                                        PCI .aux Labelled   TRUE FALSE     TRUE
# 100              PCIDSK                                     PCIDSK Database File   TRUE FALSE     TRUE
# 101            PCRaster                                     PCRaster Raster File   TRUE  TRUE     TRUE
# 102                 PDF                                           Geospatial PDF   TRUE  TRUE     TRUE
# 103                 PDS                               NASA Planetary Data System  FALSE FALSE     TRUE
# 104                PDS4                             NASA Planetary Data System 4   TRUE  TRUE     TRUE
# 105            PLMOSAIC                                  Planet Labs Mosaics API  FALSE FALSE     TRUE
# 106            PLSCENES                                   Planet Labs Scenes API  FALSE FALSE     TRUE
# 107                 PNG                                Portable Network Graphics  FALSE  TRUE     TRUE
# 108                 PNM                          Portable Pixmap Format (netpbm)   TRUE FALSE     TRUE
# 109                 PRF                                      Racurs PHOTOMOD PRF  FALSE FALSE     TRUE
# 110       PostGISRaster                                    PostGIS Raster driver  FALSE  TRUE     TRUE
# 111                   R                                      R Object Data Store  FALSE  TRUE     TRUE
# 112                 RDA                   DigitalGlobe Raster Data Access driver  FALSE FALSE     TRUE
# 113                 RIK                                  Swedish Grid RIK (.rik)  FALSE FALSE     TRUE
# 114                 RMF                                     Raster Matrix Format   TRUE FALSE     TRUE
# 115             ROI_PAC                                           ROI_PAC raster   TRUE FALSE     TRUE
# 116              RPFTOC                         Raster Product Format TOC format  FALSE FALSE     TRUE
# 117             RRASTER                                                 R Raster   TRUE  TRUE     TRUE
# 118                 RS2                                   RadarSat 2 XML Product  FALSE FALSE     TRUE
# 119                 RST                                        Idrisi Raster A.1   TRUE  TRUE     TRUE
# 120          Rasterlite                                               Rasterlite  FALSE  TRUE     TRUE
# 121                SAFE                              Sentinel-1 SAR SAFE Product  FALSE FALSE     TRUE
# 122                SAGA                  SAGA GIS Binary Grid (.sdat, .sg-grd-z)   TRUE  TRUE     TRUE
# 123            SAR_CEOS                                           CEOS SAR Image  FALSE FALSE     TRUE
# 124                SDTS                                              SDTS Raster  FALSE FALSE     TRUE
# 125           SENTINEL2                                               Sentinel 2  FALSE FALSE     TRUE
# 126                 SGI                                SGI Image File Format 1.0   TRUE FALSE     TRUE
# 127              SIGDEM                       Scaled Integer Gridded DEM .sigdem  FALSE  TRUE     TRUE
# 128              SNODAS                            Snow Data Assimilation System  FALSE FALSE     TRUE
# 129                 SRP                      Standard Raster Product (ASRP/USRP)  FALSE FALSE     TRUE
# 130             SRTMHGT                                      SRTMHGT File Format  FALSE  TRUE     TRUE
# 131                 TIL                                          EarthWatch .TIL  FALSE FALSE     TRUE
# 132                 TSX                                       TerraSAR-X Product  FALSE FALSE     TRUE
# 133            Terragen                                     Terragen heightfield   TRUE FALSE     TRUE
# 134             USGSDEM                       USGS Optional ASCII DEM (and CDED)  FALSE  TRUE     TRUE
# 135               VICAR                                          MIPL VICAR file   TRUE  TRUE     TRUE
# 136                 VRT                                           Virtual Raster   TRUE  TRUE     TRUE
# 137                 WCS                                 OGC Web Coverage Service  FALSE FALSE     TRUE
# 138                WEBP                                                     WEBP  FALSE  TRUE     TRUE
# 139                 WMS                                      OGC Web Map Service  FALSE  TRUE     TRUE
# 140                WMTS                                 OGC Web Map Tile Service  FALSE  TRUE     TRUE
# 141                 XPM                                        X11 PixMap Format  FALSE  TRUE     TRUE
# 142                 XYZ                                        ASCII Gridded XYZ  FALSE  TRUE     TRUE
# 143                ZMap                                           ZMap Plus Grid  FALSE  TRUE     TRUE
# 144              netCDF                               Network Common Data Format   TRUE  TRUE     TRUE

RasterLayer 클래스 객체를 raster() 함수를 사용해서 처음부터 직접 만들 수도 있습니다.

아래의 예는 8개의 행과 8개의 열을 가진, 총 64 개의 셀(or 픽셀)을 가진 RasterLayer 클래스를 직접 만들어본 것입니다. 이때 레스터 객체의 좌표 참조 시스템(CRS, Coordinates Reference System)은 WGS84 가 기본 설정값이며, 이는 해상도(resolution)의 단위가 도 (in degrees) 라는 의미입니다. 아래 예에서는 res = 0.5 로서 해상도를 0.5도로 설정해주었습니다. 각 셀의 값은 왼쪽 상단부터 시작하여, 행 방향(row-wise)으로 왼쪽에서 오른쪽으로 채워나가게 됩니다. 아래 예는 총 64개 셀의 각 셀의 값을 1~64 의 정수값을 좌측 상단에서 시작하여 우측으로 한 줄씩 채워나가도록 했으며, 제일 아래에 8*8 셀들의 값과 오른쪽에 plot()으로 시각화한 결과를 확인해보시기 바랍니다.

## -- Rasters can also be created from scratch using the raster() function.
## RasterLayer object, The CRS is the default of raster objects: WGS84.
## raster of 64 cells (pixels) = 8 rows * 8 cols (row-wise)
my_raster = raster(nrows = 8, ncols = 8, res = 0.5, 
                   xmn = -2.0, xmx = 2.0, ymn = -2.0, ymx = 2.0,
                   vals = 1:64)


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


## plotting
plot(my_raster, main = "my raster (64 cells = 8 rows * 8 cols)")

(2) RasterBrick Class

위 (1) 의 RasterLayer 클래스가 단지 1개 층 (only 1 layer) 으로만 구성되는 가장 간단한 형태의 레스터 클래스라고 소개했는데요, (2) RasterBrick 클래스와 (3) RasterStack 클래스는 여러개의 층(multiple layers)을 가질 수 있습니다.

특히, RasterBrick 클래스는 단일 다중 스펙트럼 위성 파일 (a single multispectral satellite file) 이나 또는 메모리의 단일 다층 객체 (a single multilayer object in memory) 의 형태로 다층의 레스터 객체를 구성합니다.

아래의 예는 raster 패키지의 brick() 함수를 사용해서 spDataLarge 패키지에 들어있는 landsat.tif 의 다층 레스터 파일을 RasterBrick 클래스 객체로 불러온 것입니다. 

nlayers() 함수로 총 몇 개의 층이 있는지 확인해 보니 4개의 층이 하나의 파일에 들어있네요.

plot() 함수로 RasterBrick 클래스 객체의 4개 층을 모두 시각화보면 아래와 같습니다.

## -- (2) RasterBrick
## : A RasterBrick consists of multiple layers, 
##  which typically correspond to a single multispectral satellite file 
##  or a single multilayer object in memory. 
multi_raster_file = system.file("raster/landsat.tif", package = "spDataLarge")
raster_brick = brick(multi_raster_file)

## 4 layers
raster_brick
# class      : RasterBrick 
# dimensions : 1428, 1128, 1610784, 4  (nrow, ncol, ncell, nlayers)
# resolution : 30, 30  (x, y)
# extent     : 301905, 335745, 4111245, 4154085  (xmin, xmax, ymin, ymax)
# crs        : +proj=utm +zone=12 +datum=WGS84 +units=m +no_defs 
# source     : /Library/Frameworks/R.framework/Versions/4.0/Resources/library/spDataLarge/raster/landsat.tif 
# names      : landsat.1, landsat.2, landsat.3, landsat.4 
# min values :      7550,      6404,      5678,      5252 
# max values :     19071,     22051,     25780,     31961


## nlayers() : number of layers in a Raster* object
nlayers(raster_brick)
# [1] 4


## plotting RasterBrick object with 4 layers
plot(raster_brick)

(3) RasterStack Class

세번째 레스터 클래스인 RasterStack 클래스도 다 층 (multi-layers) 레스터 객체로 구성이 됩니다. 같은 범위와 해상도를 가진 여러개의 RasterLayer 클래스 객체들을 리스트로 묶어서 RasterStack 클래스 객체를 만들 수 있습니다.

이때, 위 (2)번의 RasterBrick 클래스가 동일한 복수개의 RasterLayer 층으로 구성되는 반면에, 이번 (3)번의 RasterStack 클래스는 여러개의 RasterLayer와 RasterBrick 클래스 객체가 혼합되어서 구성할 수 있다는 점입니다.  연산 속도면에서 보면 일반적으로 RasterBrick 클래스가 RasterStack 클래스보다 빠릅니다.

RasterBrick 클래스와 RasterStack 객체에 대한 연산은 보통은 RasterBrack 클래스 객체를 반환합니다.

아래 예에서는 (a) raster(raster_brick, layer = 1) 함수를 사용해서 위의 (2)번에서 불러왔던 RasterBrick 클래스 객체의 1번째 층만 가져다가 raster_on_disk 라는 이름으로 레스터 객체를 하나 만들고, (b) raster() 함수로 동일한 범위와 해상도, 좌표 참조 시스템(CRS)를 가지고 난수로 셀의 값을 채운 raster_in_memory 라는 이름의 메모리에 있는 RasterLayer 클래스 객체를 만들었습니다. 다음에 stac() 함수로 raster_stack = stack(raster_in_memory, raster_on_disk) 처럼 (a) + (b) 하여 쌓아서 raster_stack 라는 이름의 RasterStack 클래스 객체를 만들었습니다. 

마지막으로 plot() 함수로 RasterStack 클래스 객체에 쌓여 있는 2개의 객체를 시각화해보았습니다. (raster_in_memory 는 난수를 발생시켜 셀 값을 채웠기 때문에 시각화했을 때 아무런 패턴이 없습니다.)

## -- (3) RasterStack
## : RasterStack allows you to connect several raster objects 
##   stored in different files or multiple objects in memory.

raster_on_disk = raster(raster_brick, layer = 1)

raster_in_memory = raster(xmn = 301905, xmx = 335745,
                          ymn = 4111245, ymx = 4154085, 
                          res = 30)
values(raster_in_memory) = sample(seq_len(ncell(raster_in_memory)))
crs(raster_in_memory) = crs(raster_on_disk)


## RasterStack is a list of RasterLayer objects with the same extent and resolution.
raster_stack = stack(raster_in_memory, raster_on_disk)
raster_stack
# dimensions : 1428, 1128, 1610784, 2  (nrow, ncol, ncell, nlayers)
# resolution : 30, 30  (x, y)
# extent     : 301905, 335745, 4111245, 4154085  (xmin, xmax, ymin, ymax)
# crs        : +proj=utm +zone=12 +datum=WGS84 +units=m +no_defs 
# names      :   layer, landsat.1 
# min values :       1,      7550 
# max values : 1610784,     19071


## plotting RasterStack (raster_in_memory + raster_on_disk)
plot(raster_stack)

(4) 언제 어떤 레스터 클래스를 사용하는 것이 좋은가?

위의 (1), (2), (3)에서 소개한 RasterLayer 클래스, RasterBrick 클래스, RasterStack 클래스의 특징과 서로 간의 차이점을 보면 알겠지만 다시 한번 정리하자면요, Raster* 클래스 중에서 무엇을 사용할 지는 투입 데이터 (input data)의 특징이 무엇이냐에 따라 달라집니다.

하나의 다층 레스터 파일이나 객체(a single multilayer file or object)를 처리하는 것이라면 RasterBrick 이 적합합니다.

반면에, 여러개의 레스터 파일들(many files)이나 여러 종류의 레스터 클래스 (many Raster*) 를 한꺼번에 연결해서 연산하고 처리해야 하는 경우라면 RasterStack Class 가 적합하겠습니다.

[ Reference ]

* Geographic data in R : geocompr.robinlovelace.net/spatial-class.html

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

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

지난번 포스팅에서는 R data.table에서 (a) 키(Key)와 빠른 이진 탐색 기반의 Subsetting 하는 방법 (rfriend.tistory.com/569), (b) 2차 인덱스 (secondary indices) 를 활용하여 data.table 의 재정렬 없이 빠른 탐색 기반 Subsetting 하는 방법을 소개하였습니다. (rfriend.tistory.com/615)

이번 포스팅에서는 R data.table에서 이진 연산자(binary operators) 인 '=='와 '%in%' 를 수행하는 과정에서 이차 인덱스(secondary indices)가 자동으로 인덱싱(Auto indexing)이 되어 빠르게 subsetting 하는 내용을 소개하겠습니다. (이 글을 쓰는 2021년 2월 현재는 '=='와 '%in%' 연산자만 자동 인덱싱이 지원되며, 향후 더 많은 연산자로 확대 전망)

(1) '==' 이진 연산자로 자동 인덱싱하고 속도 비교하기

(2) '%in%' 이진 연산자로 자동 인덱싱하고 속도 비교하기

(3) 전역으로 자동 인덱싱을 비활성화하기 (disable auto indexing globally)

자동 인덱싱의 속도 개선 효과를 확인해 보기 위해서 천만개의 행을 가진 예제 data.table을 난수를 발생시켜서 생성해 보겠습니다.  DT data.table의 크기를 object.size()로 재어보니 114.4 Mb 이네요.

## =========================
## R data.table
## : Auto indexing
## =========================

library(data.table)

## create a data.table big enough
set.seed(1L)
DT = data.table(x = sample(x = 1e5L, 
                           size = 1e7L, 
                           replace = TRUE), 
                y = runif(100L))

head(DT)
#    x         y
# 1: 24388 0.4023457
# 2: 59521 0.9142361
# 3: 43307 0.2847435
# 4: 69586 0.3440578
# 5: 11571 0.1822614
# 6: 25173 0.8130521

dim(DT)
# [1] 10000000        2

print(object.size(DT), units = "Mb")
# 114.4 Mb

(1) '==' 이진 연산자로 자동 인덱싱하고 속도 비교하기

이전 포스팅의 이차 인덱스(secondary index)에서는 setindex(DT, column) 으로 이차 인덱스를 명시적으로 설정하거나, 또는 'on' 매개변수로 subsetting을 하면 실행 중에 (on the fly) 기존 이차 인덱스가 있는지 여부를 확인해서, 없으면 바로 이차 인덱스를 설정해주다고 하였습니다.

R data.table에서 '==' 이진 연산자를 사용해서 행의 부분집합을 가져오기(subsetting)을 하면 기존 이차 인덱스가 없을 경우 자동으로 인덱싱을 해줍니다. 그래서 처음에 '=='로 subsetting 할 때는 (a) 인덱스를 생성하고 + (b) 부분집합 행 가져오기 (subsetting)를 수행하느라 시간이 오래 소요되지만, 두번째로 실행할 때는 인덱스가 생성이 되어 있으므로 속도가 무척 빨라지게 됩니다!

아래의 예에서 보면 처음으로 DT[x == 500L] 을 실행했을 때는 0.406초가 소요(elapsed time)되었습니다. names(attributes(DT)) 로 확인해 보면 애초에 없던 index 가 새로 생성되었음을 확인할 수 있고, indices(DT) 로 확인해보면 "x" 칼럼에 대해 이차 인덱스가 생성되었네요.

## -- when we use '==' or '%in%' on a single column for the first time, 
## a secondary index is created automatically, and used to perform the subset. 

## have a look at all the attribute names (no index here)
names(attributes(DT))
# [1] "names"             "row.names"         "class"             ".internal.selfref"

## run the first time
## system.time = the time to create the index + the time to subset
(t1 <- system.time(ans <- DT[x == 500L]))
#  user  system elapsed 
# 0.392   0.014   0.406

head(ans)
#    x         y
# 1: 500 0.7845248
# 2: 500 0.9612705
# 3: 500 0.4023457
# 4: 500 0.9139429
# 5: 500 0.8280599
# 6: 500 0.2847435


## secondary index is created
names(attributes(DT))
# [1] "names"             "row.names"         "class"             ".internal.selfref" 
# [5] "index"

indices(DT)
# [1] "x"

이제 위에서 수행했던 연산과 동일하게 DT[x == 500L] 을 수행해서 소요 시간(elapsed time)을 측정해보면, 연속해서 두번째 수행했을 때는 0.001 초가 걸렸습니다.

## secondary indices are extremely fast in successive subsets. 
## successive subsets
(t2 <- system.time(DT[x == 500L]))
#  user  system elapsed 
# 0.001   0.000   0.001

처음 수행했을 때는 0.406초가 걸렸던 것이, 처음 수행할 때 자동 인덱싱(auto indexing)이 된 후에 연속해서 수행했을 때 0.001초가 걸려서 400배 이상 빨라졌습니다!  와우!!!

barplot(c(0.406, 0.001), 
        horiz = TRUE, 
        xlab = "elapsed time",
        col = c("red", "blue"),
        legend.text = c("first time", "second time(auto indexing)"), 
        main = "R data.table Auto Indexing")
        

(2) '%in%' 이진 연산자로 자동 인덱싱하고 속도 비교하기

'==' 연산자와 더불어 포함되어 있는지 여부를 확인해서 블리언을 반환하는 '%in%' 연산자를 활용해서 부분집합 행을 가져올 때도 R data.table은 자동 인덱싱(auto indexing)을 하여 이차 인덱스를 생성하고, 기존에 인덱스가 생성되어 있으면 이차 인덱스를 활용하여 빠르게 탐색하고 subsetting 결과를 반환합니다.

아래 예는 x 에 1989~2912 까지의 정수가 포함되어 있는 행을 부분집합으로 가져오기(DT[ x %in% 1989:2912]) 하는 것으로서, 이때 자동으로 인덱스를 생성(auto indexing)해 줍니다.

## '%in%' operator create auto indexing as well
system.time(DT[x %in% 1989:2912])
#  user  system elapsed 
# 0.010   0.016   0.027

행을 subsetting 할 때 사용하는 조건절이 여러개의 칼럼을 대상으로 하는 경우 '&' 연산자를 사용하여 자동 인덱싱을 할 수 있습니다.

## auto indexing to expressions involving more than one column with '&' operator
(t3 <- system.time(DT[x == 500L & y >= 0.5]))
#  user  system elapsed 
# 0.070   0.025   0.097

(3) 전역으로 자동 인덱싱을 비활성화하기 (disable auto indexing globally)

지난번 포스팅에서 지역적으로 특정 칼럼의 이차 인덱스를 제거할 때 setindex(DT, NULL) 을 사용한다고 소개하였습니다. 

(a) '전역적으로 자동 인덱싱을 비활성화' 하려면 options(datatable.auto.index = FALSE) 를 설정해주면 됩니다.

(b) '전역으로 전체 인덱스를 비활성화' 하려면 options(datatable.use.index = FALSE) 를 설정해주면 됩니다.

## Auto indexing can be disabled by setting the global argument 
options(datatable.auto.index = FALSE)

## You can disable indices fully by setting global argument
options(datatable.use.index = FALSE)

 [ Reference ]

* R data.table vignettes 'Secondary indices and auto indexing'
:  cran.r-project.org/web/packages/data.table/vignettes/datatable-secondary-indices-and-auto-indexing.html

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

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

이전 포스팅에서는 R data.table에서 '키와 빠른 이진 탐색 기반의 부분집합 선택 (Key and fast binary search based subset)' 에 대해서 소개하였습니다. (rfriend.tistory.com/569)

이번 포스팅에서는 R data.table에서 '2차 인덱스 (Secondary indices)'를 사용하여 빠른 이진 탐색 기반의 부분집합 가져오기 방법을 소개하겠습니다.  이번 포스팅은 R data.table vignettes 을 참조하였습니다.

(1) 이차 인덱스 (Secondary indices) 는 무엇이, 키(Key)와는 무엇이 다른가?

(2) 이차 인덱스를 설정하고 확인하는 방법

(3) 'on' 매개변수와 이차 인덱스를 사용해서 빠르게 부분집합 가져오기

(4) Chaining 해서 정렬하기

(5) j 에 대해 계산하기 (compute or do in j)

(6) J 에 := 를 사용해서 참조하여 부분할당하기 (sub-assign by reference using := in j)

(7) by 를 사용해서 집계하기 (Aggregation using by)

(8) mult 매개변수를 사용해 첫번째 행, 마지막행 가져오기

(9) 이차 인덱스 제거하기 (remove all secondary indices)

(1) 이차 인덱스는 무엇이고, 키와는 무엇이 다른가?
    (Key vs. Secondary indices)

이차 인덱스(Secondary indices) 는 data.table의 키(Key)와 비슷하게 빠른 이진 탐색 기반의 부분집합 가져오기를 할 때 사용합니다.

하지만 키(Key)가 (a) 순서 벡터를 계산한 다음에, (b) 물리적으로 data.table을 재정렬(physically reordering)하는데 비해, 이차 인덱스(secondary indices)는 순서 벡터를 계산해서 index 를 생성해 속성으로 저장만 하고, 물리적 재정렬은 하지 않는 차이점이 있습니다. 만약 data.table의 행과 열이 매우 많은 큰 크기의 데이터셋이라면 물리적으로 재정렬하는데 많은 처리 비용과 시간이 소요될 것입니다.

또 하나 큰 차이점은, 키(Key)는 하나의 칼럼을 키로 설정했을 때 다른 칼럼을 키로 사용하려면 키를 새로운 칼럼으로 재설정하고 data.table을 물리적으로 재정렬을 해야 하는 반면에, 이차 인덱스(secondary indices)는 복수의 이차 인덱스를 설정하고 재사용할 수 있습니다.

이러한 차이점을 고려했을 때, 키(Key)는 동일 칼럼을 "반복적으로" 사용해서 빠르게 부분집합 가져오기(subsetting)를 해야 할 때 유리하며, 이차 인덱스(Secondary indices)는 복수개의 칼럼을 단발성으로 사용하면서 빠르게 부분집합 가져오기를 해야 하는 경우에 유리합니다.

[ R data.table 키(Key) vs. 이차적인 인덱스 (Secondary indices) ]

(2) 이차 인덱스를 설정하고 확인하는 방법

R data.table 패키지를 importing 하고, 예제로 사용할 데이터로는 Lahman 패키지에 들어있는 투수의 투구 통계 데이터인 "Pitching"을 참조하여 Data.Table로 불러오겠습니다.

library(data.table)

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

## coerce lists and data.frame to data.table by reference
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>

이차 인덱스는 setindex(DT, column) 함수의 구문으로 설정할 수 있습니다. 그러면 순서 벡터를 계산해서 내부에 index 라는 속성(attribute)을 생성해서 저장하며, 물리적으로 data.table을 재정렬하는 것은 하지 않습니다(no physical reordering). 

names(attributes(DT)) 으로 확인해보면 제일 마지막에 "index"라는 속성이 추가된 것을 알 수 있습니다.  indices(DT) 함수를 사용하면 모든 이차 인덱스의 리스트를 얻을 수 있습니다.  이때 만약 아무런 이차 인덱스가 설정되어 있지 않다면 NULL 을 반환합니다.

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


## alternatively we can provide character vectors to the function 'setindexv()'
# setindexv(Pitching, "teamID") # useful to program with

## 'index' attribute added
names(attributes(Pitching))
# [1] "names"             "row.names"         "class"             ".internal.selfref" 
# [5] "index"


## get all the secondary indices set
indices(Pitching)
# [1] "teamID"

(3) 'on' 매개변수와 이차 인덱스를 사용해서 빠르게 부분집합 가져오기

'on' 매개변수를 사용하면 별도로 setindex()로 매번 이차 인덱스를 설정하는 절차 없이, 바로 실행 중에(on the fly) 이차 인덱스를 계산해서 부분집합 가져오기(subsetting)을 할 수 있습니다. 

그리고 만약 기존이 이미 이차 인덱스가 설정이 되어 있다면 속성을 확인하여 존재하는 이차 인덱스를 재활용해서 부분집합 가져오기를 빠르게 할 수 있습니다 (on 매개변수는 Key에 대해서도 동일하게 작동합니다).

또 'on' 매개변수는 무슨 칼럼을 기준으로 subsetting 이 실행될지에 대해서 명확하게 코드 구문으로 확인할 수 있게 해주어 코드 가독성을 높여줍니다.

아래 예제는 Pitching data.table에서 이차 인덱스(secondary indices)를 설정한 'teamID' 칼럼의 값이 "NY2" 인 팀을 subsetting 해서 가져온 것입니다. (칼럼 개수가 너무 많아서 1~10번까지 칼럼만 가져왔습니다. [, 1:10]) 

Pirthcing["NY2", on = "teamID"], Pitching[.("NY2"), on = "teamID"], Pitching[list("NY2"), on = "teamID"] 모두 동일한 결과를 반환합니다.

## subset all rows where the teamID matches "NY2" using 'on'
Pitching["NY2", on = "teamID"][,1:10]
#    playerID yearID stint teamID lgID  W  L  G GS CG
# 1: fergubo01   1871     1    NY2   NA  0  0  1  0  0
# 2: fleetfr01   1871     1    NY2   NA  0  1  1  1  1
# 3: woltery01   1871     1    NY2   NA 16 16 32 32 31
# 4: cummica01   1872     1    NY2   NA 33 20 55 55 53
# 5: mcmuljo01   1872     1    NY2   NA  1  0  3  1  1
# 6: martiph01   1873     1    NY2   NA  0  1  6  1  1
# 7: mathebo01   1873     1    NY2   NA 29 23 52 52 47
# 8: hatfijo01   1874     1    NY2   NA  0  1  3  0  0
# 9: mathebo01   1874     1    NY2   NA 42 22 65 65 62
# 10: gedneco01   1875     1    NY2   NA  1  0  2  1  1
# 11: mathebo01   1875     1    NY2   NA 29 38 70 70 69


## or alternatively
# Pitching[.("NY2"), on = "teamID"]
# Pitching[list("NY2"), on = "teamID"]

복수개의 이차 인덱스 (multiple secondary indices)를 setindex(DT, col_1, col_2, ...) 구문 형식으로 설정할 수도 있습니다.

아래 예에서는 Pitching data.table에 "teamID", "yearID"의 2개 칼럼을 이차 인덱스로 설정하고, teamID가 "NY2", yearID가 1873 인 행을 subsetting 해본 것입니다.

## set multiple secondary indices
setindex(Pitching, teamID, yearID)
indices(Pitching)
# [1] "teamID"         "teamID__yearID"


## subset based on teamID and yearID columns.
Pitching[.("NY2", 1873), # i
         on = c("teamID", "yearID")]
#    playerID yearID stint teamID lgID  W  L  G GS CG SHO SV IPouts   H  ER HR BB SO BAOpp  ERA IBB WP HBP BK  BFP GF   R SH SF
# 1: martiph01   1873     1    NY2   NA  0  1  6  1  1   0  0    102  50  13  0  6  1    NA 3.44  NA  1  NA  0  177  5  37 NA NA
# 2: mathebo01   1873     1    NY2   NA 29 23 52 52 47   2  0   1329 489 127  5 62 79    NA 2.58  NA 23  NA  0 2008  0 348 NA NA
# GIDP
# 1:   NA
# 2:   NA

이차 인덱스도 DT[i, j, by] 의 구문 형식을 그대로 따르므로 이차 인덱스로 i 에 대해 행을 subsetting 하고, j 에 대해서 특정 칼럼들을 선택해서 가져올 수 있습니다.

아래 예에서는 이차 인덱스인 teamID가 "NY2", yearID가 1873인 행을 subsetting하고, j 부분에 .(teamID, yearID, playerID, W, L) 로 지정해줘서 칼럼은 teamID, yearID, playerID, W, L 만 선별적으로 선택해서 가져온 것입니다.

## -- select in j
## return palyerID, W, L columns as a data.table corresponding to teamID = "NY2" and yearID = 1873
Pitching[.("NY2", 1873), # i
         .(teamID, yearID, playerID, W, L), # j
         on = c("teamID", "yearID")] # secondary indices
#    teamID yearID  playerID  W  L
# 1:    NY2   1873 martiph01  0  1
# 2:    NY2   1873 mathebo01 29 23

(4) Chaining 해서 정렬하기

이차 인덱스를 사용해서 subsetting 한 후의 결과에 DT[i, j, by][order()] 처럼 chaining을 해서 특정 칼럼을 기준으로 정렬을 할 수 있습니다.

아래 예에서는 이차 인덱스 'teamID' 의 값이 "NY2" 인 행을 subsetting 하고, 칼럼은 .(teamID, yearID, playerID, W, L) 만 선별해서 가져오는데, 이 결과에 chaining을 해서 [order(-W)] 로 W (승리 회수) 를 기준으로 내림차순 정렬 (sorting in descending order) 을 해본 것입니다.  order(-W) 에서 마이너스 부호('-')는 내림차순 정렬을 하라는 의미입니다. (order()의 기본설정은 오름차순 정렬임)

## -- Chaining
## use chaining to order the W column in descending order
Pitching[.("NY2"), # i
         .(teamID, yearID, playerID, W, L), # j
         on = c("teamID")][ # secondary indices
           order(-W)] # order by W in decreasing order
#    teamID yearID  playerID  W  L
# 1:    NY2   1874 mathebo01 42 22
# 2:    NY2   1872 cummica01 33 20
# 3:    NY2   1873 mathebo01 29 23
# 4:    NY2   1875 mathebo01 29 38
# 5:    NY2   1871 woltery01 16 16
# 6:    NY2   1872 mcmuljo01  1  0
# 7:    NY2   1875 gedneco01  1  0
# 8:    NY2   1871 fergubo01  0  0
# 9:    NY2   1871 fleetfr01  0  1
# 10:    NY2   1873 martiph01  0  1
# 11:    NY2   1874 hatfijo01  0  1

(5) j 에 대해 계산하기 (compute or do in j)

이차 인덱스로 i 행을 Subsetting 한 다음에 j 열에 대해서 연산을 할 수 있습니다.

아래 예에서는 (a) 이차 인덱스 'teamID' 의 값이 "NY2"인 행을 subsetting 한 후에, 그 결과 안에서 W (승리회수) 의 최대값을 계산, (b) 복수의 이차 인덱스 'teamID', 'yearID'의 값이 각각 "NY2", 1873인 값을 subsetting 해서 W의 값의 최대값을 계산(max(W))한 것입니다.

## -- Compute or do in j

## Find the maximum W corresponding to teamID="NY2"
Pitching[.("NY2"), max(W), on = c("teamID")]
# [1] 42


Pitching[.("NY2", 1873), max(W), on = c("teamID", "yearID")]
# [1] 29

(6) j 에 := 를 사용해서 참조하여 부분할당하기
    (sub-assign by reference using := in j)

DT[i, j, by] 에서 j 부분에 := 사용해 'on'으로 이차 인덱스를 참조하여 부분 할당(sub-assign) 하면 매우 빠르게 특정 일부분의 행의 값만을 대체할 수 있습니다. 

만약 행의 개수가 매우 많은 데이터셋에서 Key() 를 사용해서 참조하여 부분할당을 하려고 한다면 data.table에 대한 물리적인 재정렬(physical reordering)이 발생하여 연산비용과 시간이 많이 소요될텐데요, 이를 이차 인덱스(secondary indices)를 사용하면 data.table에 대한 재정렬 없이 일부 행의 값을 다른 값으로 대체하는 일을 빠르게 할 수 있는 장점이 있습니다.

아래의 예는 이차 인덱스인 yearID 의 값이 '2019' 인 행의 값을 '2020' 으로 대체하는 부분할당을 해본 것입니다. (2019년을 2020년으로 바꾼 것은 별 의미는 없구요, 그냥 이차 인덱스 참조에 의한 부분할당 기능 예시를 들어본 것입니다.)

## -- sub-assign by reference using := in j
## get all yearID in Pitching
Pitching[, sort(unique(yearID))]
# [1] 1871 1872 1873 1874 1875 1876 1877 1878 1879 1880 1881 1882 1883 1884 1885 1886 1887 1888 1889 1890 1891 1892 1893 1894 1895
# [26] 1896 1897 1898 1899 1900 1901 1902 1903 1904 1905 1906 1907 1908 1909 1910 1911 1912 1913 1914 1915 1916 1917 1918 1919 1920
# [51] 1921 1922 1923 1924 1925 1926 1927 1928 1929 1930 1931 1932 1933 1934 1935 1936 1937 1938 1939 1940 1941 1942 1943 1944 1945
# [76] 1946 1947 1948 1949 1950 1951 1952 1953 1954 1955 1956 1957 1958 1959 1960 1961 1962 1963 1964 1965 1966 1967 1968 1969 1970
# [101] 1971 1972 1973 1974 1975 1976 1977 1978 1979 1980 1981 1982 1983 1984 1985 1986 1987 1988 1989 1990 1991 1992 1993 1994 1995
# [126] 1996 1997 1998 1999 2000 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 2016 2017 2018 2019

## replace 2019 with 2020 using on instead of setting keys
Pitching[.(2019L), yearID := 2020L, on = "yearID"] # no reordering
Pitching[, sort(unique(yearID))]
# [1] 1871 1872 1873 1874 1875 1876 1877 1878 1879 1880 1881 1882 1883 1884 1885 1886 1887 1888 1889 1890 1891 1892 1893 1894 1895
# [26] 1896 1897 1898 1899 1900 1901 1902 1903 1904 1905 1906 1907 1908 1909 1910 1911 1912 1913 1914 1915 1916 1917 1918 1919 1920
# [51] 1921 1922 1923 1924 1925 1926 1927 1928 1929 1930 1931 1932 1933 1934 1935 1936 1937 1938 1939 1940 1941 1942 1943 1944 1945
# [76] 1946 1947 1948 1949 1950 1951 1952 1953 1954 1955 1956 1957 1958 1959 1960 1961 1962 1963 1964 1965 1966 1967 1968 1969 1970
# [101] 1971 1972 1973 1974 1975 1976 1977 1978 1979 1980 1981 1982 1983 1984 1985 1986 1987 1988 1989 1990 1991 1992 1993 1994 1995
# [126] 1996 1997 1998 1999 2000 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 2016 2017 2018 2020

(7) by 를 사용해서 집계하기 (Aggregation using by)

만약 'on' 매개변수로 이차 인덱스를 사용해 "그룹별로 집계나 연산"을 하고 싶다면 by 를 추가해주면 됩니다.

아래 예에서는 이차 인덱스 'teamID'의 값이 "NY2"인 팀을 subsetting 해서, keyby = yearID를 사용해 연도(yearID) 그룹 별로 나누어서 승리회수(W)의 최대값을 계산한 것입니다.

## -- aggregation using by
## get the maximum W for each yearID corresponding to teamID="NY2". order the result by yearID
Pitching[.("NY2"), # i
         max(W),   # j
         keyby = yearID, # order by  
         on = "teamID"]  # secondary indices
#    yearID V1
# 1:   1871 16
# 2:   1872 33
# 3:   1873 29
# 4:   1874 42
# 5:   1875 29

(8) mult 매개변수를 사용해 첫번째 행, 마지막행 가져오기

이차 인덱스(secondary indices)로 빠르게 탐색하여 참조해 행을 subsetting을 해 온 다음에, mult = "first" 매개변수를 사용해서 첫번째 행, 또는 mult = "last"로 마지막 행만을 반환할 수 있습니다.

## -- melt argument
## subset only the first matching row where teamID matches "NY2" and "WS3"
Pitching[c("NY2", "WS3"), on = "teamID", 
         mult = "first"] # subset the first matching row
#   playerID yearID stint teamID lgID  W  L  G GS CG SHO SV IPouts   H  ER HR BB SO BAOpp  ERA IBB WP HBP BK  BFP GF   R SH SF
# 1: fergubo01   1871     1    NY2   NA  0  0  1  0  0   0  0      3   8   3  0  0  0    NA 27.0  NA  2  NA  0   14  0   9 NA NA
# 2: brainas01   1871     1    WS3   NA 12 15 30 30 30   0  0    792 361 132  4 37 13    NA  4.5  NA  7  NA  0 1291  0 292 NA NA
# GIDP
# 1:   NA
# 2:   NA

이차 인덱스로 참조할 기준이 많아지다 보면 그 조건들에 해당하는 행의 값이 존재하지 않을 때도 있습니다.  아래 예의 경우 이차 인덱스 teamID 가 "WS3" 이고 yearID가 '1873'인 행이 존재하지 않아서 mult = "last"로 마지막 을 반환하라고 했을 때 NA 가 반환되었습니다.(두번째 행)

## subset only the last matching row where teamID matches "NY2", "WS3" and yearID matches 1873
Pitching[.(c("NY2", "WS3"), 1873), on = c("teamID", "yearID"), 
         mult = "last"] # subset the last matching row
#    playerID yearID stint teamID lgID  W  L  G GS CG SHO SV IPouts   H  ER HR BB SO BAOpp  ERA IBB WP HBP BK  BFP GF   R SH SF
# 1: mathebo01   1873     1    NY2   NA 29 23 52 52 47   2  0   1329 489 127  5 62 79    NA 2.58  NA 23  NA  0 2008  0 348 NA NA
# 2:      <NA>   1873    NA    WS3 <NA> NA NA NA NA NA  NA NA     NA  NA  NA NA NA NA    NA   NA  NA NA  NA NA   NA NA  NA NA NA
# GIDP
# 1:   NA
# 2:   NA

이처럼 참조할 값이 존재하지 않을 경우 nomatch = NULL 매개변수를 추가해주면 매칭이 되는 행의 값만을 가져올 수 있습니다. (아래 예에서는 teamID "WS3" & yearID 1873 과 매칭되는 행이 존재하지 않으므로 nomatch = NULL 옵션이 추가되니 결과값에서 없어졌습니다.)

## -- the nomatch argument
## From the previous example, setset all rows only if there is a match
Pitching[.(c("NY2", "WS3"), 1873), on = c("teamID", "yearID"), 
         mult = "last", 
         nomatch = NULL] # subset only if there's a match
# playerID yearID stint teamID lgID  W  L  G GS CG SHO SV IPouts   H  ER HR BB SO BAOpp  ERA IBB WP HBP BK  BFP GF   R SH SF
# 1: mathebo01   1873     1    NY2   NA 29 23 52 52 47   2  0   1329 489 127  5 62 79    NA 2.58  NA 23  NA  0 2008  0 348 NA NA
# GIDP
# 1:   NA

(9) 이차 인덱스 제거하기 (remove all secondary indices)

이차 인덱스를 제거할 때는 setindex(DT, NULL) 처럼 해주면 기존의 모든 이차 인데스들이 모두 한꺼번에 NULL로 할당되어 제거됩니다.

## remove all secondary indices
setindex(Pitching, NULL)

indices(Pitching)
# NULL

참고로, Key를 설정, 확인, 제거하는 함수는 setkey(DT, col), key(DT), setkey(DT, NULL) 입니다.

## set Key
setkey(Pitching, teamID)

## check Key
key(Pitching)
# [1] "teamID"


## remove Key
setkey(Pitching, NULL)
key(Pitching)
# NULL

[ Reference ]

* R data.table vignettes 'Secondary indices and Auto indexing'
  : cran.r-project.org/web/packages/data.table/vignettes/datatable-secondary-indices-and-auto-indexing.html

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

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

지난번 포스팅에서는 R data.table 에서 .SD[which.max()], .SD[which.min()] 와 by 를 활용해서 그룹별로 최대값 또는 최소값을 가지는 행을 동적으로 인덱싱(dynamic indexing for the row with maximum or minimum value) 해오는 방법을 소개하였습니다. (rfriend.tistory.com/612)

이번 포스팅에서는 R data.table 에서 그룹별로 선형회귀모형을 적합하고, 적합된 모델로부터 설명변수의 추정 회귀계수를 구하는 방법을 소개하겠습니다.

(1) 선형 회귀모형 적합하고 회귀계수 가져오기 (fitting linear regression model and getting coefficients)

(2) 그룹 별로 적합된 회귀모형의 회귀계수 구하기 (regression coefficients by groups)

(3) 그룹 별로 구한 회귀계수의 히스토그램으로 분포 확인하기 (distribution of group-level coefficients)

(4) 그룹 별 회귀계수를 data.table로 저장하기 (saving coefficients as data.table, lists)

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

library(data.table)

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

## coerce lists and data.frame to data.table by reference
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) 선형 회귀모형 적합하고 회귀계수 가져오기
    (fitting linear regression model and getting coefficients)

data.table의 .SD와 by를 활용한 그룹별 회귀모형에 들어가기 전에, R 코드에 대한 이해를 돕기 위하여 먼저 R로 선형회귀모형을 적합하는 방법을 간단히 소개하겠습니다.

아래 예는 Pitching 데이터셋에 대해 반응변수(response, dependent, target variable) 인 y: ERA 와 설명변수(explanatory, independent, input variable)인 x: ERA (Eearned Run Average, 투수의 방어율 평균자책점) 와의 관계를 선형 회귀모형으로 모델링해보았습니다.  R 에서는 lm(y ~ x, data) 의 구문으로 표현합니다.

## -- fitting linear regression with W(Win) on ERA(Earned Run Average)
lm(ERA ~ W, data = Pitching)
# Call:
#   lm(formula = ERA ~ W, data = Pitching)
# 
# Coefficients:
#   (Intercept)       W  
# 6.0704        -0.2064

lm() 함수로 선형회귀모형을 적합한 결과 객체에서 coef(lm(y ~ x, data)) 로 회귀계수에 접근할 수 있습니다.

## coefficients
coef(lm(ERA ~ W, data = Pitching))
# (Intercept)            W 
# 6.0704227     -0.2064383

특정 설명변수의 회귀계수만을 가져오고 싶으면 coef(lm(y~x, data))['var_name'] 처럼 설명변수 이름(variable name) 또는 위치(position index)를 사용해서 가져올 수 있습니다. 아래 예에서는 'W' (Win) 설명변수의 회귀계수를 가져온 것입니다.

## coefficient of variable 'W'
coef(lm(ERA ~ W, data = Pitching))['W']
# W 
# -0.2064383

(2) 그룹 별로 적합된 회귀모형의 회귀계수 구하기
     (regression coefficients by groups)

R로 회귀모형을 적합하고 회귀계수에 접근하는 법을 알았으니, 이제 R data.table에서 그룹별로 선형회귀모형을 적합하는 방법을 소개하겠습니다.

아래 예에서는 팀 그룹별로 ERA(Earned Run Average, 투수 방어율 평균자책점) 와 W (승리 회수) 간의 관계 (즉, 'W'의 회귀계수)가 서로 다를 것이라는 가정 하에,

(1) 팀 그룹 별로 (by = teamID)

(2) 투수 평균자책점(ERA)에 대한 승리 회수(W) 설명변수의 회귀계수를 w_coef 라는 이름으로 저장하는데 ( .(w_coef = coef(lm(ERA ~ W))['W']),
(3) 단, 이때 팀 그룹 별로 관측치 개수가 20개 초과인 경우로 한정(if (.N > 20))해서 구하라.

는 분석 과제입니다.

## -- Grouped Regression
## use the .N > 20 filter to exclude teams with few observations
w_coef <- Pitching[ , if (.N > 20L) .(w_coef = coef(lm(ERA ~ W))['W'])
                    , by = teamID]

w_coef
#    teamID      w_coef
# 1:    CHN -0.17955149
# 2:    CN1 -0.27648701
# 3:    BSN -0.17162655
# 4:    PRO -0.07482397
# 5:    BFN -0.12261226
# 6:    CL2 -0.04856038
# 7:    DTN -0.09514190
# 8:    PT1 -0.11607060
# 9:    LS2 -0.14260380
# 10:    SL4 -0.03346271
# 11:    BL2 -0.11725059
# 12:    PH4 -0.20383108
# 13:    CN2 -0.12078548
# 14:    NY1 -0.13258517
# 15:    PHI -0.23418637
# 16:    NY4 -0.22204042
# 17:    BR3 -0.09991895
# 18:    WS8 -0.15919173
# 19:    CL3 -0.14955735
# 20:    PIT -0.21553344
# 21:    IN3 -0.45703062
# 22:    CL4 -0.16492015
# 23:    CL6 -0.22551150
# 24:    BRO -0.28905077
# 25:    CIN -0.20696370
# 26:    WAS -0.33627146
# 27:    SLN -0.19956027
# 28:    BLN -0.15588106
# 29:    LS3 -0.27273152
# 30:    CLE -0.18379506
# 31:    PHA -0.22567468
# 32:    BOS -0.19749652
# 33:    BLA -0.13577391
# 34:    CHA -0.20046931
# 35:    WS1 -0.28093311
# 36:    DET -0.22160152
# 37:    SLA -0.24721948
# 38:    NYA -0.19447885
# 39:    PTF -0.00557913
# 40:    BLF -0.17924751
# 41:    BUF -0.23175119
# 42:    BRF -0.15565687
# 43:    ML1 -0.18098399
# 44:    BAL -0.25190384
# 45:    KC1 -0.38279088
# 46:    SFN -0.17945896
# 47:    LAN -0.17251290
# 48:    MIN -0.24984747
# 49:    WS2 -0.25201226
# 50:    LAA -0.24018977
# 51:    NYN -0.21952677
# 52:    HOU -0.23061888
# 53:    CAL -0.20546834
# 54:    ATL -0.22054211
# 55:    OAK -0.19635645
# 56:    SE1 -0.43530805
# 57:    SDN -0.24318779
# 58:    KCA -0.25287613
# 59:    MON -0.33188681
# 60:    ML4 -0.20159841
# 61:    TEX -0.25846034
# 62:    SEA -0.24887196
# 63:    TOR -0.28199100
# 64:    COL -0.32371519
# 65:    FLO -0.34167152
# 66:    ANA -0.09909373
# 67:    ARI -0.31041121
# 68:    TBA -0.31435364
# 69:    MIL -0.31820497
# 70:    MIA -0.32147649
# teamID      w_coef

(3) 그룹 별로 구한 회귀계수의 히스토그램으로 분포 확인하기
    (distribution of group-level coefficients)

위의 (2)번에서 구한 팀 그룹별 설명변수 'W'에 대한 회귀계수의 분포를 히스토그램을 그려서 확인해 보겠습니다.

또 비교를 위해서 팀 그룹의 구분이 없이 전체 데이터셋을 대상으로 하나의 선형회귀모형을 적합했을 때의 'ERA'에 대한 설명변수 'W'의 회귀계수를 overall_coef 라는 이름으로 구해서 파란색 수직 점선으로 추가해보겠습니다.

## -- Overall coefficient for comparison
overall_coef <- Pitching[ , coef(lm(ERA ~ W))['W']]

overall_coef
# W 
# -0.2064383

'ERA'에 대한 설명변수 'W'의 회귀계수는 아래의 히스토그램에서 보는 것처럼 중심을 기준으로 좌우 대칭으로 퍼져있는 정규분포 형태를 띠고 있네요.  위에서 팀 그룹 구분없이 전체 데이터셋에 대해 구한 'W'의 회귀계수 overall_coef 는 중심 부근에 위치하고 있구요.

## Histogram: team-level distribution of Win coefficinets on ERA
hist(w_coef$w_coef, 20L, las = 1L
     , xlab = "Fitted Coefficient on W"
     , ylab = "Number of Teams"
     , main = "Team-Level Distribution \n Win Coefficients on ERA")

## adding vertical line
abline(v = overall_coef, lty = 2L, lwd = 3, col = "blue")

(4) 그룹 별 회귀계수를 data.table로 저장하기
    (saving coefficients as data.table, lists)

만약 여러개의 설명변수를 사용하여 그룹별 회귀모형을 적합하고, 각 그룹별 설명변수별 회귀계수를 모두 포괄하여 추정된 회귀계수들 결과를 data.table 로 저장하려면 아래 예의 Pitching[ , as.list(coef(lm(ERA ~ W + R))), by = teamID] 와 같이 as.list() 로 회귀계수를 반환해주면 됩니다.

## making regression's coefficients as lists
coef_dt <- Pitching[ , if (.N > 100L) as.list(coef(lm(ERA ~ W + R)))
                     , by = teamID]

coef_dt
#    teamID (Intercept)          W            R
# 1:    CHN    5.710833 -0.2327841  0.009008819
# 2:    BSN    6.207018 -0.1480406 -0.003578210
# 3:    NY1    5.204519 -0.1829990  0.009074176
# 4:    PHI    6.288039 -0.2852063  0.008062580
# 5:    PIT    5.816353 -0.3000605  0.014409032
# 6:    CL4    7.069498 -0.1379608 -0.003896127
# 7:    BRO    7.389586 -0.2486108 -0.006551714
# 8:    CIN    5.767821 -0.2772234  0.011879565
# 9:    WAS    6.992822 -0.4307016  0.012169679
# 10:    SLN    5.658827 -0.2652434  0.011129236
# 11:    CLE    5.603790 -0.2500250  0.012237163
# 12:    PHA    6.688209 -0.2133816 -0.002355098
# 13:    BOS    5.796617 -0.2486252  0.009668641
# 14:    CHA    5.646432 -0.2873486  0.015712714
# 15:    WS1    7.232626 -0.2110688 -0.011945155
# 16:    DET    6.277144 -0.2542801  0.005730178
# 17:    SLA    6.347031 -0.2954950  0.007275831
# 18:    NYA    5.697457 -0.2596195  0.012947058
# 19:    ML1    5.854472 -0.1546460 -0.005409552
# 20:    BAL    6.164851 -0.3403211  0.016795283
# 21:    KC1    7.266172 -0.3501301 -0.004535116
# 22:    SFN    5.198861 -0.2817773  0.019011014
# 23:    LAN    4.935047 -0.2974111  0.025402053
# 24:    MIN    6.189153 -0.3409982  0.016008137
# 25:    WS2    5.387437 -0.3462016  0.015272805
# 26:    LAA    5.789238 -0.3530398  0.020938732
# 27:    NYN    5.498020 -0.2827097  0.012036711
# 28:    HOU    5.672472 -0.3028138  0.013698871
# 29:    CAL    5.539583 -0.2999446  0.016522022
# 30:    ATL    5.656437 -0.3144744  0.017429243
# 31:    OAK    5.397168 -0.3234133  0.024444877
# 32:    SDN    5.491807 -0.3883997  0.024155811
# 33:    KCA    6.056765 -0.3851321  0.022070084
# 34:    MON    6.564910 -0.3598835  0.005035261
# 35:    ML4    5.586930 -0.3250744  0.019599030
# 36:    TEX    6.246584 -0.3922595  0.021545179
# 37:    SEA    6.111932 -0.3346116  0.014231822
# 38:    TOR    6.528287 -0.3587363  0.013046643
# 39:    COL    6.966675 -0.4478508  0.018163225
# 40:    FLO    6.648690 -0.4681611  0.020578927
# 41:    ANA    4.825658 -0.3218342  0.033327670
# 42:    ARI    6.694492 -0.3594938  0.009086247
# 43:    TBA    6.355612 -0.4242187  0.018274852
# 44:    MIL    6.340569 -0.4942962  0.027135468
# 45:    MIA    5.629828 -0.4644679  0.025092579
#     teamID (Intercept)          W            R

[ Reference ]

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

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

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

지난번 포스팅에서는 R data.table에서 .SD[]와 by를 사용해서 그룹별로 부분집합 가져오기 (Group Subsetting) 하는 방법을 소개하였습니다. (rfriend.tistory.com/611)

이번 포스팅에서는 R data.table에서 .SD[which.max()], .SD[which.min()]과 by 를 사용해서 그룹별로 최소값 행, 최대값 행을 indexing해서 가져오는 방법(Group Optima)을 소개하겠습니다.

(1) 그룹별로 특정 칼럼의 최대값인 행 가져오기 (get the minumum row for each group)

(2) 그룹별로 특정 칼럼의 최소값인 행 가져오기 (get the maximum row for each group)

먼저, data.table 패키지를 불러오고, 예제로 사용할 데이터로 Lahman 패키지에 들어있는 야구 팀들의 통계 데이터인 'Teams' 데이터셋을 Data.Table로 참조해서 불러오겠습니다.

library(data.table)

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

## coerce lists and data.frame 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>

(1) 그룹별로 특정 칼럼의 최대값인 행 가져오기

    (get the minumum row for each group)

팀 ID 그룹 별로 (by = teamID) 승리 회수가 최대인 행(which.max(W))을 indexing 해서 가져오기 (.SD[which.max(W)] 해보겠습니다.  .SD 는 data.table 그 자체를 참조하는데요, 여기에 .SD[which.max(W)]로 W 가 최대인 index 의 위치의 행 전체를 subset 해오는 것입니다.  indexing  해오는 위치가 특정 숫자로 고정된 것이 아니라 which.max() 로 최대값의 위치를 동적으로 (dynamic indexing) 가져오게 할 수 있습니다.

## (1) Get the best year for each team, as measured by 'W'(Win)
Teams[ , .SD[which.max(W)]
       , .SDcols = c('teamID', 'yearID', 'lgID', 'franchID', 'divID', 'Rank', 'W') 
       , by = teamID]
#    teamID teamID yearID lgID franchID divID Rank   W
# 1:    BS1    BS1   1875   NA      BNA  <NA>    1  71
# 2:    CH1    CH1   1871   NA      CNA  <NA>    2  19
# 3:    CL1    CL1   1871   NA      CFC  <NA>    8  10
# 4:    FW1    FW1   1871   NA      KEK  <NA>    7   7
# 5:    NY2    NY2   1872   NA      NNA  <NA>    3  34
# ---                                                  
# 145:    ANA    ANA   2000   AL      ANA     W    3  82
# 146:    ARI    ARI   1999   NL      ARI     W    1 100
# 147:    MIL    MIL   1999   NL      MIL     C    5  74
# 148:    TBA    TBA   2009   AL      TBD     E    3  84
# 149:    MIA    MIA   2017   NL      FLA     E    2  77

(2) 그룹별로 특정 칼럼의 최소값인 행 가져오기

    (get the maximum row for each group)

팀 ID 그룹별로(by = teamID) 승리 회수가 최소인 년도의 행 전체를 가져오려면 .SD[which.min(W)] 로 dynamic indexing 을 해서 그룹별 부분집합을 가져오면 됩니다.

.SDcols 는 원하는 특정 칼럼들만 선별적으로 가져올 때 사용합니다.

## (2) Get the worst year for each team, as measured by 'W'(Win)
Teams[ , .SD[which.min(W)]
       , .SDcols = c('teamID', 'yearID', 'lgID', 'franchID', 'divID', 'Rank', 'W') 
       , by = teamID]
#    teamID teamID yearID lgID franchID divID Rank  W
# 1:    BS1    BS1   1871   NA      BNA  <NA>    3 20
# 2:    CH1    CH1   1871   NA      CNA  <NA>    2 19
# 3:    CL1    CL1   1872   NA      CFC  <NA>    7  6
# 4:    FW1    FW1   1871   NA      KEK  <NA>    7  7
# 5:    NY2    NY2   1871   NA      NNA  <NA>    5 16
# ---                                                 
# 145:    ANA    ANA   1999   AL      ANA     W    4 70
# 146:    ARI    ARI   2004   NL      ARI     W    5 51
# 147:    MIL    MIL   2002   NL      MIL     C    6 56
# 148:    TBA    TBA   2002   AL      TBD     E    5 55
# 149:    MIA    MIA   2019   NL      FLA     E    5 57

참고로, .SD[which(조건, condition)] 을 해서 특정 조건을 만족하는 행을 동적으로 인덱싱 (dynamic indexing with conditions) 해서 부분집합을 가져올 수 도 있습니다. 

아래 예에서는 야구팀 그룹(by = teamID)별로 승리 회수가 100회 이상 (.SD[which(W >= 100)] 인 년도의 행들의 부분집합을 가져온 것입니다.

## Get the year over 100 Wins for each team
Teams[ , .SD[which(W >= 100)]
       , .SDcols = c('teamID', 'yearID', 'lgID', 'franchID', 'divID', 'Rank', 'W') 
       , by = teamID]
       
# ## or equivalently
# Teams[W >= 100 , .SD
#        , .SDcols = c('teamID', 'yearID', 'lgID', 'franchID', 'divID', 'Rank', 'W') 
#        , by = teamID]

#    teamID teamID yearID lgID franchID divID Rank   W
# 1:    BSN    BSN   1892   NL      ATL  <NA>    1 102
# 2:    BSN    BSN   1898   NL      ATL  <NA>    1 102
# 3:    CHN    CHN   1906   NL      CHC  <NA>    1 116
# 4:    CHN    CHN   1907   NL      CHC  <NA>    1 107
# 5:    CHN    CHN   1909   NL      CHC  <NA>    2 104
# ---                                                  
# 105:    OAK    OAK   2001   AL      OAK     W    2 102
# 106:    OAK    OAK   2002   AL      OAK     W    1 103
# 107:    KCA    KCA   1977   AL      KCR     W    1 102
# 108:    SEA    SEA   2001   AL      SEA     W    1 116
# 109:    ARI    ARI   1999   NL      ARI     W    1 100

[ Reference ]

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

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

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

지난번 포스팅에서는 조건이 있는 상태에서 데이터셋 합치기(conditional joins) 방법을 소개하였습니다. (rfriend.tistory.com/610)

이번 포스팅에서는 .SD[]를 사용해서 그룹별로 부분집합을 가져오는 방법(Group Subsetting)을 소개하겠습니다.

(1) 그룹별로 정렬 후 마지막 행 가져오기 (subsetting the last row by groups)

(2) 그룹별로 정렬 후 첫번째 행 가져오기 (subsetting the first row by groups)

(3) 그룹별로 무작위로 행 하나 추출하기 (subsetting a row randomly by groups)

먼저, data.table 패키지를 불러오고, 예제로 사용할 데이터로 Lahman 패키지에 들어있는 야구 팀들의 통계 데이터인 'Teams' 데이터셋을 DataTable로 불러오겠습니다.

## ===============================
## R data.table
## : Grouped .SD operations
## ===============================

library(data.table)

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


data("Teams")

# coerce lists and data.frame 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>

(1) 그룹별로 정렬 후 마지막 행 가져오기
     (subsetting the last row by groups)

년도를 기준으로 내림차순 정렬(order(yearID))을 한 상태에서, 'teamID' 그룹 별(by = teamID)로 마지막 행을 부분집합으로 가져오기(.SD[.N])를 해보겠습니다.

.SD는 data.table 그 자체를 참조해서 가져오는 것을 의미하며, .SD[.N] 에서 .N 은 행의 개수(Number of rows)를 의미하므로, .SD[.N] 는 (각 'teamID' 그룹별로, by = teamID) 행의 개수 위치의 값, 즉 (teamID 그룹별) 마지막 행의 값을 부분집합으로 가져오게 됩니다.

.SDcols 는 특정 칼럼만 선별해서 가져올 때 사용하는데요, 칼럼이 너무 많아서 ID들과 순위(Rank), 경기 수(G), 승리(W), 패배(L) 칼럼만 가져오라고 했습니다.

## (1) getting the most recent season of data for each team in the Lahman data.
## In the case of grouping, .SD is multiple in nature 
## – it refers to each of these sub-data.tables, one-at-a-time
library(data.table)

Teams[order(yearID) # the data is sorted by year
      , .SD[.N]       # the recent (last row) season of data for each team
      , .SDcols = c('teamID', 'yearID', 'lgID', 'franchID', 'divID', 'Rank', 'G', 'W', 'L') 
      , by = teamID]   # subsetting by teamID groups
#    teamID teamID yearID lgID franchID divID Rank   G  W   L
# 1:    BS1    BS1   1875   NA      BNA  <NA>    1  82 71   8
# 2:    CH1    CH1   1871   NA      CNA  <NA>    2  28 19   9
# 3:    CL1    CL1   1872   NA      CFC  <NA>    7  22  6  16
# 4:    FW1    FW1   1871   NA      KEK  <NA>    7  19  7  12
# 5:    NY2    NY2   1875   NA      NNA  <NA>    6  71 30  38
# ---                                                         
# 145:    ANA    ANA   2004   AL      ANA     W    1 162 92  70
# 146:    ARI    ARI   2019   NL      ARI     W    2 162 85  77
# 147:    MIL    MIL   2019   NL      MIL     C    2 162 89  73
# 148:    TBA    TBA   2019   AL      TBD     E    2 162 96  66
# 149:    MIA    MIA   2019   NL      FLA     E    5 162 57 105

위의 (1)번과 동일한 과업을 dplyr 패키지로 수행하면 아래와 같습니다. (가장 최근의 값을 가져오기 위해 dplyr에서는 내림차순으로 정렬한 후 첫번째 행을 가져왔습니다. (= 오름차순 정렬 후 마지막 행을 가져오는 것과 동일))

## -- using dplyr
## getting the last row for each teamID group
library(dplyr)

Teams %>% 
  group_by(teamID) %>% 
  arrange(desc(yearID)) %>% 
  slice(1L) %>% 
  select(teamID, yearID, lgID, franchID, divID, Rank, G, W, L)
  
# # A tibble: 149 x 9
# # Groups:   teamID [149]
# teamID yearID lgID  franchID divID  Rank     G     W     L
# <fct>   <int> <fct> <fct>    <chr> <int> <int> <int> <int>
# 1 ALT      1884 UA    ALT      NA       10    25     6    19
# 2 ANA      2004 AL    ANA      W         1   162    92    70
# 3 ARI      2019 NL    ARI      W         2   162    85    77
# 4 ATL      2019 NL    ATL      E         1   162    97    65
# 5 BAL      2019 AL    BAL      E         5   162    54   108
# 6 BFN      1885 NL    BUF      NA        7   112    38    74
# 7 BFP      1890 PL    BFB      NA        8   134    36    96
# 8 BL1      1874 NA    BLC      NA        8    47     9    38
# 9 BL2      1889 AA    BLO      NA        5   139    70    65
# 10 BL3      1891 AA    BLO      NA        4   139    71    64

(2) 그룹별로 정렬 후 첫번째 행 가져오기
     (subsetting the first row by groups)

이번에는 년도별로 내림차순으로 정렬(order(yearID)을 한 상태에서, 'teamID' 그룹별(by = teamID)로 첫번째 행의 값을 부분집합으로 가져오기(.SD[1L])를 해보겠습니다.

.SD[1L] 에서 .SD는 (teamID 그룹별로, by=teamID) data.table 그 자체를 참조하며, '[1L]' 은 첫번째 행(1st Line)의 위치의 값을 indexing해서 가져오라는 뜻입니다.

## (2) getting the first season of data for each team in the Lahman data.
Teams[order(yearID) # the data is sorted by year
      , .SD[1L]       # the first season of data for each team
      , .SDcols = c('teamID', 'yearID', 'lgID', 'franchID', 'divID', 'Rank', 'G', 'W', 'L') 
      , by = teamID]   # subsetting by teamID groups
#    teamID teamID yearID lgID franchID divID Rank   G  W  L
# 1:    BS1    BS1   1871   NA      BNA  <NA>    3  31 20 10
# 2:    CH1    CH1   1871   NA      CNA  <NA>    2  28 19  9
# 3:    CL1    CL1   1871   NA      CFC  <NA>    8  29 10 19
# 4:    FW1    FW1   1871   NA      KEK  <NA>    7  19  7 12
# 5:    NY2    NY2   1871   NA      NNA  <NA>    5  33 16 17
# ---                                                        
# 145:    ANA    ANA   1997   AL      ANA     W    2 162 84 78
# 146:    ARI    ARI   1998   NL      ARI     W    5 162 65 97
# 147:    MIL    MIL   1998   NL      MIL     C    5 162 74 88
# 148:    TBA    TBA   1998   AL      TBD     E    5 162 63 99
# 149:    MIA    MIA   2012   NL      FLA     E    5 162 69 93

(3) 그룹별로 무작위로 행 하나 추출하기
     (subsetting a row randomly by groups)

마지막으로 년도를 기준으로 내림차순 정렬한 상태(order(yearID))에서, 'teamID' 그룹별로 (by = teamID) 무작위로 1개의 행을 부분집합으로 가져오기(.SD[sample(.N, 1L)])를 해보겠습니다.

.SD 는 (여기서는 teamID 그룹별로, by = teamID) data.table 그 자체를 참조하며, .SD[sample(.N, 1L)] 에서 sample(.N, 1L) 은 (teamID 그룹별) 총 행의 개수(.N) 중에서 1개의 행(1L)을 무작위로 추출(random sampling)해서 가져오라는 의미입니다.

## (3) getting a random row for each group.
Teams[order(yearID)          # the data is sorted by year
      , .SD[sample(.N, 1L)], # one random row of data for each team
      , .SDcols = c('teamID', 'yearID', 'lgID', 'franchID', 'divID', 'Rank', 'G', 'W', 'L'), 
      , by = teamID]         # subsetting by teamID groups

# teamID teamID yearID lgID franchID divID Rank   G  W  L
# 1:    BS1    BS1   1872   NA      BNA  <NA>    1  48 39  8
# 2:    CH1    CH1   1871   NA      CNA  <NA>    2  28 19  9
# 3:    CL1    CL1   1872   NA      CFC  <NA>    7  22  6 16
# 4:    FW1    FW1   1871   NA      KEK  <NA>    7  19  7 12
# 5:    NY2    NY2   1872   NA      NNA  <NA>    3  56 34 20
# ---                                                        
#   145:    ANA    ANA   2003   AL      ANA     W    3 162 77 85
# 146:    ARI    ARI   2012   NL      ARI     W    3 162 81 81
# 147:    MIL    MIL   2007   NL      MIL     C    2 162 83 79
# 148:    TBA    TBA   2005   AL      TBD     E    5 162 67 95
# 149:    MIA    MIA   2017   NL      FLA     E    2 162 77 85

[ Reference ]

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

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

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

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