Zenn
🚗

Rによる空間分析のススメ―交通事故ヒートマップを作成してみよう―

2022/12/23に公開
Docker
R
visualization
RStudio
Leaflet
tech

はじめに

Rは統計とグラフィックスに特化した言語です。またRパッケージと呼ばれる拡張パッケージによりWebアプリケーションの開発やLaTeX形式のレポートの作成など、さまざまなアウトプットが可能です。空間分析においてもさまざまなRパッケージが開発されており、かなり幅広い仕事ができます。日本でもNipponMap、jpmeshなど国内の地理情報を扱うパッケージも開発されています。

またコミュニティも含めてとてもおもしろい言語です。r-wakalangというslackでわからないことを質問できたり、書籍もたくさん出ています。筆者は交通事故を統計分析する研究をしており、取り組みはじめたときに研究室内で「統計ならR」という声があり、Rを学び始めました。はじめは本当にささいなきっかけで学び始めたRですが、最近では一番使っている言語になりました。

この記事では、交通事故のヒートマップ作成を題材に最近の自分のRによる空間分析を紹介します。最終的に作成するヒートマップのソースはGitHubで公開しています。また、愛知県の交通事故ヒートマップを表示するデモページも用意したので確認してみてください。この記事の内容がこれから空間分析をするという人にいくらか役立てればと思います。

https://github.com/NONONOexe/traffic-accident-heatmap

Rockerによる分析環境

Rによる開発環境

Rによる開発では、RStudioという統合開発環境を使うのが最も一般的です。RおよびRStudioはWindows、Mac、LinuxいずれのOSでも利用できます。しかし、OSによりRパッケージのインストール方法が異なったり、フォントの違いでグラフを正常に表示できなかったりとOSの違いによるトラブルが発生することがあります。空間分析向けの環境作りにおいては代表的なパッケージであるsfのインストール時のトラブルを聞きます。

Rockerを用いた開発では、こうしたトラブルを避け、共通の環境で開発することができます。RockerはRStudio Serverというサーバー版のRStudioを含むDockerコンテナです。またDockerコンテナであるため、自身の環境を汚すことがなく、Dockerさえインストールしていれば、RやRStudioのインストールも不要です。

RStudio Serverでは、ブラウザによるインターフェースをクライアントに提供しています。したがってRockerではコンテナを起動するだけで、手軽にブラウザでRStudioを起動することができます。RStudio Server自体は本来クライアントサーバーモデルで計算を高性能なマシンで実行するためのものですが、ローカルで起動することでデスクトップアプリのような感覚で使いつつ、先に述べた環境に関する恩恵を受けられるようになっています。

Rockerの起動

Rockerは次の手順で起動できます。さまざまな事前設定によりカスタマイズすることもできますが、ここではまずプレーンなRockerを起動します。

  1. Dockerのインストール

    Dockerは公式サイトから手に入れたインストーラーや、wingetやbrewなどのパッケージマネージャからインストールできます。好みの方法でインストールしてください。

  2. 次のコマンドをコンソールで実行

    docker run --rm -d -e PASSWORD=yourpassword -p 8787:8787 rocker/rstudio

  3. RStudio Serverへのログイン

    http://localhost:8787へアクセスして、ユーザ名/パスワードをrstudio/yourpasswordとしてログインします。アクセスするとこのようにRStudioの画面がブラウザ上で確認できます。
    r-studio
    RStudioの起動画面

さて、Rockerを起動できました。しかし残念ながら、このプレーンな状態で日本の空間分析をおこなうには必要なRパッケージがインストールされていない、日本語フォントがなく日本語をプロットできないなどさまざまな問題があります。

Rockerはubuntuのイメージがベースとなっています。そのため、コンソールから必要なライブラリなどを揃えられます。しかし、これは大変な手間です。またこれらを揃えたとしてもDockerコンテナの性質上、コンテナを停止したときにインストールしたライブラリや作成したプログラムはなくなってしまいます。

これらの問題はdocker-composeを利用して、DockerfileやComposeファイルへ設定を施すことで解決できます。またコマンド1つで自分用にカスタマイズした分析環境をセットアップできます。

Rockerの停止

docker-composeでのRockerの設定に移る前に今立ち上がっているRockerを停止しておきます。コンテナはdocker stop <コンテナID>で停止できます。コンテナIDはコンテナを起動したときにコンソールに表示されたハッシュ値です。もしすでにコンソールになければdocker psで次のように確認できます。

「docker ps」による確認例
CONTAINER ID   IMAGE            COMMAND   CREATED             STATUS             PORTS                    NAMES
36f641b105d8   rocker/rstudio   "/init"   About an hour ago   Up About an hour   0.0.0.0:8787->8787/tcp   interesting_murdock

コンテナIDを指定するとき、文字列すべてを入力する必要はありません。コンテナを一意に特定できる桁までを入力すれば識別されます。上記の例では起動しているコンテナが1つしかないためdocker stop 3でも十分ということになります。

空間分析のためのDockerイメージの作成と起動

本題から逸れるため、あまり詳しくは説明を端折りますが、Rockerを日本における空間分析向けにカスタマイズしたDockerfileとComposeファイルを以下に示します。これらを任意のディレクトリに配置してください。以降では、このディレクトリにプログラムなどを保管します。

Dockerfile
FROM rocker/geospatial:4.2.2

# 日本語ロケールの設定
ENV LANG ja_JP.UTF-8
ENV LC_ALL ja_JP.UTF-8
RUN sed -i '$d' /etc/locale.gen \
  && echo "ja_JP.UTF-8 UTF-8" >> /etc/locale.gen \
    && locale-gen ja_JP.UTF-8 \
    && /usr/sbin/update-locale LANG=ja_JP.UTF-8 LANGUAGE="ja_JP:ja"
RUN /bin/bash -c "source /etc/default/locale"
RUN ln -sf  /usr/share/zoneinfo/Asia/Tokyo /etc/localtime

# 日本語フォントのインストール
RUN apt-get update && apt-get install -y \
  fonts-ipaexfont \
  fonts-noto-cjk

# ワーキングディレクトリの設定
RUN echo "setwd(\"/home/rstudio/workspace/\")" > /home/rstudio/.Rprofile

# rocker/geospatialに含まれないRパッケージのインストール
# CRANからのインストール
RUN install2.r -d TRUE -e -n -1 \
  celestial
# GitHubからのインストール
RUN installGithub.r \
  uribo/jpmesh \
  uribo/jpndistrict

# インストール時にダウンロードした一時ファイルの削除
RUN rm -rf /tmp/downloaded_packages/ /tmp/*.rds
docker-compose.yml
version: '3'

services:
  rstudio:
    # ビルドするDockerfileが含まれるディレクトリ
    build: .
    # 作成するイメージ名(任意: ユーザ名/プロジェクト名)
    image: nononoexe/rocker-geospatial-jp
    container_name: rstudio
    environment:
      TZ: Asia/Tokyo
      PASSWORD: yourpassword
    ports:
      - "8787:8787"
    # データの永続化
    volumes:
      - ".:/home/rstudio/workspace"

配置したら、下記のdocker-composeコマンドからRockerを起動します。最初はイメージのビルドを含むため、少し時間がかかります。なお、Dockerfileに従いイメージをビルドするために--buildオプションを付与していますが、2回目以降の起動ではすでにビルドされたイメージがあるため、これは不要になります。

カスタマイズしたRockerイメージをビルドして起動するためのコマンド
docker-compose up --build -d

起動すると、ロケールを日本へ変更したためメッセージが日本語で表示されます。あとはこの環境でプログラムを書いていくことができます。

日本語で表示されるRのコンソールメッセージ
日本語で表示されるRのコンソールメッセージ

ここで保存したソースコードやデータはDockerfileとComposeファイルを配置したディレクトリに保存されるようになっています。これはRのワーキングディレクトリにこのディレクトリをマウントするようComposeファイルに設定しているためです。

空間データを扱えるRockerイメージ

Rockerには以下の表に示すようにいくつかのイメージの種類があります。ここでは、空間分析に向けたrocker/geospatialをベースイメージにしたイメージを作成しました。このイメージを利用することで空間分析に必要なRパッケージがある程度インストールされた状態ではじめられます。

イメージ名 説明 ベースイメージ ベースイメージに追加されているもの
rocker/r-ver Rのみを含むイメージ ubuntu R
rocker/rstudio RStudioを利用できるイメージ rocker/r-ver RStudio Server
rocker/tidyverse 分析用のパッケージ群tidyverseを含むイメージ rocker/rstudio tidyverse、devtools
rocker/verse ドキュメント作成のためのパッケージを含むイメージ rocker/tidyverse \TeX、bookdownなど
rocker/geospatial 空間データを扱うパッケージを含むイメージ rocker/verse sf、terraなど

rocker/geospatialは表のとおりtidyverseや\TeXを含むイメージをベースとしています。そのため、空間データを扱うだけでなく、広くデータ分析や、レポート作成、パッケージ開発に便利なRPackgeも含まれています。

Rockerをパスワードなしで利用する方法

Rockerを起動するとユーザ名とパスワードが求められます。自分しか使用しないのに毎回入力するのは不便に思う方もいるかもしれません。そうした場合、Rockerにあるパスワードを省略するオプションが使えます。

パスワードを省略したい場合はコンテナ内の環境変数を設定するenvironmentキーに次のようにDISABLE_AUTH: trueを設定します。こうして、http://localhost:8787へアクセスするとログイン画面をとばして、RStudioの画面へ遷移することができます。

docker-compose.yml
    environment:
      DISABLE_AUTH: true

ただし、当然ながらポートが解放されている場合、同じネットワークからコンテナを経由してマウントしたファイルへまでアクセスできてしまうのでこのオプションは注意して利用ください。

交通事故のヒートマップの作成

作成するヒートマップと手順

ここでは、次の画像のようなヒートマップを作成します。この図は明るいところの事故件数が多く、暗いところの事故件数が少ないことを示しています。またこの画像は愛知県を対象にしたときの例を示していますが、好みの都道府県に変更できますので、ぜひ変更してみてください。

作成するヒートマップの例
作成するヒートマップの例

ヒートマップはよく見ると細かい格子状であり、そのマス(メッシュ)を件数に応じた色で塗り分けられています。これは次の手順で作成できます。

  1. 交通事故地点を含むデータのダウンロード
  2. Rへのデータインポート
  3. メッシュごとに交通事故件数を集計
  4. 交通事故件数によるヒートマップを作成

以降では、各手順を実行しながらヒートマップの作成について説明します。

データのダウンロード

まずは次のソースコードを実行して、交通事故の位置情報を含むデータをダウンロードします。このデータには警察庁が公開している2019年から2021年の道路交通事故のオープンデータを使います。そのうち、事故地点の緯度、経度が含まれる本票データをダウンロードします。

本票データのダウンロード
library(fs)
if (!dir_exists("data")) dir_create("data")
download.file("https://www.npa.go.jp/publications/statistics/koutsuu/opendata/2019/honhyo_2019.csv", "data/main-data-2019.csv")
download.file("https://www.npa.go.jp/publications/statistics/koutsuu/opendata/2020/honhyo_2020.csv", "data/main-data-2020.csv")
download.file("https://www.npa.go.jp/publications/statistics/koutsuu/opendata/2021/honhyo_2021.csv", "data/main-data-2021.csv")

download.file()で指定したURLのデータをダウンロードすることができます。ここではdataというディレクトリを作成し、そこにダウンロードしています。

Rへのデータインポート

次にダウンロードしたデータをRへインポートします。次のソースコードを実行します。

データインポート
library(tidyverse)
accidents <- bind_rows(
  read_csv("data/main-data-2019.csv", locale = locale(encoding = "Shift_JIS"), col_types = cols(.default = col_character())),
  read_csv("data/main-data-2020.csv", locale = locale(encoding = "Shift_JIS"), col_types = cols(.default = col_character())),
  read_csv("data/main-data-2021.csv", locale = locale(encoding = "Shift_JIS"), col_types = cols(.default = col_character()))
)

ここでは、read_csv()によりデータをインポートし、bind_rows()により各年のデータを結合しています。また、これらのデータは文字コードがShift_JISであるため、localeに読み込み時の文字コードを指定しています。

データ分析に欠かせないRパッケージ: tidyverse

ここで使用しているtidyverse(読み: タイディーバース)はデータサイエンスのためのRパッケージをまとめたRパッケージです。近年のRによる開発では、Rの標準的な機能(baseR)に加えて、このRパッケージが非常によく利用されています。

tidyverseには、さまざまな機能が含まれますが、この記事では以下のようなRパッケージの機能を使用しています。

パッケージ名 役割 使用する代表的な関数
readr データインポート read_csv(): CSVファイルの読み込み
dplyr データ操作 select(): 項目の選択、filter(): データの抽出
stringr 文字列操作 str_sub(): 部分文字列の切り出し
ggplot2 グラフのプロット geom_sf(): 空間オブジェクトの描画

インポートしたデータを確認します。View()を使うとデータを確認できます。

データの確認
View(accidents)

確認すると「資料区分」、「都道府県コード」などさまざまな項目が確認できます。たくさんの項目がありますが、終端の方に「地点　緯度（北緯）」、「地点　緯度（東経）」という項目が確認できます。この項目によりこのデータを空間データに変換します。

データの確認
データaccidentsの確認

空間データへの変換

緯度および経度のDEG形式への変換

このデータでは、空間データへの変換に緯度、経度の変換が必要です。「地点　緯度（北緯）」、「地点　緯度（東経）」という項目は小数点を省くDMS形式で登録されています。DMS形式とは、度(Degree)、分(Minute)、秒(Second)で表現された緯度、経度の表記です。一方、GISではDEG形式と呼ばれる度のみの形式を利用します。そのため、DMS形式からDEG形式への変換が必要です。

DMS形式の緯度および経度の視覚的確認

DMS形式の緯度および経度はデータをプロットすると確認できます。次のソースコードは例として、愛知県(都道府県コード: 54)のデータのうち5,000件をプロットするものです。

愛知県の事故地点のプロット
accidents |>
  filter(`都道府県コード` == "54") |>
  head(5000) |>
  ggplot() +
  geom_point(aes(x = as.integer(`地点　経度（東経）`),
                 y = as.integer(`地点　緯度（北緯）`))) +
  labs(x = "地点　経度（東経）",
       y = "地点　緯度（北緯）") +
  theme_minimal()

緯度経度の確認
愛知県が分断されたかのような結果表示

DMS形式では「60分」以降の値がないため、プロットが途切れています。この例では、北緯34.60～35.00、東経136.60～137.00の区間が不自然に途切れているのが確認できます。

DMS形式からDEG形式へ変換はcelestialというRパッケージにより実現できます。次のソースコードにより各データの緯度および経度をDEG形式に変換します。

DMS形式からDEG形式への変換
library(celestial)
accidents_converted_deg <- accidents |>
  # 項目名の変更
  rename(
    lat_dms = `地点　緯度（北緯）`,
    lon_dms = `地点　経度（東経）`,
  ) |>
  # 度・分・秒を分割
  mutate(
    lat_d = str_sub(lat_dms, 1, 2),
    lat_m = str_sub(lat_dms, 3, 4),
    lat_s = str_c(str_sub(lat_dms, 5, 6), ".", str_sub(lat_dms, 7)),
    lon_d = str_sub(lon_dms, 1, 3),
    lon_m = str_sub(lon_dms, 4, 5),
    lon_s = str_c(str_sub(lon_dms, 6, 7), ".", str_sub(lon_dms, 8))
  ) |> 
  # 変換不可能なデータを除外
  filter(
    0 <= lat_m & lat_m < 60,
    0 <= lat_s & lat_s < 60,
    0 <= lon_m & lon_m < 60,
    0 <= lon_s & lon_s < 60
  ) |>
  # DEG形式に変換
  mutate(
    lat = dms2deg(lat_d, lat_m, lat_s),
    lon = dms2deg(lon_d, lon_m, lon_s)
  )

dms2degが度、分、秒を受け取り、DEG形式の値を返してくれる関数です。ここでは、単純に緯度および経度の文字列を度、分、秒に分割して形式を変換しています。また、分、秒は0から60までという値の範囲があるため、変換不可能なものを除いています。

sfオブジェクトへの変換

緯度および経度をDEG形式にしたので、これらを座標に指定し、空間データへ変換します。次のソースコードを実行してデータを空間オブジェクトに変換します。

sfオブジェクトへの変換
library(sf)
accidents_sf <- accidents_converted_deg |>
  st_as_sf(coords = c("lon", "lat"), crs = 4326)

sfはRで空間オブジェクトを扱うためのRパッケージです。sfが提供する空間オブジェクトはSimple Featureと呼ばれる一般的な空間データの規格に基づくオブジェクトです。そのためPostGISなど他の空間データを扱うツールとのやりとりも可能です。

変換ができているかどうかはプロットから確認できます。データ件数が非常に多いため、ここでは例として、愛知県(都道府県コード: 54)のデータのうち5,000件のプロットを示します。

交通事故地点のプロット
accidents_sf |>
  filter(`都道府県コード` == "54") |>
  head(5000) |>
  ggplot() +
  geom_sf() +
  theme_minimal()

交通事故地点のプロット結果
愛知県の交通事故地点のプロット

交通事故件数の空間集計

メッシュの作成

ヒートマップは地点情報を一定の大きさの空間内に交通事故が何件あるかを集計して作成します。このような空間方向への集計を空間集計といいます。空間集計するにはまず空間を同じ大きさに分割した集計単位が必要です。

ここでは、基準地域メッシュという約1km四方の区域を集計単位にします。地域メッシュはjpmeshというRPakcageで提供されています。administration_meshは指定された都道府県コードの地域メッシュを返してくれます。すべての交通事故を扱うのは処理が大きくという理由もあり、以降の例では、都道府県コード23の愛知県を対象にします。標準地域メッシュは次のソースコードを実行して作れます。

基準地域メッシュの作成
library(jpmesh)
pref_code <- 23
grid_squares <- administration_mesh(pref_code, 1)

標準地域メッシュも作成できているかをプロットして確認します。

基準地域メッシュのプロット
grid_squares |>
  ggplot() +
  geom_sf() +
  theme_minimal()

基準地域メッシュのプロット結果
基準地域メッシュのプロット結果

交通事故件数の集計

交通事故とメッシュの両方の空間データが作成できました。これらを空間的に結合して、各メッシュに含まれる交通事故地点を数え上げます。次のソースコードを実行します。

メッシュごとの交通事故件数の算出
accident_counts_sf <-
  # 交通事故と基準地域メッシュを結合
  st_join(accidents_sf, grid_squares) |>
  # 一時的に非空間データに変換
  st_drop_geometry() |>
  # メッシュ単位に集計
  group_by(meshcode) |>
  summarise(count = n()) |>
  # 空間情報を復元
  inner_join(grid_squares, by = "meshcode") |>
  st_as_sf()

ここでは、st_join()で交通事故と基準地域メッシュの二つの空間データを結合しています。この関数はデータベースでいう、INNER JOINのため、片方の空間データが領域を持たないものは除外されます。これにしたがい、交通事故側は全国データでしたが、基準地域メッシュにより指定した都道府県のもののみとなります。

次に、集計前にst_drop_geometry()により一時的に非空間データにしています。これは空間データをそのまま集計するとデータサイズが大きく時間がかかるためです。データが持つ空間情報は基準地域メッシュであるため、集計後に再度grid_squaresと結合して空間情報を復元しています。

R Consoleにaccident_count_sfと入力して、集計された値を確認すると次のようにcountという項目で交通事故件数が集計されていることが確認できます。

集計結果
#> Simple feature collection with 2993 features and 2 fields
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 136.6625 ymin: 34.575 xmax: 137.7625 ymax: 35.425
#> Geodetic CRS:  WGS 84
#> # A tibble: 2,993 × 3
#>    meshcode count                                                       geometry
#>    <chr>    <int>                                                  <POLYGON [°]>
#>  1 51376091     1 ((137.0125 34.575, 137.025 34.575, 137.025 34.58333, 137.0125…
#>  2 51376092     1 ((137.025 34.575, 137.0375 34.575, 137.0375 34.58333, 137.025…
#>  3 51377002     2 ((137.025 34.58333, 137.0375 34.58333, 137.0375 34.59167, 137…
#>  4 51377003     1 ((137.0375 34.58333, 137.05 34.58333, 137.05 34.59167, 137.03…
#>  5 51377004     1 ((137.05 34.58333, 137.0625 34.58333, 137.0625 34.59167, 137.…

交通事故件数の可視化

メッシュごとに集計できたので、これを事故件数ごとに塗分けて、ヒートマップを作成します。メッシュをプロットしたときとほとんど同じようにプロットすることができます。次のソースコードを実行して確認します。

基準地域メッシュの交通事故件数による塗分け
accident_counts_sf |>
  ggplot() +
  geom_sf(aes(fill = count)) +
  scale_fill_continuous("事故件数(件)") +
  theme_minimal()

交通事故件数により塗分けられた基準地域メッシュ
交通事故件数により塗分けられた基準地域メッシュ

交通事故件数が多い部分が明るくなるようなプロットができました。プロットから最も多いメッシュではこの3年間で500件以上の事故が発生していることが確認できます。ただ、このプロットではどのメッシュが地理的にどの位置であるかは分かりにくいです。そこで、次にこのヒートマップを地図に重ねてみます。

Web地図への可視化

インタラクティブなヒートマップの作成

tmapというR Pacakgeを使うと先ほどのような通常のプロットだけでなく、Web地図によるインタラクティブな可視化が簡単に実現できます。次のソースコードを実行により、先ほど空間集計した交通事故件数のデータをWeb地図へ出力します。

tmapによるヒートマップ作成
library(tmap)
tmap_mode("view")
tm_shape(accident_counts_sf) +
  tm_polygons(
    col        = "count",
    alpha      = 0.8,
    title      = "事故件数 (件)",
    border.col = "gray",
    id         = "meshcode",
    popup.vars = c("事故件数" = "count"))

Web地図上へプロットされたヒートマップ
Web地図上へプロットされたヒートマップ

tmapでは、空間オブジェクトの情報をそのオブジェクトをクリックしたときのポップアップにより表示できます。また、このポップアップの表示項目はpopup.varsに指定できます。ここでは、事故件数のみが表示されるように設定しています。

事故件数のポップアップ表示
事故件数のポップアップ表示

tmapのplotモード

tmapにはggplot2を用いた通常のプロットのモードと、leafletを用いたインタラクティブな可視化のモードがあります。先ほどのコードのtmap_mode()がこれを制御しています。先ほどのコードもtmap_mode("plot")とすると、以下の画像のように通常のプロットを出力できます。

によるプロット
tmapによる通常のプロット

複数のレイヤーによる可視化

行政区域レイヤーの追加

tmapのLeafletをインターフェースとした可視化はさらに複数のレイヤーを重ねて表示することができます。先ほどのヒートマップに行政区域を追加します。行政区域の境界データはjpndistrictというRパッケージにより提供されています。

行政区域レイヤーの追加
# 交通事故件数のレイヤーを作成
accident_counts_tm <-
  tm_shape(
    accident_counts_sf,
    name = "交通事故件数") +
  tm_polygons(
    col        = "count",
    alpha      = 0.8,
    title      = "事故件数 (件)",
    border.col = "gray",
    id         = "meshcode",
    popup.vars = c("事故件数" = "count"))

# 行政区域のレイヤーを作成
library(jpndistrict)
cities_tm <-
  tm_shape(
    jpn_cities(pref_code),
    name = "行政区域") +
  tm_polygons(
    col   = "white",
    alpha = 0.5,
    id    = "city")

# レイヤーを重ねて表示
cities_tm + accident_counts_tm

行政区域を追加した地図
行政区域を追加した地図

ここでは、先ほどヒートマップを出力したコードを変数に代入しています。この変数がヒートマップのレイヤーを表します。tmapでは、この各レイヤーを定義したのち、+演算子によりそれぞれのレイヤーを重ねることができます。このとき後に追加した方が上に重なります。

重ねたレイヤーは左上のアイコンから表示または非表示に切り替えることができます。またこのアイコンからベースマップも変更できます。これにより可視化したい場所について好きな組み合わせのレイヤーで見られます。

レイヤーの切り替え

死亡事故地点レイヤーの追加

さらに注目する地域の事故状況を具体的に知るために、事故地点のレイヤーを加えてみます。ただし、すべての事故地点の追加は数が多くかえって見づらくなってしまうため、ここでは危険度の高い死亡事故地点に限定したレイヤーを加えます。以下のようにデータの抽出は必要ですが、レイヤーの追加は行政区域レイヤーと同様に実現できます。

死亡事故地点レイヤーの追加
# 死亡事故地点のレイヤーを作成
fatal_accidents_tm <-
  accidents_sf |>
  st_filter(grid_squares) |>
  filter(0 < as.integer(死者数)) |>
  select(!c(
    lat_dms, lon_dms,
    lat_d, lat_m, lat_s,
    lon_d, lon_m, lon_s)) |>
  tm_shape(name = "交通死亡事故地点") +
  tm_dots(
    col        = "royalblue",
    border.col = "transparent",
    id         = "fatal_accidents")

# レイヤーを重ねて表示
cities_tm + accident_counts_tm + fatal_accidents_tm

死亡事故地点を追加した地図
死亡事故地点を追加した地図

事故情報のポップアップ表示
事故情報のポップアップ表示

このようにtmapを使えば比較的短いソースコードで空間データのレイヤーを追加できます。さらに独自のレイヤーを作って拡張してみてもらえたらと思います。

おわりに

この記事では交通事故ヒートマップの作成を例にRでの空間分析について紹介しました。空間分析のためのRパッケージはここで紹介したもの以外にもたくさんあります。調べてみると自分のやりたいことに応じたRパッケージが開発されていることがあるので、これから空間分析を始めようという人はぜひ一つの手段として調べてみてもらえたらと思います。

Rによる空間分析についてはいくらか書籍も出版されています。特に「Geocomputation with R」は日本語化もされ、公開されているため勉強したい人へオススメです。

📚オススメのRによる空間分析に関する書籍

Rは勉強にあたりさまざまなRパッケージやそれらが提供する関数があり、はじめは覚えることが多く大変かもしれません。ただ、それぞれの命名や設計はわかりやすく、使い慣れると実現したい可視化や分析が少ない手数でできるようになり、便利だと思っています。

ぜひ分からないことや、もっとよい方法の提案など、この記事についてのコメントまたは質問があれば気軽に書いてください。

PR: Road Traffic Accident Database

筆者は研究でこの記事のように交通事故のデータを扱うためのデータベースをRで開発しています。現在このデータベースはクローズドなデータが必須となってしまっていますが、オープンデータのみで広いユーザが扱えるように改修しています。またドキュメントも執筆しています。こちらもRでの空間データプログラミングの参考にしていただいたり、使ってみたいという要望やコメントなどいただけたら嬉しいです。

https://github.com/NONONOexe/rtadb

Discussion