지난번 포스팅에서는 좌표계 (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])

South Korea (from R spData package's 'world' dataset)

 

 

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

 

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

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

 

Posted by R Friend Rfriend

댓글을 달아 주세요

지난번 포스팅에서는 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) ]

 

Latitude and Longitude (*Credit : Illinois State University), https://journeynorth.org/tm/LongitudeIntro.html

 

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

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

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

 

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

 

Geoid, Ellipsoid & Sphere (*credit: https://www.slideserve.com/raisie/where-am-i)

 

 

타원면(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) ]

image source: http://faculty.kutztown.edu/courtney/blackboard/Physical/05Project/project.html

 

 

 

(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; 오른쪽) ]

image source: https://geocompr.robinlovelace.net/spatial-class.html

 

 

(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

 

 

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

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

 

Posted by R Friend Rfriend

댓글을 달아 주세요

지난번 포스팅에서는 레스터 데이터의 유형으로 (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

 

 

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

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

 

Posted by R Friend Rfriend

댓글을 달아 주세요

지난번 포스팅에서는 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

 

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

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

 

 

Posted by R Friend Rfriend

댓글을 달아 주세요

이전 포스팅에서는 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

 

 

 

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

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

 

Posted by R Friend Rfriend

댓글을 달아 주세요

지난번 포스팅에서는 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

 

 

 

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

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

 

 

 

 

 

Posted by R Friend Rfriend

댓글을 달아 주세요

지난번 포스팅에서는 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

 

 

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

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

 

Posted by R Friend Rfriend

댓글을 달아 주세요

지난번 포스팅에서는 조건이 있는 상태에서 데이터셋 합치기(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

 

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

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

 

Posted by R Friend Rfriend

댓글을 달아 주세요

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

 

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

 

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

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

 

 

 

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

 

library(data.table)

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

data("Pitching")

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


## data.table
setDT(Pitching)

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


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

 

 

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

 

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

 

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

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

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

 

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

 

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

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

 

 

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

 

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

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

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

 

 

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


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

 

 

 

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

 

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

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

 

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

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

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

 

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

 

 

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

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


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

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


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

 

 

 

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

 

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

 

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

 

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

 

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

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

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

 

 

 

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

 

구분

data.table

dplyr

(1)

conditional

left join

by key

w/ a list

library(data.table)

Pitching[G > 5, # condition

         rank_in_team := frank(ERA), 

         by = .(teamID, yearID)]

library(dplyr)

rank_in_team_df <- Pitching_df %>% 

  filter(G > 5) %>% 

  group_by(teamID, yearID) %>% 

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

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

## left outer join {dplyr}

Pitching_df <- left_join(

  Pitching_df, 

  rank_in_team_df, 

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

(2)

conditional

left join

by key

w/ a data.table

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

## -- merging team performance

Pitching_df <- left_join(

  Pitching_df, 

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

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

## condition: rank_in_team == 1

Pitching_df$team_performance <- ifelse(

  Pitching_df$rank_in_team == 1, 

  Pitching_df$Rank, # if TRUE

  NA) # if FALSE

 

 

[ Reference ]

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

 

 

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

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

 

Posted by R Friend Rfriend

댓글을 달아 주세요

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

 

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

 

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

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

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

 

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

 

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

 

library(data.table)

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

 

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

 

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

 

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

 

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


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

 

 

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

 

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

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

 

 

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

 

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

 

 

 

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

 

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

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

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



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

 

 

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

 

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

 

 

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

 

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

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

 

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

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

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

 

 

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

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

 

 

[ Reference ]

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

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

 

 

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

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

 

 

 

Posted by R Friend Rfriend

댓글을 달아 주세요