5.3 空間データ

ロンドンにおける1854年のコレラ発生

  • イギリスの内科医ジョン・スノーは、著書『コレラの伝染形態』(On the Mode of Communication of Cholera) において、コレラ感染による死亡事例の空間的な分布をビジュアル化し、地図の効果的な使い方を示した。

  • スノーは、1854年のコレラ発生時にロンドンのソーホー近郊における感染による死亡事例についての空間ポイント・データを収集し、その情報を地図上にプロットした。

ジョン・スノーによるロンドンにおけるコレラ感染による死亡事例の地図を再現したもの

  • コレラ感染による死亡事例を示す黒い長方形の領域は、地図の中央に位置しているブロード通り (Broad Street) の揚水ポンプの周りに集中している。

    • 地図上にはすべての揚水ポンプが「PUMP」というラベル付きの黒丸で示されている。

  • この地図からスノーは、コレラ感染による死亡事例はブロード通り周辺に集中していることを発見した。

  • 彼は、コレラが下水で汚染された水によって広がったのではないかと推測した。

  • 水道水の詳しい検査や地域住民へのインタビューを含む大規模な調査の後、スノーはブロード通りとケンブリッジ通りとの交差点にある揚水ポンプがコレラ発生源であると結論づけた。

  • また、スノーはサザーク&ボクソール社 (Southwark and Vauxhall Company) の水道水供給にロンドンのコレラ拡大の責任があることを示すために「壮大な自然実験 (grand natural experiment)」を採用した。

  • 自然実験 (natural experiment)とは、研究者の介入はないが実験と似た現実世界の状況を意味する。

  • 以下の図は、スノーが自然実験を行った地域をビジュアル化するために用いた空間ポリゴン地図を再現したものである。

ジョン・スノーによる自然実験の地図を再現したもの

  • ランベス社 (Lambeth Company) は南部地域（赤色で示されている）へ、サザーク&ボクソール社はテムズ川沿岸地域（青色で示されている）へ水道水を供給していた。

  • スノーは、ロンドンのテムズ川周辺の地域で、一部の人々はランベス社から（清潔な）水道水供給を受け、別の人々はサザーク&ボクソール社から（汚染された）水道水供給を受けていたことに注目した。

  • 十分な調査の後、スノーはいかなる交絡因子 (confounder) もこの自然実験には影響を与えなかったと結論づけた。

    • この文脈における交絡因子とは、水道会社とある地域でのコレラ発生率の両方と関連をもつ変数を意味する。

  • つまり、コレラ感染率における違いは水道会社の選択に起因する。

  • コレラにより死亡した人の住所と、彼らに水道水を供給していた会社とを照合することで、スノーは大多数の死亡事例はサザーク&ボクソール社によって水道水を供給されていた家庭で起きたことを明らかにした。

  • スノーの著書は「地図を用いて空間データをビジュアル化することで研究者は過去に知ることのできなかったパターンを発見し、それらを納得のいく形で説明できる」という空間データ分析の威力を示している。

Rにおける空間データ

  • 以下の図は、2008年のアメリカ大統領選挙の結果を効果的にビジュアル化している。

2008年アメリカ大統領選挙での選挙人団の投票先マップ

  • これは、空間ポリゴン・データの例で、各州がいくつかの点をつないで境界を作るポリゴン（多角形）を描く。

  • こうしておくと、バラク・オバマ（ジョン・マケイン）がある州で相対多数の票を得れば、そのポリゴンすなわち州を青（赤）に色づけすることができる。

    • アメリカ特有の選挙人団 (electoral college) 制度については Wikipedia や CNN を参照。

  • R では、maps パッケージが様々な地図作成ツールと多くの空間データベースを提供しており、世界の様々な都市の空間データベースが含まれている。

  • 例えば、us.cities と呼ばれるアメリカ合衆国の都市のデータフレームなどがある。

  • パッケージに収録されているどのデータフレームも、data() 関数を使って読み込むことができる。

library(maps)
data(us.cities)
head(us.cities)
##         name country.etc    pop   lat    long capital
## 1 Abilene TX          TX 113888 32.45  -99.74       0
## 2   Akron OH          OH 206634 41.08  -81.52       0
## 3 Alameda CA          CA  70069 37.77 -122.26       0
## 4  Albany GA          GA  75510 31.58  -84.18       0
## 5  Albany NY          NY  93576 42.67  -73.80       2
## 6  Albany OR          OR  45535 44.62 -123.09       0

  • このデータセットには、以下の変数が含まれている。

    • name：名前

    • country.etc：州

    • pop：人口

    • lat：緯度

    • long：経度

    • capital：首都であれば capital = 1、州都であれば capital = 2、どちらでもなければ capital = 0

  • map() 関数を用いれば、空間データベースにアクセスし、その中のデータをビジュアル化することができる。

  • 例えば、アメリカ合衆国を図示するには、引数 database を "usa" に設定する。

  • 空間上の点は、経度と緯度の情報を x 座標と y 座標の入力値として points() 関数を使って簡単に地図上へ加えることができる。

  • 各州都は実線の円で描かれ、その大きさを cex = capitals$pop / 500000 のように州都の人口に比例させることができる。

  • 地図が描かれた後で title 関数を用いてタイトルを加えることもできる。

map(database = "usa") 
capitals <- subset(us.cities, capital == 2) # 州都を部分集合化

## 経度と緯度を用いて人口に比例した点を追加する
points(x = capitals$long, y = capitals$lat, col = "blue", cex = capitals$pop / 500000, pch = 19)
par(family = "HiraginoSans-W3") #日本語が文字化けしないようにおまじない
title("アメリカ合衆国の州都") # タイトルを追加する

  • 別の例として、カリフォルニア州をプロットしてみる。

  • アメリカ合衆国の空間ポリゴン・データを含む state データベースを用い、引数 regions を "California" に指定する。

map(database = "state", regions = "California")

  • このカリフォルニア州の地図に、州内の人口の多い7大都市を加える。

  • これらの都市をデータから抽出するには、順番が付された指標ベクトルを返す order 関数を用いる。

    • order() 関数では、引数 decreasing を TRUE や FALSE に指定することで降順や昇順で並べ替え（ソート）することができる。

    • sort() 関数は、順序づけられたベクトルそのものを返す。

cal.cities <- subset(us.cities, subset = (country.etc == "CA"))
head(cal.cities)
##              name country.etc    pop   lat    long capital
## 3      Alameda CA          CA  70069 37.77 -122.26       0
## 10    Alhambra CA          CA  88857 34.08 -118.13       0
## 11 Aliso Viejo CA          CA  41975 33.57 -117.73       0
## 15    Altadena CA          CA  43280 34.19 -118.13       0
## 20     Anaheim CA          CA 334909 33.84 -117.87       0
## 27     Antioch CA          CA 109485 37.99 -121.80       0
sind <- order(cal.cities$pop, decreasing = TRUE) # 人口を降順に並べる
top7 <- sind[1:7] # 人口の多い上位7都市
cal.cities[top7, ] # 人口の多い上位7都市のデータを表示
##                 name country.etc     pop   lat    long capital
## 521   Los Angeles CA          CA 3911500 34.11 -118.41       0
## 801     San Diego CA          CA 1299352 32.81 -117.14       0
## 804      San Jose CA          CA  897883 37.30 -121.85       0
## 802 San Francisco CA          CA  723724 37.77 -122.45       0
## 517    Long Beach CA          CA  486571 33.79 -118.16       0
## 778    Sacramento CA          CA  480392 38.57 -121.47       2
## 337        Fresno CA          CA  472517 36.78 -119.79       0
sort(cal.cities$pop, decreasing = TRUE)[1:7] # sort関数は順序づけられた人口を返す
## [1] 3911500 1299352  897883  723724  486571  480392  472517

  • points() 関数を用いてこれらの都市を地図に加え、それらの名前も text() 関数を用いて追加する。

    • cex = 0.6 として文字サイズを小さく（0.6倍に）している。
map(database = "state", regions = "California")
points(x = cal.cities$long[top7], y = cal.cities$lat[top7], pch = 19, col = "blue")

## 円と重なることを避けるため、経度に定数を加える
text(x = cal.cities$long[top7] + 2, y = cal.cities$lat[top7], label = cal.cities$name[top7],
     cex = 0.6, col = "blue")
par(family = "HiraginoSans-W3") #日本語が文字化けしないようにおまじない
title("カリフォルニア州の大都市")

  • ここで、R で空間ポリゴン・データがどのようになっているのか検討する。

  • プロット作成が行われないようにするため、map 関数の引数 plot を FALSE に設定すると、x（\(x\)座標すなわち経度）と y（\(y\)座標すなわち緯度）として保存された座標の数列を伴うリスト・オブジェクトを返す。

  • このリストの中では名前が names として保存されている異なるポリゴンは NA によって分割されている。

usa <- map(database = "usa", plot = FALSE) # 地図を保存
names(usa)  # 要素の一覧
## [1] "x"     "y"     "range" "names"

  • ベクトル x の長さを計算して、アメリカ合衆国の地図を作成するために用いられた座標の数を確認することができる。

  • cbind() 関数を使って\(x\)座標と\(y\)座標を行列にまとめた後、最初の数件の観察値を表示する。

length(usa$x) 
## [1] 7252
head(cbind(usa$x, usa$y)) # ポリゴンの最初の6つの座標
##           [,1]     [,2]
## [1,] -101.4078 29.74224
## [2,] -101.3906 29.74224
## [3,] -101.3620 29.65056
## [4,] -101.3505 29.63911
## [5,] -101.3219 29.63338
## [6,] -101.3047 29.64484

  • アメリカ合衆国の地図は、7252件の座標ペアで構成されていることがわかる。

  • map 関数はこれらの点を連結し、地図を作成する。

  • 空間データは、空間に広がるパターンについての情報を含み、地図を用いてビジュアル化できる。

  • 空間ポイント・データでは事象の地点を地図上の点として表すのに対し、空間ポリゴン・データは地理的領域を地図上の点を連結することで表す。

日本の空間データ

  • 日本の情報を確認したい場合は mapdata パッケージを用いる。
library(mapdata)
map('japan')
map.axes() # 緯度と経度の枠を表示

  • 緯度と経度の範囲を xlim と ylim で指定して拡大することもできる。

  • オプション interior = FALSE で都道府県の境界を非表示にして、boundary = 以降で境界を波線で表示するように指定する。

map("japan", interior = FALSE, xlim = c(138,141), ylim = c(34,36))
map("japan", boundary = FALSE, lty = 2, add = TRUE)
map.axes()

法政大学（BT）の緯度と経度を調べて表示する。

  • 街区レベル（xx町x丁目x番）の緯度と経度は国土交通省からダウンロードできる。

  • より細かい緯度と経度を調べるには、例えばココを使う。

map("japan", interior = FALSE, xlim = c(138.5,140.5), ylim = c(35,36))
map("japan", boundary = FALSE, lty = 2, add = TRUE)
points(139.740631, 35.695374, col = "blue", pch=20)
par(family = "HiraginoSans-W3") #日本語が文字化けしないようにおまじない
text(139.740631, 35.695374 + .05, "法政大学", col = "blue")
map.axes()

Rにおける色

  • 地図に色づけする方法を学ぶ。

  • 一般に、色はビジュアル化に重要である。

  • ここまでは、色は "red" や "blue" のような名前を用いて指定してきた。

  • R で使える657色の名前をすべて見るために、colors() 関数の出力を見る。

allcolors <- colors()

head(allcolors) # いくつかの色
## [1] "white"         "aliceblue"     "antiquewhite"  "antiquewhite1"
## [5] "antiquewhite2" "antiquewhite3"
length(allcolors) # 色の名前の数
## [1] 657

  • R はこれよりもずっと多くの色を生成することができる。

  • 生成可能な色の全範囲からある色を生成するには、16進数カラーコード (hexadecimal color code) を用いることができる。

    • 16進数は基数が16で、0から15までの値を表す整数0-9と英文字A-Fによる記数法。

    • 16進数カラーコードはハッシュマーク (#) で始まる6文字の並び。

  • 2桁ずつで三原色 (Red, Green, Blue: RGB) の各色の強さ（彩度）を表し、それぞれが0から255（すなわち\(2^8\)段階のいずれか）の値をとる。

    • 例えば、紫色は半分の強さの赤と青 RGB = (127, 0, 127) によって生み出される。

    • 127は16進法では7Fに等しいので、#7F007F という16進数カラーコードが導かれる。

  • R では、rgb() 関数によって数値から16進数カラーコードを作成することができる。

  • 3つの引数 red, green, blue は各色の彩度で0から1の値をとり、0から255までの整数値に変換され、16進数で表示される。

    • 16進数のカラー表記を見つけるための情報源がインターネット上にはたくさんある。

  • まずは原色から始める。

red <- rgb(red = 1, green = 0, blue = 0) # 赤
green <- rgb(red = 0, green = 1, blue = 0) # 緑
blue <- rgb(red = 0, green = 0, blue = 1) # 青
c(red, green, blue) # 結果
## [1] "#FF0000" "#00FF00" "#0000FF"
  • 黒と白は各原色をそれぞれ0%と100%にすることで表される。
black <- rgb(red = 0, green = 0, blue = 0) # 黒
white <- rgb(red = 1, green = 1, blue = 1) # 白
c(black, white) # 結果
## [1] "#000000" "#FFFFFF"

  • rgb() 関数の各引数は2以上の長さのベクトルを撮ることができ、一度に2つ以上のカラーコードを作成できる。

  • 紫（50%の赤と50%の青）と黄色（100%の赤と100%の緑）は次のように作ることができる。

rgb(red = c(0.5, 1), green = c(0, 1), blue = c(0.5, 0)) # 紫と黄色
## [1] "#800080" "#FFFF00"

  • 16進数カラーコードの最後に00からFFまでの2桁を足すことで透過色にでき、透過度もコントロールできる。

  • rgb() 関数の4番目の引数 alpha に0から1の彩度の尺度を指定し、16進数カラーコードに変換する。

## 透過色の青
blue.trans <- rgb(red = 0, green = 0, blue = 1, alpha = 0.5) 

## 透過色の黒
black.trans <- rgb(red = 0, green = 0, blue = 0, alpha = 0.5)

  • 16進数カラーがわかれば、"red" や "blue" のようなすでに名前のつけられている色を用いた時と同じように、それらを（文字オブジェクトとして）作図に用いることができる。

  • 次の図では、透過色の円は重なり合ったときでも容易に識別できるが、不透過色の円同士は識別しにくいことがわかる。

## 不透過色の円は識別しにくい
plot(x = c(1, 1), y = c(1, 1.2), xlim = c(0.5, 4.5), ylim = c(0.5, 4.5), 
     pch = 16, cex = 5, ann = FALSE, col = black)
points(x = c(3, 3), y = c(3, 3.2), pch = 16, cex = 5, col = blue)

## 透過色にすると識別しやすい
points(x = c(2, 2), y = c(2, 2.2), pch = 16, cex = 5, col = black.trans) 
points(x = c(4, 4), y = c(4, 4.2), pch = 16, cex = 5, col = blue.trans)

  • この図では視認性向上のため、plot() 関数の引数 ann を FALSE に設定して、デフォルトで表示される軸ラベルを消している。

アメリカ大統領選挙

  • R でどのように色が表示されるかわかったので、地図を色づけする。

  • 例として、2008年の大統領選挙の結果を用いて、アメリカ合衆国の地図を色づけする。

  • 選挙データファイル pres08.csv に含まれる変数の名前と説明は次の通りである。

変数 説明
state 州の略称
state.name 省略されていない州名
Obama オバマの得票率（パーセンテージ）
McCain マケインの得票率（パーセンテージ）
EV その州の選挙人票の数
  • 各州を2つのやり方で色づけする。

  1. オバマが勝利した州に青を、マケインが勝利した州に赤を用いる。

    • こうすると「青の州と赤の州」に色分けされた地図が作成される。

  2. RGBのカラースキームでは青と赤の混合比によって様々な色調の紫が作成できることを活用して、2大政党の得票率を計算し、青の彩度を民主党の2大政党間の得票率として、赤の彩度を共和党の得票率として設定する。

    • こうすることで、州の色が民主党候補や共和党候補への支持の程度を反映したものになる。
  • まず、データを読み込み、2大政党間の得票率を計算し、2大政党の得票率に基づいてカリフォルニア州におけるRGBカラースキームを設定する。
library(maps)
pres08 <- read.csv("pres08.csv")

## 2大政党の得票率
pres08$Dem <- pres08$Obama / (pres08$Obama + pres08$McCain)
pres08$Rep <- pres08$McCain / (pres08$Obama + pres08$McCain)

## カリフォルニア州への彩色
cal.color <- rgb(red = pres08$Rep[pres08$state == "CA"], 
                 blue = pres08$Dem[pres08$state == "CA"], 
                 green = 0)

  • 次に、カリフォルニア州の地図を2つのやり方で色づけする。

  • 地図に色を加えるには引数 col を指定する。

  • さらに各州を指定した色で塗りつぶすために引数 fill を TRUE に設定する。

  1. オバマは2008年の選挙でカリフォルニア州で勝利したので、この州を青く彩色する。
## 青い州としてのカリフォルニア
map(database = "state", regions = "California", col = "blue", 
    fill = TRUE)

  1. RGBカラースキームを使って、2大政党の得票率に基づいてカリフォルニア州を彩色する。
## 紫色の州としてのカリフォルニア
map(database = "state", regions = "California", col = cal.color,
    fill = TRUE)

  • 最後に、この作業を、ループを用いてすべての州について行う。

  • 地図はハワイ、アラスカ、ワシントンD.C.は含んでいないため、これらの州については省略する。

  • 各州に色をつけるため引数 add を TRUE にする点に注意する。

  • 2つの地図を作成するためのコードはほとんど同じで、各州を彩色するのにどの色を用いるかという点が唯一の違いである。

  1. オバマが勝利した州は青、マケインが勝利した州は赤という2択のカラースキームを用いる。
## 赤い州と青い州としてのアメリカ
map(database = "state") # 地図を作成
for (i in 1:nrow(pres08)) {
    if ((pres08$state[i] != "HI") & (pres08$state[i] != "AK") &
        (pres08$state[i] != "DC")) {
        maps::map(database = "state", regions = pres08$state.name[i],
            col = ifelse(pres08$Rep[i] > pres08$Dem[i], "red", "blue"),
            fill = TRUE, add = TRUE)
    }
}

  1. 各州の2大政党の得票率に基づいたRGBカラースキームを用いる。
## 紫色の州としてのアメリカ
map(database = "state") # 地図を作成
for (i in 1:nrow(pres08)) {
    if ((pres08$state[i] != "HI") & (pres08$state[i] != "AK") &
        (pres08$state[i] != "DC")) {
        map(database = "state", regions = pres08$state.name[i], 
            col = rgb(red = pres08$Rep[i], blue = pres08$Dem[i],
                green = 0), fill = TRUE, add = TRUE)
    }
}

  • 赤と青の地図から、オバマが西海岸と東海岸の多くの州で勝利している一方で、マケインは中西部でとりわけ強いことがわかる。

  • しかし、紫色の地図からは民主党か共和党で完全に占められている州はないことがわかる。

  • 各州には両方の支持者がいて、勝者総取り方式の選挙制度のせいで各州が民主党支持の州か共和党支持の州として見えるだけである。

ウォルマートの拡大

  • 小売ディスカウント百貨店・倉庫型店舗として商業的成功を収めているアメリカ発祥の多国籍チェーンのウォルマートの拡大について検討する。

  • ウォルマートは、1962年アーカンソー州ロジャーズ (Rogers, Arkansas) に最初の店舗をオープンさせた。

  • その後の数十年間で、アメリカ国内外において数多くの店鋪をオープンさせ、世界最大の多国籍小売業者の1つになった。

  • ウォルマート店舗の新規開店に関するデータ walmart.csv には、最初の開店である1962年3月1日から2006年8月1日までのウォルマートの店舗の開店に関する空間と時間の情報が含まれている。

  • このデータセットに含まれる変数の名前と説明は以下の表の通りである。

変数 説明
opendate 店舗の開店日
st.address 店舗所在地の住所
city 店舗所在地の市
state 店舗所在地の州
type 店舗のタイプ (Wal-MartStore, SuperCenter, DistributionCenter)
long 店舗所在地の経度
lat 店舗所在地の緯度

  • まず、全店舗の位置を地図上にプロットする。

  • データセットは3つの異なるタイプの店舗の情報を含み、タイプの違いは変数 type によって表されている。

    • Wal-MartStore：通常のウォルマートの店舗

    • SuperCenter：通常のウォルマートにスーパーマーケットが併設された店舗（スーパーセンター）

      • 薬局、園芸店、カー用品店、その他の専門店もしばしば併設されている。

    • DistributionCenter：食品や商品を通常のウォルマートやスーパーセンターへ配送する店舗（流通センター）

  • 3タイプの店舗を区別するために異なる色を用いる。

    • 通常のウォルマート：赤

    • スーパーセンター：緑

    • 流通センター：青

  • 各店舗を表す円が重なっても識別できるように透過色にする。

    • rgb() 関数で alpha = 1 / 3 と設定する。

  • 他の2タイプの店舗よりも数が少ない流通センターは目立つよう大きな円で示す。

walmart <- read.csv("walmart.csv")

walmart$storecolors <- NA # 空のベクトルを作成

## 赤 = 通常のウォルマート, 緑 = スーパーセンター, 青 = 流通センター
walmart$storecolors[walmart$type == "Wal-MartStore"] <- 
    rgb(red = 1, green = 0, blue = 0, alpha = 1/3)
walmart$storecolors[walmart$type == "SuperCenter"] <-
    rgb(red = 0, green = 1, blue = 0, alpha = 1/3)
walmart$storecolors[walmart$type == "DistributionCenter"] <-
    rgb(red = 0, green = 0, blue = 1, alpha = 1/3)

## 流通センターは大きめの円にする
walmart$storesize <- ifelse(walmart$type == "DistributionCenter", 1, 0.5)

  • 地図を作成しウォルマートの店舗の位置をそれに追加する。

  • また、legend() 関数を用いて凡例を追加する。

    • x 座標と y 座標を設定して凡例の場所を指定し、凡例で表示する文章のベクトルを引数 legend として示す。

    • 引数 bty を設定しないデフォルトの状態では凡例は枠線で囲まれるが、bty = "n" と設定すれば枠線は消える。

  • プロットするオブジェクトのタイプを指定するには引数 pch を用いる。

  • ここでは色つきの円を用い、そのサイズは引数 pt.cex でコントロールする。

## 凡例付き地図
map(database = "state")

points(walmart$long, walmart$lat, col = walmart$storecolors, 
       pch = 19, cex = walmart$storesize)

par(family = "HiraginoSans-W3") #日本語が文字化けしないようにおまじない

legend(x = -120, y = 32, bty = "n", 
       legend = c("ウォルマート", "スーパーセンター", "流通センター"),
       col = c("red", "green", "blue"), pch = 19, # 色付きの実線の円
       pt.cex = c(0.5, 0.5, 1)) # 円の大きさ

  • この地図はウォルマートの経営戦略をはっきりと示している。

  • スーパーセンターは中西部と南部全体に広く散らばっている。

  • 一方、北東部や西海岸および都市部ではそれほど広がっておらず、通常のディスカウントストア型を超える拡大は行っていない。

Rにおけるアニメーション

  • 前回のウォルマート店舗の分析では、以下のデータセットを用いた。
library(maps)

walmart <- read.csv("walmart.csv")
head(walmart)
##     opendate                    st.address              city state      long
## 1 1962-03-01 5801 SW Regional Airport Blvd       Bentonville    AR -94.23982
## 2 1962-07-01              2110 WEST WALNUT            Rogers    AR -94.07141
## 3 1964-08-01              1417 HWY 62/65 N          Harrison    AR -93.09345
## 4 1965-08-01             2901 HWY 412 EAST    Siloam Springs    AR -94.50208
## 5 1967-10-01        3801 CAMP ROBINSON RD. North Little Rock    AR -92.30229
## 6 1967-10-01         1621 NORTH BUSINESS 9         Morrilton    AR -92.75858
##        lat               type
## 1 36.35088 DistributionCenter
## 2 36.34224        SuperCenter
## 3 36.23698        SuperCenter
## 4 36.17990        SuperCenter
## 5 34.81327      Wal-MartStore
## 6 35.15649        SuperCenter

  • データセットは開店日の情報 opendate を含んでいるにもかかわらず、前回の分析では時間の次元を無視していた。

  • 空間パターンのみではなく時空間パターンを調べることで、ウォルマートが時間とともにどのように拡大していったかをより良く理解することができる。

  • ここでは様々な時点での全ての営業店舗を表す一連の地図を作成する。

  • そのために、特定の日付のデータを部分集合化し、前回行ったようにウォルマート店舗の地図を作成する関数を定義する。

  • 前回の一連のコードを walmart.map() という以下の2つの入力値をとる関数として作成する。

    1. data：店舗の開店日を示す opendate という名の変数 (Date クラスでなければならない) を含むデータフレーム

    2. date：地図が作成される時点を定義する別の Date オブジェクト

  • この関数は指定した日付以前に開店したすべての店舗を部分集合化し、それらの場所を地図上にプロットする。

walmart.map <- function(data, date) {
    walmart <- subset(data, subset = (opendate <= date))
    map(database = "state")
    points(walmart$long, walmart$lat, col = walmart$storecolors, 
           pch = 19, cex = walmart$storesize)
}

  • この関数を使えば、ある時点における地図を容易に作成できる。

  • ここでは地図を10年ごとに作成する。

walmart$storecolors <- NA # 空のベクトルを作成

## 赤 = 通常のウォルマート, 緑 = スーパーセンター, 青 = 流通センター
walmart$storecolors[walmart$type == "Wal-MartStore"] <- 
    rgb(red = 1, green = 0, blue = 0, alpha = 1/3)
walmart$storecolors[walmart$type == "SuperCenter"] <-
    rgb(red = 0, green = 1, blue = 0, alpha = 1/3)
walmart$storecolors[walmart$type == "DistributionCenter"] <-
    rgb(red = 0, green = 0, blue = 1, alpha = 1/3)

## 流通センターは大きめの円にする
walmart$storesize <- ifelse(walmart$type == "DistributionCenter", 1, 0.5)

walmart$opendate <- as.Date(walmart$opendate)

walmart.map(walmart, as.Date("1974-12-31"))
title("1975")

walmart.map(walmart, as.Date("1984-12-31"))
title("1985")

walmart.map(walmart, as.Date("1994-12-31"))
title("1995")

walmart.map(walmart, as.Date("2004-12-31"))
title("2005")

  • 上のような時空間データをビジュアル化するその他の方法として、地理的なパターンが時間とともにどのように変化するかを動的に見せる アニメーション (annimation) がある。

  • animation パッケージで様々な場所と時間において、ウォルマートがどのようにその店舗を拡大していったかを示すことができる。

  • アニメーションの準備として以下を行う。

    1. アニメーション化される地図の数を設定する。

    2. データセットの最初から最後まで等間隔の日付で構成されたベクトルを作成する。

n <- 25 # アニメーション化する地図の数

dates <- seq(from = min(walmart$opendate), 
             to = max(walmart$opendate), length.out = n)

  • 基本的には、animation パッケージを使うことと、時間の経過にしたがって一連の地図を作成するループを作ることとの間に大差はない。

  • 実際、もう1つ saveHTML() 関数を追加して、ループ全体をそれでくるむだけでよい。

  • この関数は R のコードをメインの入力として波括弧 {} の中にとり、ループ中に作成されたすべての図をHTMLファイルに挿入する。

  • できあがったHTMLファイルは、ウェブブラウザでアニメーション表示できる。

  • saveHTML() 関数には次のような便利な引数がある。

    • htmlfile：HTMLファイルの名前を指定

    • title：アニメーションのタイトルを指定

    • outdir：完成したファイルが保存されるディレクトリ (フォルダ) 名を指定

    • autobrouwse：出力が自動でブラウザに表示されるかを管理

  • 次の一連のコードは、アニメーションを作成し、walmart.html という名のHTMLファイルと他のすべてのファイルを作業ディレクトリに保存する。

    • saveHTML() 関数はループ中に先ほど作成した walmart.map() 関数を繰り返し呼び出す。

    • 作業ディレクトリへのパスを返す getwd() 関数を引数 outdir への入力として指定することで、すべての出力ファイルがそのディレクトリへ保存される。

# install.packages("animation") # インストールしていない場合
library("animation")
saveHTML({
  for (i in 1:length(dates)) {
    walmart.map(walmart, dates[i])
    title(dates[i])
  }
}, title = "Wal-Mart Expansion", htmlfile = "walmart.html",
  outdir = getwd(), autobrowse = FALSE)
## HTML file created at: walmart.html

  • アニメーションは出来上がった walmart.html ファイルをウェブブラウザで開くことで再生できる。

  • HTMLに加えて利用可能なフォーマットにはビデオファイルのための saveVideo() 関数がある。

  • 次の一連のコードで walmart.mp4 というMP4ファイルが作業ディレクトリに作成される。

    • 上記コードの saveHTML を saveVideo に変更

    • video.name でビデオファイルの名前を指定

    • ffmpeg という動画を変換、編集できるフリーソフトをインストールする必要がある。

## library("animation")
## saveVideo({
##   for (i in 1:length(dates)) {
##     walmart.map(walmart, dates[i])
##     title(dates[i])
##   }
## }, video.name = "walmart.mp4", 
##   outdir = getwd(), autobrowse = FALSE)

  • アニメーションを見ると、ウォルマートのフランチャイズについて、次のようなことがはっきりとわかる。

    • 南部発祥

    • 1970年代と1980年代に次第に地域全体に広がっていく様子

    • 1990年代中頃からの国内の残りの地域への急速な出店拡大

    • 出店拡大を見込んだ新しい流通センター設立のタイミングと場所

日本におけるスターバックスの拡大

  • コーヒーチェーン Starbucks の店舗拡大を時空間データを用いて分析する。

  • 店舗データは このサイト の 開店日順データ を用いる。

  • 1996年から2020年までのデータをコピーしてスプレッド・シートに貼り付け、不要なデータを削除して変数名を変更したファイル starbucks20.csv を読み込む。

starbucks20 <- read.csv("starbucks20.csv")
head(starbucks20)
##   id   opendate
## 1  1 1996-08-02
## 2  2 1996-09-05
## 3  3 1996-11-12
## 4  4 1997-02-28
## 5  5 1997-03-18
## 6  6 1997-07-18
##                                                                        address
## 1                                       東京都 中央区 銀座3-7-14 ＥＳＫビル 1F
## 2                            東京都 千代田区 神田駿河台2-5 村田ビルディング 1F
## 3                                   東京都 中央区 八重洲2-1 八重洲地下街 中3号
## 4                 東京都 町田市 原町田6-4-1 町田東急ツインズ ウエスト３階 W306
## 5                   東京都 千代田区 内幸町2-2-3 日比谷国際ビルヂング B1F(B132)
## 6 神奈川県 横浜市西区 みなとみらい2-3-9 クイーンズスクエア横浜［アット！］１Ｆ
##   prefecture
## 1     東京都
## 2     東京都
## 3     東京都
## 4     東京都
## 5     東京都
## 6   神奈川県
##                                                                        note
## 1                                                                 日本1号店
## 2                                                                 日本2号店
## 3                                                                 日本3号店
## 4    日本4号店 旧称:まちだ東急店､TOKYUまちだandYOU店 2015/06/26リニューアル
## 5                                             日本5号店 2018/05リニューアル
## 6 旧称:横浜［アット！］店(-2017/10/26) 神奈川県1号店 2005/09/30リニューアル

  • starbucks20.csv の address をコピーして谷謙二研究室（埼玉大学教育学部人文地理学）が提供している KTGIX.net の 地名・施設名からジオコーディング・地図化 を利用して緯度と経度を取得する。

  • 取得結果をコピーしてスプレッド・シートに貼り付け、不要なデータを削除して変数名を変更したファイル starbucks20longlat.csv を読み込む。

longlat20 <- read.csv("starbucks20longlat.csv")
head(longlat20)
##   id     long      lat
## 1  1 139.7667 35.67164
## 2  2 139.7608 35.70039
## 3  3 139.7696 35.67987
## 4  4 139.4465 35.54240
## 5  5 139.7544 35.67037
## 6  6 139.6338 35.45609
##                                                                        address
## 1                                       東京都 中央区 銀座3-7-14 ＥＳＫビル 1F
## 2                            東京都 千代田区 神田駿河台2-5 村田ビルディング 1F
## 3                                   東京都 中央区 八重洲2-1 八重洲地下街 中3号
## 4                 東京都 町田市 原町田6-4-1 町田東急ツインズ ウエスト３階 W306
## 5                   東京都 千代田区 内幸町2-2-3 日比谷国際ビルヂング B1F(B132)
## 6 神奈川県 横浜市西区 みなとみらい2-3-9 クイーンズスクエア横浜［アット！］１Ｆ
##                            address.match
## 1              東京都中央区銀座3丁目7-14
## 2         東京都千代田区神田駿河台2丁目5
## 3               東京都中央区八重洲2丁目1
## 4             東京都町田市原町田6丁目4-1
## 5           東京都千代田区内幸町2丁目2-3
## 6 神奈川県横浜市西区みなとみらい2丁目3-9
##                                                                         note
## 1                                                                       <NA>
## 2 他候補；東京都千代田区神田駿河台2丁目5-4；東京都千代田区神田駿河台2丁目5-5
## 3            他候補；東京都中央区八重洲2丁目1-10；東京都中央区八重洲2丁目1-4
## 4                                                                       <NA>
## 5                                                                       <NA>
## 6                                                                       <NA>

  • 住所によってはエラーになり緯度と経度が取得できない場合がある。

    • starbucks20longlat.csv では、エラーの場合に long と lat が -1 になるように修正している。
error <- longlat20$long == -1 # エラー (-1) のときTRUE、そうでなければFALSEとなるベクトル
head(error)
## [1] FALSE FALSE FALSE FALSE FALSE FALSE
head(longlat20[error, ]) # エラーとなった店舗を表示
##      id long lat
## 50   50   -1  -1
## 52   52   -1  -1
## 68   68   -1  -1
## 82   82   -1  -1
## 116 116   -1  -1
## 159 159   -1  -1
##                                                                   address
## 50    京都府 京都市下京区 四条通柳馬場西入ル 立売中之町106 ヤサカ四条ビル
## 52          大阪府泉佐野市泉州空港北一番地 関西国際空港旅客ターミナルビル
## 68          京都府 京都市中京区 三条通河原町東入ル中島町113 近江屋ビル 1F
## 82  京都府 京都市中京区 三条通烏丸西入ル 御倉町85-1 ＫＤＸ京都烏丸ビル 1F
## 116                                               大阪府枚方市市岡東町8-6
## 159                     大分県 大分市 玉沢楠本755-1 トキハわさだタウン 1F
##     address.match note
## 50           <NA> <NA>
## 52           <NA> <NA>
## 68           <NA> <NA>
## 82           <NA> <NA>
## 116          <NA> <NA>
## 159          <NA> <NA>

  • 店舗データ starbucks20 と緯度経度データ longlat20 をmerge 関数を用いて結合する。

    • 2つのデータフレームに含まれる同じ変数を用いて結合する場合は、その変数名（例えば id）をオプション by を用いて（例えば by = "id" のように）指定する

    • 変数名が異なる場合は、by.x と by.yという引数を用いて、それぞれのデータフレームで結合用の変数名を指定することもできる。

    • デフォルトでは結合されたデータフレームには引数 by.x で指定されたデータフレーム x の変数名が保存される．

  • ここでは、2つのデータフレームに共通の id を用いて次のようにデータを結合する。

starbucks <- merge(starbucks20, longlat20, by = "id")
head(starbucks)
##   id   opendate
## 1  1 1996-08-02
## 2  2 1996-09-05
## 3  3 1996-11-12
## 4  4 1997-02-28
## 5  5 1997-03-18
## 6  6 1997-07-18
##                                                                      address.x
## 1                                       東京都 中央区 銀座3-7-14 ＥＳＫビル 1F
## 2                            東京都 千代田区 神田駿河台2-5 村田ビルディング 1F
## 3                                   東京都 中央区 八重洲2-1 八重洲地下街 中3号
## 4                 東京都 町田市 原町田6-4-1 町田東急ツインズ ウエスト３階 W306
## 5                   東京都 千代田区 内幸町2-2-3 日比谷国際ビルヂング B1F(B132)
## 6 神奈川県 横浜市西区 みなとみらい2-3-9 クイーンズスクエア横浜［アット！］１Ｆ
##   prefecture
## 1     東京都
## 2     東京都
## 3     東京都
## 4     東京都
## 5     東京都
## 6   神奈川県
##                                                                      note.x
## 1                                                                 日本1号店
## 2                                                                 日本2号店
## 3                                                                 日本3号店
## 4    日本4号店 旧称:まちだ東急店､TOKYUまちだandYOU店 2015/06/26リニューアル
## 5                                             日本5号店 2018/05リニューアル
## 6 旧称:横浜［アット！］店(-2017/10/26) 神奈川県1号店 2005/09/30リニューアル
##       long      lat
## 1 139.7667 35.67164
## 2 139.7608 35.70039
## 3 139.7696 35.67987
## 4 139.4465 35.54240
## 5 139.7544 35.67037
## 6 139.6338 35.45609
##                                                                      address.y
## 1                                       東京都 中央区 銀座3-7-14 ＥＳＫビル 1F
## 2                            東京都 千代田区 神田駿河台2-5 村田ビルディング 1F
## 3                                   東京都 中央区 八重洲2-1 八重洲地下街 中3号
## 4                 東京都 町田市 原町田6-4-1 町田東急ツインズ ウエスト３階 W306
## 5                   東京都 千代田区 内幸町2-2-3 日比谷国際ビルヂング B1F(B132)
## 6 神奈川県 横浜市西区 みなとみらい2-3-9 クイーンズスクエア横浜［アット！］１Ｆ
##                            address.match
## 1              東京都中央区銀座3丁目7-14
## 2         東京都千代田区神田駿河台2丁目5
## 3               東京都中央区八重洲2丁目1
## 4             東京都町田市原町田6丁目4-1
## 5           東京都千代田区内幸町2丁目2-3
## 6 神奈川県横浜市西区みなとみらい2丁目3-9
##                                                                       note.y
## 1                                                                       <NA>
## 2 他候補；東京都千代田区神田駿河台2丁目5-4；東京都千代田区神田駿河台2丁目5-5
## 3            他候補；東京都中央区八重洲2丁目1-10；東京都中央区八重洲2丁目1-4
## 4                                                                       <NA>
## 5                                                                       <NA>
## 6                                                                       <NA>

  • データフレーム同士が同じ名前の変数を含んでいれば、結合されたデータフレームの当該変数には末尾に.xや.yが加えられ、それぞれの変数がどちらのデータフレームのものかが示される。

  • subset 関数を用いて、エラー (long と lat の値が -1) ではないデータを抽出し、必要な変数のみを残す。

starbucks <- subset(starbucks, lat >=0)
starbucks <- subset(starbucks, select = c(id, opendate, address.match, prefecture, long, lat))
head(starbucks)
##   id   opendate                          address.match prefecture     long
## 1  1 1996-08-02              東京都中央区銀座3丁目7-14     東京都 139.7667
## 2  2 1996-09-05         東京都千代田区神田駿河台2丁目5     東京都 139.7608
## 3  3 1996-11-12               東京都中央区八重洲2丁目1     東京都 139.7696
## 4  4 1997-02-28             東京都町田市原町田6丁目4-1     東京都 139.4465
## 5  5 1997-03-18           東京都千代田区内幸町2丁目2-3     東京都 139.7544
## 6  6 1997-07-18 神奈川県横浜市西区みなとみらい2丁目3-9   神奈川県 139.6338
##        lat
## 1 35.67164
## 2 35.70039
## 3 35.67987
## 4 35.54240
## 5 35.67037
## 6 35.45609

  • mapdata パッケージを用いて、日本地図に店舗を表示する。

    • ここ を参考に Starbucks ロゴの色を指定する。

    • 最小値と最大値のベクトルを返す range 関数を用いて、店舗がある緯度と経度の範囲±1度に地図を拡大する。

library(mapdata)

## starbucksロゴの色
starbucks$color <- rgb(red = 0.0118, green = 0.4, blue = 0.2078, alpha = 1/3)
starbucks$color[1]
## [1] "#03663555"
## 地図の範囲
x.range <- range(starbucks$long) + c(-1, 1)
y.range <- range(starbucks$lat) + c(-1, 1)

## 地図を表示
map("japan", interior = FALSE, xlim = x.range, ylim = y.range)
map("japan", boundary = FALSE, lty = 2, add = TRUE)
points(starbucks$long, starbucks$lat, col = starbucks$color, pch = 19, cex = 0.5)
map.axes()

  • 上で作成した starbuccks データを starbucks.csv として保存しておく。

  • 以下では、ウォルマートの例を参考に、starbucs.csv を用いて年ごとの地図やアニメーションを作成し、日本においてスターバックスがどのように店舗を拡大していったかを分析する。

  • データ starbucs.csv を読み込む。

    • 文字化けしないように、fileEncoding = "utf8" のようにエンコーディングを指定する。
starbucks <- read.csv("starbucks.csv", fileEncoding = "utf8")

  • 日付の変数 opendate を R の Date オブジェクトに変更する。

    • 開店日がないデータは削除する。
starbucks$opendate <- as.Date(starbucks$opendate)
starbucks <- subset(starbucks, is.na(opendate) == FALSE) # 開店日がないデータを削除
  • mapdata パッケージを読み込む。
library(mapdata)
  • 指定した日付以前に開店したすべての店舗を部分集合化し、それらの場所を地図上にプロットする関数 starbucks.map を作成する。
starbucks.map <- function(data, date) {
  starbucks <- subset(data, subset = (opendate <= date))
  map("japan", interior = FALSE)
  map("japan", boundary = FALSE, lty = 2, add = TRUE)
  points(starbucks$long, starbucks$lat, col = starbucks$color, pch = 19, cex = 0.5)
  map.axes()
}

  • Starbucksロゴの色を指定する。

  • 作成した starbucks.map を用いて、日本1号店がオープンした1996年から2000年までの各年と、それ以降の5年ごとに地図を作成する。

    • 店舗を示す点には、Starbucksロゴの色を指定する。
## starbucksロゴの色
starbucks$color <- rgb(red = 0.0118, green = 0.4, blue = 0.2078, alpha = 1/3)

starbucks.map(starbucks, as.Date("1996-12-31"))
title("1996")

starbucks.map(starbucks, as.Date("1997-12-31"))
title("1997")

starbucks.map(starbucks, as.Date("1998-12-31"))
title("1998")

starbucks.map(starbucks, as.Date("1999-12-31"))
title("1999")

starbucks.map(starbucks, as.Date("2000-12-31"))
title("2000")

starbucks.map(starbucks, as.Date("2005-12-31"))
title("2005")

starbucks.map(starbucks, as.Date("2010-12-31"))
title("2010")

starbucks.map(starbucks, as.Date("2015-12-31"))
title("2015")

starbucks.map(starbucks, as.Date("2020-12-31"))
title("2020")

  • 上記の地図から次のようなことがわかる。

    • 初店舗から数年は東京近郊で店舗を拡大した。

    • 1990年後半から大阪や神戸に出店し、2000年までに名古屋、福岡、仙台などの大都市に拡大した。

    • 2005年までに日本海側の都市や北海道および沖縄に店舗を拡大した。

    • 2006年以降も全国に店舗を拡大し続けている。

  • また、以下のようにアニメーションを作成することもできる。

n <- 25 # アニメーション化する地図の数

dates <- seq(from = min(starbucks$opendate), 
             to = max(starbucks$opendate), length.out = n)

# install.packages("animation") # インストールしていない場合
library("animation")
saveHTML({
  for (i in 1:length(dates)) {
    starbucks.map(starbucks, dates[i])
    title(dates[i])
  }
}, title = "Starbucks Japan Expansion", htmlfile = "starbucks.html",
outdir = getwd(), autobrowse = FALSE)
## HTML file created at: starbucks.html
  • アニメーションからも、上記の点がはっきりわかる。