• Home

국토교통부 실거래가 데이터 크롤링 코드

국토교통부 실거래가 데이터 스크래핑 코드를 공유한다. 블로그에 썼던 데이터 분석에 사용한 데이터는 친분이 있는 분으로 부터 받은 데이터인데, 새로운 매매 데이터가 올라가면서 매번 요청하기 힘들것 같아서 코드를 작성해 봤고, 아래와 같이 간단하게 스크래핑 코드를 만들 수 있었다.

데이터 스크래핑 코드는 항상 그렇듯이 임시방편적이고, 코드가 깨끗하지 않다. 그래서 좀 색다르게 magrittr 방식으로 코딩해 봤으나 그다지 나아 보이지 않는다. ;;

수백만건의 데이터를 아래와 같은 간단한 코드로 수집할 수 있다는게 데이터 크롤링의 재미가 아닐까 생각하며, 다시한번 rvest, stringi, XLConnect, data.table과 같은 유용한 패키지의 존재에 감사한다.

library(rvest)
library(httr)
library(stringi)
library(XLConnect)
library(data.table)

## 엑셀 데이터 다운로딩 부분
# 데이터는 working dir아래 data 디렉토리에 쌓인다. 따리서 data 디렉토리를 만들어 주자!

su <- file("succ.txt", "w")

agent_nm <- 
  "Mozilla/5.0 (Macintosh; Intel Mac OS X 10.10; rv:35.0) Gecko/20100101 Firefox/35.0"

maxidx <- html('http://rt.molit.go.kr/rtFile.do?cmd=list') %>% 
  html_nodes('.notiWhite1 .notiBorad01') %>%
  html_text %>% as.numeric %>% max

for(i in maxidx:1){
  urls_view <- sprintf("http://rt.molit.go.kr/rtFile.do?cmd=view&seqNo=%d", i)
  r <- GET(urls_view,
           user_agent(agent_nm))
  htxt <- html(r, "text")

  if((html_nodes(htxt, "td.notiBorad14")[[2]]  %>% 
        html_text()  %>% 
        stri_trim_both) == '첨부파일이 존재하지 않습니다.') next

  download_tags <- html_nodes(htxt, "td.notiBorad14")[[2]] %>%
    html_nodes('a[href^="javascript:jsDown"]')  

  for(dtag in download_tags){

    dtag %>% html_attr("href") %>% 
      stri_match_all_regex(pattern="javascript:jsDown\\('([0-9]+)','([0-9]+)'\\);") %>%
      .[[1]]  %>%  
          { 
            f_idx <<- .[2] %>% as.numeric
            s_idx <<- .[3] %>% as.numeric
          }

    f_nm <- dtag %>% html_text

    urls <- sprintf("http://rt.molit.go.kr/rtFile.do?cmd=fileDownload&seq_no=%d&file_seq_no=%d", f_idx,s_idx)
    r <- GET(urls,
             user_agent(agent_nm))
    bin <- content(r, "raw")

    #1kb 미만의 데이터는 버림(에러 페이지?) 
    if(length(bin) < 1000) next
    writeBin(bin, sprintf("data/%s",f_nm))
    cat(sprintf("%d, %d\n", f_idx,s_idx), file = su)
    print(sprintf("%d, %d", f_idx,s_idx))
  }
}

close(su)




## 엑셀 데이터에서 테이블을 추출해 하나의 아파트 매매 데이터로 통합하는 코드 
## 연립 다세대, 단독 다가 데이터도 간단한 코드 변환으로 통합할 수 있다. 


f_list <- list.files('data') %>%  
  stri_trans_nfc  %>% 
  .[stri_detect_fixed(.,'매매아파트')]

total_list <- list()

#파일의 각 시트(지역)를 data.table로 변환하고 필요 필드 추가함 
for(xlsf in f_list){
  wb <- loadWorkbook(paste0('data/',xlsf))
  sells <- list()
  fname <- stri_replace_last_fixed(xlsf, '.xls',"")
  yyyymm <- substring(fname, 1, 6)
  typenm <- substring(fname, 9)
  for(nm in getSheets(wb)) {
    df <- data.table(readWorksheet(wb, sheet = nm, header = TRUE))
    df[,`:=`(region=nm, yyyymm=yyyymm, typenm=typenm)]
    df[,`거래금액.만원.`:= stri_replace_all_fixed(`거래금액.만원.`, ',', '')]
    sells[[nm]] <- df
  }
  total_list[[paste0(yyyymm,typenm)]] <- rbindlist(sells)
}

result_sales_dt <- rbindlist(total_list)

#결과 데이터 저장 
save(result_sales_dt, file='result_sales_dt.RData')

실거래가 페이지가 바뀌지 않는 이상 동작을 할 것이고, result_sales_dt객체에 2015년 1월 기준 4,964,004건의 매매데이터가 쌓여 있는 것을 확인해 볼 수 있을 것이다.

사족이나 그동안 스크래핑은 Python과 같은 다른 언어를 사용하거나 webkit 라이브러리를 이용해 브라우저 객체를 통해 스크래핑 하곤 했는데, 이젠 여러 패키지 덕분에 R에서도 매우 간단해졌다. 10년 전만해도 HTTP 프로토콜 책을 옆에 끼고 크롤러 코딩을 했었는데, 격세지감을 느낄 따름이다.

윈도우에서 동작여부가 매우 궁금하나, 이런 부분에 대해서는 case by case라 시간 소모하고 싶지 않다.

그러나 참고를 위해서 코드가 동작하는 환경 정보는 추가한다.

## R version 3.1.2 (2014-10-31)
## Platform: x86_64-apple-darwin13.4.0 (64-bit)
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] RWordPress_0.2-3    knitr_1.9           XLConnect_0.2-10   
## [4] XLConnectJars_0.2-9 stringi_0.4-1       rvest_0.2.0        
## [7] data.table_1.9.4   
## 
## loaded via a namespace (and not attached):
##  [1] bitops_1.0-6    chron_2.3-45    digest_0.6.8    evaluate_0.5.5 
##  [5] formatR_1.0     htmltools_0.2.6 httr_0.6.1      magrittr_1.5   
##  [9] markdown_0.7.4  plyr_1.8.1      Rcpp_0.11.4     RCurl_1.95-4.5 
## [13] reshape2_1.4.1  rJava_0.9-6     rmarkdown_0.5.1 stringr_0.6.2  
## [17] tools_3.1.2     XML_3.98-1.1    XMLRPC_0.3-0    yaml_2.1.13

내가 사는 아파트 매매 가격 분석

아직 무주택자인 관계로 국내 부동산 매매 시장에 매우 관심이 많다. 그러다 몇번의 국토교통부 실거래가를 수집하려 틈틈히 노력을 기울였으나 대부분 전수 데이터를 얻는덴 실패하고 그러다 잊고 지내기를 몇번을 거쳐왔는데 최근 실거래가 조회 페이지에 자료제공항목이 생기면서 이런 저런 과정을 통해 원천 전수 데이터를 얻을 수 있게 되었다.

데이터를 얻어서 할 수 있는 일은 역시 분석을 통한 정보의 습득이다. 그리고 그 분석이 내가 지금 살고 있는 아파트와 관련이 있으면 또한 재미가 있을 것인데 이유는 지금 전세를 살고 있는 아파트의 환경이 꽤 맘에 들기 때문이다.

부동산 데이터는 멀티레벨 분석에 최적화된 데이터 특징을 가지고 있다. 이는 지역별 고정효과가 있고 지역의 아파트 브랜드에 따른 랜덤효과가 분명 있을 것이기 때문이다. 더 나아가서는 층별효과 등이 존재할 것이다.

궁극적으로 하고 싶은 분석은 멀티레벨 분석이나 멀티레벨의 시작점은 역시 회귀모형이니 작은 단위로 모델링이 가능한지 확인하고 간단한 추론 과정도 거쳐보도록 하겠다.

2006년부터 2014년 12월까지 국토교통부에 등록된 전국 아파트 매매 데이터는 약 500만건이다. 이 중에서 현재 살고 있는 신림푸르지오1차 아파트의 전체 매매 데이터 423건을 가져와 분석을 시작해 보겠다.

suppressPackageStartupMessages({
  library(sjPlot)
  library(dplyr)
})

load('myapt.RData')

#대략 어떤 포맷인지?
glimpse(myaptdf)
## Observations: 423
## Variables:
## $ id          (dbl) 1693, 43839, 43840, 43841, 105791, 105792, 105793,...
## $ yyyymm      (chr) "200601", "200602", "200602", "200602", "200603", ...
## $ region      (chr) "서울", "서울", "서울", "서울", "서울", "서울", "서울", "서울", "서...
## $ district    (chr) "서울특별시 관악구 신림동", "서울특별시 관악구 신림동", "서울특별시 관악구 신림동",...
## $ street_main (chr) "1730", "1730", "1730", "1730", "1730", "1730", "1...
## $ street_sub  (chr) "0", "0", "0", "0", "0", "0", "0", "0", "0", "0", ...
## $ name        (chr) "신림푸르지오", "신림푸르지오", "신림푸르지오", "신림푸르지오", "신림푸르지오", ...
## $ dealtype    (chr) "BUY", "BUY", "BUY", "BUY", "BUY", "BUY", "BUY", "...
## $ renttype    (chr) "NA", "NA", "NA", "NA", "NA", "NA", "NA", "NA", "N...
## $ area_excl   (dbl) 138.74, 84.79, 114.84, 114.84, 84.79, 84.79, 84.79...
## $ date        (chr) "11~20", "21~28", "1~10", "21~28", "1~10", "11~20"...
## $ property    (dbl) 63500, 41200, 49000, 49700, 44800, 39000, 40500, 4...
## $ deposit     (dbl) 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ fee_month   (dbl) 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ floor       (dbl) 22, 16, 6, 14, 11, 24, 16, 5, 10, 9, 19, 1, 5, 14,...
## $ build_year  (chr) "2005", "2005", "2005", "2005", "2005", "2005", "2...
## $ roadname    (chr) "남부순환로", "남부순환로", "남부순환로", "남부순환로", "남부순환로", "남부순환...
## $ yyyy        (chr) "2006", "2006", "2006", "2006", "2006", "2006", "2...

myaptdf 데이터의 yyyy 필드는 년도별 효과를 보기위해 필자가 생성한 변수이며, 원본의 yyyymm에서 추출한 필드이다.

모델 설명변수로는 층수, 매매년도, 면적을 넣는다.

fit <- lm(property ~ floor + yyyy + area_excl, data=myaptdf)

summary(fit)$adj.r.squared
## [1] 0.7752066
plot(fit, 1)

plot of chunk _mdl

plot(fit, 2)

plot of chunk _mdl

shapiro.test(resid(fit))
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(fit)
## W = 0.974, p-value = 7.509e-07

오차의 등분산성과 오차의 정규성 그래프를 확인해보면 세 레코드의 아웃라이어 데이터가 존재하는 것을 볼 수 있다. 원인이 뭔지 모르겠지만 예측값에 비해 매우 높은 가격을 주고 매매가 이뤄진 경우라 볼 수 있는데, 이 데이터는 제거하고 다시 모델링을 해본다(가장 고점이였던 2007년 2008년을 고려해 보더라도 이 세 매매건은 말도 안되는 가격에 거래가 된 것이라 볼 수 있다. 예측 가격 대비 적게는 1억 많게는 2.5억 더 주고 거래가 된 케이스니..). 게다가 모델은 오차 정규성 검정도 실패를 했다는 것을 알 수 있다.

fit2 <- lm(property ~ floor + yyyy + area_excl, data=myaptdf[-c(54,67,69),])

summary(fit2)$adj.r.squared
## [1] 0.7923766
AIC(fit2)
## [1] 8191.741
plot(fit2, 1)

plot of chunk unnamed-chunk-2

plot(fit2, 2)

plot of chunk unnamed-chunk-2

shapiro.test(resid(fit2))
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(fit2)
## W = 0.9962, p-value = 0.4254

아웃라이어 데이터를 제거한 이후 다소 바람직한 모델이 만들어 졌다는 것을 알 수 있다. 오차의 분산을 보니 다소 비선형적인 패턴이 존재하나 아래에서 다시 확인해 보도록하겠다. adj. \(R^2\) 값도 올라가서 모델이 데이터의 분산을 설명하는 정도가 다소 높아졌다. 게다가 오차 정규성 검정도 통과했다.

혹시나 특정 변수와 비선형 관계가 있을지 각 단변수대비 예측변수의 관계를 플로팅해봤다.

p <- sjp.lm(fit2, type='pred',printPlot=FALSE)
print(p$plot.list[[1]])

plot of chunk _pred

print(p$plot.list[[3]])

plot of chunk _pred

그래프로 보면 층수변수와 매매가격은 다소 비선형 관계로 보이기도 한다. 5층 미만과 최고층이 중간층보다 낮은 매매 가격대를 보이는 것을 알 수 있는데 상식적으로도 매우 맞는 상황이라 생각한다. 따라서 위로 볼록한 층수 변수의 특징을 반영하기 위해 변수 변환 및 추가를 해본다.

fit3 <- lm(property ~ I(floor^2) + floor + yyyy + area_excl, data=myaptdf[-c(54,67,69),])

summary(fit3)$adj.r.squared
## [1] 0.8294764
AIC(fit3)
## [1] 8110.035
#plot(fit3, 1)

#plot(fit3, 2)

#shapiro.test(resid(fit3))

층변수에 2차 함수 특징을 부여하고 모델링을 해보니 adj. \(R^2\) 값이 이전 모델보다 올랐으며 변수가 추가되었음에도 AIC가 낮아지는 현상이 일어났다. 반면 변수 해석은 다소 어려워진 측면이 있다.

그럼 이제 추론을 해보도록 하자!

library(ggplot2)

res <- predict(fit3, newdata = data.frame(floor=1:24, yyyy='2014', area_excl=114.84),level = 0.80, interval = "confidence")

#층수에 따른 가격 변동 예측 
ggplot(data.frame(floor=1:24, estimated_property=res[,1], estimated_property_l=res[,2], estimated_property_h=res[,3]), aes(floor, estimated_property)) + geom_line() + geom_point() + geom_errorbar(aes(ymax=estimated_property_h, ymin=estimated_property_l),linetype=5)

plot of chunk _param

sjp.lm(fit3,sort.est = FALSE)

plot of chunk _param

다른 변수를 고정시키고 층수 변화를 보면서 가격 변동을 보니 8층에서 20층 사이가 가장 매매가가 높은 층인 것으로 알 수 있으며 꼭데기 층은 5층과 비슷한 가격이며, 1층이 최저가격을 보이는 것을 알 수 있다.

모델의 forest plot을 출력해본 결과 2008년에 매매 가격이 정점이였다가 점차 하락을 하고 2013년을 끝으로 다시 매매 가격이 올라가는 양상을 보이고 있다. 2014년의 가격은 2012년의 가격을 약간 넘어서는 결과를 보여준다. 당연한 이야기지만 1 평방미터가 증가할 때마다 평균 366만원 매매가격 증가 효과를 보이고 있다.

그렇다면 내가 살고 있는 아파트의 현재 매매가격은 얼마가 적당할까?

predict(fit3,newdata = data.frame(floor=23, yyyy='2014', area_excl=114.84), interval = "prediction")
fit lwr upr
56675.66 49270.2 64081.12

5.7억정도로 예상이 되며 급매와 같은 기회를 잡는다면 5억 초반도 가능하다는 생각이 든다. 근데 문제는 그럴만한 총알이 없다는게 문제! 재밋게도 작년 전세 계약때 공인중개사 왈 거래가 최근 많이 줄었지만 층수 고려해서 5억 9천정도에 실거래가 이뤄진다는 말이 생각난다. 데이터로도 봐도 그 말이 크게 틀린말이 아니였음을 보여준다.

2015년 1월에 많은 매매가 이뤄졌다는 뉴스가 올라와서 1월 매매 데이터가 궁금하지려 하고 있다. 과연 가격이 더 떨어져서 많은 매매가 이뤄졌는지 아니면 다른 원인이 있는지 매우 궁금하다.


Ps. 마지막으로 데이터를 모두 수집하고 정리해준 회사 모 매니저님에게 감사의 마음을 전한다. ;

이젠 ggmap으로 네이버지도 기반 시각화를 즐기자!

새해 첫 포스팅을 네이버 지도 API 연동 코드가 ggmap 패키지(개발버전)에 통합된 것을 알리는 것으로 시작한다.

제작년(2013) 11월에 ggmap의 구글 API의 오류를 살펴보다가 한국의 PoI(point of interest)가 포함된 네이버 맵 API를 통합하면 어떨까 해서 약 6시간동안 뚝딱 만들어서 pull request를 보냈던 기억이 난다. 그러나 pull request에 대한 어떠한 피드백도 받지 못하고 있다가 22일전 쯤에 ggmap의 제작자로부터 코드 커밋에 대한 커멘트가 올라왔던 것을 오늘 확인했고, 실제 개발 버전을 설치해 테스트를 해보니 잘 동작하는 것을 보고 이렇게 글을 쓴다.

hw

너무 오랜 기다림이고 잊고있었던 상태라 다소 예상치 못한 즐거움을 준것도 사실임…. ㅎㅎ

이 코드가 오늘도 통합이 안되어 있었다면 ggnavermap과 같은 패키지를 만들려고까지 했는데, 그런 코드 브랜치는 만들 필요가 없어져서 매우 다행이라는 생각을 해본다. 다만 네이버의 지도 static API지원이 앞으로도 쾌적하게 지원되기만을 바랄 뿐이다.

그리고 다음 CRAN 릴리즈가 언제가 될지 모르지만 아래와 같이 개발버전을 사용해도 큰 문제는 없다고 생각한다.

앞으로 네이버맵 static API관련 ggmap의 코드 서포트를 해야될거 같아 보이니 현재 ggmap의 구조도 다시 파악해야 될 것 같다.

모두들 네이버 지도 API로 즐거운 지도 시각화가 되길 바란다. :)

## devtools::install_github("dkahle/ggmap")
library(ggmap)
library(lubridate)

#한국 지진 데이터
#https://dl.dropboxusercontent.com/u/8686172/eq2.tsv
eq <- read.table('eq2.tsv', sep = "\t", header = T, stringsAsFactors = F)
eq$latitude <- unlist(strsplit(eq$latitude, " "))[seq(from = 1, to = nrow(eq),
    by = 2)]
eq$longitude <- unlist(strsplit(eq$longitude, " "))[seq(from = 1, to = nrow(eq),
    by = 2)]


eq$longitude <- as.double(eq$longitude)
eq$latitude <- as.double(eq$latitude)
eq$year <- as.factor(substr(eq$date, 1, 4))
eq$date <- ymd_hm(eq$date)

cent <- c(126.96136, 37.52962)
bmap <- ggmap(get_navermap(center = cent, level = 4, baselayer = "default",
    overlayers = c("anno_satellite", "traffic"), marker = data.frame(cent[1],
        cent[2]), key = "c75a09166a38196955adee04d3a51bf8", uri = "www.r-project.org"),
    extent = "device", base_layer = ggplot(eq, aes(x = longitude, y = latitude)))

bmap + geom_point(aes(size = power, colour = date), data = eq, alpha = 0.7) +
    geom_density2d()