[R 지리공간 데이터 분석] sf 패키지를 활용한 벡터 데이터 세계지도 시각화 (vector data map visualization with sf package)

지난번 포스팅에서는 지리공간 데이터의 포맷으로서 벡터 데이터(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") library(sf) library(raster) library(spData) library(spDataLarge)

참고로, 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/최대값/결측값의 요약통계량을 계산해 보겠습니다.

##-- columns of spatial object names(world) # [1] "iso_a2"   "name_long" "continent" "region_un" "subregion" "type"  "area_km2"  "pop"       # [9] "lifeExp"   "gdpPercap" "geom" ##-- summary statistics summary(world["lifeExp"]) # lifeExp                 geom     # Min.   :50.62   MULTIPOLYGON :177   # 1st Qu.:64.96   epsg:4326    :  0   # Median :72.87   +proj=long...:  0   # Mean   :70.85                       # 3rd Qu.:76.78                       # Max.   :83.59                       # NA's   :10 summary(world["pop"]) # pop                       geom     # Min.   :5.630e+04   MULTIPOLYGON :177   # 1st Qu.:3.755e+06   epsg:4326    :  0   # Median :1.040e+07   +proj=long...:  0   # Mean   :4.282e+07                       # 3rd Qu.:3.075e+07                       # Max.   :1.364e+09                       # NA's   :10

world 데이터셋의 제일 마지막 'geom' 칼럼에는 177개 국가에 대한 WGS 84 좌표 시스템의MultiPolygon 정보 (좌표 리스트)가 들어있습니다.

##-- geom world$geom # Geometry set for 177 features # geometry type:  MULTIPOLYGON # dimension:      XY # bbox:           xmin: -180 ymin: -90 xmax: 180 ymax: 83.64513 # geographic CRS: WGS 84 # First 5 geometries: #   MULTIPOLYGON (((180 -16.06713, 180 -16.55522, 1... #   MULTIPOLYGON (((33.90371 -0.95, 34.07262 -1.059... #   MULTIPOLYGON (((-8.66559 27.65643, -8.665124 27... #   MULTIPOLYGON (((-122.84 49, -122.9742 49.00254,... #   MULTIPOLYGON (((-122.84 49, -120 49, -117.0312 ...

행과 열의 부분집합을 subset 하는 방법도 기본 R 과 동일합니다. 아래 예에서는 1행~2행, 1열~3열까지의 데이터셋을 subset 해보겠습니다.

##-- subset world_mini = world[1:2, 1:3] world_mini # Simple feature collection with 2 features and 3 fields # geometry type:  MULTIPOLYGON # dimension:      XY # bbox:           xmin: -180 ymin: -18.28799 xmax: 180 ymax: -0.95 # geographic CRS: WGS 84 # # A tibble: 2 x 4 # iso_a2 name_long continent                                                       geom # <chr>  <chr>     <chr>                                                           <MULTIPOLYGON [arc_degree]> #   1 FJ     Fiji      Oceania   (((180 -16.06713, 180 -16.55522, 179.3641 -16.80135, 178.7251 -17.01204, 178.596~ #   2 TZ     Tanzania  Africa    (((33.90371 -0.95, 34.07262 -1.05982, 37.69869 -3.09699, 37.7669 -3.67712, 39.20~

(3) sf 패키지의 plot()으로 기본 지도 시각화하기

 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 으로 나누어준 것은 시각화한 지도 크기에 맞추어서 너무 크지도, 작지도 않은 원을 덮어쓰우기 위해 숫자를 조정한 것입니다.)

##-- overlaying circles with cex argument plot(world["continent"], reset = FALSE) pop_sqrt = sqrt(world$pop) / 10000 world_centroid = st_centroid(world, of_largest = TRUE) plot(st_geometry(world_centroid), add = TRUE, cex = pop_sqrt)

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

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