FOSS4Gを活用した国土数値情報のベクトルタイル化とデータPNG化の試行について

:::note info
この記事は [FOSS4G Advent Calendar 2024](https://qiita.com/advent-calendar/2024/foss4g) 9日目の記事です。
:::
# はじめに
2024年2月、国土交通省政策統括官付地理空間情報課のX公式アカウント[@GIS_MLIT](https://x.com/GIS_MLIT)が開設されました。また、2024年10月からは[第1回 国土交通省 地理空間情報データチャレンジ ～国土数値情報編～](https://signate.jp/competitions/1484)が開催されており、さらに2024年11月には[地理空間情報課ラボの公式ホームページ](https://www.mlit-gis-lab.jp/)が開設されるなど、地理空間情報の活用に対する機運が高まっています。また、地理空間情報課ラボの公式ホームページで紹介されている[実験的取り組み](https://www.mlit-gis-lab.jp/lab/)の一つとして、[国土数値情報のベクトルタイル化に向けた検証](https://www.mlit-gis-lab.jp/lab/lab04/)の取り組みが紹介されています。ベクトルタイルは、データ容量が軽く、機械判読可能であることから、ベクトルタイル化により国土数値情報の利活用シーンがさらに拡大することが期待されています。
本記事は、国土数値情報のメッシュデータについてベクトルタイル化およびデータPNG化を試行し、ベクトルタイルが不得意とするメッシュデータのデータPNGの適用可能性を検証したものです。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/205045/2ca394b0-4baf-1199-0f71-fec974ff0b9b.png)

# ベクトルタイル
ベクトルタイルは、地理空間データを地図として効率的に配信・表示するための技術となっています。通常のラスタータイル（画像タイル）とは異なり、ベクトルタイルは地図要素（道路、建物、境界線など）のジオメトリと属性情報を持つベクトル形式で表現されます。この形式により、データを動的にスタイル設定したり、ズームレベルに応じて詳細な情報を効率的に表示したりすることが可能となっています。
ベクトルタイルのつくりかたについては、こちらの記事に詳しく書かれています。

https://qiita.com/Kanahiro/items/ceeb20c158b4c70b62b6

# データPNG
データPNG(Data PNG)は、画像ファイルフォーマットである PNG (Portable Network Graphics) を使ってデータを表現します。PNGは標準で圧縮されウェブブラウザでも利用しやすく、インターネットを使って大量データを高速に扱えます。地図タイルでの利用を想定して、産業技術総合研究所の西岡氏によって考案され、CC BY 4.0 互換のライセンスで利用可能です。データPNGは、[国土地理院 標高タイル](https://maps.gsi.go.jp/development/ichiran.html#dem)や[国土交通省 重ねるハザードマップ](https://disaportal.gsi.go.jp/)等で利用されています。
データPNGは、例えば、256 x 256 = 65,536 個のピクセルを持つ正方形画像の一つのピクセルに、RGB各8ビット（合計24ビット）のデータを格納できます。これを利用して、標高、ハザードマップ、気象、人流のデータなど、メッシュ（グリッド）に分割されるデータの表現や転送に利用します。
データPNGの仕様については、こちらの産総研ホームページに詳しく書かれています。

https://gsj-seamless.jp/labs/datapng/

# 国土数値情報をベクトルタイル化する際の課題
国土数値情報はさまざまなジオメトリタイプ（ポイント、ライン、ポリゴン、メッシュ等）で提供されていますが、ベクトルタイルは細かいメッシュデータ（今のところの筆者所感では100mメッシュ以下など）を不得意としている印象です。100mメッシュよりも細かいメッシュ（10mメッシュや5mメッシュ等）を元データとする、国土地理院 標高タイルや国土交通省 重ねるハザードマップではデータPNG化で軽量化を行い、ラスタータイルとして配信されています。
例えば、こちらは[国土数値情報の土地利用細分メッシュデータ（100mメッシュ）](https://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-L03-b-2021.html)をベクトルタイル化（PMTiles）して[PMTiles Viewer](https://pmtiles.io/)で表示したものですが、Webマップでのメッシュの表示に時間がかかる状況です。
 - 前提条件
     - ベクトルタイル化前:15.6GB(GeoJSON)、ベクトルタイル化後:1.76GB(PMTiles)
 - 国土数値情報 土地利用細分メッシュデータ（100mメッシュ）のベクトルタイル（PMTiles）
https://pmtiles.io/?url=https://shiworks.xsrv.jp/pmtiles-data/mlit-ksj/L03-b-21.pmtiles

メッシュデータをベクトルタイル化する際にはズームレベルごとの地物・頂点間引き、簡素化処理が適用されず、ファイルサイズが大きくなりがちです。
 - [（参考）Mapbox ベクトルタイルの仕様](https://docs.mapbox.com/data/tilesets/guides/vector-tiles-standards/)

# 国土数値情報をベクトルタイル化する際の課題への対応
上記の課題への対応として、以下の2つの手法を用いた、国土数値情報のベクトルタイル化とデータPNG化を試行して、比較を行いました。
 - 手法1
     - 元データをディゾルブ（統合）して、ベクトルタイル化する手法
 - 手法2
     - 元データをディゾルブ（統合）し、データPNGに変換して、ラスタータイル化する手法

# 前提条件
- 対象データ
    - [国土数値情報 土地利用細分メッシュデータ（2021年度（令和3年度））](https://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-L03-b-2021.html)
    - 100mメッシュ単位のGeoJSONを日本全国でマージして作成
    - 日本全国約4,900万メッシュ
- 使用ツール
    - [GDAL/OGR](https://gdal.org/en/latest/)
        ```gdalinfo --version```
        ```GDAL 3.9.3, released 2024/10/07``` 

    - [felt/tippecanoe](https://github.com/felt/tippecanoe)
        ```tippecanoe -version```
        ```tippecanoe v2.70.1```

# 手法１：元データをディゾルブ（統合）して、ベクトルタイル化する手法
- 土地利用細分メッシュデータを土地利用種別でディゾルブ（融合）
```
ogr2ogr -f Parquet L03-b-21_dissolved.parquet L03-b-21.parquet -dialect sqlite -sql "SELECT ST_Union(geometry) AS geometry, \"土地利用種別\" FROM \"L03-b-21\" GROUP BY \"土地利用種別\"" --config OGR_SQLITE_CACHE 32768 --config CPL_DEBUG ON
```

- 投影法を変換
```
ogr2ogr -t_srs EPSG:4326 L03-b-21_dissolved_reprojected.parquet L03-b-21_dissolved.parquet
```

- GeoJSONの生成
```
ogr2ogr -f "GeoJSON" L03-b-21_dissolved_reprojected.geojson L03-b-21_dissolved_reprojected.parquet
```

- PMTilesの生成
```
tippecanoe -o L03-b-21_dissolved.pmtiles -Z0 -z11 L03-b-21_dissolved_reprojected.geojson -pf -pk
```

- 作成したベクトルタイルはこちらから閲覧可能です
https://pmtiles.io/?url=https://shiworks.xsrv.jp/pmtiles-data/mlit-ksj/L03-b-21_dissolved.pmtiles
- ファイルサイズ（PMTiles形式）はディゾルブ（統合）前と比較して概ね94%減（1.76GB⇒95.7MB）と軽量化
- ただし、「土地利用種別」をもとにディゾルブ（統合）を行ったため、元データに存在する「細分メッシュコード」や「衛星写真撮影年月日」の属性情報は失われている

# 手法2：元データをディゾルブ（統合）し、データPNGに変換して、ラスタータイル化する手法
 - 手法1で作成したディゾルブ（統合）後のデータを使用
 - ラスタライズ
     - もとの100mメッシュのサイズ(0.00125秒×0.0008333333秒)に合わせる
```
gdal_rasterize -a "土地利用種別" -init 0 -ot Int32 -tr (4.5/3600) (3.0/3600) -tap -co COMPRESS=DEFLATE -a_srs EPSG:4326 L03-b-21_dissolved_reprojected.parquet L03-b-21_dissolved_reprojected_gray.tif
```

 - 1バンド→3バンドに変換
     - 土地利用細分メッシュ（ラスタ版）データの土地利用種別とRGB値の対応表を参考にカラーテーブルを用意しておく
     - データなしのエリアは（0,0,0,0）とする
     https://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-L03-b_r.html
     ![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/205045/e598063e-e21c-e349-abde-d3a3768a4de3.png)
     - [color_table.txt](https://shiworks.xsrv.jp/pmtiles-data/mlit-ksj/color_table.txt)
     ![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/205045/4ce2c5ed-51dd-75d7-0ebf-d7797479d963.png)


```
gdaldem color-relief L03-b-21_dissolved_reprojected_gray.tif color_table.txt L03-b-21_dissolved_reprojected_rgba.vrt -of VRT -alpha
```

 - 高解像度化
     - gdal2tilesのリサンプリングはあまり厳密ではないので、処理するtifを高解像度にして誤差を減らす

```
gdalwarp -tr (4.5/3600/20) (3.0/3600/20) -r near L03-b-21_dissolved_reprojected_rgba.vrt L03-b-21_dissolved_reprojected_rgba_high_resolution.vrt
```

 - タイル化
```
gdal2tiles.py --xyz -r near -z 4-12 --processes=16 L03-b-21_dissolved_reprojected_rgba_high_resolution.vrt L03-b-21_dissolved
```

- 作成したラスタータイルはこちらから閲覧可能です
https://shiworks.xsrv.jp/raster-tiles/mlit-ksj/L03-b-21_dissolved/openlayers.html
- ラスタータイルのファイルサイズはPMTiles形式（95.3MB）と比較して同程度（86.0MB）
- ただし、データPNGで保持できるのは一つの属性情報のみ（土地利用種別のみ）
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/205045/73ef001e-3f1d-200d-2f65-2ed5db69b242.png)

引き続き検証を行います。