Chapter 13 geographic visualization

지리정보를 시각화하는 방법을 살펴보겠습니다. 제가 참고한 파일은 아래와 같습니다. 정말 감사드립니다.

저자	URL
lee minho	https://lumiamitie.github.io/r/visualization/leaflet-in-R/
토리	https://m.blog.naver.com/PostView.nhn?blogId=hss2864&logNo=221645854030&proxyReferer=https:%2F%2Fwww.google.com%2F

13.1 데이터 다운 로드

데이터는 아래의 사이트에서 다운로드 합니다. 원데이터를 사이트에 가서 받아 보시는 것을 추천드립니다. 바쁘신 분은 2020년 5월 업데이트된 파일을 제 구글드라이브에 링크해놓았으니 참고하시기 바랍니다. 즉, 둘중 하나만 받으시면 되겠습니다.

데이터 소스	url
원데이터	http://data.nsdi.go.kr/dataset/12942
jinhaslab	https://drive.google.com/file/d/1LR0KwE6Y3DwGMpqw2TIXsw5s0GCUft2Q/view?usp=sharing

원하는 곳에 저장하면 되겠습니다. 그리고 경로를 잘 기억해야 합니다. 저는 아래와 같이 data/gis에 저장했습니다.

shp file

library(rgdal)
library(tidyverse)
library(ggplot2)
library(raster)
library(DT)
dat <- shapefile('data/gis/Z_NGII_N3A_G0010000.shp')
dat_d <- fortify(dat)

데이터를 살펴보겠습니다.long와 lat가 있습니다. 위도 경도입니다. 지도의 좌표입니다.

dat_d %>% 
  slice(1:100) %>%
  DT::datatable()

## Warning in instance$preRenderHook(instance): It seems your data is too big
## for client-side DataTables. You may consider server-side processing: https://
## rstudio.github.io/DT/server.html

다운로드 받고 불러온 파일을 가지고 시각화 해보도록 하겠습니다.

dat %>%
ggplot() + 
  geom_polygon(aes(x=long, y=lat, group=group), 
               fill='white', color='grey')

그런데 아래의 표를 보면 0~16까지 번호가 있습니다. 우리나라 17개 도를 표시한 것 같은데 0이 무엇인지 1이 무엇인지 나와있지 않습니다. 어떻게 알아보면 될까요?

dat_d %>%
  group_by(id) %>%
  count() %>%
  arrange(as.numeric(id)) %>%
  datatable()

우선 지역 데이터와 열지수 데이터를 가져와서 어떻게 지역을 맵핑해보도록 하겠습니다. 아래는 기상 데이터를 저장해 놓은 것을 가져오기 위한 데이터베이스 연결입니다. 상기 과정 및 서버 설정은 R and database에서 다루게 됩니다.

# db connect -------
con <-connect(dbms = "postgresql", 
              connectionString = 
                    paste0("jdbc:postgresql://", 
                    rstudioapi::askForPassword("Server address"), 
                    "/ocdm"),
              user=rstudioapi::askForPassword("Database user"), 
              password=rstudioapi::askForPassword("Database password"))

## Connecting using PostgreSQL driver

저는 이와 같은 DB가 연결되어 있습니다.

location_lkup <- dbGetQuery(con, 
                            "SELECT  region_id, region_name FROM heat.location_lkup") %>%
  unique()
location_lkup

##    region_id region_name
## 1         11        서울
## 3         26        부산
## 4         27        대구
## 6         28        인천
## 9         29        광주
## 10        30        대전
## 11        31        울산
## 12        36        세종
## 13        41        경기
## 18        42        강원
## 33        43        충북
## 38        44        충남
## 44        45        전북
## 54        46        전남
## 70        47        경북
## 84        48        경남
## 98   49 / 50        제주

그렇다면 0-16까지 번호로 되어 있는 id 와 기상정보에 있는 region_id 를 어떻게 정리할 지를 고민해 보도록 하겠습니다. 우선 0번은 어느 지역일까요? 아래와 같이하면 0 번은 제주도라는 것을 알 수 있습니다. 여기서 id == '0' 를 통해 알아본 것입니다. id가 0이면 TRUE, 아니면 FALSE라고 인식하게 되므로, TRUE인 파란색이 id ==0 으로 저장되어 있음을 알 수 있습니다.

dat_d %>%
ggplot() + 
  geom_polygon(aes(x=long, y=lat, 
                   group=group, fill=id == '0'))

이에 location_lkup을 업데이트 할 수 있을 것 같습니다.

location_lkup <- location_lkup %>%
  mutate(id = case_when(
    region_name == '제주' ~ 0
  ))

그러면 나머지도 알아보도록 하겠습니다.

즉, 아래와 같이 id 를 매치하였습니다. 마지막으로 제주를 49/50을 4950으로 변경하였습니다.

	region_id	region_name	id
1	4950	제주	0
2	48	경남	1
3	47	경북	2
4	46	전남	3
5	45	전북	4
6	44	충남	5
7	43	충북	6
8	42	강원	7
9	41	경기	8
10	36	세종	9
11	31	울산	10
12	30	대전	11
13	29	광주	12
14	28	인천	13
15	27	대구	14
16	26	부산	15
17	11	서울	16

그림을 그릴 자료인 dat_d 에

여기서 제공드리지는 못하지만 특정 데이터에서 열사병으로 병원에 입원한 횟수를 지역별, 일별로 계산한 자료가 heat_dz 입니다. 이를 년도별 자료로 변경하여 map data 에 병합하겠습니다.

heat_dz <- heat_dz %>%
  mutate(date = as.Date(date, origin = "1970-01-01") )

heat_summary <- heat_dz %>%
  mutate(months = month(date), year = year(date)) %>%
  group_by(region_id, year) %>%
  summarize(heat = sum(heat)) %>%
  mutate(date =
           make_datetime(year  = year, 
                         month = 8,
                         day   =1))

이를 이용해 mm데이터를 만들었습니다.

mm <- heat_summary %>%
  left_join(location_lkup, by = 'region_id') %>%
  left_join(dat_d %>% mutate(id = as.numeric(id)), by = 'id')

13.2 visualizaton 1

즉 dat_d 파일은 제공드렸기 때문에 갖고 계실 것이고, 여기에 지역별 정보를 합친데이터가 mm 이라고 판단하시면 될 것 입니다. 각자의 데이터를 지역별로 만들어 dat_d 합친다면 mm과 동일한 자료 구조가 될 것 입니다. mm을 갖고 그려보도록 하겠습니다.

mm %>% 
  filter(year %in% c(2002)) %>%
  ggplot() +
  geom_polygon(aes(x=long, y=lat, 
                   group = group, 
                   fill = heat), 
                   color = 'grey') +
  theme_minimal() 

2002년 자료를 만들었습니다.

2020 korea heat

다른 년도도 반복적으로 만들어 볼수 있겠습니다.

heat_map <- function(yrs){
p<- mm %>% 
  filter(year %in% c(yrs)) %>%
  ggplot() +
  geom_polygon(aes(x=long, y=lat, 
                   group = group, 
                   fill = heat), 
                   color = 'grey') +
  theme_minimal() +
  ggtitle(paste0('Year:', yrs)) +
  scale_fill_gradientn(colours=rev(topo.colors(50)),
                       na.value = "transparent",                                                 breaks=c(0, 50, 100),
                       labels=c(0, 50,100),
                       limits=c(0, 100)) 

p
ggsave(plot = p,
       filename = paste0('animation/', yrs, ".png"), 
       device = 'png'
       )  
}

seq(from = 2002, to=2015, by=1) %>% 
  map_df(heat_map)

library(magick)
list.files(path = "animation", pattern = "*.png", full.names = T) %>% 
  map(image_read) %>% # reads each path file
  image_join() %>% # joins image
  image_animate(fps=2) %>% # animates, can opt for number of loops
  image_write("animation/heat_first.gif") # write to current dir

heat korea

13.3 visualization 2

shinyapp 을 통해 시각화를 할 수 있습니다. 간단한 shinyapp을 이용해 보겠습니다. shinyapp은 shinyapp 쳅터에 설명이 있습니다.

library(shiny)
yrs =c(2002:2015)
shinyApp(
  
  ui = fluidPage(
    selectInput("yrs", "Year:",
                choices = unique(yrs)),
    plotOutput("map_korea")
  ),
  
  server = function(input, output) {
    output$map_korea = renderPlot({
      hh <- image_read(paste0('animation/', input$yrs, '.png'))
      plot(hh)
      
    })
  },
  
  options = list(height = 500)
)