서울 생활인구 데이터를 격자로 재할당하기 ver2

이글에서는 서울시에서 KT와 협약을 맺어 공개하는 데이터인 '생활인구'데이터를 격자로 재할당하는 방법을 다룬다.

 

이 글은 전에 올렸던 두 개의 글과 연관이 있다.

 

첫번째로, 2년전에 올렸던 <서울 생활인구 데이터의 격자 재할당>

 

서울 생활인구 데이터를 격자 단위로 재할당하기

생활인구를 건물 연면적에 가중치를 두어 50m 격자로 재할당했는데, 이 방식은 몇 가지 문제가 있다.

1. 건물 연면적을 가중치로 두기 때문에, 공원 처럼 건물이 없는 곳은 인구를 배분하지 않는다. '생활인구'의 개념 자체가 집에 거주하는 인구가 아닌, 움직이는 인구이므로 가중치를 연면적으로만 부여하는 것은 인구가 잘못 배분될 가능성이 크다.

2. 현실적으로 연면적이 제대로 들어있는 데이터를 구하기 어렵다.

 

물론, 생활인구가 워낙 자세한 집계구 단위로 공개되기 때문에 그냥 연면적으로 배분해도 분포를 보는 정도로는 크게 문제는 없다. 오히려 현실적인 두번째 문제, 연면적 데이터를 구하기 어렵다는 점이 문제가 된다.

 

그리고 2년여간 데이터를 종종 볼 일이 있었는데, 자세히 배분하는 것 보다는 250m 격자 수준 정도에서 분석하는 것이 좀 더 현실적이라는 생각이 들었다.

 

 

두 번째 연관되는 글은 며칠 전에 올린 통계청 집계구 인구의 격자 할당

 

통계청 집계구 인구를 격자로 재할당하기

이 글에서는 위의 방식을 이용해서 생활인구를 재할당한다. 

 

집계구 데이터를 사용한다고 할 때, 250m 나 500m 정도의 격자로 할당하면 연면적 가중치 없이 배분해도 어느 정도 분포를 보는데는 문제가 없을 것 같다. 결국 자세한 지역이 궁금할때는 집계구 단위로 되돌아가서 보면 되기 때문이다.

 

지난 글의 코드를 최대한 재사용하여 500m 격자로 생활인구를 배분해보겠다.

 

 

 

 

생활인구 데이터 준비

 

 

생활인구는 서울 열린데이터 광장에서 다운받을 수 있다.

위의 메뉴에서 서울생활인구(내국인)과 통계지역경계(집계구.shp, 2016)을 내려받으면 된다.

이 글에서는 2월과 3월 데이터를 사용했다. 

 

 

다운을 받아서 월별로 폴더를 지정한 후 압축을 풀어두자. 

3월의 데이터가 위의 그림처럼 일별로 들어가 있다.

 

 

 

 

 

본격 격자 배분

 

 

격자 할당에 대한 내용은 지난 글에서 충분히 설명했다. 겹치는 내용은 생략하고 넘어가므로 궁금한 점이 있다면 지난 글을 참고하면 된다.

 

 

시작은 지난번과 거의 같다. 

library(foreach)
library(doParallel)
library(tidyverse)
library(data.table)
library(reshape2)
library(scales)
######### 병렬 처리 ##########
# 코어 개수 획득
numCores <- parallel::detectCores() - 1

R 상에서 병렬 처리가 없으므로 Rcpp 로 보낼 코어 개수만 취한다.

 

 

library(Rcpp)
library(sf)

folderWrite <- "d:/지난번에 만든 건물 250m 그리드 파일이 있는 폴더/" 

#Rcpp 파일 로드(별도로 설명한다)
sourceCpp("distributeValueToGrid.cpp")

#건물이 존재하는 250m 그리드를 읽는다. 기존에 저장한 파일
bldgGrid <- fread(paste0(folderWrite,"bldgGrid_parallel.tsv"),
                            sep = "\t", header = TRUE, stringsAsFactors = FALSE)

#Rcpp 변수에 입력한다. 추후 인구를 할당할 때 사용할 준비작업
bldgGridToSet(as.integer(bldgGrid$gridx), as.integer(bldgGrid$gridy))
rm(bldgGrid)

지난번에 만든 파일을 불러온다. 건물이 존재하는 250m 그리드를 저장해놓은 파일이다. 파일이 없다면, 지난 글(www.vw-lab.com/85)을 참고하자.

 

 

##집계구 경계를 읽는다.

folderAdm <- "d:/서울시 생활인구/통계지역경계(2016년+기준)/"
fileName <- paste0(folderAdm,"집계구.shp")
system.time( admBndry <- fileName %>% read_sf())


#앞에서 읽은 집계구별 총 인구에 집계구를 join한다.
#집계구 도형 순서와 일치하는 집계구~인구 데이터프레임을 만드는 작업이다.
admCD <- as.data.frame(admBndry$TOT_REG_CD)
colnames(admCD) <- c("TOT_REG_CD")


#집계구 경계를 입력한다. 한번만 입력하면 된다.
admcoord <- admBndry %>% st_coordinates()

#서울시 생활인구에서 제공하는 집계구는 MultiPolygon이 아닌 Ploygon 형식이므로 기존 함수에 보내는 변수가 약간 다르다.
##4번째 변수와 5번째 변수를 같게 한다. 기존 코드를 재사용하기 위함
putAdmBoundary(admcoord[,1], admcoord[,2],
                  as.integer(admcoord[,3]), as.integer(admcoord[,4]), 
                  as.integer(admcoord[,4]))
rm(admcoord)
rm(admBndry)

위에서 다운받은 서울시 집계구 경계를 읽는다.

지난번과 약간 차이가 있는데, 통계청에서 제공하는 집계구는 MultiPolygon 형식이고, 서울시 집계구는 Polygon 형식이다. 도형 데이터의 집합 위계에서 한 단계가 없으므로, 기존 함수를 사용할 때 변수를 약간 바꿔주어야 한다. 위의 putAdmBoundary 함수 호출에서 네번째와 다섯번째 변수가 동일한데, 이번에는 그렇게 보내야 제대로 작동한다. 

 

 

#재분배 함수

library(lubridate)

#50m 그리드를 상위 그리드로 재집계 하는 함수
aggregatePopuGrid <- function(gridPopu, fromGrid, toGrid) {
  
  if (fromGrid >= toGrid) {
    print("그리드 설정 오류")
    return(NULL)
  }
  
  temp <- gridPopu %>% mutate(x = as.integer(as.integer(x/toGrid)*toGrid) + (toGrid/2) ,
                            y = as.integer(as.integer(y/toGrid)*toGrid) + (toGrid/2)) %>%
    group_by(x,y) %>%
    summarise(.groups="keep", value = sum(value)) %>%
    ungroup()
  
  return(temp)
  
}

지난번에 사용했던 함수와 동일하다.

이번에는 읽고 -> 50m로 배분한 즉시 -> 다시 500m로 재취합 할 것이므로 함수를 미리 정의한다.

 

 

distributePopulation <- function(fileName, admCD_) {
   
      
  # 통계청 집계구별 인구를 읽는다.
  system.time(popu <- fread(fileName,
                              sep = ",", header = TRUE, stringsAsFactors = FALSE,
                              quote="\"", encoding="unknown",
                              select = c(1,2,4,5), col.names = c("yyyymmdd","hh","admcd","popu"),
                              colClasses = list(character=1:4, numeric=5)))
  
  
  
  hours <- c("00","01","02","03","04","05",
             "06","07","08","09","10","11",
             "12","13","14","15","16","17",
             "18","19","20","21","22","23") 
  ymdchr <- popu[1]$yyyymmdd
  
  numCores <- parallel::detectCores()-1
  
  if (exists("result_")) rm(result_)
  
  for (hour in hours) {
  
    #Rcpp 변수에 입력한다. 집계구 polygon의 일련번호와 일치하는 인구다.
    #뒤에서 그리드에 할당할 때 사용한다.
    #시간별로 데이터가 주어지므로, 한번에 하나의 시간대 데이터를 필터링해서 격자배분
    admJoined <- admCD_ %>% left_join(., popu %>% filter(hh==hour),
                           by=c("TOT_REG_CD"="admcd"))
    putAdmPopu(admJoined$popu)
    
    #격자 배분
    resultHourly <- distributeValue(numCores) 
    
    #곧바로 500m 격자로 재취합한다.
    system.time(resultHourly <- aggregatePopuGrid(resultHourly, 50, 500) %>% mutate(hh = hour))
    
    print(paste0(hour,"시 생활인구 :",sum(resultHourly$value, na.rm = T)))
    
    #누적
    if (exists("result_"))  { 
      result_ <- result_ %>% rbind(resultHourly) 
    }  else {
      result_ <- resultHourly
    }
    
  }
  
  #시각변수 추가
  result_ <- result_ %>% mutate(date = ymd_hms(paste0(ymdchr," ",hh,"0000"))) %>%
    mutate(popu = replace_na(value, 0)) %>% 
    select(date,x,y,popu)
  print(paste0("평균 인구 집계(파일 총계) :",sum(result_$popu)/24))
  
  rm(resultHourly)
  rm(admJoined)
  rm(popu)
  return(result_)
}

 

재분배하는 함수는 지난번과 약간 달라졌다.

우선, 데이터 형식이 다르므로 읽어오는 내용이 달라졌다.

 

생활인구는 다양한 성연령 정보를 담고 있는데, 이 글에서는 4번 열인 총생활인구수만 가져와서 이용한다.

위에서 fread 함수를 보면 기준일ID, 시간대구분, 집계구코드, 총생활인구수의 4개 열만 읽어오도록 한 부분을 확인할 수 있다.

성연령으로 구분된 데이터를 사용하려면 적절한 열 번호를 읽고 결측치 처리를 해주면 된다. 나머지 부분은 여기서 소개하는 코드를 그대로 이용할 수 있다.

 

파일 하나가 24시간의 정보를 모두 담고 있다. 즉, 각 집계구마다 24개의 인구가 있으므로 시간대 구분을 보존하여 격자배분하기 위해 다음과 같은 과정을 거쳤다. 24시간을 각각 처리해야 한다.

 

1. 파일 전체에서 00시 추출(필터링)

2. 00시의 데이터를 격자로 배분

3. 곧바로 500m 격자로 재취합

3. 배분된 결과에 다시 00시 붙임

4. rbind로 누적취합

5. 시간을 바꿔서 반복

 

물론, 이러한 처리 방식은 기존 Rcpp 코드를 그대로 재사용하기 위함이다. 불필요하게 50m로 잘게 나눴다가 500m로 재취합 하는 과정을 거치는데다가, 서울 인구만으로 Rcpp에서 병렬 처리를 하다보니 map-reduce 하는 과정이 쓸데없이 낭비스러워서 효율이 떨어진다. 코드를 수정하고 싶다면, 하나의 파일을 하나의 코어에서 처리하도록 하는 방법이 효율적일 것 같다.

 

 

folderPopu <- "d:/서울시 생활인구/내국인/LOCAL_PEOPLE_202002/"
src_files <- list.files(folderPopu)

if (exists("resultFinal")) rm(resultFinal)

#날짜별로 읽어서 보낸다.
for (file in src_files) {
  
  fileName <- paste0(folderPopu, file)
  resultDaily <- distributePopulation(fileName, admCD)
  
  if (exists("resultFinal"))  { 
     resultFinal <- resultFinal %>% rbind(resultDaily) 
  }  else {
     resultFinal <- resultDaily
  }
    
  print(fileName)
}

rm(resultDaily)

위에서 함수를 정의했으므로, 이제 파일을 읽고 함수를 호출한다.

2월 데이터 전체가 있는 폴더를 읽은 후, 파일 각각을 loop를 돌면서 위에서 정의한 함수를 호출한다.

일별 정보이므로 받는 족족 누적해서 파일을 병합한다.

 

한 달 정보를 처리하는데 몇 분의 시간이 걸린다. 로그가 지저분하게 많이 출력되므로, 계속 보고 있으면 기다리는 시간이 그렇게 지루하지는 않다.

 

folderPopu <- "d:/서울시 생활인구/내국인/LOCAL_PEOPLE_202003/"
src_files2 <- list.files(folderPopu)

#날짜별로 읽어서 보낸다.
for (file in src_files2) {
  
  fileName <- paste0(folderPopu, file)
  resultDaily <- distributePopulation(fileName, admCD)
  
  if (exists("resultFinal"))  { 
     resultFinal <- resultFinal %>% rbind(resultDaily) 
  }  else {
     resultFinal <- resultDaily
  }
    
  print(fileName)
}

rm(resultDaily)

이번에는 3월 정보를 읽어서 변환한다. 같은 변수에 누적시킨다.

 

변환이 모두 끝났다. 저장하자.

 

folderWork <- "d:/work/"


fwrite(resultFinal %>% mutate(popu = round(popu, digits=2)), file=paste0(folderWork,"resultDistributed_500m_Feb_Mar.tsv"),
         quote=FALSE, sep = "\t", row.names = FALSE, col.names = TRUE)

##gis 용으로는 dcast로 변환해서 저장한다.
resultFinalStr <- resultFinal %>% filter(date <= ymd(20200210)) %>%
  mutate(date = as.character(date)) %>%
  mutate(date = paste0("D",substr(date,3,4), substr(date,6,7), substr(date,9,10), substr(date,12,13)))

resultFinalCasted <- resultFinalStr %>%
    reshape2::dcast(x+y~date, value.var = "popu")

popuShp <- st_as_sf(resultFinalCasted, coords = c("x", "y"), crs = 5179) 
st_write(popuShp, paste0(folderWork,"popu500_Feb20200210.shp"), driver="ESRI Shapefile")

처리가 끝난 파일은 인구 부분에서 소수점 2자리 정도만 취해서 저장한다.

 

그리고 shp를 만들기 위해서 reshape2 라이브러리의 dcast 함수를 이용하여 같은 좌표에 여러 시각의 정보를 모두 담는다. 처음에 그냥 저장하려고 하니 dbf에 256열이 한계일 수도 있다고 해서 일단, 2월 10일까지만 필터링했다.

 

 

 

 

 

 

데이터 탐색을 위해 그려보기

 

 

이제 R에서 간단히 확인해보자.

#하루 평균치를 구해서 그려보자.
resultDaily <- resultFinal %>% mutate(date = date(date)) %>%
  group_by(date, x,y) %>%
  summarise(.groups="keep",popu = sum(popu)/24) %>%
  ungroup()


## 그래프 그리기
ggplot(data = resultDaily %>% filter(x>=953000 & x<=957000 & y>=1949000 & y<=1953000),
            aes(x=date, y=popu),
            alpha=0.7, size = 0.7)+
  geom_line()+ 
  geom_area()+  
  facet_grid( -y~x)+
  coord_cartesian(ylim = c(0,14000))+
  theme_classic() +
  theme(axis.title.y.right = element_blank(),                # hide right axis title
        axis.text.y.right = element_blank(),                 # hide right axis labels
        axis.ticks.y = element_blank(),                      # hide left/right axis ticks
        axis.text.y = element_text(margin = margin(r = 0)),  # move left axis labels closer to axis
        panel.spacing = unit(0, "mm"),                       # remove spacing between facets
        strip.background = element_rect(size = 0.5))  +      # match default line size of theme_classic
  labs(x="time", y="popu")

ggsave(paste0(folderWork,"서울전체 생활인구 격자_일부.png"),
       antialias = "default", width = 400, height = 400, unit = c("mm"), dpi = 300)

 

 

두 달치의 모든 시간대를 다 보면 정신없을 것 같아서 다시 일별데이터로 만들었다. 24시간을 모두 더한 뒤 평균값을 구했다.

 

ggplot 에서는 facet_grid 방식을 이용해서 500m 격자당 두 달치의 변화를 보는 형식으로 만들어보았다.

 

 

 

 

 

맨 위와 오른쪽에 있는 숫자들은 각 500m 그리드의 중심점을 가리킨다. 각 그리드에서는 해당 그리드의 일별 변동을 볼 수 있다. 종로와 남산이 있는 영역인데, 종로 부근에서 대구 코로나 확진자가 발생한 2월 말 이후 생활인구가 줄어드는 것을 볼 수 있다. 3월말에는 다시 약간 증가한다.

위의 지역들은 대부분 상업지역이므로 주중에는 인구가 많고 주말에는 인구가 적다. 울퉁불퉁 톱니처럼 표현되었다.

 

저 그래프가 무엇인지 잘 와닿지 않을 것 같아서 서울 전체를 그려봤다.

 

이번에는 그리는데 몇 분 걸린다. 인내심을 가지고 기다려보자.

 

 

ploted <- ggplot(data = resultDaily, aes(x=date, y=popu),
                 alpha=0.7, size = 0.7)+
  geom_line()+ 
  geom_area()+
  facet_grid( -y~x)+
  coord_cartesian(ylim = c(0,14000))+
  theme_classic() +
  theme(axis.title.y.right = element_blank(),                # hide right axis title
        axis.text.y.right = element_blank(),                 # hide right axis labels
        axis.ticks.y = element_blank(),                      # hide left/right axis ticks
        axis.text.y = element_text(margin = margin(r = 0)),  # move left axis labels closer to axis
        panel.spacing = unit(0, "mm"),                       # remove spacing between facets
        strip.background = element_rect(size = 0.5))  +      # match default line size of theme_classic
  labs(x="time", y="popu")


ggsave(plot = ploted, paste0(folderWork,"서울전체 생활인구 격자.png"),
       antialias = "default", width = 440, height = 400, unit = c("mm"), dpi = 300)

 

 

공간과 시간을 한번에 담는 방식이다. 

서울 전체를 보면, 각 그리드에서의 두 달치 일별 변동이 한눈에 들어온다. 큰 모니터에서 보면 좀 더 잘 보인다.

facet_grid 방식은 공간 데이터의 EDA에도 나름 유용하게 사용할 수 있는 것 같다.

 

 

 

 

 

 

 

 

코드 전체는 지난번과 같은 github에 올려두었다.

 

vuski/populationDistribution