Modern Data Science with R_Cpt14

14. Working with spatial data

: 공간데이터를 관리하고 시각화 하기

14.1. Motivation: What’s so great about spatial data?

가장 유명한 초기의 공간데이터 분석은 1954년 의사 John Snow의 런던 콜레라 데이터 분석으로 알려져있다. mdsr 패키지의 CholeraDeaths 데이터를 이용하여 콜레라로 사망한 사람의 위치를 플로팅 해보자.

knitr::opts_chunk$set(
    echo = TRUE,
    message = FALSE,
    warning = FALSE
)
library(mdsr)
library(sp)

plot(CholeraDeaths["Count"])

이를 통해 의사 존 스노우는 콜레라 사망이 도시의 물펌프 주변에서 발생한다는 패턴을 발견했다. 적절한 플롯은 통계 모델보다 더 설득력 있게 다가갈 수 있다.

14.2. Spatial data structures

공간 데이터는 데이터프레임이 아닌 다른 형식으로 저장된다. 본 교재에서는 Shapefile을 사용하였으며, 이는 데이터프레임보다 더 복잡한 데이터 컨테이너다. 데이터 컨테이너로 지리적 컨텐스트를 제공할 수 있으며, 확장자는 .shp, .shx, .dbf 등이 있다.

EX: John Snow의 콜레라 맵 재현하기 다음 링크를 클릭하여 데이터 파일을 다운로드 받고 압축을 푼 뒤 예제를 시작한다.
http://rtwilson.com/downloads/SnowGIS_SHP.zip

library(rgdal)

dsn <- fs::path_wd("SnowGIS_SHP", "")

list.files(dsn)

 [1] "Cholera_Deaths.dbf"          "Cholera_Deaths.prj"          "Cholera_Deaths.sbn"         
 [4] "Cholera_Deaths.sbx"          "Cholera_Deaths.shp"          "Cholera_Deaths.shx"         
 [7] "OSMap.tfw"                   "OSMap.tif"                   "OSMap_Grayscale.tfw"        
[10] "OSMap_Grayscale.tif"         "OSMap_Grayscale.tif.aux.xml" "OSMap_Grayscale.tif.ovr"    
[13] "Pumps.dbf"                   "Pumps.prj"                   "Pumps.sbx"                  
[16] "Pumps.shp"                   "Pumps.shx"                   "README.txt"                 
[19] "SnowMap.tfw"                 "SnowMap.tif"                 "SnowMap.tif.aux.xml"        
[22] "SnowMap.tif.ovr"            

경로를 현재 작업폴더로 설정한 뒤, 다운로드 받은 SnowGIS_SHP폴더 안의 모든 파일을 살펴보았다. Cholera_Deaths 파일 6개와 Pumps 파일 5개가 있는데, 이는 layers 라고 하는 두 개의 shapefile 집합이다. 다음 코드에서 확인할 수 있다.

ogrListLayers(dsn)

[1] "Cholera_Deaths" "Pumps"         
attr(,"driver")
[1] "ESRI Shapefile"
attr(,"nlayers")
[1] 2

Cholera_Deaths 레이어를 로드해보자.

ogrInfo(dsn, layer = "Cholera_Deaths")

Source: "C:\Users\user\Desktop\Lab_drive\R_project\[2023]Modern Data Science With R\SnowGIS_SHP", layer: "Cholera_Deaths"
Driver: ESRI Shapefile; number of rows: 250 
Feature type: wkbPoint with 2 dimensions
Extent: (529160.3 180857.9) - (529655.9 181306.2)
CRS: +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +units=m +no_defs 
LDID: 87 
Number of fields: 2 

결과를 보면 ESRI 형식의 shapefile이며 250개의 데이터 행으로 구성된 것을 확인할 수 있다. R로 로드하기 위해 readOGR() 함수를 사용한다.

CholeraDeaths <- readOGR(dsn, layer = "Cholera_Deaths")

OGR data source with driver: ESRI Shapefile 
Source: "C:\Users\user\Desktop\Lab_drive\R_project\[2023]Modern Data Science With R\SnowGIS_SHP", layer: "Cholera_Deaths"
with 250 features
It has 2 fields

summary(CholeraDeaths)

Object of class SpatialPointsDataFrame
Coordinates:
               min      max
coords.x1 529160.3 529655.9
coords.x2 180857.9 181306.2
Is projected: TRUE 
proj4string :
[+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +units=m +no_defs]
Number of points: 250
Data attributes:
       Id        Count       
 Min.   :0   Min.   : 1.000  
 1st Qu.:0   1st Qu.: 1.000  
 Median :0   Median : 1.000  
 Mean   :0   Mean   : 1.956  
 3rd Qu.:0   3rd Qu.: 2.000  
 Max.   :0   Max.   :15.000  

클래스를 보면 SpatialPointDataFrame이라고 나타나 있는데, 이는 S4 클래스이므로 @표기법을 통해 데이터 슬롯에 액세스할 수 있다.

str(CholeraDeaths@data)

'data.frame':   250 obs. of  2 variables:
 $ Id   : int  0 0 0 0 0 0 0 0 0 0 ...
 $ Count: int  3 2 1 1 4 2 2 2 3 2 ...

각 포인트(지점)에 대한 ID와 사망 횟수가 출력되었다.

14.3. Making maps

14.3.1. Static maps with ggmap

앞서 콜레라 데이터를 ggplot2로 플로팅 한다면 다음과 같다.

cholera_coords <- as.data.frame(coordinates(CholeraDeaths))

ggplot(cholera_coords) +
  geom_point(aes(x = coords.x1, y = coords.x2)) + coord_quickmap()

하지만 위 결과는 컨텍스트가 없다는 점에서 단순 plot() 결과보다 낫진 않다. 우리는 런던 지도위에 포인트들이 오버레이되길 원하고 이는 ggmap을 통해서 가능하다. ggmap은 ggplot2 오브젝트로서 Google지도의 쿼리 결과를 반환한다. 모든 ggmap 오브젝트는 ggplot2 오브젝트이므로 유사한 문법을 사용할 수 있다. 런던 지도위에 포인트를 찍어보자. 주의 사항으로 ggmap 패키지는 library()로 불러온 뒤 register_google(key = “키를 입력하세요”)로 키를 인증받아야만 아래의 get_map 등의 기능을 사용할 수 있다.

library(ggmap)

m <- get_map("John Snow, London, England", zoom = 17, maptype = "roadmap")

Error in `get_googlemap()`:
! Google now requires an API key; see `]8;;ide:help:ggmap::register_googleggmap::register_google]8;;()`.
Backtrace:
 1. ggmap::get_map("John Snow, London, England", zoom = 17, maptype = "roadmap")
 2. ggmap::get_googlemap(...)

플롯 위에 포인트가 나타나지 않았다. 그 이유는 CholeraDeaths 오브젝트의 좌표와 ggmap의 지도 오브젝트 ’m’의 좌표 단위가 다르기 때문이다.

14.3.2. Projections

: 3차원 좌표계를 2차원 표현으로 변환하는 프로세스를 의미한다. 변환 과정에서 모양(각도)과 면적을 보존할 수 있다.

다양한 투영 방법 중 일반적인 것은 Lambert 등각원뿔투영과 Albers 동일면적원뿔투영이다. 전자는 각도를 보존하지만 후자는 둘다 보존하지 못한다. 하지만 심한 왜곡은 최소화된다는 장점이 있다.

library(mapproj)
library(maps)

map("state", projection = "lambert",
    parameters = c(lat0 = 20, lat1 = 50), wrap = TRUE)

map("state", projection = "albers",
    parameters = c(lat0 = 20, lat1 = 50), wrap = TRUE)

Lambert와 Albers 지도

지리적 위치 추적에는 좌표참조 시스템, CRS가 필요하다. R에서는 PROJ.4 라이브러리를 사용하며, proj4string() 함수로 사용할 수 있다.

proj4string(CholeraDeaths) %>% strwrap()

[1] "+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy"
[2] "+units=m +no_defs"                                                                   

1. +proj=tmerc 는 데이터가 Transverse Mercator 프로젝선으로 인코딩 되었음을 나타낸다.
   
2. +ellps=airy 는 Airy 타원체가 사용되었음을 나타낸다.
   
3. +units=m 은 단위는 미터(m)를 사용하였음을 나타낸다.
   
   
   

EPSG코드는 proj4 문자열에 대한 속기(shorthand)를 제공한다. 종류는 다음과 같다.

1. EPSG:4326 - GPS시스템 및 구글 어스의 표준이다.
2. EPSG:3857 - Google Maps, Open Street Maps 등의 tiles4 지도에 쓰이는 Mercator 투영법이다.
3. EPSG:27700 - 영국에서 일반적으로 쓰이는 표준이다.

CRS() 함수는 EPSF코드를 PROJ.4 문자열로 변환해준다.

CRS("+init=epsg:4326")

Coordinate Reference System:
Deprecated Proj.4 representation: +proj=longlat +datum=WGS84 +no_defs 
WKT2 2019 representation:
GEOGCRS["WGS 84",
    ENSEMBLE["World Geodetic System 1984 ensemble",
        MEMBER["World Geodetic System 1984 (Transit)",
            ID["EPSG",1166]],
        MEMBER["World Geodetic System 1984 (G730)",
            ID["EPSG",1152]],
        MEMBER["World Geodetic System 1984 (G873)",
            ID["EPSG",1153]],
        MEMBER["World Geodetic System 1984 (G1150)",
            ID["EPSG",1154]],
        MEMBER["World Geodetic System 1984 (G1674)",
            ID["EPSG",1155]],
        MEMBER["World Geodetic System 1984 (G1762)",
            ID["EPSG",1156]],
        MEMBER["World Geodetic System 1984 (G2139)",
            ID["EPSG",1309]],
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1],
            ID["EPSG",7030]],
        ENSEMBLEACCURACY[2.0],
        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]]],
    USAGE[
        SCOPE["unknown"],
        AREA["World."],
        BBOX[-90,-180,90,180]]] 

이를 이용하여 앞의 콜레라 문제에 적용시켜 보자.

cholera_latlong <- CholeraDeaths %>% spTransform(CRS("+init=epsg:4326"))

bbox(cholera_latlong)

                 min        max
coords.x1 -0.1400738 -0.1329335
coords.x2 51.5118557 51.5158345

드디어 우리가 보기에 익숙한 위도와 경도 값을 얻었기에 포인트를 지도에 투영할 수 있다.

ggmap(m) + geom_point(aes(x = coords.x1, y = coords.x2,
                          size = Count), data = as.data.frame(cholera_latlong))

하지만 여전히 클러스터의 중심 위치도 다르고 포인트도 물가가 아닌 다른 곳에 잘못 찍힌 듯해 보인다. 그 이유는 CholeraDeaths 데이터 셋이 EPSG형식이 아니기 때문이다. 그나마 EPSG:27700과 유사하지만, 해당 형식에서 필요로 하는 +datum과 +towgs84 태그가 없다. 따라서 EPSG:27700으로 먼저 사용하고 첫 의도대로 EPSG:4327으로 투영해보자.

proj4string(CholeraDeaths) <- CRS("+init=epsg:27700")

cholera_latlong <- CholeraDeaths %>%
  spTransform(CRS("+init=epsg:4326"))

snow <- ggmap(m) +
  geom_point(aes(x = coords.x1, y = coords.x2,
                 size = Count), data = as.data.frame(cholera_latlong))

#pumps
pumps <- readOGR(dsn, layer = "Pumps")

OGR data source with driver: ESRI Shapefile 
Source: "C:\Users\user\Desktop\Lab_drive\R_project\[2023]Modern Data Science With R\SnowGIS_SHP", layer: "Pumps"
with 8 features
It has 1 fields

proj4string(pumps) <- CRS("+init=epsg:27700")
pumps_latlong <- pumps %>% spTransform(CRS("+init=epsg:4326"))

snow + geom_point(data = as.data.frame(pumps_latlong),
                  aes(x = coords.x1, y = coords.x2), size = 3, color = "red")

제대로 포인트가 찍히고 있다.

14.3.3. Geocoding, routes, and distances

Geocoding이란 읽을 수 있는 주소를 지리적 좌표로 변환하는 것을 의미한다. ggmap의 geocode() 함수를 통해 실행할 수 있다.

smith <- "Smith College, Northampton, MA 01063"
geocode(smith)

Northampton의 스미스 대학교 검색하기(예시)

앞서 사용한 구글의 geocode는 하루에 검색할 수 있는 쿼리 제한이 있지만, R의 RGoogleMaps 패키지는 제한이 없다. getGeoCode() 기능을 사용해보자.

library(RgoogleMaps)

amherst <- "Amherst College, Amherst, MA"
getGeoCode(amherst)

      lat       lon 
 42.37038 -72.51607 

mapdist() 함수를 사용하여 거리를 검색할 수도 있고, route() 함수를 사용하여 여러 장소들 간의 경로를 찾을 수도 있다. 마지막으로 qmap()은 ggmap과 get_map을 같이 나타낼 수 있게한다.

mapdist(from = smith, to = amherst, mode = "bicycling")

legs_df <- route(smith, amherst, alternatives = TRUE)

head(legs_df) %>%
  select(m, km, miles, miles, seconds, minutes, hours, start_lon, start_lat)

qmap("The Quarters, Hadley, MA", zoom = 12, maptype = 'roadmap') +
  geom_leg(aes(x = start_lon, y = start_lat, xend = end_lon, yend = end_lat),
           alpha = 3/4, size = 2, color = "blue", data = legs_df)

14.3.4. Dynamic maps with leaflet

: Leaflet은 HTML로 대화형 지도를 구축하기 위한 오픈소스 JavaScript 라이브러리다. tidyverse 패키지에 포함되어 있으므로 dplyr과 ggplot2와 동일한 관점에서 이해할 수 있다.

library(leaflet)

white_house <- geocode("The White House, Washington, DC")

map <- leaflet() %>%
  addTiles() %>%
  addMarkers(lng = ~lon, lat = ~lat, data = white_house)

white_house <- white_house %>%
  mutate(title = "The White House", address = "2600 Pennsylvania Ave")

map %>%
  addPopups(lng = ~lon, lat = ~lat, data = white_house,
            popup = ~paste0("<b>", title, "</b></br>", address))

leaflet()으로 맵 위젯을 생성하고 레이어를 추가하였다. 위 예제에서는 OpenStreetsMap의 정적 지도를 포함하는 타일 레이어를 추가하였고, 두 번째 레이어로 포인트 위치를 지정하는 마커를 추가했다.