위도, 경도 좌표의 기준이 어떻게 되는지, 누가 언제 정의한 좌표 체계인지를 정의해 놓은 것입니다.
오픈소스 PostGIS 에서는 Spatial Reference System Table에 SRID 별로 'auth_name', 'auth_srid', 'srtext', 'proj4text' 등의 정보가 저장되어 있는데요, R에서도 좌표계, 좌표 참조 시스템(CRS)이 있어서 좌표계 간에 변환을 할 수 있게 해줍니다.
참고로, 전세계를 한번에 나타내야 할 때 많이 사용하는 전지구 좌표계로는 WGS84(EPSG:4326), Bessel 1841(EPSG:4004, EPSG:4162), GRS80(EPSG)4019, EPGS:4737), Google Mercator(EPSG:3857)등이 있습니다.
[전지구 좌표계] *WGS84 경위도: GPS가 사용하는 좌표계EPSG:4326, EPSG:4166 (Korean 1995)+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs *Bessel 1841 경위도: 한국과 일본에 잘 맞는 지역타원체를 사용한 좌표계EPSG:4004, EPSG:4162 (Korean 1985)+proj=longlat +ellps=bessel +no_defs +towgs84=-115.80,474.99,674.11,1.16,-2.31,-1.63,6.43 *GRS80 경위도: WGS84와 거의 유사EPSG:4019, EPSG:4737 (Korean 2000)+proj=longlat +ellps=GRS80 +no_defs *Google Mercator: 구글지도/빙지도/야후지도/OSM 등 에서 사용중인 좌표계EPSG:3857(공식),
(1) 특정 좌표 참조 시스템(CRS, Coordinate Reference Systems) 지정해주기
기본 지리특성 칼럼(sfc, Simple feature columns)의 좌표 참조 시스템의 기본 값은 NA (Not Available) 입니다.
library(sf)
## Coordinate Reference Systems (CRS) ## default CRS is NA (Not Available) point_sfg_1 = st_point(c(1, 2)) point_sfg_2 = st_point(c(3, 4)) points_sfc = st_sfc(point_sfg_1, point_sfg_2)
st_crs(points_sfc) # Coordinate Reference System: NA
st_sfc() 함수의 crs 매개변수를 사용해서 특정 좌표 참조 시스템(CRS)을 설정해줄 수 있습니다.
가령, 아래 예에서는 EPSG 코드 4326 번으로 좌표 참조 시스템을 설정하여 보았습니다.
* EPSG 코드의 EPSG는 European Petroleum Survey Group 의 앞글자를 딴 줄임말입니다. EPSG 그룹은 좌표 시스템 정보 데이터베이스와 지도 투영법과 데이터에 대한 표준 문서 자료를 발행하는 곳입니다.
* WGS84 는 World Geodetic System 1984 의 앞글자를 딴 줄임말입니다. WGS84 는 전지구 좌표계로서, 우리가 매일 사용하는 GPS(Global Positioning System), 지도 제작, 측지학 등에서 사용합니다.
## EPSG: European Petroleum Survey Group ## : They publish a database of coordinate system information ## plus some very good related documents on map projections and datums
raw proj4string 을 직접 입력해서 좌표 참조 시스템을 지정해주려면 st_sfc() 함수의 crs 매개변수에 직접 따옴표 안에 proj, datum 을 입력해주면 됩니다. 아래 예는 WGS84를 raw proj4string로 해서 좌표 참조 시스템을 등록한 것입니다.
## to specify CRS using a raw proj4string point_sfc_wgs_2 = st_sfc(point_sfg_1, point_sfg_2, crs = "+proj=longlat +datum=WGS84 +no_defs")
지난번 포스팅에서는 R의 sf 패키지를 사용해서 지리공간 벡터 데이터의 기하 유형(geometry types of vector data)로서 단순 지리특성 기하(Simple feature geometry)인 점(Point), 선(LineString), 면(다각형) (Polygon), 다중점(MultiPoint), 다중선(MultiLineString), 다중면(MultiPolygon), 기하집합(GeometryCollection)에 대해서 소개하였습니다 .
이번 포스팅에서는 앞서 소개했던 벡터 데이터의 '단순 지리특성 기하 (Simple feature geometry, sfg)' 들을 sf 패키지의 st_sfc() 함수를 사용하여 '단순 지리특성 칼럼 (Simple feature column, sfc)' 으로 합치는 방법을 소개하겠습니다.
여기서 단순 지리특성 칼럼(Simple feature column, sfc)은 단순 지리특성 기하들의 리스트(a list of sfg) 로서, 두 개의 지리특성(features)를 하나의 칼럼 객체로 합쳐놓은 것입니다. (한국말로 feature, sfc, sfg 용어를 번역하기가 쉽지 않네요.)
-- 동일한 단순 지리특성 기하 유형 합치기
(1) 두개의 단순 지리특성 기하 점(2 sfg points)를 한개의 단순 지리특성 칼럼(1 sfc) 객체로 합치기
(2) 두개의 단순 지리특성 기하 면(2 sfg polygons)를 한개의 단순 지리특성 칼럼(1 sfc) 객체로 합치기
-- 서로 다른 단순 지리특성 기하 유형 합치기
(3) 단순 지리특성 기하 점과 면을 합쳐서 한개의 단순 지리특성 칼럼(1 sfc) 객체로 만들기
(1) 두개의 단순 지리특성 기하 점(2 sfg points)를 st_sfc() 함수로 한개의 단순 지리특성 칼럼(1 sfc) 객체로 합치기
R의 sf 패키지를 사용할 것이므로, 먼저 library(sf) 로 패키지를 불러옵니다.
다음으로, 예제로 사용할 점(Point) 두개를 st_point() 함수를 사용해서 만들었습니다.
그리고, st_sfc() 함수를 사용하여, 단순 지리특성 기하(Simple feature geometry, sfg)인 두개의 점을
--> 단순 지리특성 칼럼(Simple feature columns, sfc) 의 하나의 객체로 합쳤습니다.
library(sf)
## Simple feature geometry (sfg) ## Point point_sfg_1 = st_point(c(1, 2)) point_sfg_1 # POINT (1 2)
point_sfg_2 = st_point(c(3, 4)) point_sfg_2 # POINT (3 4)
## Simple feature columns (sfc) ## st_sfc() function combine two simple features into one object with two features. ## sfc represents the geometry column in sf data frames. points_sfc = st_sfc(point_sfg_1, point_sfg_2)
points_sfc # Geometry set for 2 features # geometry type: POINT # dimension: XY # bbox: xmin: 1 ymin: 2 xmax: 3 ymax: 4 # CRS: NA # POINT (1 2) # POINT (3 4)
st_geometry_type() 함수를 사용하면 기하 유형을 확인해 볼 수 있습니다. 동일한 기하 유형인 두 개 점을 합쳤으므로 하나의 객체 집합에 두 개의 점(Points)가 들어있습니다.
## check the geometry type of sfc
st_geometry_type(points_sfc) # [1] POINT POINT # 18 Levels: GEOMETRY POINT LINESTRING POLYGON MULTIPOINT MULTILINESTRING MULTIPOLYGON ... TRIANGLE
(2) 두개의 단순 지리특성 기하 면(2 sfg polygons)를 st_sfc() 함수로
한개의 단순 지리특성 칼럼(1 sfc) 객체로 합치기
이번에는 동일한 기하 형태인 두개의 단순 지리특성 기하 면 (2 Simple feature geometry polygons) 을 st_polygon() 함수로 만들어보겠습니다.
그리고 이를 st_sfc() 함수를 사용해서 하나의 단순 지리특성 칼럼(Simple feature column) 객체로 합쳐보겠습니다.
## Polygon using list of matrices of points polygon_list_1 = list(rbind(c(1, 1), c(1, 3), c(3, 5), c(3, 1), c(1, 1))) polygon_sfg_1 = st_polygon(polygon_list_1) polygon_sfg_1 # POLYGON ((1 1, 1 3, 3 5, 3 1, 1 1))
이번 포스팅에서는 지리공간 데이터 중에서 실제 세계를 도형으로 표현(graphical representation) 하는 벡터 데이터 모델 (Vector data model)에서 사용하는 기하 유형 (geometry types)에 대해서 소개하겠습니다.
(* 참고로 레스터 데이터는 실제 세계를 픽셀로 이루어진 그리드 안에 지리특성 값을 가지는 형태라고 하였습니다.)
벡터 데이터의 기하 유형에는 총 17개가 있는데요, 이번 포스팅에서는 그중에서도 많이 사용되는 아래의 7개의 대표적인 유형에 대해서만 설명을 하겠습니다.
(1) 점 (Point)
(2) 선 (LineString)
(3) 면 (다각형) (Polygon)
(4) 다중점 (MultiPoint)
(5) 다중선 (MultiLineString)
(6) 다중면 (MultiPolygon)
(7) 기하집합 (GeometryCollection)
총 17개의 기하 유형 중에서 위에서 설명한 7개 이외에 나머지 유형에 대해서는 PostGIS 사이트를 참고하시기 바랍니다. ( http://postgis.net/docs/using_postgis_dbmanagement.html )
[ 벡터 데이터의 기하 유형 (Geometry Types of Vector data) ]
(1) 점 (Point), 다중점 (MultiPoint)
점은 2D (위도, 경도), 3D (위도, 경도, 해발고도), 또는 4D (위도, 경도, 고도, 측정 정확도) 공간 상의 하나의 위치 좌표에 해당합니다.
예: Point (3 3)
점 (Point)은 0 차원으로서, 각 기하도형의 가장 기본이 되는 형태입니다. 여러개의 점을 이어주면 선이 되고, 선을 닫아주면 면이 되는 식이니깐요. 이렇게 점 --> 선 --> 면 혹은 점 --> 다중점, 선 --> 다중선, 면 --> 다중면 방향으로 기하도형을 변환하는 것을 구성(composition) 한다고 하며, 그 반대 방향으로의 변환을 해체(Decompostion) 한다고 합니다.
점은 우물, 랜드마크, 관심지역(Point of Interest) 등의 지리특성을 정의할 때 사용합니다.
R에서 점(point)을 만들 때 st_point() 함수를 사용합니다. st_point(dim="XYM") 에서 dim 은 차원(dimension)을 의미하고, M은 보통 Measurement 의 M (측정 정확도, measurement accuracy)을 의미합니다. 점의 2D, 3D, 4D 좌표는 c() 로 묶어줍니다.(concatenate)
참고로, 접두사 'ST_' 는 PostGIS의 함수에서도 동일하게 사용되며, 'standard Spatial Type (ST)' 의 앞글자를 딴 것입니다.
install.packages("sf")
library(sf)
##---------------- ## point
## XY point st_point(c(3, 3)) # POINT (3 3)
## XYZ point st_point(c(3, 3, 2)) # POINT Z (3 3 2)
## XYM point st_point(c(3, 3, 1), dim = "XYM") # POINT M (3 3 1)
## XYZM point st_point(c(3, 3, 2, 1)) # POINT ZM (3 3 2 1)
(2) 선 (LineString)
선(LineString) 은 점들을 직선으로 연결한, 여러개 점들의 연속적인 집합입니다. 선은 아래의 예처럼 몇 개 점들의 집합으로 표현할 수 있습니다(이들 점을 순서대로 직선으로 연결).
예: LineString (1 1, 1 3, 3 5)
선은 1차원(dimension of 1)이며 길이(length)가 있습니다. (반면, 점은 0 차원이고 길이가 없음)
선은 길, 철도, 강, 전선 등과 같이 선형 특성을 가진 지형지물을 정의하는데 사용합니다.
R에서 선(LineString)은 st_linestring() 함수를 사용하여 만들 수 있습니다. 선을 이루는 여러개 점들의 좌표는 rbind()로 행렬을 만들어서 사용합니다.
면(Polygon)은 2차원의 표면을 가지는 형태로 점들이 연결된 것으로서, 안에 구멍(hole)이 있을 수도 있고 또는 없을 수도 있습니다. (위의 그림 예에서는 면 안에 구멍이 없음). 면이 되기 위해서는 시작점과 끝점이 동일해야 하며, 이를 닫혀있어야 한다(closed, not open)고 말합니다. 아래 예의 경우 면의 시작점과 끝점이 (x=1, y=1) 로서 동일하며, 즉 닫혀있습니다.
예: Polygon ((1 1, 1 3, 3 5, 3 1, 1 1))
면(Polygon)은 육지, 저수지, 건물 등과 같이 공간적인 외곽 (spatial extent)과 면적(area)을 가진 지리지형을 나타내는데 사용할 수 있습니다.
R의 st_polygon() 함수를 사용해서 면(Polygon)을 만들 수 있습니다. 면을 만들 때는 점 행렬의 리스트(list of matrices of points)를 입력값으로 사용합니다.
## Polygon using list of matrices of points polygon_list = list(rbind(c(1, 1), c(1, 3), c(3, 5), c(3, 1), c(1, 1)))
다중점(MultiPoint)은 점을 여러개 모아놓은 것(a collection of points)입니다.
예: MultiPoint (1 1, 3 3, 3 5)
R에서 점은 st_point() 함수로 만들었다면, 다중점(MultiPoint)은 st_multipoint() 로 만들 수 있습니다. 이때 st_point() 에서는 점의 좌표를 c()로 묶어주기만 했다면, 다중점(MultiPoint)을 만들 때는 rbind()로 묶어서 점의 행렬(matrices)을 만들어서 st_multipoint() 함수를 사용합니다.
## MultiPoint using matrices of points multipoint_matrix = rbind(c(1, 1), c(3, 3), c(3, 5))
R에서 면(Polygon)은 st_polygon() 함수로 만들었다면, 다중면(MultiPolygon)은 st_multipolygon() 함수로 만들 수 있습니다. 면(Polygon)을 만들 때 점 행렬의 리스트를 st_polygon() 함수에 넣어주었다면, 다중면(MultiPolygon)을 만들 때는 여러개의 면이 모인 것이므로 점 행렬 리스트들의 리스트 (list of lists)를 st_multipolygon() 함수에 넣어줍니다.
기하(도형)집합(GeometryCollection)은 위에서 설명한 점(Point), 선(LineString), 면(다각형)(Polygon), 다중점(MultiPoint), 다중선(MultiLineString), 다중면(MultiPolygon) 중에서 여러개의 기하 유형을 모아놓은 것이며, 어떤 조합(any combination of geometries)도 가능합니다.
아래 예에서는 점(Point), 다중선(MultiLineString), 면(Polygon)이 혼합된 조합의 기하집합(GeometryCollection) 입니다.
R의 st_geometrycollection() 함수를 사용해서 기하집합(GeometryCollection)을 만들 수 있습니다. 이때 각 기하 유형을 리스트(list)로 묶어서 넣어주면 됩니다.
## GeometryCollection using list of combination of geometries point_coord = c(3, 3) multilinestring_list = list(rbind(c(1, 1), c(1, 3), c(3, 5)), rbind(c(2, 2), c(5, 2))) polygon_list = list(rbind(c(4, 3), c(4, 5), c(5, 5), c(5, 3), c(4, 3)))
이번 포스팅에서는 R의 leaflet 패키지와 Open Street Map을 활용하여 RStudio 또는 Jupyter Notebook R kernel IDE 에서 Interactive Map 을 그리는데 있어,
(1) 기본 상호작용하는 지도 생성하기 (create a basic interactive map)
(2) 지도 스타일 추가하기 (add map styles)
(3) 지도 좌표에 표시 풍선 추가하기 (add marker popup at latitude/longitude coordinates)
(4) 상호작용하는 동적 지도를 웹 페이지로 저장하기 (save an interactive map as web page)
하는 방법을 소개하겠습니다.
아래의 포스팅 예시는 RStudio IDE 에서 진행하였습니다. (Jupyter Notebook R kernel 에서도 가능)
(1) 기본 상호작용하는 지도 생성하기 (create a basic interactive map)
먼저 R의 leaflet 패키지를 설치하고 불러오기를 해보겠습니다.
# install and importing leaflet package install.packages("leaflet") library(leaflet)
다음으로 leaflet 패키지의 leaflet() 함수를 사용하여, 기본 OpenStreetMap 타일을 추가하고 (add default OpenStreetMap map tiles), 대한민국의 위도(latitude), 경도(longitude) 위치에 zoom = 6 의 축소/확대 비율로 설정해서 동적인 지도(interactive map)를 그려보겠습니다.
마우스 커서로 좌측 상단의 '+'를 클릭하면 지도가 'zoom in' 확대되어 더 상세하게 볼 수 있으며, 반대로 좌측 상단의 '-' 를 클릭하면 지도가 'zoom out' 축소되어 더 넓은 지역을 볼 수 있게 됩니다. (마우스 휠을 사용해서 zoom in, zoom out 할 수도 있습니다.)
애초에 기본 지도를 생성할 때 지도의 zoom 수준을 setView(zoom = xx) 에서 숫자로 설정해줄 수도 있습니다. 바로 위에서 생성했던 지도가 zoom = 6 이었다면, 아래 지도는 zoom = 8, zoom = 16 으로 기본 지도를 생성하게 해보겠습니다. 아래에 보면, 지도가 좀더 zoom in 으로 확대되어 더 상세하게 볼 수 있습니다.
(물론 interactive map 이기 때문에 마우스를 이용해서 언제든지 자유롭게 지도를 동적으로 zoom in, zoom out 할 수 있습니다.)
addProviderTiles("Tile Name Here") 를 이용하여 외부 지도 타일을 추가할 수 있습니다.
아래는 NASA에서 제공하는 우주의 인공위성에서 바라본 2012년 밤의 지도 타일("NASAGIBS.ViirsEarthAtNight2012") 입니다. 별도의 파일을 미리 다운로드 하거나 할 필요는 없구요, 그냥 타일 이름을 addProviderTiles() 메소드 안에 큰 따옴표로 해서 적어주기만 하면 됩니다. (물론, 외부와 통신할 수 있는 인터넷망에 접속해 있어야 합니다.)
한국의 경우 남한은 밤에 환한 반면에, 북한은 평양 빼고는 어둡네요. 38선 경계도 불빛으로 환하게 구분이 되는 모습을 보노라니 약간 서글프네요.
위의 두 예 외에 leaflet() 에 다른 외부 타일 제공자의 타일 지도를 추가하고 싶으면 아래의 url 의 Leaflet-prodivers preview 에 방문해서 원하는 형태의 지도 스타일을 찾았다면 그 이름을 addProviderTiles("Tile Name Here") 메소드 안에 타일 이름을 큰 따옴표 안에 써주기만 하면 됩니다.
(add marker popup at latitude/longitude coordinates)
addMarkers() 메소드를 사용하면 위도(latitude), 경도(longitude) 좌표 위치에 풍선 모양의 표식과 커서를 클릭했을 때 팝업으로 나타나는 설명을 추가할 수 있습니다.
아래 예제에서는 대한민국 삼성동의 COEX Mall (코엑스 몰), GBC (현대자동차 Global Business Center), Samsil-Stadium (올림픽 주 경기장) 의 3개 주요 POI (Point Of Interest) 지점에 말풍선 표식을 추가해보았습니다. 마우스 커서를 말풍선 표식에 가져가서 클릭을 해보면 popup 설명이 튀어나옵니다.
(아래 예제는 스크린 캡쳐를 한 이미지 파일이어서 interactive map 기능이 안됩니다. RStudio 나 Jupyter Notebook R kernel 에서 interactive map 기능 사용 가능합니다.)
지난번 포스팅에서는 지리공간 데이터의 포맷으로서 벡터 데이터(vector data)와 레스터 데이터(raster data) 모델에 대해서 간단하게 소개(https://rfriend.tistory.com/589)를 하였습니다.
이번 포스팅에서는 R의 sf 패키지를 활용해서 spData 에 내장되어 있는 세계지도(world map)를 시각화해보겠습니다.
(1) 지리공간 데이터 처리 및 분석을 위한 R 패키지 설치
(2) 요약통계량 summary() 및 부분 집합 subset 하기
(3) sf 패키지의 plot()으로 기본 지도 시각화하기
(4) 다른 지도 층을 추가하기 (add plots as layers)
(5) 원과 텍스트를 덮어쓰우기 (overlaying circles and texts over the plot)
(1) 지리공간 데이터 처리 및 분석을 위한 R 패키지 설치
R의 sf 패키지는 지리공간 벡터 데이터(vector data) 분석을 위한 패키지로서, 이전의 sp 패키지를 대체합니다. 그리고 GEOS와 GDAL과 R의 인터페이스를 제공하며 rgeos, rgdal 패키지를 대신합니다.
R의 raster 패키지는 지리공간 레스터 데이터(raster data)를 처리 및 분석하는데 사용합니다.
그리고 spData, spDataLarge 패키지에는 지리공간 데이터 샘플을 내장하고 있습니다. spDataLarge 패키지는 내장 데이터의 크기가 너무 크기 때문에 R CRAN 사이트에 패키지가 올라가 있지 않습니다. 따라서 CRAN에서 설치하려고 하면 에러가 나며, 아래처럼 Github 에서 type = "source" 로 해서 설치를 해주어야 합니다.
(만약 spData 패키지의 개발자 버전을 Github에서 바로 받아서 설치를 하려면 devtools::install_github("nowosad/spData") 로 하면 됩니다.)
# install packages
install.packages("sf") # for vector data, ‘spatial data frame’ install.packages("raster") # for raster data install.packages("spData") install.packages("spDataLarge", repos = "https://nowosad.github.io/drat/", type = "source")
참고로, spData 패키지에는 37개의 지리공간 데이터셋이 내장되어 있으며, https://nowosad.github.io/spData/ 깃헙 페이지에 가면 전체 리스트를 확인할 수 있습니다. 아래의 데이터셋 리스트 중에서 'world' (World country polygons) 데이터셋을 가지고 sf 패키지를 사용해서 요약통계량 및 시각화를 예를 들어보겠습니다.
(2) 요약통계량 summary() 및 부분 집합 subset 하기
spData 패키지에 내장된 'world' 벡터 데이터에는 세계 대륙의 국가별로 10개의 집계데이터 칼럼과 마지막에 geom (지리공간 Polygon 리스트) 으로 구성되어 있습니다.
sf 패키지는 base R의 기본 함수를 그래도 사용할 수 있습니다. names()로 지리공간 데이터의 칼럼를 리스트업 할 수도 있고, summary() 함수로 요약통계량을 계산할 수 있습니다.
summary() 함수로 전세계 국가별 인구("pop")와 예상수명("lifeExp")에 대해 177개 Multipolygon 별로 최소/Q1/중앙값/Q3/최대값/결측값의 요약통계량을 계산해 보겠습니다.
world 세계지도 데이터 칼럼 중 7열~10열까지 4개 칼럼 별로 177 개 MultiPolygon별 값에 따라 색깔을 달리해서 R의 plot() 함수로 시각화를 해보겠습니다. ('area_km2', 'pop', 'lifeExp', 'gdpPercap' 칼럼의 값에 따라서 'geom' 의 177개 MultiPolygon 별 색깔이 달라짐.)
단 plot() 함수로 시각화를 하면 매우 간단하게 지도를 시각화할 수 있지만, 상호작용하는 동적인 지도(interactive map)가 아니라 그냥 정적인 지도만 그려지게 됩니다. (interactive map 은 다음번 포스팅에서 소개 예정)
##-- making a basic map using sf's plot() # not interactive plot plot(world[7:10])
위처럼 여러개의 칼럼에 대해 한꺼번에 여러개의 지도를 시각화하면 색깔이 무엇을 의미하는지 범례가 표시되지 않습니다.
반면에, 아래처럼 특정 한개의 변수에 대해서만 plot() 으로 시각화를 하면 색깔이 의미하는 바에 대한 범례(legend)를 색깔 막대 형태로 표시해줍니다.
## plot of 'lifeExp' per countries plot(world["lifeExp"])
plot() 으로 지리공간 데이터를 시각화할 때 col 매개변수로 특정 색깔을 지정해줄 수도 있습니다. 대신, 특정 색깔을 col 로 지정해주면 위의 지도 시각화처럼 오른쪽에 색깔 막대 범례는 필요가 없으므로, 비록 1개 변수만 시각화한다고 해도 범례는 없습니다.
## no continuous palette legend plot(world["lifeExp"], col="yellow")
(4) 다른 지도 층을 추가하기 (add plots as layers)
위의 (3)번 처럼 지도를 먼저 한번 시각화를 하고, add = TRUE 매개변수를 사용하면 나중에 그 위에 다른 지도를 겹쳐서, 즉 층을 추가하여 지도를 덮어쓰기로 그릴 수 있습니다.
단, 첫번째 지도 그래프에 키(key)가 있을 경우에는 reset = FALSE 매개변수를 꼭 설정해준 다음에, 이후에 다음번 plot(add = TRUE)를 사용해주어야 합니다.
# add plots as layers world_asia = world[world$continent == "Asia", ] asia = st_union(world_asia) # combining countries in Asia
# if the first plot has a key, reset=FALSE must be used. plot(world["lifeExp"], reset = FALSE) plot(asia, add = TRUE, col = "red")
(5) 원을 덮어쓰우기 (overlaying circles over the plot)
cex 매개변수를 사용하면 시각화한 지도 위에 원을 덮어쓰울 수 있습니다.
아래 예에서는 세계지도 위에 국가별 인구에 비례해서 원을 덮어쓰워보겠습니다. (pop 변수에 제곱근을 취하고 10000 으로 나누어준 것은 시각화한 지도 크기에 맞추어서 너무 크지도, 작지도 않은 원을 덮어쓰우기 위해 숫자를 조정한 것입니다.)
지리공간 데이터 (GeoSpatial data)를 처리하고 분석하는데 있어서 첫번째 관문이자 큰 도전사항 중에 하나가 지리공간 데이터 포맷이 매우 다양하다는 것입니다.
아래에 다양한 지리공간 데이터(various GeoSpatial data foramts)의 리스트를 소개하고, 특히 이중에서 점, 선, 다각형으로 구성된 벡터 데이터 포맷의 이미지 시각화를 예시로 보였습니다. 지리공간 데이터 포맷이 상당히 많지요?
이들 지리공간 데이터 포맷별로 데이터를 DB나 R로 불러오기 (importing)할 때 사용하는 DB utility tools 이나 R의 package가 달라지다 보니 번거롭고 또 어려운 점이 있습니다.
[ 다양한 지리공간 데이터 포맷 (various GeoSpatial data formats) ]
R을 활용한 지리공간 데이터의 처리 및 분석, 시각화를 본격적으로 들어가기 전에 먼저, 이들 지리공간 데이터 포맷들 중에서 특히 벡터 데이터(Vector data)와 레스트 데이터 (Raster data) 모델에 대해서 이들이 무엇이고, 어떻게 활용이 되며, 무슨 R 패키지를 사용해서 분석할 수 있는지에 대해서 알아보겠습니다.
[ 지리공간벡터 데이터(Vector data) vs. 레스터 데이터 (Raster data) ]
(1) 지리공간 벡터 데이터 (Vector data)
벡터 데이터에는 KML(.kml or .kmz), GML, GeoJSON, Shapefile (.shp), WKT 등의 데이터 포맷이 있습니다.
GeoJSON 데이터 포맷은 이름에서 짐작할 수 있듯이 JSON 기반으로 간단한 지리공간 데이터와 그 외 일반 데이터를 저장합니다. GeoJSON 데이터는 인터넷으로 지리공간 & 일반 데이터를 공유하는데 역시 많이 사용됩니다.
Shapefile 데이터 포맷은 GIS (Geographic Information System) 소프트웨어를 위한 지리공간 벡터 데이터입니다. Shapefile 은 GIS 의 국제적인 제공사인 Esri(Environmental Systems Research Institute)에서 개발하고 관리하며, GIS 소프트웨어 간 상호운용성(interoperability)를 보장합니다.
WKT 데이터 포맷은 Well-Known Text 의 약자로서, 벡터 지리공간 데이터를 표현하는데 텍스트 마크업 언어(Text Markup Language)를 사용합니다. WKB (Well-Known Bianry)는 WKT와 같은 정보를 저장하는데 있어 이진(binary) 포맷을 사용해 보다 간소하고 컴퓨터가 처리하기에 편리하도록 하며, 대신 사람이 읽을 수는 없습니다.
벡터 데이터는 실제 세상을 그래픽으로 재표현(graphical representation of the real world)한 것으로서, 점, 선, 다각형(points, lines, polygons) 유형의 그래픽을 이용합니다. 벡터 데이터는 지구 표면의 객체나 특징을 일반화하여 표현하는데 사용됩니다.
벡터 데이터는 별개로 분리되고, 경계가 잘 정의되어 있어서 보통 높은 수준의 정밀도 (high level of precision) 을 가지고 있습니다. 이런 이유로 벡터 데이터는 사회 과학 (social sciences) 분야에서 많이 사용됩니다.
R 의 sf 패키지 (spatial data frame) 를 사용하여 벡터 데이터를 불러오고, 처리 및 분석, 시각화를 할 수 있습니다. (다음 포스팅에서 소개) sf 패키지는 이전의 sp 패키지, rgeos, rgdal 패키지를 모두 아우르고 있고, GEOS, GDAL, PROJ 와 R 의 interface를 제공해주어서, R로 지리공간 벡터 데이터를 다루는데 있어 매우 편리하고 강력합니다.
[ 강과 도심 지역을 나타낸 벡터 데이터(vector data)와 레스터 데이터(raster data) 비교 ]
레스터 데이터(Raster data)에는 ESRI Grid, GeoTIFF, JPEG 2000, NITF 등이 데이터 포맷이 있습니다.
레스터 데이터는 픽셀의 격자(grid of pixels) 로 지구의 표면을 표현합니다. 각 픽셀 안에는 색, 측정 단위 등과 같이 질문의 요소에 대한 정보를 전달하는 값이 있습니다.
레스터 데이터는 인공위성이나 항공장비에서 지구 표면을 향해 위에서 아래로 수직으로 찍은 사진으로 생각하면 이해하기가 쉽습니다.(예: NASA에서 제공하는, 인공위성에서 찍은 지구의 야간 사진 등) 이 지구표면을 수직으로 찍은 사진을 픽셀의 격자로 나누어서 각 픽셀(pixel, cell)에 지리특성정보 값을 가지고 있는 것입니다.
레스터 데이터의 품질은 사진을 찍었던 장비의 해상도의 한계나, 활용하고자 하는 분야의 목적에 따라서 다양합니다. 레스터 데이터는 많은 환경관련 과학 분야 (environmental sciences)에 많이 사용되고 있습니다.
R의 raster 패키지를 사용하면 R에서 레스터 데이터를 처리할 수 있습니다.
위에서 각각 소개한 벡터 데이터와 레스터 데이터는 상호 간에 변환(converting from vector to raster data, from raster to vector data)이 가능하며, 하나의 분석 목적에 두 유형의 데이터 포맷이 동시에 사용되기도 합니다.
다음번 포스팅에서는 R의 spData 패키지에 내장되어 있는 지리공간 벡터 데이터 모델(Vector data model)을 가지고 sf 패키지로 시각화하는 간단한 예를 소개하겠습니다.
k-means 군집분석의 경우 입력 모수(input parameter)로서 '군집의 개수 k'를 결정하는 것이 어려움이라면, DBSCAN 은 입력 모수로서 '(a) 점으로 부터의 반경 Eps (Epsilon)'와 '(b) Eps 내 최소 점의 개수 기준인 MinPts' 를 결정하는 것이 중요하고도 어려운 문제 중에 하나입니다.
그런데 무슨 연립방정식 풀듯이 이론적으로 증명된 DBSCAN의 입력모수 MinPts 와 Eps 를 구할 수 있는 공식 같은 것, 객관적인 통계량 같은 것은 없습니다. 다만 MinPts와 Eps를 결정하는데 도움을 받을 수 있는 주관적인 Heuristic method 가 있을 뿐입니다.
이번 포스팅에서는 밀도 기반 군집분석 DBSCAN 의 입력 모수 Eps, MinPts 를 결정하는 Heuristic 방법 (Determining the Parameters Eps and MinPts) 을 소개하겠습니다.
(1) MinPts 결정하는 Heuristic 방법: ln(n)
(2) Eps 결정하는 Heuristic 방법: elbow (knee) method using sorted k-dist plot
(1) DBSCAN에서 MinPts 결정하는 Heuristic 방법
MinPts 는 한 점으로부터 반경 Eps 인 원을 그렸을 때 그 점이 코어 점 (core points), 군집이 되기 위해 Eps 안에 필요한 최소한의 점 개수를 말합니다.
MinPts 가 만약 너무 작은 수이면 잡음(noise)으로 구분되어야 할 점들 마저도 코어 점(core points)나 또는 경계점(border points)로 잘못 구분이 되어 원래 데이터셋 내의 군집 개수보다 더 많은 수의 군집이 형성될 수가 있으므로 주의가 필요합니다.
MinPts 를 결정할 때는 데이터 특성과 구조에 대해서 잘 알고 있는 업 전문가 (domain expert) 의견을 반영할 필요가 있습니다. 그런데 현실은 데이터 분석을 할 때 업 전문가가 없을 수도 있고, 있더라도 MinPts를 잘 결정할 수 없을 수도 있으므로 Heuristic 방법을 알아 둘 필요가 있습니다.
DBSCAN의 원 논문(참조 [1])에서는 2차원 데이터에 대해 실험을 해보니 MinPts 가 4개와 5개 이상 간의 k-dist plot (아래 설명 예정) 의 큰 변동이 없는 반면에 MinPts 가 점점 커질 수록 연산량(computation)이 상당히 커지므로 2차원 데이터에서는 MinPts = 4 개로 하는 것을 권장하고 있습니다.
2차원보다 많은 변수를 가지고 있는 데이터셋의 경우 MinPts = 2 * dim 을 추천하는 논문(참조 [2])도 있습니다.
데이터셋별로 데이터의 구조나 객체의 개수 n이 서로 다를 수 있으므로, 데이터셋별 객체 개수 n 특성을 감안해서MinPts를 결정하는 Heuristic 방법으로 ln(n)을 사용할 수 있습니다.(참조 [3]) 여기서 n 은 데이터 개수 (number of points in database) 를 말합니다.
(2) DBSCAN에서 Eps 결정하는 Heuristic 방법
: Elbow (Knee) method using sorted k-dist plot
DBSCAN의 원 논문(참조 [1])에서는 아래의 sorted k-dist graph 를 그린 후 Elbow method 를 사용해서 첫번째 계곡(first "valley") 지점의 점이 구분 기준점(threshold point)이 되고, 이 기준점의 왼쪽은 잡음(noist), 기준점의 오른쪽은 군집으로 구분하고, 꺽이는 부분의 k-dist 를 Eps 로 결정하는 Heuristic 방법을 소개합니다.
sorted k-dist graph 를 그리는 방법은, 먼저 MinPts 를 k개라고 했을 때, 하나의 점으로부터 k개의 가장 가까운 점들 간의 거리, 즉 k_NN (k-Nearest Neighbor)의 거리 k-dist 를 구해고, k-dist 를 내림차순으로 정렬하여, X축은 정렬된 점들 별로 Y축에는 k-dist 를 그려줍니다.
(단, R의 dbscan 패키지에서는 k-dist 가 오름차순으로 정렬이 된 sorted k-dist plot 을 그려줘서 원 논문과는 그래프의 좌우가 반대입니다.)
R의 factoextra 패키지에 내장되어 있는 multishapes 데이터셋을 가지고 R의 dbscan 패키지를 사용해서 최적의 Eps 모수 값을 결정해보겠습니다.
예제 데이터셋 multishapes 은 아래의 산점도처럼 크기가 다른 원 고리형 2개, 선형 2개, 원형 1개의 5개 군집과 잡음들로 구성이 되어있습니다.
## factoextra for visualizing clusters if(!require(devtools)) install.packages("devtools") devtools::install_github("kassambara/factoextra")
다음으로 위에서 설명했던 sorted k-dist plot 을 R의 dbscan 패키지의 kNNdistplot(data, k) 함수를 사용해서 k=5로 지정하고 그려보겠습니다. (k-dist 의 오름차순 기준으로 정렬이 되어서 원 논문과는 반대 모양임.)
sorted k-dist plot 에서 꺽이는 팔꿈치(elbow) 부분의 점이 threshold point 가 되겠습니다. 아래의 예처럼 비교적 눈에 띄게 구분이 되는 경우도 있고, 어디를 꺽인 팔꿈치 (혹은 무릎 knee) 부분인지 콕 집기가 애매한 데이터셋도 있습니다. (객관적이라기 보다는 다분히 주관적입니다.)
이 기준점을 기준으로 왼쪽까지는 군집(clusters)에 속하는 점들이 되겠고, 오른쪽 부터는 잡음 점(noise points)로 간주합니다. 그리고 이 기준점(threshold point)의 k-NN distance 를 최적의 Eps로 결정하면 됩니다. 이렇게 하면 k-dist(p) 와 같거나 작은 값을 가지는 점들은 모두 코어 점(core points)가 됩니다.
아래 예에서는 꺽이는 팔꿈치 부분의 점의 5-NN distance (k-dist(p)) 가 0.15 이므로 Eps 를 0.15로 하면 되겠네요. (k 를 4 ~ 7 까지 바꾸어 가면서 sorted k-dist plot 을 그려보니 k가 커질수록 elbow 지점의 Eps 가 조금씩 커지기는 하는데요, 큰 차이는 없네요.)
## Determining the optimal Eps value using k-dist plot & elbow method
* [1] Martin Ester, Hans-Peter Kriegel, Jorg Sander, Xiaowei Xu, 1996, "A Density-Based Algorithm for Discovering Clusters in Large Spatial Databases with Noise", KDD-96
* [2] Erich Schubert, Jorg Sander, Martin Ester, Hans Peter Kriegel, Xiaowei Xu, 2017, "DBSCAN Revisited: Why and How You Should (Still) Use DBSCAN
* [3] Chossing eps and minpts for DBSCAN (R)? : https://stackoverflow.com/questions/12893492/choosing-eps-and-minpts-for-dbscan-r
* [4] DBSCAN in R : http://www.sthda.com/english/wiki/wiki.php?id_contents=7940
이번 DBSCAN(Density-Based Spatial Clustering of Applications with Noise) 알고리즘에 대한 포스팅은 Martin Ester, el.al, "A Density-Based Algorithm for Discovering Clusters in Large Spatial Databases with Noise"[1]논문을 참조하여 작성하였습니다.
군집분석은 공간 데이터(spatial data)의 그룹, 구조, 구성요소 등을 식별 (class identification)하는 과업에 활용될 수 있습니다. 하지만 대용량의 공간 데이터에 대한 군집화는 다음의 3가지 요건을 충족시킬 수 있어야 합니다.
(1) 대용량 데이터를 다룰 때 적당한 입력 모수에 대한 선험적 지식이 종종 알려져있지 않으므로, 입력 모수(input parameter)를 결정하기 위해 필요한 업 지식은 최소화되어야 하고 (minimal requirements of domain knowledge),
(2) 공간데이터의 형태는 구형, 옆으로 퍼진 형태, 선형, 가늘고 긴 형태 등 다양할 수 있으므로, 임의의 형태의 군집 탐색 가능 (discovery of clusters with arbitrary shape)해야 하고,
(3) 단지 수천개의 객체를 가진 작은 데이터셋 뿐만이 아니라, 대용량 공간데이터에 대해서도 효율적으로 군집화 연산이 가능해야 함.
이번에 소개하는 DBSCAN 알고리즘은 위의 3가지 대용량 공간 데이터에 대한 군집화 조건을 모두 만족합니다. DBSCAN 알고리즘은 이론과 현실 문제 적용에서의 우수성을 인정받아서 2014년에 ACM SIGKDD 데이터마이닝 컨퍼런스에서 상을 받기도 했습니다.
아래의 시각화는 여러가지 모양과 형태의 데이터셋들에 대해 다양한 군집분석 알고리즘을 적용하여 그 결과를 비교한 것입니다. 아래 10개의 군집화 알고리즘 중에서 DBSCAN의 군집화 결과가 데이터셋의 형태에 상관없이 매우 좋은 군집화 결과를 보여주고 있습니다. DBSCAN은 군집에 속하지 않는 잡음(noise), 이상치(outlier) 도 탐지를 할 수 있어서 anomaly detection 에도 쓸 수가 있습니다.
(정의 1) 반경 Eps 이내 이웃점 (Eps-neighborhood of a point)
: 공간 데이터베이스에 속한 점들 중에서 두 점 p와 q의 거리가 반경 Eps() 이내인 점
* D : 데이터베이스(Database), 데이터셋
* dist(p, q) :점 p와 q 의 거리(distance)
이때, 분석가가 입력해줘야 하는 모수로서,
* Eps (Epsilon, 엡실론): 점 p 로 부터의 반경
* MinPts : 최소 기준 점 개수 (minimum number of points)
DBSCAN 알고리즘은 밀도 기반(Density-Based)라고 했는데요, 이때 데이터 점의 밀도(density)는 하나의 점으로 부터 반경 Eps 이내에 점이 몇 개나 있는지로 측정합니다.
(정의 2) 코어 점(core points), 경계 점(border points)
* 코어 점 (core points) : 군집 내 점 (points inside of the cluster), 한 점으로 부터 반경이 Eps 인 원을 그렸을 때 그 원 안에 이웃점(Eps-neighborhood of a point)이 MinPts 이상의 점이 있는 점.
* 경계 점 (border points) : 군집의 경계에 있는 점 (points on the border of the cluster)
* 출처: [1] Martin Ester, el.al
(정의 3) 직접적으로 밀도(기반)-도달가능한 (directly density-reachable)
점 p가 점 q의 반경 Eps 이내에 있는 이웃점에 속하고, 점 q의 반경 Eps 이내 이웃점의 개수가 MinPts 이상일 때 점 p는 점 q 로 부터 직접적으로 밀도(기반)-도달가능하다고 합니다.
(a)
(b)
(core point condition)
위의 figure 2 에서 점 p는 점 q로 부터 직접적으로 밀도(기반)-도달가능합니다. 하지만 점 q는 점 p로 부터 직접적으로 밀도(기반)-도달가능하지 않습니다. 왜냐하면 점 p는 코어 점(core point)의 조건을 충족시키지 못하기 때문입니다.
(정의 4) 밀도(기반)-도달가능한 (density-reachable)
만약 연쇄적인 점들
들이 있고, 점
이 점
로 부터 직접적으로 밀도(기반)-도달가능하다면(directly density-reachable), 점 p는 점 q로 부터 반경 Eps 내 MinPts 조건 하에 밀도(기반)-도달가능(density-reachable)하다고 합니다.
(정의 3)의 '직접적으로 밀도(기반)-도달가능'은 대칭적(symmetric)인 반면에, (정의 4)의 '밀도(기반)-도달가능'은 비대칭적(asymmetric)입니다.
(정의 5) 밀도(기반)-연결된 (density-connected)
만약 두 점 p와 q가 모두 어떤 점 o 로 부터 반경 Eps 내 MinPts 조건 하에 밀도(기반)-도달가능(density-reachable)하다면 점 p는 점 q와 반경 Eps 내 최소 점 개수 MinPts 조건 하에 밀도(기반)-연결되었다고 합니다.
(다르게 말하면, p 가 o의 친구이고, q도 o의 친구이면, p와 q는 친구 o를 통해 서로 연결되었다고 보면 됩니다.)
(정의 6) 군집 (cluster)
드디어 이제 DBSCAN 알고리즘이 군집(cluster)를 어떻게 정의하는지 말해볼 때가 왔군요. 위의 정의(1)~(5)까지의 용어와 개념을 사용하여 군집을 정의해보면, "군집은 밀도(기반)-도달가능한 최대치의 밀도(기반)-연결된 점들의 집합이다 (A cluster is defined to be a set of density-connected points which is maximal with respect to density-reachablility.)" 라고 할 수 있겠습니다.
* 출처: [1] Martin Ester, el.al
D를 공간 점들의 데이터베이스(데이터셋) 라고 하고, C 를 군집(Cluster) 라고 했을 때 군집 C는 반경 Eps와 최소 점 개수 MinPts 조건이 주어졌을 때, 아래의 (a) 최대의 밀도(기반)-도달가능조건과 (b) 연결성 조건을 만족하는 D의 비어있지 않은 부분집합이라고 할 수 있습니다.
(a) Maximality wrt. Density-reachability
(b) Connectivity
(정의 7) (잡음, noise)
잡음 점은 군집에 속하지 못하는 점. 즉, 코어 점도 아니고 경계 점도 아닌 점을 말합니다.
분석의 목적이 군집화라면 잡음점을 무시하거나 제거하면 되구요, 만약 분석의 목적이 군집화가 아니라 anomaly detection, outlier detection 이라면 잡음 점(noises)들이 주요 관심사가 되겠습니다.
DBSCAN 알고리즘의 군집화 절차는 아래와 같이 정리할 수 있습니다.
입력 모수(input parameters)로서, 점으로 부터의 (a) 반경 Eps와 (b) Eps 반경 내 최소 점 개수 기준인 MinPts 조건이 주어졌을 때,
(1) 공간 데이터셋으로부터 초기값(seed)으로서 코어 점(core points)의 조건을 만족하는 임의의 점을 선택합니다.
(2) 초기값으로 부터 밀도(기반)-도달가능한 점들을 뽑아서 코어 점(core points)과 경계 점(border point)을 구분하고, 이에 속하지 않은 점들을 잡음(noises)으로 구분합니다.
(3) 반경 Eps 인 원 주위에 있는 코어 점들을 서로 연결합니다.
(4) 연결된 코어 점들을 하나의 군집으로 정의합니다.
(5) 모든 경계점들을 어느 하나의 군집으로 할당합니다. (만약 경계점 중에 여러 군집에 걸쳐있는 경우는 반복 과정에서 먼저 할당된 군집으로 할당함.)
왼쪽의 그림은 Wikipedia에 소개된 내용인데요, 원 논문의 그림보다 좀더 이해하기 쉬울거 같아서 한번 더 소개합니다.
점으로 부터의 반경이 Eps 이고 MinPts = 4 라고 했을 때,
점 A를 포함해서 가운데의 빨간점 6개는 코어 점(core points) 입니다.
그리고 코어점은 아니지만 코어점과 연결 가능한 노란색 점 B, C 는 경계점(border points) 입니다.
그리고 코어 점도 아니고 경계 점도 아닌 파란색 점 N 은 잡음(noise) 점이 되겠습니다.
군집화의 과정은 각 점들을 순회하면서 재귀적(recursive)으로 도달가능, 연결가능을 평가하면서 진행이 됩니다.
(2) k-means Clustering vs. DBSCAN 군집화 알고리즘 비교
k-means 와 DBSCAN 군집화 알고리즘을 아래의 표에 비교해서 정리해보았습니다. 유사성(혹은 비유사성 거리) 기반의 k-means 대비, 밀도 기반의 DBSCAN 의 경우 분석가가 미리 군집의 개수 (k)를 입력해주지 않아도 되고, 잡음/이상치에도 견고하며, 계산 복잡도도 상대적으로 작습니다. 게다가 군집으로 찾아낼 수 있는 모양도 구형, 원형, 길게 늘어선 형태, 선형 등 임의의 모양에 대해서 비교적 잘 군집화를 하고, 잡음/이상치도 별도로 구분을 해낼 수 있습니다. 이래저래 DBSCAN이 k-means 대비 우수한 점이 많습니다.
(3) R을 이용한 DBSCAN 군집화 (예시)
DBSCAN 분석을 위한 R코드는 참조 [3] 사이트의 코드를 거의 그대로 사용하였습니다.
예제로 사용할 데이터셋으로 중앙이 비어있는 원형, 선형, 구형 등 다양한 형태의 군집과 잡음으로 구성된, factoextra 패키지에 내장되어 있는 multishapes 데이터셋을 사용하겠습니다.
## factoextra for visualizing clusters if(!require(devtools)) install.packages("devtools") devtools::install_github("kassambara/factoextra")
## multishapes dataset str(multishapes) # 'data.frame': 1100 obs. of 3 variables: # $ x : num -0.804 0.853 0.927 -0.753 0.707 ... # $ y : num -0.853 0.368 -0.275 -0.512 0.811 ... # $ shape: num 1 1 1 1 1 1 1 1 1 1 ...
plot(multishapes$x, multishapes$y)
k-means clustering 과 DBSCAN clustering 알고리즘을 비교해보기 위해, 위의 multishapes 데이터셋에 대해 군집의 개수 k=5로 해서 k-means 군집화를 해보겠습니다.
k-means 군집화의 경우, (1) 상단에 위치한 중앙이 비어있는 원형 군집 2개를 제대로 구분하지 못하고 있고, (2) 좌측 하단에 위치한 선형 2개 군집도 제대로 구분하지 못하고 있으며, (3) 우측 하단의 물방울 형태 군집도 잡음/이상치가 섞여서 군집화가 되었습니다. 전혀 만족스럽지 않은 군집화 결과네요.
이번에는 DBSCAN 알고리즘으로 입력 모수로서 Eps = 0.15, MinPts = 5 로 하여 multishapes 데이터셋에 대해 군집화를 해보겠습니다.
R의 fpc 패키지나 dbscan 패키지를 사용하여 DBSCAN 알고리즘으로 군집화를 할 수 있습니다. 아래 예시에서는 fpc 패키지를 사용하였으며, 패키지명::함수명() 형태로서 fpc::dbscan(data, eps, MinPts) 으로 패키지 이름을 명시적으로 입력해주었습니다.
아래 표(참조 [4])는 프로그래밍 언어별 DBSCAN 을 할 수 있는 패키지들을 비교한 것입니다. R 의 dbscan 패키지가 지원하는 기능 면에서는 가장 강력하네요. R fpc 패키지는 dbscan 에만 특화된 simple 한 패키지이구요.
##-- DBSCAN using fpc or dbscan package install.packages("fpc") library(fpc)
## Compute DBSCAN using fpc package set.seed(1004) db <- fpc::dbscan(df, eps = 0.15, MinPts = 5)
## or # install.packages("dbscan") # library(dbscan)
마지막으로 factoextra 패키지의 fviz_cluster() 함수로 DBSCAN 군집화 결과를 시각화해보았습니다.
DBSCAN 군집화 결과를 보면, (1) 상단에 위치한 중앙이 비어있는 원형 군집 2개를 잘 구분하였고, (2) 좌측 하단에 위치한 선형 2개 군집도 제대로 구분하였으며, (3) 우측 하단의 물방울 형태 군집도 잘 군집화하였을 뿐만 아니라, (4) 잡음/이상치도 잘 구분하였습니다. 매우 만족스러운 군집화 결과네요!
각 클러스터별로 요약통계량 (가령, 평균)을 계산해서 profiling 해보면 클러스터의 특징을 파악하는데 많은 도움이 될 것입니다. 아래 코드는 dplyr 패키지의 group_by() 와 summarise() 함수를 사용해서 각 DBSCAN Cluster 별로 x, y 변수의 평균을 구해본 것입니다.
* Martin Ester, Hans-Peter Kriegel, Jorg Sander, Xiaowei Xu, 1996, "A Density-Based Algorithm for Discovering Clusters in Large Spatial Databases with Noise", KDD-96
* Martin Ester, Hans-Peter Kriegel, Jorg Sander, Xiaowei Xu, 1996, "A Density-Based Algorithm for Discovering Clusters in Large Spatial Databases with Noise", KDD-96
* Martin Ester, Hans-Peter Kriegel, Jorg Sander, Xiaowei Xu, 1996, "A Density-Based Algorithm for Discovering Clusters in Large Spatial Databases with Noise", KDD-96
* Martin Ester, Hans-Peter Kriegel, Jorg Sander, Xiaowei Xu, 1996, "A Density-Based Algorithm for Discovering Clusters in Large Spatial Databases with Noise", KDD-96
* Martin Ester, Hans-Peter Kriegel, Jorg Sander, Xiaowei Xu, 1996, "A Density-Based Algorithm for Discovering Clusters in Large Spatial Databases with Noise", KDD-96
* Martin Ester, Hans-Peter Kriegel, Jorg Sander, Xiaowei Xu, 1996, "A Density-Based Algorithm for Discovering Clusters in Large Spatial Databases with Noise", KDD-96
* Martin Ester, Hans-Peter Kriegel, Jorg Sander, Xiaowei Xu, 1996, "A Density-Based Algorithm for Discovering Clusters in Large Spatial Databases with Noise", KDD-96
* [1] Martin Ester, Hans-Peter Kriegel, Jorg Sander, Xiaowei Xu, 1996, "A Density-Based Algorithm for Discovering Clusters in Large Spatial Databases with Noise"
군집분석 결과를 현장에 적용하려면 각 군집의 특성을 이해하고, 군집별로 차별화된 대응 전략을 도출할 수 있어야만 합니다. 이번 포스팅에서는 군집분석 결과를 해석하는 3가지 방법을 소개하겠습니다.
(1) 군집별 변수별 중심 좌표 (Center Points by Clusters)
(2) 군집별 변수별 평행 좌표 그림 (Parallel Coordinates Plot by Clusters)
(3) (차원축소 후) 군집별 산점도 (Scatter Plot by Clusters using Reduced Dimension)
(0) 샘플 데이터 준비 : USArrests {datasets}
예제로 사용할 샘플 데이터는 R datasets 패키지에 내장되어있는 USArrests 데이터셋을 사용하겠습니다. 미국 내 50개 주별 범죄율 관련된 살인(murder), 폭행(assault), 도시인구(UrbanPop), 강간(rape) 의 4개 변수 통계 데이터가 들어있습니다.
군집분석은 변수별 척도에 민감하므로 scale() 함수를 사용해서 표준화(standardization)를 해주었습니다.
## importing USArrests dataset from {datasets} data(USArrests)
그런데 kmeans_k4$centers 로 조회한 4개 각 군집의 변수별 중심 좌표 (평균) 은 여기서는 표준화한 후의 데이터에 대한 변수별 평균입니다. 표준화하기 전의 원래 데이터의 군집별 변수별 중심 좌표 (평균)을 알고 싶으면, 아래의 예시처럼 원래의 데이터에 군집화 결과 변수를 추가해주고(미국 주 이름별로 이미 정렬이 된 상태이기 때문에 merge가 아니라 그냥 칼럼을 추가해주었습니다.), 각 군집별 변수별 평균을 구해주면 됩니다.
(2) 군집별 변수별 평행 좌표 그림 (Parallel Coordinates Plot by Clusters)
다변량 데이터셋으로 군집분석을 했을 때 다수 변수들에 대해 관측치 개개의 값들이 군집별로 어떻게 다른지를 보려면 평행 좌표 그림을 사용할 수 있습니다.
GGally 패키지의 ggparcoord() 함수를 이용해서 ggplot2 기반의 군집별 변수별 평행좌표그림 (parallel coordinates plot)을 그려보겠습니다.
표준화하기 전의 원본 데이터를 사용하고, scale = "std" 를 해주면 알아서 각 변수들을 표준화해줍니다.
(혹은 표준화한 데이터를 바로 사용해도 됩니다).
단, 관측치 개수가 너무 많으면 하나의 평행좌표그림에 모든 관측치를 표현할 경우 서로 겹쳐보여서 군집 특성을 파악하는 것이 힘들 수도 있습니다.
ggplotly 패키지의 ggplotly() 함수를 사용하면 interactive plot 을 그려서 군집별로 하이라이트해서 좀더 수월하게 그래프를 볼 수도 있습니다. (저는 R 4.02 버전을 쓰고 있는데요, 이 버전에서는 ggplotly 가 지원이 안된다고 나오네요 -_-;)
p <- ggparcoord(data = USArrests, columns = c(1:4), groupColumn = "cluster", scale = "std") + labs(x = "Violent Crime Rates by US State", y = "value in scaled", title = "Parallel Coordinates Plot by Cluster")
k-means 군집분석한 결과를 R의 Factoextra 패키지의 fviz_cluster() 함수를 사용하면 주성분분석(PCA) 으로 다변량을 2차원으로 축소해주고, 이를 군집별로 색깔과 모양을 달리하여 R ggplot2 기반의 세련된 시각화를 해줍니다.
아래 예의 경우 Dim1 (62%), Dim2 (24.7%) 라고 되어있는데요, 이는 전체 분산 중 가로 축 Dim1 이 62%를 설명하고, 세로 축 Dim2 가 24.7% 를 차지한다는 뜻입니다. USArrests 데이터셋의 경우 '살인(murder)', '폭행(assault)', '강간(rape)'의 3개 변수가 "범죄"와 관련이 있어서 Dim1에 의해 많이 설명이 되고, 나머지 변수인 '도시인구(UrbanPop)'가 Dim2에 의해 설명이 되겠네요.
혹은 다변량 데이터셋에서 군집을 구분해서 잘 설명해줄 수 있는 대표적인 변수 2개만 선택해서 ggplot2로 직접 2차원 산점도를 그릴 수도 있습니다.
이번 USArrests 데이터셋의 경우 '살인(murder)', '폭행(assault)', '강간(rape)'의 3개 변수가 "범죄"와 관련이 있으므로 이중에서 '살인(murder)' 변수를 하나 뽑고, 나머지 변수인 '도시인구(UrbanPop)'를 나머지 하나로 선택해서 2차원 산점도를 그려보겠습니다.
## Scatter plot by Murder and UrbanPop by clusters library(ggplot2) library(repr) options(repr.plot.width=12, repr.plot.height=12)
군집분석 중에서 k-means clustering, k-median clustering, k-medoid clustering, fuzzy k-means clustering, Gaussian mixture model clustering 등의 알고리즘은 군집의 개수 k를 미리 정해주어야 합니다.
그런데 군집분석이 정답 Y가 없는 상태에서 데이터 구조/군집을 탐색하고 마이닝하는 비지도학습(unsupervised learning)이다 보니 최적의 군집의 개수 k를 결정하고 군집 결과를 평가하는 것이 다분히 주관적인 편입니다. 해당 업과 데이터에 대해서 잘 아는 SME(Subject Matter Expert)의 지식과 경험, 그리고 활용 목적을 고려한 주관적인 의견을 반영하여 최종결정하는 것이 필요하고 중요합니다.
단, 이때 주관적인 경험과 지식에만 의존하는 것이 아니라, 군집의 개수 k를 선택하는데 있어 계량적인 데이터 분석 수치와 시각화의 도움을 받을 수 있다면 좀더 쉽고 또 보다 데이터 구조와 특성을 잘 반영한 의사결정을 할 수 있을 것입니다. (물론, 이마저도 정답이 없다보니 분석가의 주관적인 평가와 판단이 들어가게 되고, 명확한 판단을 내리기가 애매한 경우도 생길 수 있습니다.)
이번 포스팅에서는 군집의 개수 k 를 미리 정해주어야 하는 알고리즘에서 군집의 개수 k를 결정하고 선택하는데 도움을 받을 수 있는 3가지 방법(Determining the number of clusters)에 대해서 알아보겠습니다.
(1) 계층적 군집분석의 덴드로그램 시각화를 이용한 군집의 개수 결정 (Hierarchical Clustering: Dendrogram)
(2) 팔꿈치 방법을 이용한 군집의 개수 결정 (The Elbow Method)
(3) 실루엣 방법을 이용한 군집의 개수 결정 (The Silhouette Method)
0. 샘플 데이터 USArrests {datasets} 불러오기 및 전처리
예제로 사용할 데이터는 R의 datasets 패키지에 내장되어 있는, 미국 주별 폭력적인 범죄 비율("Violent Crime Rates by US State)의 통계를 포함한 USArrests 데이터셋입니다. 이 데이터셋은 1973년 미국의 50개 주 별로 10만 명당 살인(Murder), 폭행(Assault), 강간(Rape) 체포 건수와 도시 인구 비율 (Percent Urban Population)의 네 개의 변수로 구성되어 있습니다.
## importing USArrests dataset from {datasets} data(USArrests)
다음으로, 표준화한 데이터에 대해서 50개의 각 주별 쌍의 유사도 행렬 (proximity matrix)로서 유클리드 거리 행렬 (Euclidean distance matrix) 을 dist(method = "euclidean") 함수를 사용해서 계산하였습니다. 자, 이제 군집분석을 위한 데이터셋 준비가 다 되었네요.
(1) 계층적 군집분석의 덴드로그램 시각화를 이용한 군집의 개수 결정 (Hierarchical Clustering: Dendrogram)
계층적 군집분석(hierarchical clustering)에는 군집과 군집 간의 거리를 측정하는 다양한 연결법(linkage method)으로서, 단일연결법, 완전연결법, 평균연결법, 중심연결법, Ward연결법 등이 있습니다. 그중에서 Ward 연결법(Ward Linkage Method)에 의한 계층적 군집분석을 해보고, 이 결과를 덴드로그램(Dendrogram)으로 시각화를 해보겠습니다.
덴드로그램을 시각화 한 후에 왼쪽의 높이(Height) 축을 기준으로, 위에서 아래로 높이를 이동하면 군집의 개수가 1개, 2개, 4개, 6개, 8개, ..., 50개, 이런 식으로 변하게 됩니다.
이때 군집 간 높이의 차이가 큰 군집의 개수를 선택하면 군집 내 응집력은 높고, 군집간 이질성이 큰 적절한 군집을 구할 수 있습니다. 아래의 예의 경우 군집의 개수가 2개 또는 4개가 적당해 보이네요.
다만, 이 방법은 (a) 데이터 개수가 많을 경우 (가령, 수십만개, 수백만개) 군집화 시간이 오래걸리고, 덴드로그램으로 시각화가 거의 불가능할 수 있습니다. (b) 통계량을 산출하는게 아니고 덴드로그램 시각화 결과를 보고 분석가의 판단에 의지하다보니 다분히 주관적입니다.
set.seed(1004) # for reproducibility hc_ward <- hclust(dist_eucl, method="ward.D")
# 덴드로그램 plot(hc_ward)
(2) 팔꿈치 방법을 이용한 군집의 개수 결정 (The Elbow Method)
다음으로는 비계층적 군집분석(nonhierarchical clustering, partitioning clustering) 중에서 k-means 군집분석(또는 k-median, k-medoids 등)의 k를 다르게 바꾸어가면서 군집화를 했을 때, 군집의 수 k별 "군집 내 군집과 개체 간 거리 제곱합의 총합 (tot.withinss: Total within-cluster sum of squares)" 의 그래프가 팔꿈치 모양으로 꺽이는 지점의 k 를 최적 군집의 개수를 선택하는 방법이 팔꿈치 방법(The Elbow Method) 입니다.
군집 Gi의 중심 측도가 평균인 경우,
Total Within-Cluster Sum of Squares
tot.withinss =
아래의 USArrests 데이터셋 예의 경우 군집의 개수가 4개일 때 Total within-cluster sum of squares 가 팔꿈치 모양으로 꺽이고, 군집 k가 5개 이상일 때부터는 tot.withinss의 변화가 매우 작으므로, 군집의 개수를 4개로 하는 것이 가장 좋아보이네요. (분석가의 주관적인 요소가 들어가긴 합니다. ^^;)
# find the optimal number of clusters using Total within-cluster sum of squares tot_withinss <- c()
for (i in 1:20){ set.seed(1004) # for reproducibility kmeans_cluster <- kmeans(USArrests_scaled, centers = i, iter.max = 1000) tot_withinss[i] <- kmeans_cluster$tot.withinss }
plot(c(1:20), tot_withinss, type="b", main="Optimal number of clusters", xlab="Number of clusters", ylab="Total within-cluster sum of squares")
또는, 군집의 개수 k 별로 전체 분산(Total Variance) 중에서 그룹 간 분산(Between-Group Variance)의 비율로 계산(F-test 라고도 함)하는 "설명된 분산의 비율(the Percentage of Variance Explained)"을 계산하여 시각화한 후에, 팔꿈치 모양으로 꺽이는 지점의 군집의 개수 k를 선택하는 방법을 사용할 수 도 있습니다. (이 경우 위의 그룹 내 분산(Within-cluster Sum of Squares)을 사용했을 때와는 그래프 모양이 반대가 됩니다.)
## the percentage of variance explained as a function of the number of clusters ## Percentage of variance explained is the ratio of the between-group variance to the total variance, ## also known as an F-test. r2 <- c()
for (i in 1:20){ set.seed(1004) # for reproducibility kmeans_cluster <- kmeans(USArrests_scaled, centers = i, iter.max = 1000) r2[i] <- kmeans_cluster$betweenss / kmeans_cluster$totss }
plot(c(1:20), r2, type="b", main="The Elbow Method - Percentage of Variance Explained", xlab="Number of clusters", ylab="Percentage of Variance Explained")
(3) 실루엣 방법을 이용한 군집의 개수 결정 (The Silhouette Method)
Peter J. Rousseeuw 는 "Silhouettes: a grapical aid to the interpretation and validation of cluster analysis" (1987) 라는 논문에서 분할적 군집분석(partitioning clustering)의 결과를 시각화하여 해석하고 평가할 수 있는 지표로서 실루엣(Silhouette)을 제안하였습니다.
실루엣을 계산하기 위해서 먼저 아래의 용어를 정의합니다. (아래의 논문 예시 그림 참조)
군집 A가 객체 i와 다른 객체들을 포함하고 있을 때,
a(i) = 군집 A 내의 객체 i 와 군집 A 내의 다른 객체들 간의 평균 비유사성 (즉, 거리 (distance))
(average dissimilarity of i to all other objects of A)
d(i, C) = 군집 A 내의 객체 i와 다른 군집 C에 속한 모든 객체와의 평균 비유사성 (즉, 거리)
(average dissimilarity of i to all objects of C)
b(i) = d(i, C) 의 최소값 =
* source: Reter J. Rousseeuw (1987)
실루엣은 아래 공식을 보면 알 수 있는 것처럼, 밀집성(tightness)와 분리성(separation)의 비교에 기반한 지표로서,실루엣은 어떤 객체가 군집 안에 잘 놓여 있는지, 그리고 어떤 객체는 단지 군집들 사이의 어딘가에 놓여있는지를 보여줍니다.
전체 군집분석의 실루엣들을 하나의 그래프에 모아서 보여줌으로써, 데이터 구성과 형상의 개요(overview of the data configuration)와 군집들의 상대적인 질(the relative quality of the clusters)을 평가할 수 있습니다.
평균 실루엣 폭(the average silhouette width)은 군집화 결과에 대한 타당성(validity)을 평가하고, 또 적절한 군집 개수를 선택(to select an 'appropritate' number of clusters)하는데 활용할 수 있습니다.
실루엣은 -1 <= s(i) <= 1 사이의 값을 가집니다. 위 공식에 의하면"실루엣 값이 1에 가까우면 다른 군집과의 거리보다 동일 군집 내 객체 간 거리가 가깝다는 뜻이므로 객체 i가 잘 군집화가 되었다는 의미"입니다. (이와 반대로 실루엣 값이 -1에 가까우면 객체 i의 군집화가 잘못되었다는 뜻입니다.)
만약 대부분의 객체가 높은 실루엣 값을 가진다면 군집의 구성이 적절하게 잘 되었다고 평가할 수 있습니다. 반면 많은 객체들이 낮거나 혹은 음(-)의 실루엣 값을 가진다면 군집 구성이 너무 많거나 또는 너무 적은 수의 군집을 가졌기 때문이라고 판단할 수 있습니다.
(아래 예에서는 군집4 에 정렬된 객체들의 뒷부분에 낮은 실루엣 값, 또는 음의 실루엣 값을 가진 객체들이 여러개 보이네요.)
## fviz_silhouette: Visualize Silhouette Information from Clustering install.packages("factoextra") library(factoextra)
