Rで福岡温泉オープンデータのボロノイ図を描く


このエントリーをはてなブックマークに追加

Rでオープンデータを扱ってみる

ぱっと検索したところ、下記が見つかったのでまずやってみたところ

オープンデータを解析する- ggplot2を用いたボロノイ分割で厚木市のコンビニ出店を見てみよう!- - Data Science by R and Python

やや関数の使い方などが変わってきてたようなので(ggmapとか)

練習がてら、九州の温泉のオープンデータを用いてボロノイ図の作成をやってみた

九州の温泉一覧|オープンデータ共有&ダウンロード|LinkData

(用いたデータはkyushu_onsen_kihon_listというやつ)

ただし全部表示すると数が多かったので福岡県に絞って表示した


コード

は下記になる、APIを使って読み込んでもいいのだが練習がてらExcel形式のやつをダウンロードして読み込んだ

packageは入っていない場合は適宜 install.package() などでインストールしておくこと

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
library(readxl)
library(ggplot2)
library(deldir)
library(reshape2)
library(plyr)
library(ggmap)
library(scales)
### 最初の9行目までは不要なので skip=8 で読みとばす, ダウンロードしたファイルのディレクトリは適当に自分のを使う
data <- read_excel("/path/to/kyushu_onsen_kihon_list.xls", skip = 8)
### 緯度経度のデータは数値(double型)で読みたかったので as.numeric() で変換
data[,c(4,5)] <- as.numeric(as.matrix(data[,c(4,5)]))
### 温泉の名前(使わなかったけど一応読み込み), 緯度, 経度, 所在県の列だけ読み込み
data2 <- data[,c(2,4,5,8)]
### headerに名前をつける
names(data2) <- c("name", "latitude", "longitude", "prefecture")
### 全部のデータを表示させると地図が見にくかったので福岡県のデータだけ抜粋
data2 <- subset(data2, data2[,4] == "福岡県")
### mapの中心となる場所を、データの緯度経度それぞれの中央値とした(これがしたかったのでas.numeric()した)
loc=c(median(data2$longitude), median(data2$latitude))
### ggmapを用いてgoogle mapのデータを取得
M=ggmap(get_map(location = loc, zoom = 9, source = "google")) + xlab("") + ylab("")
### 緯度経度情報だけ取得
dd <- data2[,2:3]
### 点データを凸囲いする
hull = chull(dd)
hullPointsIndices <- c(hull, hull[1])
hullPoly_df <- data2[hullPointsIndices,]
### ここまでのデータのプロット、必要ならコメントアウトを外してプロットしてみる
# qplot(longitude, latitude, data = dd) +
# geom_polygon(data = hullPoly_df, colour = 'red', fill = 'red', alpha = I(1/10), size = 1)
### Voronoi division
xrng <- expand_range(range(dd$longitude), .05)
yrng <- expand_range(range(dd$latitude), .05)
dd.deldir <- deldir(x=dd$longitude, y=dd$latitude, rw = c(xrng, yrng))
### plot Delaunay diagram
# qplot(longitude, latitude, data = dd) +
# geom_segment(aes(x = x1, y = y1, xend = x2, yend = y2), size = .5, data = dd.deldir$delsgs, colour = 'red', alpha = .5)
### triangulate with polygons
trilist <- unclass(triang.list(dd.deldir))
tridf <- melt(trilist, id.vars = c('x', 'y'))[, c(3, 1:2)]
names(tridf) <- c('tri', 'x', 'y')
### ここまでのデータのプロット、必要ならコメントアウトを外してプロットしてみる
# qplot(longitude, latitude, data = dd) +
# geom_polygon(aes(x=x, y=y, fill=factor(tri)), data = tridf, colour = 'black', alpha = .5) + scale_fill_discrete(guide = FALSE)
### ここまでのデータのプロット、google mapも利用する。必要ならコメントアウトを外してプロットしてみる
# M +
# geom_point(aes(x=longitude, y=latitude, colour=factor(prefecture)), size=3, data = data2) +
# geom_segment(aes(x = x1, y = y1, xend = x2, yend = y2), size = .25, data = dd.deldir$dirsgs, linetype = 2) +
# theme_bw(base_family = "HiraMaruProN-W4")
##### voronoi with polygons
tilelist <- unclass(tile.list(dd.deldir))
tilelist <- lapply(tilelist, function(l) {
data.frame(x = l$x, y = l$y)
})
tiledf <- melt(tilelist, id.var = c("x", "y"))
names(tiledf) <- c('x', 'y', 'tile')
### 最終的なボロノイ図のプロット
M +
geom_polygon(aes(x = x, y = y, fill = factor(tile)), data = tiledf, colour = 'black', alpha = .3) +
geom_point(aes(x = longitude, y = latitude, colour=factor(prefecture)), size=2, data = data2) +
theme_set(theme_bw(base_size = 12, base_family = "HiraMaruProN-W4"))

下記のような画像が表示されたら成功


その他メモ

たとえばggmap使用でエラーなど出たのでそのあたりメモ

Excelデータなのでreadxl packageを使って read.xls 関数で読み込む(比較的新しいExcel読み込みpackageらしい)

  • readxlパッケージ 1.0.0の主要な変更点 - まだ厨二病
  • 最初の9行目まではいらないデータなので skip=8 で飛ばす
  • headerの名前が登録されてないので names() で登録する
  • 必要なデータは “latitude”, “longitude”, “prefecture”, 一応温泉の名前もとってきてみたが使わなかった


ボロノイ図を書くためにいろいろなパッケージを使う, 無い場合はインストールしておくこと


ggmapを使ってるとエラーが出た


参考リンク先のやり方(get_map()関数で緯度経度のmin maxによる範囲指定)ができない

1
2
loc=c(min(data2$longitude), min(data2$latitude), max(data2$longitude), max(data2$latitude))
M=ggmap(get_map(location = loc, zoom = 12, source = "google")) + xlab("") + ylab("")
  • 上記部分で、エラー: improper location specification, see ?get_map.と出てきてしまう
  • 2015年ごろはできていたみたい、今はできない?
  • stackoverflowが出てきた
  • そのまま読み込むと、全て型がcharacterになってしまうのでなんとかしてfloatとかに変えたい
    • as.numeric()でdouble型になら変換可能


chull() データの点集合の一番外側の点を抜粋する


ボロノイ分割, ボロノイ図

  • ある距離空間上の任意の位置に配置された複数個の点に対して、同一距離空間上の他の点がどの母点に近いかによって領域分けされた図のこと
  • とくに二次元ユークリッド平面の場合、領域の境界線は、各々の母点の二等分線の一部になる
  • ボロノイ分割とは? | 用語集とGISの使い方 | 株式会社パスコ
  • なんか細胞と核みたいな見た目だなぁ
  • triang.list
    • class “deldir”
    • return a list of the Delaunay triangles
    • とりあえず、ドロネー図で使われるドロネートライアングル(ドロネー三角形)を形成するポイントを返す?のかな
  • unclass()

    • classを一時的にとっぱらう関数
    • listとかに分解される?
      1
      2
      3
      4
      5
      6
      7
      8
      9
      10
      11
      12
      13
      14
      15
      16
      17
      18
      19
      > y <- data.frame(cbind(x=1, y=1:10))
      > y
      x y
      1 1 1
      2 1 2
      3 1 3
      4 1 4
      5 1 5
      6 1 6
      7 1 7
      8 1 8
      9 1 9
      10 1 10
      > class(y)
      [1] "data.frame"
      > class(unclass(y))
      [1] "list"
  • melt()

    • class reshape2
    • data.frameの複数列の値を、カテゴリ変数1列と値1列の組に変換する
    • これにより、data.frameがやや冗長化するがggplot2での作図にさまざまな操作が行えて便利
    • 本当は、reshape2は古くて今だとdplyr+tidyrでやったほうがいいらしい?
    • tidyr: シンプルなデータ変形ツール - Heavy Watal

ひとまず以上


このエントリーをはてなブックマークに追加