R 지도 시각화 3

움직이는 지도로 표현하기

Dr.Kevin 8/1/2018

지난 포스팅에서 행정경계구역 데이터를 전처리하고 단계구분도를 그려보았습니다. 이번에는 움직이는 지도로 표현하는 방법을 알아보겠습니다. leaflet 패키지를 이용하면 가능합니다.

이번 포스팅에서도 행정경계구역 데이터 전처리를 해야 하는데, 지난 포스팅에서 다룬 것과 다른 내용입니다. 우리나라 정부 기관은 국제 표준(?)인 WGS84 좌표계 대신 GRS80 좌표계로 된 데이터를 공공데이터로 제공하고 있습니다. 따라서 지난 포스팅에서 소개해드린 단계구분도는 별 상관 없이 그릴 수 있지만 실제 지도 위에 시각화하는 것은 좌표계를 바꿔주기 전에는 불가능합니다.

그러므로 GRS80 위경도 좌표를 WGS84 위경도 좌표로 바꿔주는 작업을 먼저 해야 합니다. 그리고 나서 데이터의 행 수가 매우 많으므로 이를 축소해주어야 하는데 지난 포스팅에서 사용한 방법을 쓰면 leaflet() 함수가 제대로 작동하지 않습니다. 따라서 지난 포스팅에서 언급했던 두 번째 방법을 이용하여 shp 파일을 크게 줄이는 방법을 소개해드리겠습니다.

# 필요 패키지를 불러옵니다. 
library(ggplot2)
library(rgdal)
library(sp)

# 단계구분도를 그릴 때 사용할 시/군/구 행정경계구역 데이터를 불러옵니다.
sigg <- 
  readOGR(
    dsn = 'bnd_sigungu_00_2016',
    layer = 'bnd_sigungu_00_2016',
    encoding = 'CP949')

shp 파일을 불러왔으면 WGS84 좌표계로 변환하기 전에 위경도 좌표가 어떻게 생겼는지 육안으로 먼저 확인해보겠습니다. sigg는 s4 객체라서 @를 이용하여 하위 요소를 지정할 수 있습니다. sigg@polygons은 250개의 지방자치단체 경계구역 데이터를 담고 있는 리스트입니다. 첫 번째 지방자치단체의 행정경계구역 데이터가 몇 건인지 확인해보겠습니다.

# 첫 번째 행정경계구역 데이터의 개수를 확인합니다. 
sigg@polygons[[1]]@Polygons[[1]]@coords %>% nrow()

## [1] 2343

첫 번째 요소의 경계데이터는 총 2,343개임을 알 수 있습니다. 그럼 처음 10행만 확인해볼까요?

# 첫 번째 행정경계구역 데이터를 확인합니다. 
sigg@polygons[[1]]@Polygons[[1]]@coords %>% head(n = 10L)

##           [,1]    [,2]
##  [1,] 953683.8 1959210
##  [2,] 953647.3 1959057
##  [3,] 953672.0 1958971
##  [4,] 953780.2 1958977
##  [5,] 953976.8 1958998
##  [6,] 953977.1 1958998
##  [7,] 953994.1 1958990
##  [8,] 954025.5 1958974
##  [9,] 954040.2 1958968
## [10,] 954074.5 1958954

좌표계가 10만 단위의 경도와 100만 단위의 위도로 되어 있습니다. 이제 데이터를 WGS84 좌표계로 바꿔보겠습니다. sp 패키지의 spTransform() 함수를 사용하고 CRSobj 인자에 적당한 값을 할당하면 됩니다.

# GRS80 좌표계를 WGS84 좌표계로 변환합니다.
sigg <- spTransform(x = sigg, CRSobj = CRS('+proj=longlat +datum=WGS84'))

# 같은 행정경계구역 데이터를 다시 확인합니다. 
sigg@polygons[[1]]@Polygons[[1]]@coords %>% head(n = 10L)

##           [,1]     [,2]
##  [1,] 126.9751 37.63118
##  [2,] 126.9747 37.62981
##  [3,] 126.9750 37.62903
##  [4,] 126.9762 37.62909
##  [5,] 126.9784 37.62929
##  [6,] 126.9784 37.62929
##  [7,] 126.9786 37.62921
##  [8,] 126.9790 37.62908
##  [9,] 126.9791 37.62902
## [10,] 126.9795 37.62889

결과를 출력해보니 좌표계가 바뀌어 있습니다. 우리나라 휴전선이 38선으로 불렸다는 것을 떠올린다면 서울의 위도가 38도 아래라는 것을 짐작할 수 있을 것입니다.

이제 sigg 객체를 shp 파일로 저장합니다. 이번에는 writeOGR() 함수를 이용하면 됩니다. dsn 인자에는 shp 파일을 저장할 폴더명을 지정하고 layer 인자에는 파일명을 할당합니다. driver에는 orgDrivers에 의해 반환되는 driver 이름을 문자열로 할당하라고 되어 있지만, 솔직히 잘 모르겠습니다. 그냥 온라인 상에서 확인한 것으로 지정했습니다.

# save sigg as shapefiles
writeOGR(
  obj = sigg,
  dsn = '../shp/simple',
  layer = 'sigungu',
  driver = 'ESRI Shapefile')

dns 폴더에 파일명이 sigungu이고 확장자가 서로 다른 4개의 파일이 저장되어 있는지 확인하시기 바랍니다. 확장자명은 bdf, prj, shp, shx입니다. 제대로 저장되었다면 이제 mapShaper로 이동하여 shp 파일을 축소할 시간입니다. mapShaper를 클릭하면 웹브라우저에 아래와 같은 웹페이지가 열립니다.

화면 상단 중앙에 있는 select 버튼을 클릭하면 Windows 탐색기 또는 MacOS Finder가 열립니다. 방금 저장한 sigungu.shp 파일을 선택하면 아래와 같은 화면으로 바뀝니다.

이 상태에서 그대로 import 버튼을 클릭하면 아래 그림에서와 같이 웹브라우저에 shp 파일로 그린 지도가 출력됩니다. 그리고 오른쪽 상단에 있는 Simplyfy 버튼을 클릭하면 화면 중앙에서 보이는 것처럼 팝업창이 하나 열립니다. 이 때 맨 위에 있는 Simplification menu의 prevent shape removal 항목을 선택합니다. Method는 Visvalingam / weighted area가 기본으로 선택되어 있습니다. 변경하지 말고 이 상태로 Apply 버튼을 클릭하면 됩니다.

새로 열리는 화면 상단에 있는 메뉴를 조절하여 축소할 비율을 선택하면 됩니다. 문건웅 교수님은 웹에서 클릭만으로 하는 통계분석 책에서 1.0%로 설정해도 충분하다고 하셨는데요. 그 가이드를 따라 1.0%로 설정한 후 오른쪽 상단에서 맨 끝에 있는 Export 버튼을 클릭합니다.

마지막 화면입니다. 화면 상단 중앙에 열린 팝업에서 저장할 파일의 형태를 결정합니다. shapefile로 설정하는 것이 가장 무난합니다. 이 상태에서 Export 버튼을 클릭하면 파일이 압축된 상태로 자동 다운로드 됩니다.

압축 파일을 풀면 맨 처음 우리가 저장했던 shp 파일과 동일한 이름의 4개 파일이 저장되어 있음을 확인할 수 있습니다. 본인의 기호에 따라 그대로 쓰셔도 무방하겠지만 저는 편의상 파일명 앞에 simple_을 추가하도록 하겠습니다. 이제 저장한 파일을 불러오도록 하겠습니다.

여기서 한 가지 주의해야 할 점은 dsn에 폴더명 대신 shp 확장자를 포함하여 지정해야 한다는 것입니다. 아직은 그 이유를 파악하지 못 했습니다만, 다음에 기회가 된다면 보완하도록 하겠습니다.

# import simple shp file
spSigg <- 
  readOGR(
    dsn = '../shp/simple/simple_sigungu.shp', 
    layer = 'simple_sigungu',
    encoding = 'CP949')

우리는 shp 파일을 축소하기 전에 첫 번째 지방자치단체의 행정경계구역 데이터 수가 2,343개였다는 것을 알고 있습니다. 얼마나 줄었는지 한 번 확인해 보겠습니다.

# 첫 번째 행정경계구역 데이터의 개수를 확인합니다. 
spSigg@polygons[[1]]@Polygons[[1]]@coords %>% nrow()

## [1] 210

spSigg@polygons 첫 번째 요소의 행 개수가 210으로 출력되었습니다. 거의 1/20 수준으로 줄었음을 알 수 있습니다. 이제 나머지 전처리 과정을 진행하도록 하겠습니다. 한 가지 확인드리고 싶은 것은, 두 번째 포스팅에서는 fortify() 함수를 이용하여 sigg 객체에서 데이터 프레임을 추출했었는데요. 이번 포스팅에서는 spSigg 객체 그대로 사용합니다. 대신 spSigg@data에 모든 정보를 집중하고 leaflet() 함수를 이용하여 시각화할 때 필요한 정보를 꺼내 쓰도록 하겠습니다.

spSigg@data를 보니 mapShaper를 통해 데이터를 축소하니 기존 sigg 객체의 data 요소에 있던 정보들이 모두 사라지고 FID 하나만 남아 있습니다. 따라서 기존 sigg@data의 컬럼 일부를 spSigg@data에 cbind()한 다음 새로 붙인 컬럼의 이름을 편의상 siggCd와 siggNm으로 변경하겠습니다. 그리고 나서 siggCd의 첫 두 글자만 추출하여 sidoCd 컬럼을 새로 만듭니다.

# sigg@data의 모든 컬럼을 spSigg@data에 열기준으로 붙입니다. 
spSigg@data <- cbind(spSigg@data, sigg@data[, 2:4])

# 새로 붙인 컬럼의 이름을 변경합니다. 
colnames(x = spSigg@data)[2:4] <- c('year', 'siggCd', 'siggNm')

# siggCd에서 앞 두 글자만 추출하여 광역시도명을 붙입니다.
spSigg@data$sidoCd <- str_sub(string = spSigg@data$siggCd, start = 1, end = 2)

# 처음 10행만 미리보기 합니다. 
head(x = spSigg@data, n = 10L)

##   FID year siggCd   siggNm sidoCd
## 0   0 2016  11010   종로구     11
## 1   1 2016  11020     중구     11
## 2   2 2016  11030   용산구     11
## 3   3 2016  11040   성동구     11
## 4   4 2016  11050   광진구     11
## 5   5 2016  11060 동대문구     11
## 6   6 2016  11070   중랑구     11
## 7   7 2016  11080   성북구     11
## 8   8 2016  11090   강북구     11
## 9   9 2016  11100   도봉구     11

다음으로 선거 결과 데이터인 result 객체와 병합할 때 기준 컬럼을 만들어야 합니다. 우선 광역시도 행정경계구역 데이터를 불러온 후, SIDO_CD와 SIDO_NM 컬럼만으로 sidoGb 객체를 하나 만들고 컬럼명을 sidoCd와 sidoNm으로 변경한 다음 spSigg@data와 병합합니다. 그리고 나서 sidoNm과 siggNm을 붙인 areaNm 컬럼을 새로 만듭니다. 이 부분은 지난 포스팅에서 소개해드린 것과 동일합니다.

# 광역시도명 코드와 이름을 얻기 위해 시도 행정경계구역 데이터를 불러옵니다.
sido <- 
  readOGR(
    dsn = 'bnd_sido_00_2016',
    layer = 'bnd_sido_00_2016',
    encoding = 'CP949')

# 광역시도명만 추출합니다.
sidoGb <- sido@data[, c('SIDO_CD', 'SIDO_NM')]
colnames(x = sidoGb) <- c('sidoCd', 'sidoNm')

# spSigg@data에 sidoGb를 병합합니다.
spSigg@data <- 
  merge(x = spSigg@data, 
        y = sidoGb, 
        by = 'sidoCd', 
        all.x = TRUE)

# 컬럼 순서를 변경합니다. 
spSigg@data <- spSigg@data[, c(2:5, 1, 6)]

# sidoNm과 siggNm을 합친 areaNm 컬럼을 생성합니다.
spSigg@data$areaNm <- str_c(spSigg@data$sidoNm, spSigg@data$siggNm, sep = ' ')

이제 해야 할 일은 첫 번째 포스팅에서 저장했던 선거 결과 객체가 담긴 RDS 파일을 불러오는 일입니다.

# 데이터가 저장된 폴더로 작업경로를 변경합니다. 
setwd(dir = 'RDS가 저장된 폴더를 지정하세요')

# RDS 파일을 읽습니다.
result <- readRDS(file = '2017_Presidential_Election_Result.RDS')

선거 결과 데이터와 지금까지 작업한 spSigg@data$areaNm과 서로 같은지 교집합 원소의 개수를 확인해보겠습니다. 250이면 모든 원소가 다 겹치는 것입니다.

# 선거결과 데이터의 지역명과 spSigg@data$areaNm이 서로 같은지 확인합니다.
intersect(x = result$지역명, 
          y = spSigg@data$areaNm %>% unique() %>% sort()) %>% length()

## [1] 250

다행하게도 250이 출력되었습니다.

드디어 행정경계구역 데이터 전처리 마지막 단계에 도달했습니다. result 객체에서 후보자별 득표율, 최대R, 색상 및 시청/군청/구청명으로 받아놓은 위경도 좌표 lon과 lat 컬럼을 spSigg@data에 병합합니다.

# spSigg@data에 선거 결과 데이터를 붙입니다. 
spSigg@data <- 
  merge(x = spSigg@data, 
        y = result[, c(1, 22:30)], 
        by.x = 'areaNm', 
        by.y = '지역명', 
        all.x = TRUE)

# 컬럼 순서를 변경합니다. 
spSigg@data <- spSigg@data[, c(2:7, 1, 8:16)]

# FID 순서로 내림차순 정렬합니다. 
spSigg@data <- spSigg@data %>% arrange(FID)

이로써 전처리 과정이 모두 완료되었습니다. 이제 spSigg 객체를 RDS 파일로 저장합니다.

# RDS로 저장합니다.
saveRDS(object = spSigg, file = 'spSigg for leaflet.RDS')

leaflet 패키지를 활용한 지도 시각화

Environment 패널이 텅 비었다고 가정하고 방금 저장한 RDS 파일을 불러오도록 하겠습니다.

# 데이터가 저장된 폴더로 작업경로를 변경합니다. 
setwd(dir = 'RDS가 저장된 폴더를 지정하세요')

# RDS 파일을 읽습니다.
spSigg <- readRDS(file = 'spSigg for leaflet.RDS')

leaflet 패키지를 이용하면 움직이는 지도 시각화를 할 수 있습니다. 함수는 아주 간단합니다. 자세한 내용이 필요하신 분은 별도로 찾아보시기 바랍니다. 지도 시각화는 두 가지를 선보이도록 하겠습니다.

시/군/구청 위치에 작은 원을 출력하고, 그 원을 출력하면 팝업이 생겨 후보별 득표율이 보이도록 합니다.
단계구분도를 그리고 최다 득표 정당 색으로 채운 뒤, 다각형을 클릭하면 후보별 득표율이 저장된 팝업을 띄웁니다.

먼저 leaflet 패키지를 불러온 뒤 첫 번째 코드를 실행하겠습니다.

# 필요 패키지를 불러옵니다.
library(leaflet)

# 시/군/구청 위치에 작은 원을 그리고 클릭하면 팝업이 생기도록 설정합니다.
spSigg %>% 
  leaflet() %>% 
  setView(
    lng = median(x = spSigg@data$lon),
    lat = median(x = spSigg@data$lat),
    zoom = 7) %>% 
  addProviderTiles(provider = 'CartoDB.Positron') %>% 
  addCircles(
    lng = ~ lon,
    lat = ~ lat,
    radius = ~ exp(x = 7),
    weight = 1,
    color = 'black',
    opacity = 1.0, 
    popup = ~ str_c('[', areaNm, '] ', '문재인 : ', 문재인R, ', ', '홍준표 : ', 홍준표R),
    fillColor = ~ 색상,
    fillOpacity = 0.5)

마지막으로 두 번째 코드를 실행하겠습니다.

# 지방자치단체 경계구역 데이터로 다각형을 만들어 단계구분도를 그리고 
# 다각형을 클릭하면 팝업이 생기도록 설정합니다.
spSigg %>% 
  leaflet() %>% 
  setView(
    lng = median(x = spSigg@data$lon),
    lat = median(x = spSigg@data$lat),
    zoom = 7) %>% 
  addProviderTiles(provider = 'CartoDB.Positron') %>% 
  addPolygons(
    weight = 1,
    smoothFactor = 0.5,
    color = 'black',
    opacity = 1.0, 
    popup = ~ str_c('[', areaNm, '] ', '문재인 : ', 문재인R, ', ', '홍준표 : ', 홍준표R),
    fillColor = ~ 색상, 
    fillOpacity = 0.5)

위에서 보이는 지도 이미지는 shinyapps.io에 별도로 올린 뒤 <iframe>으로 삽입한 것입니다. shinyapps.io에 올린 이미지는 작게 등록되기 때문에 답답한 면이 있는데요. 이를 해소하는 방법으로 RPubs로 게시하는 것입니다. 여기를 클릭하면 화면 전체 크기의 지도를 확인할 수 있습니다.

이로써 선거 결과 지도 시각화를 3회에 걸쳐 모두 마무리하였습니다. 번외로 선거관리위원회 홈페이지에서 제공되고 있는 선거 결과 데이터를 크롤링하는 포스팅을 소개하는 것으로 끝내겠습니다.

Written on August 1, 2018