voronoi

Voronoi分割してコンビニ名で塗り分ける。

下記記事で、geom_polygonを使ってうまく塗り分けられないとのことなので、一つの解決法を提示している。

• http://tomoshige-n.hatenablog.com/entry/2014/08/15/002345
• http://wafdata.hatenablog.com/entry/2015/09/22/114711

library(dplyr)
library(deldir)
library(scales)
library(ggmap)
quartzFonts(maru=quartzFont(rep("HiraMaruProN-W4",4)),kaku =quartzFont(rep("HiraKakuPro-W6",4)))

df <- read.csv("http://linkdata.org/api/1/rdf1s973i/konbini_kari.csv", as.is=TRUE)
head(df)

##   konbini_kari                              address      lat     long
## 1            1 神奈川県厚木市三田南２丁目２１番９号 35.47149 139.3483
## 2            2 神奈川県厚木市妻田南１丁目１７−２５ 35.45299 139.3600
## 3            3     神奈川県厚木市下荻野１０４９−１ 35.48385 139.3415
## 4            4         神奈川県厚木市金田１１９−１ 35.47002 139.3683
## 5            5     神奈川県厚木市妻田東３丁目２−７ 35.45966 139.3597
## 6            6             神奈川県厚木市七沢２３８ 35.44230 139.2991
##   X..purl.org.dc.terms.1.1.isPartOf                      label
## 1                  セブン－イレブン       セブンイレブン三田店
## 2                  セブン－イレブン     セブンイレブン妻田南店
## 3                  セブン－イレブン セブンイレブン厚木下荻野店
## 4                  セブン－イレブン   セブンイレブン厚木金田店
## 5                  セブン－イレブン   セブンイレブン厚木妻田店
## 6                  セブン－イレブン   セブンイレブン厚木七沢店
##         label2
## 1 046-242-1171
## 2 046-223-8455
## 3 046-241-5830
## 4 046-225-5572
## 5 046-221-0676
## 6 046-248-0531

df  <-  df[,c(3,4,5)]
colnames(df) <- c("lat","lon","store")
xrng <- scales::expand_range(range(df$lon), .05)
yrng <- scales::expand_range(range(df$lat), .05)

# Voronoi分割
res_deldir <- deldir::deldir(x=df$lon,y=df$lat, rw = c(xrng, yrng))
tilelist <- unclass(deldir::tile.list(res_deldir))

## ここまでは下記記事とほぼ同じコードなので詳細はそちら参照のこと
## http://tomoshige-n.hatenablog.com/entry/2014/08/15/002345

tiledf <- dplyr::bind_rows(lapply(tilelist, function(l){
  data.frame(tile = l$ptNum, x = l$x, y = l$y)
  }))

# ptに店の緯度経度とほぼ同じ緯度経度が入っている
tiledf2 <- dplyr::bind_rows(lapply(tilelist, function(l){
  data.frame(tile = l$ptNum, lon = l$pt["x"], lat = l$pt["y"])
  }))

# tiledf2の緯度経度とdfの緯度経度を比較して最も近い店の名前を取得
# 近傍点を返す関数を定義
getNearPoint <- function(x, y){
  tmp <- dplyr::bind_rows(x, y)[,c("lon", "lat")]
  res <- as.matrix(dist(tmp))[-1,1]
  result <- y[which.min(res), "store"]
  return(result)
}

# 近傍点のコンビニ名を入れる
for(i in seq(nrow(tiledf2))){
  tiledf2$store[i] <- getNearPoint(tiledf2[i,], df)
}

# タイル番号でtiledfとtiledf2を結合
tiledf <- tiledf %>% dplyr::inner_join(tiledf2)

##############
## 描画 ######
##############

# ggmapで描画する地図の範囲指定
loc <- c(min(df$lon), min(df$lat), max(df$lon), max(df$lat))
library("ggmap")
ggmap(get_map(location = loc, zoom = 12, source = "google")) + xlab("") + ylab("") + 
  geom_polygon(aes(x = x, y = y, group = factor(tile), fill = factor(store)) ,data = tiledf, color = "black", alpha = 0.3) +
  geom_point(aes(x = lon, y = lat, color = factor(store)), size = 2, data = tiledf)+
  labs(color="", fill="")　+ theme_bw(base_family='HiraKakuProN-W3')