エシェロン解析で交通事故発生地点を分類する

はじめに

タクシーアプリ『GO』でデータアナリストとしてデータ分析をしている赤池です。

タクシーアプリのデータ分析をしていると、タクシー需要の局所的な広がりや、ユーザーの乗降地点の地域属性を知りたくなることがあります。 そういった時に、単純に地図上に点をプロットするだけでも概況はつかめますが、密度が高すぎると集中の度合いが判別しづらくなるという問題があります。 一方で、メッシュ単位で集計して可視化したとしても、一部の突出した大きな値に全体が引っ張られ、中小規模の重要な「山」が埋もれてしまいがちです。

そこで今回は、空間に分布する1変数を用いて「エシェロン解析」によって客観的に空間構造を抽出する方法をご紹介します。

1. エシェロン解析って何？

1-1. エシェロン解析とは

エシェロン解析とは、1変数の空間データを「標高」に見立てて地形の起伏を階層的に捉え、領域分割やホットスポット抽出を行う手法です。 エシェロンとは階層のことであり、値の大小や隣接関係をもとに直線や平面のデータを立体的に捉えることで、分析の基本単位である「エシェロン」を抽出します。それらのエシェロンは、複数のエシェロンで構成される「圏」の分類や、集積性の高い重要領域の検出に活用します。

出典：小田・石岡・正木・栗原（2012） p.21

1-2. エシェロンの分類と抽出方法

エシェロンは、データの極大値を含む「ピーク」とそれ以外の「ファウンデーション」に分類されます。 ピーク → ファウンデーションの順で、隣接する格子との大小比較を繰り返して抽出します。

エシェロンの抽出手順

抽出順序	抽出対象	処理
1	ピーク	極大値に隣接する格子との大小関係の比較を繰り返し、ピークのエシェロンを抽出する
2	ファウンデーション	ピークのエシェロンを除いた格子の中で、隣接する格子との大小比較を繰り返し、ファウンデーションのエシェロンを抽出する

エシェロン解析では主に格子データが用いられます。これは、空間において範囲が区切られた離散的なデータをさし、メッシュのような規則的な格子だけでなく行政区域のような非規則的なものも含みます。

1-3. エシェロン解析で分かること

エシェロン解析を行うことで、以下のようなことが客観的に分かります。

• ホットスポットの特定: 単に数値が高いだけでなく、統計的に意味のある「集積地（ホットスポット）」がどこにあるかを特定できます。
• 階層構造の把握: 「大きな山（メインの流行や集積）」の中に「小さな山（局所的な流行）」が含まれているのか、それとも独立した山なのか、という親子関係や独立性（エシェロン木と呼ばれる樹形図で表現されます）を可視化できます。
• ノイズの除去: たまたま数値が高いだけの「ノイズ」と、意味のある「ピーク」を区別しやすくなります。

他の分析手法との違いを、具体的なデータを使って見てみましょう。

2. 交通事故発生地点を分類する

2-1. 分析対象の整理

今回は警察庁が公開している「交通事故統計情報のオープンデータ」^1を用いて、交通事故の発生件数の広がりを分類してみます。

このデータは発生した事故ごとに発生時点や緯度経度、事故の状況などが記録されている空間点パターンデータです。今回は2024年の「人対車両」「車両相互」の交通事故に絞ります。 分析対象地域は、現在私が住んでいる宮城県仙台市とします。

2-2. 可視化による集積の確認

まず、分析対象とした仙台市は以下の地図中で黄色に塗りつぶした場所です。 東京から新幹線で約1時間半。首都圏からのアクセスもよく、食べ物もおいしいいい街です。

【事故発生地点をポイントで可視化】

交通事故の発生地点をポイントで可視化すると以下のようになります。 仙台市の西側は自然が豊かなので、街が栄えている東側に事故が集中しています。 ここで、東側のどこで事故が発生しやすいかを分類しようとすると、ポイントの集積を評価しなければいけませんが、基準が無いため分析者によって異なる結果となりそうです。

【事故発生地点をクラスターで可視化】

次に、ポイントをクラスターで可視化してみましょう。 ポイントよりも、事故が集中している場所が分かるようになりました。しかし、同じサイズのクラスターも多く、そのクラスターがどの範囲まで広がっているのかは判然としません。 また、クラスター数などのパラメータの設定によっては、異なる可視化結果となることもあり、結果の解釈が安定しません。

【事故発生地点をヒートマップで可視化】

ヒートマップによる可視化もポイントと比べれば事故が集まる場所が分かりますが、どこにピークがあって、それがどの程度広がっているのかはよく分かりません。

【事故発生地点をHexでビニングして可視化】

最後に六角形のHexでデータをビニングし、事故件数に応じて色付けしてみます。こちらは、局所的に事故が多いピークになっていそうなHexはわかるものの、その広がりや範囲はいまいち掴めません。

【事故発生地点をHexでビニングしてフィルターした後に可視化】

事故件数が3件以上のHexに絞り込んでみます。これにより、おおよその事故の集中エリアが抽出できるようになりましたが、中心地の事故件数が多すぎるために「周囲と比べると事故件数が多いものの、分析対象における事故件数ランキングをつけた場合に中位程度」の集積を見逃してしまう可能性があります。

2-3. エシェロン解析を試す 〜実装〜

さて、ここでエシェロン解析を適用してみます。 エシェロン解析を行うためのパッケージがPythonにないため、Rのechelonを使用します。 以下、私の「慣れ」の都合でエシェロンの分類のみRで行い、その前後の処理はPythonで行いました。

まず、交通事故統計情報のオープンデータをダウンロードし、以下の形に加工しておきます。（コードは省略）

• 変数 df_hex （上位5件）

hex	事故件数
"892e2b45b1bffff"	1
"892e2b47113ffff"	4
"892e2b445c3ffff"	2
"892e2b40877ffff"	1
"892e2b45413ffff"	2

エシェロン解析の実行時に格子間の隣接行列が必要となるため、これも出力しておきます。 ちなみに、ポリゴンデータやメッシュデータであればR側でspdepパッケージのpoly2nbやcell2nbで隣接行列を取得できますが、今回はHex粒度で集計するため、以下の処理で出力しました。 Hexをポリゴンに変換してからpoly2nbで隣接行列を取得することもできますが、Pythonのh3パッケージのare_neighbor_cells関数でHex同士が隣接しているかを判定しています。

import polars as pl
import h3


# Hexの隣接行列を取得
def are_neighbor_cells(hexes: list) -> list:
    h1, h2 = hexes
    return h3.are_neighbor_cells(h1, h2)


df_nb = (
    df_hex.select("hex")
    .join(df_hex.select("hex"), how="cross")
    .with_columns(
        is_neighbor = (
            pl.concat_list("hex", "hex_right")
            .map_elements(are_neighbor_cells, return_dtype=pl.Boolean)
        )
    )
    .with_columns(neighbor_flg = pl.sql_expr("if(is_neighbor, 1, 0)"))
    .pivot(
        index          = "hex",
        on             = "hex_right",
        values         = "neighbor_flg",
        maintain_order = True
    )
)

次に、R側でエシェロン解析を実行します。実行自体は関数一つで完了できます。

library(tidyverse)
library(echelon)


# (注) 以下、Python側で作成したデータを以下のデータフレームに読み込んでいるとします。
# - `df`   : Hexごとの事故件数データ
# - `df_nb`: 隣接行列


# 関数用にマトリクスに変換
nb = df_nb |>
    select(-hex) |>
    as.matrix()

# エシェロン解析の実行
result <- echelon(x = df$事故件数, nb = nb)

上記のresultからエシェロンの番号とそれに紐づくHexを取り出せます。後続の圏の分類を今回Pythonで行う都合上、表形式に変換した上でcsvファイルに保存しておきます。

# エシェロンのIDとそれに紐づくHexのインデックス値と事故件数（文字列）のリストを
# 圏の統合用に表形式に加工して保存する。
eche_idx = c()
eche_val = c()

i <- 1
for (eche in result$Echelons) {
    for (idx_val in eche) {
        eche_idx <- eche_idx |> append(i)
        eche_val <- eche_val |> append(idx_val)
    }
    i <- i + 1
}

result.echelon <- tibble(
    echelon   = eche_idx,
    index_val = eche_val
)

この結果をもとに圏を抽出したいのですがRのechelonパッケージには関数がなかったため、（行ったり来たりで恐縮ですが）Python側で圏の抽出を行います。 ピークのエシェロンをそれぞれの圏の起点とし、隣接するエシェロンを併合していくことで圏を抽出します。便宜上、複数の圏に接するエシェロンはランダムに割り当てました。

圏の抽出用データ作成

import polars as pl
import polars_st as pst
import polars.selectors as cs
import h3
from shapely import to_wkt
from shapely.geometry import Polygon


# (注) 以下、Python、Rで作成したデータを以下のデータフレームに読み込んでいることとします。
# - エシェロン解析用に準備したデータ
#     - `df_hex`: Hexごとの事故件数データ
#     - `df_nb` : Hex粒度の隣接行列
# - エシェロン解析の結果データ
#     - `df_echelon_table`    : エシェロン別の情報データ（R側で`result$Table`で得られるデータ。トップ、親子関係など。）
#     - `df_echelon_index_val`: エシェロンのIDとそれに紐づくHexのインデックス値と事故件数（文字列）のリスト

# エシェロン別の情報
df_echelon_table = (
    df_echelon_table
    .with_row_index(offset=1)
    .rename({"index": "echelon"})
)

# Hex（のインデックス）ごとのエシェロンの判定結果
df_echelon = (
    df_echelon_index_val
    .with_columns(pl.col("index_val").str.replace("\)", "").str.split("("))
    .with_columns((pl.col("index_val").list.get(0).cast(pl.Int64) - 1).alias("index"))
    .select("index", "echelon")
    .sort("index")
)

# 判定用データ作成
df_prep = (
    df_hex
    .join(df_echelon,       on="index",   how="left")
    .join(df_echelon_table, on="echelon", how="left")
    .rename({"Order": "eche_order"})
    .select("index", "hex", "事故件数", "echelon", "eche_order")
)

# Hex別の隣接Hexを抽出
df_hex_and_nb_hex = (
    df_nb
    .unpivot(index="hex", variable_name="h2", value_name="neighbor_flg")
    .filter(pl.col("neighbor_flg").eq(1))  # 隣接しているHexに絞る
    .join(df_prep, on="hex", how="left")
    .join(
        df_prep.select(cs.all().name.prefix("nb_")),
        left_on  = "h2",
        right_on = "nb_hex",
        how      = "left"
    )
    .rename({"h2": "nb_hex"})
    .select(
        "index", "hex", "echelon", "eche_order", "事故件数",
        "nb_index", "nb_hex", "nb_echelon", "nb_eche_order"
    )
    .sort(by=["index", "nb_index"])
)

圏の抽出処理

# Topのエシェロン
df_eche_top = (
    df_prep
    .filter(pl.col("eche_order").eq(1))
    .rename({"echelon": "top_echelon"})
    .group_by("index", "hex", "事故件数", "top_echelon")
    .agg(
        pl.col("top_echelon").alias("top_echelons"),
        pl.n_unique("top_echelon").alias("n_top_echelon")
    )
    .select("index", "hex", "top_echelons", "n_top_echelon", "事故件数")
)

# 最終成果物の初期化
df_result = df_eche_top.clone()
n_hex     = df_prep["hex"].n_unique()

while len(df_result) < n_hex:
    # Topのエシェロンに接するHexを抽出
    _df_hex_and_nb_eche = (
        df_hex_and_nb_hex
        .join(df_result, on="hex", how="anti")  # 判定済みのHexを除外
        .join(df_result, left_on="nb_hex", right_on="hex", how="inner")  # 判定済みのHexに隣接するHexに絞り込み
        .explode("top_echelons")
        .select("index", "hex", "top_echelons", "事故件数")
        .unique()
        .group_by("index", "hex", "事故件数")
        .agg(
            pl.col("top_echelons"),
            pl.n_unique("top_echelons").alias("n_top_echelon"),
        )
        .select("index", "hex", "top_echelons", "n_top_echelon", "事故件数")
    )

    # 判定済みHexの追加
    df_result = pl.concat([df_result, _df_hex_and_nb_eche]).sort("index")

df_result = (
    df_result
    .with_columns(
        # リストの並び番が上記処理の実行の度にランダムにソートされるため、特に乱数を生成せずに1つ目の要素を取得
        pl.col("top_echelons").list.get(0).alias("top_echelon")
    )
    .select("index", "hex", "top_echelons", "top_echelon", "n_top_echelon", "事故件数")
)

圏ごとに集計・ポリゴン化

def cell_to_hex_polygon(hex: str):
    boundary = h3.cell_to_boundary(hex)
    polygon = Polygon([[latlon[1], latlon[0]] for latlon in boundary])
    return to_wkt(polygon)

# トップのエシェロンごとにHexの輪郭をポリゴン化
df_top_echelon_polygon = (
    df_result
    .with_columns(
        hex_lonlat = pl.col("hex").map_elements(cell_to_hex_polygon, return_dtype=pl.String)
    )
    .with_columns(
        pst.from_wkt(pl.col("hex_lonlat")).st.set_srid(4326).alias("geometry")
    )
    .group_by("top_echelon")
    .agg(
        pl.mean("事故件数").round(1).alias("平均事故件数"),
        pl.col("geometry").st.union_all().st.to_wkt()
    )
    .sort("top_echelon")
)

2-4. エシェロン解析を試す 〜結果確認〜

上記処理で抽出した圏を平均事故件数に応じて色分けして可視化すると、以下のようになります（周辺に事故が発生したHexがない箇所については無視）。

【圏の抽出結果】

年間の事故件数の平均が1.5件以上の圏に絞り込むと、特に歓楽街を中心とした地域や駅周辺、バイパス、大きい道路が交差する地域などに事故が集中している傾向が抽出されています。

Hexを用いて格子データを作成した都合で圏の境界はなめらかではありませんが、エシェロン解析によって局所的な値の広がりを捉えることができました。

【圏の抽出結果（平均事故件数1.5以上）】

【圏の抽出結果（平均事故件数1.5以上・例：中心地）】

3. おわりに

今回はエシェロン解析を用いて空間上のデータの広がりを取得・可視化する方法についてご紹介しました。空間データ分析の一助となれば幸いです。

4. Reference

• 栗原考次・石岡文生（2021）『エシェロン解析　階層化して視る次空間データ』, 共立出版
• 小田牧子・石岡文生・正木隆・栗原考次（2012）「エシェロン階層構造を利用した森林の分割」, 『データ分析の理論と応用 Vol. 2, No. 1』, pp.17–31