MapLibre GL JS＋ベクトルタイルでウェブ地図に動的な集計円グラフを描く

作成例

警察庁オープンデータ「交通事故統計情報」（参照リンク）より、2019-2022年の4年間に起きた交通事故について、場所と属性値を動的に集計して、マップ上に円グラフで表示しています。

デモサイトとなるウェブ地図のリンクはこちら（参照リンク）
円グラフで示すと、事故が集中している箇所について、どのような性質の事故が多いのかという考察に繋がります。

参考にしたもの

MapLibre GL JSの公式ページにあるExamplesのうち、下記のものをベースに作成しました。（参照リンク）

ただし、こちらはGeoJSONデータを読み込む形のため、データサイズが大きい場合などでベクトルタイル化して読み込む場合は、タイル生成の段階から若干のアレンジが必要になります。
この記事では、そのアレンジの要点について説明します。

作業工程

工程は概ね以下の３段階です。

1. 前処理：オープンデータをGeoJSONに加工して属性値を追加付与する
2. ベクトルタイル生成（Felt／Tippecanoeを利用）
3. MapLibre GL JSで可視化

前処理について

「警察庁交通事故統計情報」（参照リンク）から2019年〜2022年に起きた約130万件の交通事故データをCSV形式で取得し、緯度経度を加工した上で、QGISで読み込んでGeoJSON形式に変換しました。
QGISでシンプルに表示するとこういう状態です。

この段階での重要なポイントは、円グラフの集計対象としたい項目に「１」のフラグを立てておくことです。
この例では、交通事故のうち「車と歩行者の衝突事故」について'pedestrian_flag'という属性項目を追加した上で、該当する場合は「1」とフラグ立てしています。
QGISのフィールド計算機から下記の式で、事故当事者に歩行者が含まれているものに「1」を振りました。

case 
when "当事者種別（当事者A）" = 61 then 1
when "当事者種別（当事者B）" = 61 then 1
else 0
end

属性テーブルは下記のような形です。「歩行者関連」以外にも「夜間事故か否か」などいくつかフラグ項目を作成しました。

これをGeoJSONで出力して前処理は完了です。

ベクトルタイルの生成

さて、100万件を超えるポイントデータをウェブ地図でスムーズに読み込むには、ベクトルタイル化するのが最も手軽な手法だと思います。
ベクトルタイルの生成には、Felt／Tippecanoe（参照リンク） を使いました。
Dockerが使える環境であれば、上記の参照リンク先にあるFelt／TippecanoeのGitHubリポジトリを丸ごとダウンロードして、下記のコマンドでDockerイメージを生成すればセットアップできます。

docker compose up -d --build

なお、私の場合は作業フォルダをDockerコンテナと繋げる（マウントする）ため、Dockerイメージのセットアップ前に下記の「docker-compose.yml」ファイルをあらかじめtippecanoeと同じディレクトリに入れてから上記コマンドを実行しています。

docker-compose.yml

version: '1'
services:
  map-tippecanoe:
    restart: always
    build: .
    container_name: 'tippecanoe_felt'
    working_dir: '/root/'
    tty: true
    volumes:
      - /Users/作業ディレクトリのパス/opt:/root/opt

Felt／Tippecanoeの実行環境が整ったら、作成したGeoJSONデータを対象に下記のTippecanoeコマンドでタイル化を実行します。

tippecanoe -zg -o ta_jp_flags_clustered.pmtiles --no-tile-compression -r1 --cluster-distance=50 --cluster-densest-as-needed --accumulate-attribute='{"pedestrian_flag":"sum"}' ta_jp_flags.geojson

Tippecanoeのオプションは多様なので、解説します。

コマンド	概要
-zg	ズームレベルを自動設定する（ここはお好みで指定してもOK）
-o 出力ファイル名	出力ファイル名と形式を指定。今回はPMTilesで出力しています。
--no-tile-compression	出力結果を圧縮しない。（圧縮するとうまく読み込めないと思う）
-r1	低ズームレベルの時にポイントを間引かない（代わりにクラスタ化するので）
--cluster-distance=50	クラスタリングするピクセル距離（公式は10を推奨しているが、今回の目的ではもう少し粗いほうがいいので50に設定）
--cluster-densest-as-needed	密集箇所をクラスタ化する。クラスタにまとめられたポイントは「"clustered": true, "point_count":件数」という属性項目が自動付与される。
--accumulate-attribute='{"属性名":"sum"}'	指定した属性名の各点の値がクラスタ代表点に加算(sum)されて格納される。集計処理したい属性項目が複数ある場合は -accumulate-attribute='{"num_accident": "sum", "daynight_flag":"sum"}'という形で追加する。
-	最後の項目は入力ファイルとなるGeoJSONを指定しています。

こんな感じで処理が走ります。

For layer 0, using name "ta_jp_flags"
1296282 features, 147619798 bytes of geometry, 10556017 bytes of string pool
Choosing a maxzoom of -z12 for features typically 209 feet (64 meters) apart, and at least 18 feet (6 meters) apart
Choosing a maxzoom of -z17 to keep most features distinct with cluster distance 50 and cluster maxzoom 24
  99.9%  17/114260/52100  
# 100.0%  17/110297/56351  

注意点として、約130万件のベクトルタイル化＋クラスタリング処理を回すのに確か10時間くらいかかった気がします。
「--cluster-distance=50」の数値が大きくなると、より長く処理時間がかかる様子でした。

MapLibre GL JSで可視化

記事冒頭で載せたデモサイトの全コードはGitHub（参照リンク）でご覧いただけます。
この記事では、円グラフ描画に関連する箇所のみ抜粋して解説していきます。

//必要なモジュールの読み込み
import * as maplibregl from "maplibre-gl";
import * as pmtiles from 'pmtiles'; //この例ではPMTilesを使うのでこれが必要
import 'maplibre-gl/dist/maplibre-gl.css';
import './style.css';

//PMTiles用のプロトコル設定
const protocol = new pmtiles.Protocol();
maplibregl.addProtocol("pmtiles",protocol.tile);

//カラーセットの設定（後述のSVG生成コードに直接記述しても良い）
const colors = ['#4169e1', '#87cefa'];//Colors for pie charts. Can be changed by selected categories, but so far only one pattern is set. 

//マップの基本設定（お好みに応じて変更）
const map = new maplibregl.Map({
    container: 'map',
    style: 'https://tile.openstreetmap.jp/styles/osm-bright-ja/style.json',//下地図にはこれを使わせてもらってます。
    center: [140.000, 36.000],
    interactive: true,
    zoom: 7.5,
    minZoom: 2,
    maxZoom: 21,
    maxPitch: 60,
});

map.on('load', () => {
    map.addSource('ta_cluster', {
        'type': 'vector',//ポイントデータなのでベクタータイプとして読み込む
        'url': 'pmtiles://app/pmtiles/ta_jp_flags_clustered.pmtiles',//PMTilesファイルのパス
    });
    map.addLayer({
        'id': 'ta_pseudo',//円グラフ設置箇所としてマップ中に透明表示しておく擬似的レイヤ
        'source': 'ta_cluster',
        'source-layer':'ta_jp_flags',//レイヤ名はベクトルタイル生成に使ったGeoJSONのファイル名が該当する
        "minzoom": 2,
        "maxzoom": 16,
        'layout': {
            'visibility': 'visible',
        },
        'type': 'circle',
        'paint': {
            'circle-color': 'transparent',//透明にしておく
            'circle-stroke-color':'transparent',//透明にしておく
            'circle-radius': 1
        },
    });
    map.addLayer({
        'id': 'ta_single',
        'source': 'ta_cluster',
        'source-layer':'ta_jp_flags',
        "minzoom": 2,
        "maxzoom": 16,
        'filter': ['!=', 'clustered', true],//クラスタ化されていない単独の点は円グラフで表現せず、別レイヤで描画する
        'layout': {
            'visibility': 'visible',
        },
        'type': 'circle',
        'paint': {
            'circle-stroke-width':2,
            'circle-color': 'transparent',
            'circle-stroke-color':['step',['get','pedestrian_flag'],colors[1],1,colors[0]],//クラスタ化されていない単独の点は、歩行者関連事故か否かを枠線の色で区別する
            'circle-stroke-opacity': 0.9,
            'circle-radius': 8
        },
    });
    map.addLayer({
        'id': 'ta_cluster_label',
        'type': 'symbol',
        'source': 'ta_cluster',
        'source-layer':'ta_jp_flags',
        'minzoom': 2,
        'maxzoom': 16,
        'filter': ['!=', 'clustered', true],//クラスタ化されていない単独の点は円グラフで表現せず、別レイヤで描画する
        'layout': {
            'text-field': '1',//クラスタ化されていない、つまり１件の点なので「１」とラベル表示
            'text-font': ['Open Sans Semibold', 'Arial Unicode MS Bold'],
            'text-size': 12
        },
        'paint': {
            'text-color': '#111',
        }
    });
    
    //ここからが円グラフマーカーの生成部分 (see -> https://maplibre.org/maplibre-gl-js/docs/examples/cluster-html/)
    const markers = {};
    let markersOnScreen = {};   
    //マーカーを新たに更新する処理
    function updateMarkers() {
        const newMarkers = {};//このリストを初期化する
        const features = map.queryRenderedFeatures({layers: ['ta_pseudo']});
         ////マップ表示中の各featureに対して繰り返し処理
        for (let i = 0; i < features.length; i++) {
            const coords = features[i].geometry.coordinates;//座標値を格納
            const props = features[i].properties;//属性プロパティを格納
            if (!props.clustered) continue;//当該フィーチャがクラスタ化されていない、つまり単独の点の場合は処理をスキップする
            const id = props.fid + '_' + props.point_count;///クラスタのユニークIDとして、点のfidとポイント数の組み合わせを格納。fidだけだと、ズームイン・アウトした際に集約件数が違うのに同一クラスタ扱いされて情報が更新されないケースが出てくる。

            let marker = markers[id];//マーカーリストのクラスタidをmarkerに格納
            if (!marker) {//マーカーが既存の表示範囲内に存在しない場合
                const el = createDonutChart(props);//円グラフSVG生成を実行
                marker = markers[id] = new maplibregl.Marker({
                    element: el
                }).setLngLat(coords);//指定の座標値に円グラフを設置して当該markerをmarkeridのリストに含める。つまり設置済みのmarkerはidリストに記録された状態となる。
            }
            newMarkers[id] = marker;//newMarkersのリストに現在のマーカーのクラスタidを格納
            if (!markersOnScreen[id]) marker.addTo(map);//画面上に当該マーカーが存在しない場合はマップに追加
        }
        //画面外のマーカーがリストに入っている場合はそれらを除去する
        for (let id in markersOnScreen) {
            if (!newMarkers[id]) markersOnScreen[id].remove();
        }
        markersOnScreen = newMarkers;//現在表示範囲内のクラスタidが格納されたリストをmarkersOnScreenに格納
    }

    //円グラフの対象データが読み込まれた後や画面移動するたびに画面上のマーカーを更新する
    map.on('data', (e) => {
        if (e.sourceId !== 'ta_cluster' || !e.isSourceLoaded) return;//対象データが読み込まれていない状態や対象となるソース以外が読み込まれた場合はスルー
        map.on('move', updateMarkers);//マップを動かしている最中でも更新（どの動作の時に情報更新するかはお好みで調整）
        map.on('moveend', updateMarkers);//マップを動かし終わったら更新
        updateMarkers();//最初の読み込み時点で表示
    });
});

//フィーチャの属性情報に基づいて円グラフSVGを生成する処理
function createDonutChart(props) {
    const offsets = [];
    const counts = [props.pedestrian_flag, props.point_count - props.pedestrian_flag];//円グラフの各項目に対応する件数を格納（この例では、各クラスタにおける歩行者関連事故と、それ以外の事故の件数をリスト格納）
    let total = 0;//合計値を初期化
    for (let i = 0; i < counts.length; i++) {
        offsets.push(total);//カウントアップされた積み上げ数値がここにリスト格納される。これは円グラフ描画の際に必要となる。
        total += counts[i];//リストに入れた各項目の数値分カウントアップして最終的にクラスタ内の総件数になる。
    }
    //円グラフ中央に表示するクラスタ内総件数のフォント設定（お好みで調整）
    const fontColor = total >= 30000 ? "red" : "black";
    const fontSize = total >= 10000 ? 18 : total >= 1000 ? 16 : total >= 100 ? 14 : 11;
    //円グラフの件数に応じた大きさ設定（お好みで調整）
    const r = total >= 30000 ? 40 : total >= 10000 ? 32 : total >= 5000 ? 28 : total >= 1000 ? 24 : total >= 100 ? 21 : total >= 10 ? 18 : 12;
    const r0 = Math.round(r * 0.6);
    const w = r * 2;
    //設定に応じてSVGで図形を記述
    let html =
        `<div><svg width="${
            w
        }" height="${
            w
        }" viewbox="0 0 ${
            w
        } ${
            w
        }" text-anchor="middle" style="font: ${
            fontSize
        }px sans-serif; fill: ${fontColor}; display: block">`;

    for (let i = 0; i < counts.length; i++) {//円グラフの各項目（今回の場合は歩行者関連事故と、それ以外の事故の件数）について、donutSegment（後述の円グラフ各項目描画処理）を実施
        html += donutSegment(
            offsets[i] / total,
            (offsets[i] + counts[i]) / total,
            r,
            r0,
            colors[i]
        );
    }
    //円グラフの中央部分に白抜きでクラスタ内総件数をテキスト表示する
    html +=
        `<circle cx="${
            r
        }" cy="${
            r
        }" r="${
            r0
        }" fill="white" /><text dominant-baseline="central" transform="translate(${
            r
        }, ${
            r
        })">${
            total.toLocaleString()
        }</text></svg></div>`;
    //div要素にSVGを入れて出力
    const el = document.createElement('div');
    el.innerHTML = html;
    return el.firstChild;
}

//円グラフの各項目部分を描画する処理
function donutSegment(start, end, r, r0, color) {//ここのstartとendには、offsetsにリストで格納した積み上げ件数比率が入る
    if (end - start === 1) end -= 0.00001;
    const a0 = 2 * Math.PI * (start - 0.25);
    const a1 = 2 * Math.PI * (end - 0.25);
    const x0 = Math.cos(a0),
        y0 = Math.sin(a0);
    const x1 = Math.cos(a1),
        y1 = Math.sin(a1);
    const largeArc = end - start > 0.5 ? 1 : 0;

    return [
        '<path d="M',
        r + r0 * x0,
        r + r0 * y0,
        'L',
        r + r * x0,
        r + r * y0,
        'A',
        r,
        r,
        0,
        largeArc,
        1,
        r + r * x1,
        r + r * y1,
        'L',
        r + r0 * x1,
        r + r0 * y1,
        'A',
        r0,
        r0,
        0,
        largeArc,
        0,
        r + r0 * x0,
        r + r0 * y0,
        `" fill="${color}" fill-opacity="0.8"/>`//円グラフの透明度は好みで調整
    ].join(' ');
}

長くなりましたが、JavaScriptコードは以上です。
あとはMapLibreのmap要素を読み込んで表示するためのHTMLやCSSを記述すればOK。

MapLibre GL JSの円グラフ生成コードにおける要点

基本的なアプローチは、queryRenderedFeatures 関数によって表示範囲内のクラスタ化されているポイントデータ情報を収集し、その情報に基づきSVG形式の円グラフを生成し、マーカーとしてクラスタ中心点に配置する流れとなっています。
ここで重要なのは、queryRenderedFeatures関数はその性質上、現在の画面内に表示されている範囲のフィーチャ情報しか取得できない ことです。
つまり、画面を動かすたびに新たに表示された／画面外となったポイントデータの差分を把握して、新しい部分は円グラフを追加描画する仕組みが必要となります。その一連の機能が「updateMarkers」の部分です。

また、差分を把握するために、「id」には円グラフ用クラスタの代表点のIDとポイント数（何件の点がその代表点に集約されているか）を組み合わせたユニークIDを振っています。
公式サイトのサンプルコードでは"id"のみを格納していますが、クラスタリングありでベクトルタイル化した際には、同じIDを持つ代表点であっても集約している件数がズームレベルによって異なるため、ID+ポイント数の組み合わせでユニーク判定しています。

補足：追加アレンジについて

冒頭のデモサイト（参照リンク）では、ここからさらにアレンジを加えて、下記のような作り込みを行なっています。

円グラフ関連の追加アレンジ

• 円グラフ集計する内容を複数用意してセレクトボックスで選択可能にする
• 十分にズームした状態では円グラフではなく個別の事故情報表示に切り替える
• その表示状態に応じて凡例も動的に変更する

さらに余談となりますが、円グラフと関係ない部分では下記の工夫を加えています。

円グラフ関連以外で工夫したこと

• 表示中の位置とズームレベルをURLのハッシュ値で保持し、読み込み時のURLにハッシュ値が付いていればその情報を初期表示に反映させる（特定の場所の事故状況を他人にシェアしたい、というケースに対応するため）
• 現在地表示ボタンを設置。MapLibreのGeolocateControl （参照リンク） を使っても良いのですが、個人的にアイコンの意味がわかりにくい気がしたので別で実装しました。
• 場所検索機能を設置。これはOpenStreetMap Nominatimを使ったMapLibre公式サンプル（参照リンク）そのままです。
• 拡大時に個別の事故情報をポップアップで表示する際、表形式ではなく、文章形式で表現した（好みの問題ですが、ナラティブな方が個人的に内容が頭に入るため）

ハッシュ値付きURLやポップアップ表現などの具体的イメージは下図の通りです。

そのコードはこちらのGitHub（参照リンク）でご覧いただけます。

このアプローチの課題

この手法は、ベクトルタイルの恩恵によって100万件以上のポイントデータをスムーズに表示できることや、表示中の位置とズームレベルに応じて自動的に円グラフを描画できるのは良いのですが、クラスタリングの粒度や集約範囲、代表点の位置（中心点）を細かくコントロールできない のが難点です。

それが要因となって、直感的でない、ミスリードを招きかねないマップとなるケースがあります。
一例として、下図の中央部を見ると、まるで公園の中で11件の事故が起こっているように見えてしまいます。

拡大すると、実際には11件の事故は公園の周囲や駐車場で起きていることがわかります。
これらの事故を1つに自動的にまとめたので、公園内に11件分のポイント数を持つクラスタがプロットされてしまったわけです。

そのため、分析の際にはマクロな視点で地域の事故特性を掴むだけでなく、気になる部分は拡大して詳細な事故内容を確認することも同時に重要となる でしょう。

おわりに

この記事では、MapLibre GL JS＋ベクトルタイルでウェブ地図に動的な集計円グラフを描く一連の流れを解説しました。

前処理にて属性値にフラグ立てをしておくこと、ベクトルタイル生成時にクラスタリング設定をしておくこと、MapLibre GL JSのコードは公式サンプルからアレンジが必要なこと、それぞれの要点を説明しました。

コードの部分は冗長でわかりづらい部分もあるかと思いますが、コメントを振っているので参考になればと思います。

また、似たような可視化アプローチを使って、QGIS上で分析を行うケースについても解説しております。（関連Qiita記事リンク）。
もし興味を持っていただけましたら、こちらもあわせてご覧ください。