MapLibre でヒートマップやクラスタリング #GIS

はじめに

これは MapLibre AdventCalendar 2023 23日目の記事です。

アドベントカレンダーに初めて参加させていただきます。

以前、小縮尺地図におけるポイントデータのデザインという記事を書いてみました。Web 地図においてポイントデータを表示する際に、データが大量にある場合、パフォーマンス上や見た目の観点により、小縮尺に全てのポイントを表示することはできない場合があります。その場合の対応についてアイデアを書き連ねてみたのですが、今回はこれらのアイデアについて、MapLibre GL JS での実装を確認しようという趣旨の記事になります。

表示方法を考える際の観点

MapLibre GL JS で多量のポイントデータの表示方法を考える際の観点を簡単に整理してみます。

手法について

詳細は、以前の記事に記載しておりますが、多量のポイントデータを表示するテクニックとして、以下のようなものが考えられます。これらの手法を組み合わせていくことになります。

省略：そもそも一部のデータを表示しない（間引き）。
クラスタリング：複数のポイントをまとめて１つのポイントして表示する。
ヒートマップ：ポイントデータを密度として表示する。

用意するデータについて

データについて、クライアントのブラウザ上で処理するか、データの段階でクラスタリングなどの処理を適用しておくかも判断する必要があります。

MapLibre GL JS で利用するベクタ型の主なデータ形式は、GeoJSON かベクトルタイル（Mapbox Vector Tile）かと思いますが、以下のように対応が異なります。

GeoJSON：ブラウザ上での処理も可能な場合が多いと思います。ただし、あまりに大きいデータサイズの場合は取り扱いが難しくなりそうです。
ベクトルタイル：大きいデータの場合、タイル化することで表示パフォーマンスを向上させられます。ただし、ブラウザ上での加工は難しいので、ベクトルタイル生成の段階でクラスタリングなどの処理を行っておきたいところです。

（参考）今回の例示用データについて

国土地理院が提供している指定緊急避難場所のデータを例に、実装例を紹介できればと思います。

指定緊急避難場所は、災害対策基本法により、市町村長が災害種別（異常な現象の種類）ごとに指定することになっています。

データの属性値としては、以下の内容が含まれています。

属性名	説明
NO	市町村別に1から連番
施設・場所名	指定緊急避難場所の名称
住所	指定緊急避難場所の住所、「都道府県／市町村／丁／番地」の順
洪水	対象とする異常な現象の種類が「洪水」の指定緊急避難場所、該当は「1」、非該当は無記入
崖崩れ、土石流及び地滑り	対象とする異常な現象の種類が「崖崩れ、土石流及び地滑り」の指定緊急避難場所、該当は「1」非該当は無記入
高潮	対象とする異常な現象の種類が「高潮」の指定緊急避難場所、該当は「1」、非該当は無記入
地震	対象とする異常な現象の種類が「地震」の指定緊急避難場所、該当は「1」、非該当は無記入
津波	対象とする異常な現象の種類が「津波」の指定緊急避難場所、該当は「1」、非該当は無記入
大規模な火事	対象とする異常な現象の種類が「大規模な火事」の指定緊急避難場所、該当は「1」、非該当は無記入
内水氾濫	対象とする異常な現象の種類が「内水氾濫」の指定緊急避難場所、該当は「1」、非該当は無記入
火山現象	対象とする異常な現象の種類が「火山現象」の指定緊急避難場所、該当は「1」、非該当は無記入
指定避難所との住所同一	指定避難所と住所が一致するする指定緊急避難場所、該当は「1」、非該当は無記入 ※2023年6月1日からしばらくの間、「指定避難所との重複」データが「指定避難所との住所同一」データとして公開されているようです。
緯度	緯度の値、 10進法
経度	経度の値、 10進法
備考	補足事項がある場合に記載

GeoJSON の場合

データの省略（フィルタ）

個々の指定緊急避難場所の優先度をつけるのは難しいですが、災害種別ごとに表示するという方法が考えられます。スタイル設定の filter を用いることで実装可能です。

「洪水」と「内水氾濫」を兼ねている指定緊急避難場所のみ表示する Style 設定例

map.addLayer({
  id: "overlay",
  type: "circle",
  source: "overlay",
  filter: ["all",
    ["==", ["get", "洪水"], "1"],
    ["==", ["get", "内水氾濫"], "1"],
  ],
  paint: {
    // 省略
  }
});

施設・場所名に「小学校」か「中学校」を含む指定緊急避難場所のみ表示する Style 設定例

map.addLayer({
  id: "overlay",
  type: "circle",
  source: "overlay",
  filter: ["any",
    ["in", ["get", "施設・場所名"], "小学校"],
    ["in", ["get", "施設・場所名"], "中学校"],
  ],
  paint: {
    // 省略
  }
});

また、そもそも MapLibre GL JS へ渡す前に、GeoJSON を加工してしまうという方法も考えられます。

オンザフライでクラスタリング

GoeJSON の場合、MapLibre GL JS ではオンザフライでクラスタリングをしてくれます。以下のように、MapLibre の source に cluster や clusterMaxZoom を設定することでクラスタリングが可能です。（なお、ベクトルタイルの場合はこのクラスタリング機能はありません。）

GeoJSON Source のクラスタリング向け設定例

map.addSource("overlay-cluster", {
  "type": "geojson",
  "data": "./path/to/your.geojson",
  "cluster": true,
  "clusterMaxZoom": 10, // クラスタリングを行う最大ズームレベル
  "clusterProperties": {
    "合計-洪水": ["+", ["to-number",["get", "洪水"]]],
    "合計-崖崩れ、土石流及び地滑り": ["+", ["to-number",["get", "崖崩れ、土石流及び地滑り"]]],
    "合計-高潮": ["+", ["to-number",["get", "高潮"]]],
    "合計-地震": ["+", ["to-number",["get", "地震"]]],
    "合計-津波": ["+", ["to-number",["get", "津波"]]],
    "合計-大規模な火事": ["+", ["to-number",["get", "大規模な火事"]]],
    "合計-内水氾濫": ["+", ["to-number",["get", "内水氾濫"]]],
    "合計-火山現象": ["+", ["to-number",["get", "火山現象"]]]
  }
});

まず、clusterMaxZoom は、クラスタリングを行う最大ズームレベルです。データの量や密度と相談し、どこまでクラスタリングを行うか、どこから個々のポイントデータを表示するようにするかも腕の見せどころかと思います。（簡単に設定できるので、私は何度も試行錯誤するタイプです。）

生成されたクラスタリング後のポイントデータには、以下のような属性が追加されます。

cluster：クラスターの場合 true となる。
cluster_id：クラスタリングの際に付与されたユニークな ID。
point_count：クラスタリングでグループ化されたポイント数。
point_count_abbreviated：省略表記のポイント数（例えば、7392 は 7.4k となる）。

さらに、clusterProperties で指定した条件で合計または最大値を計算して、指定した属性名で追加してくれるようです。上記の例では、各災害種別で指定されている指定緊急避難場所の数を合計するように設定しています。（例えば、クラスタリングされた「洪水」の指定緊急避難場所の数は、「合計-洪水」の属性値に格納されます。）

他にもクラスタリングされる際の半径 clusterRadius など、調整できるパラメータがあります。

クラスタリングされたデータの表示

上記の GeoJSON Source の設定でクラスタリングされたデータについては、すでにポイントデータですので、通常のポイントデータと同じように設定すれば表示ができます。

クラスタリングでよく見るのは、クラスタリングされたデータの数を反映した大きさの「円」とその「数字」を表示したものとなりますので、例として以下のような設定をお示しいたします。

クラスタリング用の Style 設定例

// 「円」用のスタイル
map.addLayer({
  id: "overlay-cluster",
  type: "circle",
  source: "overlay-cluster",
  filter: ["has", "point_count"],
  paint: {
    "circle-blur": 0.2,
    // クラスタリングされたデータ数によって色を変更
    "circle-color": [
      "interpolate", ["linear"],
      ["get", "point_count"],
      1,
      "#51bbd6",
      100,
      "#f1f075",
      750,
      "#f28cb1"
    ],
    // クラスタリングされたデータ数によって大きさを変更
    "circle-radius": [
      "interpolate", ["linear"],
      ["get", "point_count"],
      1,
      10,
      100,
      20,
      750,
      30
    ]
  }
});
// 「数字」用のスタイル
map.addLayer({
  id: "overlay-cluster-count",
  type: "symbol",
  source: "overlay-cluster",
  filter: ["has", "point_count"],
  layout: {
    "text-field": ["get", "point_count_abbreviated"],
    "text-font": ["NotoSansJP-Regular"],
    "text-size": 12
  },
  paint: {
    "text-halo-width": 1,
    "text-halo-color": "#FFFFFF"
  }
});

今回のアドベントカレンダーでは、@sanskruthiya 様がさらに凝った表示方法を紹介されていました。（@sanskruthiya 様は、ベクトルタイルを利用されております。クラスタリング自体はベクトルタイルの生成時に行っており、MapLibre GL JS 上では、ベクトルタイルを用いた描画結果を GeoJSON として取得した上で、Marker としてリッチな装飾を実現されています。）

ヒートマップ

ヒートマップは、ポイントデータの密度を表現できる方法です。MapLibre GL JS でもヒートマップ用のスタイル設定（heatmap タイプ）があります。

特にたくさんのデータがある場合、heatmap-intensity の設定が重要です。これは個々のポイントデータの重みづけの設定ですが、密度や数によって調整しないと見づらくなってしまいます。また、Web 地図の場合は、ズームレベルによって、点密度が変化しますので、ズームレベルごとにどのような調整にするかは結構難題かと思います。あまり良い考えは持っておりませんが、ズームレベルが１上がるごとに、ある対象領域が画面上に表示される面積は４倍になりますから、そのような観点をうまく使うのが良いかもしれません。（とはいえ、以下のサンプルは単純に比例関係で調整しておりますが……。）

ヒートマップ用の Style 設定例

map.addLayer({
  "id": "overlay-heat",
  "type": "heatmap",
  "source": "overlay",
  "maxzoom": 11,
  "paint": {
    "heatmap-color": [
      "interpolate",
      ["linear"],
      ["heatmap-density"],
      0,
      "rgba(33,102,172,0)",
      0.2,
      "rgb(103,169,207)",
      0.4,
      "rgb(209,229,240)",
      0.6,
      "rgb(253,219,199)",
      0.8,
      "rgb(239,138,98)",
      1,
      "rgb(178,24,43)"
    ],
    "heatmap-intensity": [
      "interpolate",
      ["linear"],
      ["zoom"],
      6,
      0.02,
      9,
      1
    ],
    "heatmap-opacity": 0.5
  }
});

ヒートマップ・クラスタリングのパフォーマンス

さて、ここで、ヒートマップとクラスタリングを行うことでどの程度パフォーマンスが向上するのか、描画速度を比較してみたいと思います。

クラスタリングの場合、点の描画数を抑えることができるからか、「クラスタリングを行う」という作業が入るとしても、大幅にパフォーマンスを向上させることができます。以下の実験では、クラスタリングの方は clusterProperties を用いて、８つの災害種別ごとの集計も行わせていますが、それでも圧倒的にクラスタリングを行った方が早いです。

ヒートマップはパフォーマンス改善よりも、「見た目の問題」を解決する手法という印象を持っています。実際に、今回の指定緊急避難場所データを全て（約11万点）をそのままポイント（circle レイヤ）とヒートマップ（heatmap レイヤ）を比較しましたが、有意差は見られませんでした。

実験条件

MapLibre GL JS (v3.6.2) で Map オブジェクトの作成から、idle イベントが起きるまでの時間を測定。
指定緊急避難場所のみを表示（背景地図などは表示しない）。
Chrome のデベロッパーツールで、キャッシュを無効化して実施。
GeoJSON データ、Maplibre GL JS のライブラリ（js、css、map）は、localhost （python による簡易サーバでホスト）から取得。
表示のたびに、処理をランダムに変更。（ランダムに表示させている関係上、各処理数ごとのサンプル数はばらついている。）
比較は、Welch の t 検定を使用（R v4.2.1）。
その他、実験に用いたコードはこちら。

クラスタリングとヒートマップの組み合わせ

ここまでを踏まえ、個人的に、クラスタリングとヒートマップを組み合わせる方法の採用が有効なのではないかと考えています。

パフォーマンスを向上させるためには、クラスタリングによりポイントデータの数を減らすことが有効です。しかし、クラスタリングを行うことによって、誤解を招くような位置に表示されてしまうような事例（例えば、クラスタリング後の代表点が他の市町村に表示されてしまう等）があります。（同様の指摘は、先にご紹介した@sanskruthiya 様の記事でも指摘されています。）

そこで、クラスタリング後のデータに対して、さらにヒートマップ表示を行うことで、あえてポイントの存在を曖昧にさせる、というのが、この方法の意図です。規模感をつかむために、クラスタリングの数字データのみは表示させるような工夫も有効だと思います。

デモサイト

国土地理院が提供している指定緊急避難場所の CSV データを利用者がアップロードし、検索や地図表示、結果のダウンロード等を行えるサイトです。（お試しの際は「サンプルデータ」ボタンからどうぞ。）

内部的に CSV を GeoJSON へ変換し、上述したようなテクニックを利用して地図表示を行っています。データ容量的にはベクトルタイル化したいところですが、利用者のデータを分析するという観点から、ベクトルタイルではなく、（利用者によってアップロードされた CSV を変換した）GeoJSON を利用しています。

レポジトリはこちら

ベクトルタイルの場合

ベクトルタイルの場合、MapLibre GL JS ではオンザフライのクラスタリングはできません。そのため、ベクトルタイル生成時に、クラスタリングを行っておく必要があります。

属性値に合わせたグループ化

通常のクラスタリングについては、@sanskruthiya 様の記事で紹介されておりますので、今回の場合、別のアプローチをとってみたいと思います。

指定緊急避難場所は市町村長が災害種別ごとに指定するという趣旨から、単純に距離が近いものをまとめるのではなく、「市町村」や「都道府県」ごとに集約するのも自然かと考えられます。小縮尺では、指定者である「市町村」ごとのポイントデータとして表示し、さらに小縮尺になれば、各市町村が所属する「都道府県」ごとのポイントデータとして表示するという形にすれば、スムーズにズームレベル間の移動ができるのではないかと期待できます。

今回の元データは CSV なので、集計は Excel のピボットテーブルで行い、都道府県ごと・市町村ごとの結果を CSV で出力しました。CSV は、元データも合わせて GeoJSON へ変換し、felt/tippecanoe でベクトルタイル（PMTiles 形式）へ変換しました。

ZL0～6: 都道府県ごとのデータ（集計は Excel のピボットテーブル）
ZL6～8: 市町村ごとのデータ（集計は Excel のピボットテーブル）
ZL9～10: 距離をもとにクラスタリング（tippecanoe の機能を利用）
ZL11: 個々の指定緊急避難場所のデータ

tippecanoe のコマンド例

# 都道府県ごとのデータを変換
tippecanoe -l group \
           都道府県SKHB.csv.geojson \
           -o pref.pmtiles \
           --force --no-tile-size-limit --no-tile-compression --no-feature-limit \
           --minimum-zoom=0 --maximum-zoom=6 --base-zoom=0

# 市町村ごとのデータを変換
tippecanoe -l group \
           市町村SKHB.csv.geojson \
           -o muni.pmtiles \
           --force --no-tile-size-limit --no-tile-compression --no-feature-limit \
           --minimum-zoom=7 --maximum-zoom=8 --base-zoom=7

# 個々の指定緊急避難場所のデータを変換
# （クラスタリング処理も行う）
tippecanoe -l cluster \
           全国指定緊急避難場所データ.csv.geojson \
           市町村SKHB.csv.geojson \
           -o clustered.pmtiles \
           --force --no-tile-size-limit --no-tile-compression --no-feature-limit \
           --minimum-zoom=9 --maximum-zoom=11 --base-zoom=9 \
           --cluster-distance=50 --cluster-maxzoom=10 \
           --accumulate-attribute=洪水:sum \
           --accumulate-attribute=崖崩れ、土石流及び地滑り:sum \
           --accumulate-attribute=高潮:sum \
           --accumulate-attribute=地震:sum \
           --accumulate-attribute=津波:sum \
           --accumulate-attribute=内水氾濫:sum \
           --accumulate-attribute=火山現象:sum \

# 最後に上記３つのベクトルタイルを統合
tile-join  -o skhb.pmtiles \
           pref.pmtiles muni.pmtiles clustered.pmtiles \
           --force --no-tile-size-limit --no-tile-compression

特に、--accumulate-attribute を活用することで、GeoJSON Source 設定の clusterProperties と同じようなことができます。

作成したベクトルタイルを表示すると以下のようになります。MapLibre GL JS におけるスタイル設定については、GeoJSON の場合と変わりません。

スタイルなどは要調整ですが、都道府県→市町村→クラスタリング→個々のポイントデータとフォーカスして行けるかと思います。

都道府県ごと

市町村毎

距離ベースでのクラスタリング

個別のポイント表示

デモサイト

もう少し調整が必要な様子がわかると思います……。

終わりに

小縮尺地図におけるポイントデータのデザインで考えた内容を Maplibre GL JS で実際に実装してみました。スタイル設定だけでも色々なことができる一方、見やすさにこだわろうとすると、さらに細かい調整をクライアントサイドでもデータ生成段階でも行う必要があろうかと思います。

うれしいことに、MapLibre GL JS（及びその派生元の Mapbox GL JS）や tippecanoe 等、Web 地図を取り巻くエコシステムは充実してきており、Web 地図の構築やベクトルタイルの生成はずいぶん簡単にできるようになってきたと感じています。だからこそ、今回のような見た目やパフォーマンスの問題は、最初から詳細に設計していくよりも、作って見て壊しての試行錯誤が求められるものであると認識しています。