ラスタタイルと addProtocol でメッシュ統計

## はじめに

:::note
これは [MapLibre AdventCalendar 2024](https://qiita.com/advent-calendar/2024/maplibre) 6日目の記事です。
:::

昨年に引き続き、今年も MapLibre アドベントカレンダーに参加させていただきます。

最近、[MapLibre でタイルの差分確認](https://qiita.com/mg_kudo/items/d934b7bd51a2605284d3)という記事を書きました。MapLibre GL JS の `addProtocol()` 機能を用いて、２種類のラスタタイルの演算を行うアイデアを試してみた記事です。

今回は、この記事の内容を応用して、ラスタタイルを「メッシュ統計データ」と考え、統計処理（という程、大それたものでもありませんが）をいろいろと試してみるという記事になります。

## メッシュ統計としてのラスタタイル

ラスタタイルは、地図データを一定の大きさ（一般的に 256 px 又は 512 px 四方）の画像データに分割したものとなります。

それぞれの画像データ（タイル）は、さらに個々の画素（ピクセル）に分割されたものとなりますので、この１つ１つのピクセルをデータと考えると、ラスタタイルは **「（タイル状に分割された）メッシュデータ」** とみなすことができます。

産総研では、「[グリッド PNG タイル](https://gsj-seamless.jp/labs/datapng/gridpngtile.html)」という概念を提案していますが、これは、数値データ等をピクセルの RGB 値にエンコードして格納するというアイデアに基づいています。

例えば、国土地理院では、標高値を RGB へエンコードした「[標高タイル](https://maps.gsi.go.jp/development/demtile.html)」が提供されています。実際、[地理院地図](https://maps.gsi.go.jp/)では、標高タイルを用いて、任意の位置の標高を表示できる他、標高タイルを面的に用いた「自分で作る色別標高図」や「地理院地図3D」といった機能も提供されています。

:::note
「タイルをデータベースとして扱う」という観点については、「[地図タイルへの入門知識を考える](https://qiita.com/mg_kudo/items/9f0b7a4f108ad96d9c25)」という記事でも紹介しています。
:::

さらに、１種類のタイルだけでなく、複数のタイルを組み合わせて利用するという発想も考えられます。

複数のメッシュデータを利用する際に重要となるのが、それらの **メッシュが適合するかどうか** です。メッシュがずれている場合、それぞれのメッシュを揃える処理が必要となります。

ラスタタイルに関して言えば、元データの測地系や投影法によって多少誤差は出ますが、基本的には、同じ仕様のタイルであれば、同じ大きさ・領域のタイル（及びピクセル）に分割されていると考えることができます。

つまり、一般的なラスタタイルを作成すれば統一されたメッシュが得られるため、結果的に「どのラスタタイルも標準化されたメッシュデータとして扱える」ということになります。

地図タイルはすっかりデファクトスタンダードとなり、様々な機関からタイルデータが提供されています。これらを全てメッシュ統計として活用できるという視点で捉えると、その応用範囲は大きく広がる可能性があります。もちろん意図しない使い方には注意しないといけませんが、より手軽にメッシュ統計の解析ができる環境が広がっているのかもしれません。


## Web Cookable Data と `addProtocol()`

[グリッド PNG タイル](https://gsj-seamless.jp/labs/datapng/gridpngtile.html) の設計コンセプトとして **「Web Cookable Data」(Webで調理可能なデータ）** というものが示されています。先ほど紹介した地理院地図の「自分で作る色別標高図」や「地理院地図3D」も、全てブラウザ上の処理で完結しています。

MapLibre GL JS では、`addProtocol()` という機能を利用することでタイルデータを加工することができます。この機能では、１種類のタイルのみだけでなく、複数のタイルデータを組み合わせて加工することも可能です。これもブラウザ上で完結する処理となっていますから、 MapLibre GL JS の `addProtocol()` 機能は「Web Cookable Data」を扱うのにうってつけと言えるでしょう。


## 実例

本記事では、あまり自分で試してみなかった組み合わせでラスタタイルを応用してみたいと思います。

:::note
オープンデータとなっている[地理院タイル](https://maps.gsi.go.jp/development/ichiran.html)と[ハザードマップポータルサイト](https://disaportal.gsi.go.jp/hazardmapportal/hazardmap/copyright/opendata.html)を利用しました。
:::

:::note warn
なお、提供側からすれば想定外の利用方法をしている場合がありますので、ご注意ください。
:::

### 内水・洪水・津波・高潮の最大浸水深

ハザードマップポータルサイトでは、内水氾濫、洪水、津波、高潮の各浸水ランクのタイルデータが提供されています。[重ねるハザードマップ](https://disaportal.gsi.go.jp/hazardmap/maps/index.html)では、特定の地点において、これらの災害の想定最大浸水深は調べることができますが、面的に把握することはできません。

そこで、これら４種のタイルデータと `addProtocol()` を用いて、最大の浸水ランクを表示する地図を作成してみました。４つのタイルを読み込んで演算するので、結構重くなります。

:::note warn
なお、各災害種別ごとに未整備の範囲がありますのでご注意ください。
:::
![サンプル1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2842832/d631f455-8e6b-5ee5-1958-c5700a6817e7.png)

https://mghs15.github.io/raster-tile-calc/maxdepth.html

### 高潮と津波の想定浸水深の差

海から水がやってくるという点で似たように思える高潮と津波ですが、場所によって想定される浸水深が異なることがあります。

そこで、浸水ランクの差を表示した地図を作成してみました。

青が高潮の方が高く、赤が津波の方が高いと想定される場所です。高潮は洪水もセットでシミュレーションされると思われますので、河川の周りが着色されています。どちらが高いかについては、全国まちまちのような気がします。

:::note warn
こちらも、高潮と津波で未整備の範囲がありますのでご注意ください。
:::
![サンプル2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2842832/37428e8c-093b-e8b3-bc29-cdb9c53ed535.png)

https://mghs15.github.io/raster-tile-calc/tsunami-takashio.html

### 傾斜量図と洪水の想定浸水深

[傾斜量図](https://www.gsi.go.jp/bousaichiri/slopemap.html)は、地表面の傾きの量を算出し、その大きさを白黒の濃淡で表現したものです。

洪水の際は、基本的に地形の緩やかなところに流れていくと思われますが、傾斜の激しいところも浸水するのか気になりましたので比較してみます。

浸水ランクが大きい（深い）ほど青く、傾斜量が大きいほど不透明度を高く表示しています。予想通り、堤防以外の部分で特に傾斜量が大きそうな場所はあまり見当たりません。ただし、急峻な地形ほど、勾配があっても浸水する様子が見られます。また、傾斜量が大きいところで浸水は止まるため、浸水範囲の辺縁部の傾斜量は大きくなる傾向です。

![サンプル3.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2842832/add6e8ba-4546-9132-d261-8918f9b50a76.png)

https://mghs15.github.io/raster-tile-calc/flood-slope.html

### 空中写真の差分比較

空中写真を比較する際には差分比較をしてみたいときが多いでしょう。本来は、機械学習等を用いるべき分野なのでしょうが、今回は非常に簡易的に RGB の合計値がどの程度変化したかで見てみます。

平成30年（2018年）の北海道胆振東部地震では、土砂災害で大きな被害が出ました。この土砂災害の際には国土地理院が空中写真を撮影しているので、地震前の写真（2011年）と比較することで、どの程度土砂災害が抽出できるか試してみました。

なお、国土地理院が当時土砂災害のエリアを判読しているので、これを正解として以下に並べてみます。

結果は……非常に微妙です。

やはり、もう少し工夫は必要かもしれません（が、正射画像さえあれば、ブラウザ上で簡単に試行錯誤して見られるのは利用可能性があるのではないかと思います。）

![サンプル4-1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2842832/8ef5d3c3-a3a5-b880-dd70-33d4bf9334b3.png)
![サンプル4-2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2842832/b70f6cdc-3dad-45e1-3615-eb74099f9a60.png)

https://mghs15.github.io/raster-tile-calc/photosdiff.html#12.5/42.76242/141.86785

### 植生指数と傾斜量図

日本においては、森林は山岳部に残っているというイメージですが、傾斜の少ないところに残っている森林があるのか気になります。（実用的には、森林管理のしやすさという指標につながりそうです。）

植物の量や活性を示すのには、[植生指数](https://www.gsi.go.jp/kankyochiri/ndvi.html)というものがよく使われます。国土地理院からは NDVI という植生指数のタイルが提供されていますが、今回は空中写真から、GRVI という、赤と緑だけを用いて計算できる指標を用いて、植物の量を示してみます。

```JavaScript:GRVI の算出
GRVI = (G - R) / (G + R)
```
以下の地図は、次のような条件で傾斜量が小さく、植生が高い場所を示した地図となります。

* GRVI が 0.15 以上のものに限定して表示
* 傾斜量図が一定の値よりも緩やかな場所のみを抽出（傾斜量図の濃さが傾斜量と厳密な比例関係なのか検証していませんが、傾斜15度をイメージして、RGBの平均値が `(90-15)/90*255` よりも緩やかな部分でフィルタしています。）

建物の影や水域も GRVI が高く表示されてしまう他、森林以外に農地も GRVI が高くなったり、そもそも写真の色合いで大きく変わったりするなど、課題はたくさんありますが、なんとなく状況はつかめそうです。

![サンプル5-1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2842832/366668c7-5437-279b-a8cf-b92dbe4e8c2f.png)
![サンプル5-2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2842832/36e1f739-6c6f-5ea9-28e6-6c3c3cc2ae70.png)

https://mghs15.github.io/raster-tile-calc/green.html

### 建物と災害リスク

地理院タイルの「標準地図」の建物データは橙色で目立ちます。この橙色の部分（建物に該当するピクセル/メッシュ）を抽出し、他の情報と組み合わせることで、特定の条件の建物を抽出できるのではないかという発想になりました。

今回は、重ねるハザードマップの洪水浸水想定や傾斜量図（土砂災害等のリスク）を組み合わせてみることにしました。

以下は、傾斜量に応じて建物部分を着色したり、浸水想定区域内の建物を抽出したりした例となります。単純に重ねるよりも明瞭に目的のピクセルを強調できているのではないでしょうか？

:::note warn
標準地図は **基本測量成果** です。２次利用される際は、測量法に抵触しないようにご注意ください。
:::

:::note warn
建物を抽出しようとして、国道もご抽出されていますのでご注意ください。
:::

![サンプル6-1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2842832/a8167ed4-d610-0ee5-8614-486c268e4cdb.png)

https://mghs15.github.io/raster-tile-calc/building-slope.html

![サンプル6-2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2842832/68a03a9a-94c5-ab87-ae6f-da3c8c42bc19.png)

https://mghs15.github.io/raster-tile-calc/building-flood.html

### レポジトリ

https://github.com/mghs15/raster-tile-calc

## 課題

今回、試行錯誤してみた結果、以下のような課題を感じました。

* 統計処理といっても単にピクセル毎に RGB 値を計算するだけとなるため、算出できる指標や応用方法には限界がある
* 適用可能なモデル式があるならともかく、読み込んだ範囲を標本としてパラメータを調整しながら結果を出力するという処理は、`addProtocol()` で行うにはそぐわない

ただし、これらはブラウザ上で処理を行う制約に由来する課題であり、GIS ソフトや統計ソフトを使えば解決するものです。

## おわりに

実例は少し微妙だったかもしれませんが、ラスタタイルの特性を活かしていくつかの分析をしてみました。

最近はラスタタイルからベクトルタイルへの移行も進んでおりますが、メッシュ統計としてタイルをとらえると、ラスタタイルにはベクトルタイルには代え難い魅力を感じています。

しかも、MapLibre GL JS の力を借りることで、データを用意することなく、新しい指標を地図に表現できます。ウェブ地図特有のシームレスなタイル読込・表示管理を気にすることなく、データ処理ができる `addProtocol()` 機能の強力さを改めて実感しました。

タイルデータは今や様々な機関から提供されており、統計データの蓄積という点で、その他のメッシュ統計セットと比べても有望なリソースになっているのではないかと思います。日本では、「地域メッシュ」によるメッシュ統計に歴史がありますが、最近は[空間 ID の仕様](https://www.ipa.go.jp/digital/architecture/guidelines/4dspatio-temporal-guideline.html)でタイルの仕組みが登場する等、メッシュ単位としてタイルを利用する動きも広がってきていると思われます。

実は、今回は人口のデータを使いたかったのですが、タイル化されている手ごろなものが見つけられず、実例には加えませんでした。もしも、タイルをメッシュ統計として利用するトレンドが盛り上がっていくとしたら、各種統計データがタイルデータとして提供されるようになるのかもしれません。

## 参考文献

記事中で引用したもの以外の参考文献を記載します。また、本記事は一部 ChatGPT を用いて添削しています。

* 佐藤彰洋, 2019, 『メッシュ統計』, 共立出版
* 黄川田智洋 他, 2022, 『飼料作物育種におけるGreen Red Vegetation Index（GRVI）と他のRGB植生指標との比較』, 育種学研究 2022 24 巻 2 号