Go to Qiita Advent Calendar Top
LoginSignup
5
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

Day 7

FOSS4G Advent Calendar 2025

@shi-works(shi works)

数値PNGタイルの作成とMapLibre GL JSでの表示について

Last updated at Posted at 2025-12-06

この記事は FOSS4G Advent Calendar 2025 7日目の記事です。

はじめに

  • 洪水浸水想定区域図データ等のハザードマップ系データは、通常、パレットPNG で配信されている。しかし、パレットPNGは、カテゴライズされた「浸水深ランク」を扱うため、浸水深そのものを実数値として保持できない。
  • 一方で、 数値PNG は、浸水深を実数値として保持することが可能であり、フロントエンド側で「RGB → 浸水深 → パレット変換」を行うことで、例えば、 浸水深ランクのレンジやカラーの変更を柔軟に行うことができる。
  • 本記事では、香川県 洪水浸水想定区域図 をもとに、数値PNGタイルの作成と、作成した数値PNGタイルを MapLibre GL JS で表示する方法について検証する。
  • 具体には、香川県 土器川水系の洪水浸水想定区域図（想定最大規模）の 5m メッシュ（属性：浸水深[m]） から数値PNGタイルを作成し、MapLibre GL JS の addProtocol を用いて RGB → 浸水深 → パレット変換 を行いながら表示する。

数値PNGの概要

  • 数値PNGは、産業技術総合研究所の西岡氏により考案された データPNG の一種である。

image.png
出典: データPNG https://gsj-seamless.jp/labs/datapng/

  • 数値PNGは、平面上の格子データを格納するためのデータPNG であり、標高など 連続的に変化するスカラー量 を表現する仕様になっている。
  • 実数値を RGB 値にエンコードして各ピクセルに格納することで、ブラウザなどのクライアント側で画像として描画できる。
  • 可視化時には、クライアント側で任意のカラーマップ（例：標高タイルに対する陰影段彩図など）を適用できる。

参考文献

数値PNGタイルの作成方法

以下では、土器川水系の洪水浸水想定区域図（5m メッシュ）から、数値PNGタイルを作成する一連の処理フローを示す。

  1. シェープファイルの内容確認
  2. シェープファイルからGeoTIFFの生成（ラスタライズ）
  3. 実数値（m）ラスターを整数値（mm）ラスターに変換
  4. 32bit 整数値を RGB 各バンドに分解（24bit 化）
  5. α バンド（透明度）の作成
  6. 各バンドの NoData メタ情報をクリア
  7. 4 バンド（R, G, B, A）を結合して GeoTIFF を作成
  8. バンドごとの色解釈（Color Interpretation）の明示
  9. GeoTIFFにCRSを付与
  10. 高解像度化
  11. ラスタータイルの作成

処理 PC 環境

  • OS：Windows 11 Pro（WSL2 / Ubuntu）
  • CPU：AMD Ryzen 7 5700X
  • メモリ：64GB
  • SSD：4TB（空き容量約1TB）

使用ツール

  • GDAL 3.11.3
    • ogrinfo, gdal_rasterize, gdal_calc.py, gdalwarp, gdal2tiles.py, gdal_translate, gdal_edit.py, gdalbuildvrt

使用データ

香川県のオープンデータ：
https://opendata.pref.kagawa.lg.jp/dataset/572.html

対象：

  • 土器川水系 5m メッシュ（想定最大規模）
  • シェープファイル（約300MB）
├─241226_公表図面shp(土器川)
    006_土器川水系_5mメッシュ.dbf
    006_土器川水系_5mメッシュ.prj
    006_土器川水系_5mメッシュ.shp
    006_土器川水系_5mメッシュ.shx

シェープファイルの内容確認

  • まず、シェープファイルの座標参照系（"JGD2011"）や対象範囲（Extent: (133.749125, 34.093542) - (134.017125, 34.323417)）をogrinfoで確認する。
$ ogrinfo -al -so 006_土器川水系_5mメッシュ.shp
INFO: Open of `006_土器川水系_5mメッシュ.shp'
      using driver `ESRI Shapefile' successful.

Layer name: 006_土器川水系_5mメッシュ
Metadata:
  DBF_DATE_LAST_UPDATE=2024-12-26
Geometry: Polygon
Feature Count: 1700181
Extent: (133.749125, 34.093542) - (134.017125, 34.323417)
Layer SRS WKT:
GEOGCRS["JGD2011",
    DATUM["Japanese Geodetic Datum 2011",
        ELLIPSOID["GRS 1980",6378137,298.257222101,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    USAGE[
        SCOPE["Horizontal component of 3D system."],
        AREA["Japan - onshore and offshore."],
        BBOX[17.09,122.38,46.05,157.65]],
    ID["EPSG",6668]]
Data axis to CRS axis mapping: 2,1
MESH: Integer64 (15.0)
標高: Real (12.2)
浸水深: Real (19.11)

シェープファイルからGeoTIFFの生成（ラスタライズ）

  • 浸水深（m）をもとに、5mメッシュのラスターデータを作成する。
  • 細分メッシュの1辺の長さと緯度経度の対応表は以下のとおりである。
メッシュ1辺の長さ（約m） 緯度（秒） 経度（秒） 緯度（度） 経度（度）
100.00 3.00000 4.50000 0.0008333333333333330 0.001250000000
50.00 1.50000 2.25000 0.0004166666666666670 0.000625000000
40.00 1.20000 1.80000 0.0003333333333333330 0.000500000000
25.00 0.75000 1.12500 0.0002083333333333330 0.000312500000
20.00 0.60000 0.90000 0.0001666666666666670 0.000250000000
12.50 0.37500 0.56250 0.0001041666666666670 0.000156250000
10.00 0.30000 0.45000 0.0000833333333333333 0.000125000000
5.00 0.15000 0.22500 0.0000416666666666667 0.000062500000
4.00 0.12000 0.18000 0.0000333333333333333 0.000050000000
2.50 0.07500 0.11250 0.0000208333333333333 0.000031250000
2.00 0.06000 0.09000 0.0000166666666666667 0.000025000000
1.25 0.03750 0.05625 0.0000104166666666667 0.000015625000
1.00 0.03000 0.04500 0.0000083333333333333 0.000012500000

出典: 浸水想定区域図データ電子化ガイドライン（第 5 版）
https://www.mlit.go.jp/river/shishin_guideline/bousai/saigai/tisiki/syozaiti/pdf/e-guideline_5.pdf#page=32

  • gdal_rasterizeの主要オプションは以下のとおりである。
オプション 内容
-a suisin 浸水深属性を値として使用
-a_nodata -9999 無効値を-9999に設定
-ot Float32 出力型を32bit浮動小数点
-a_srs EPSG:6668 JGD2011地理座標系
-tr 0.0000625 0.0000416666666666667 ピクセルサイズを5m×5m
-co COMPRESS=DEFLATE -co TILED=YES 圧縮・タイル化設定
  • gdal_rasterizeを実行して、シェープファイルからGeoTIFFを生成する（ラスタライズ）。
gdal_rasterize -a "浸水深" -a_nodata -9999 -ot Float32 -a_srs EPSG:6668 \
  -tr 0.0000625 0.0000416666666666667 -tap \
  -co COMPRESS=DEFLATE  -co TILED=YES \
  "006_土器川水系_5mメッシュ.shp" flood_l1_shinsuishin_gray.tif

実数値（m）ラスターを整数値（mm）ラスターに変換

gdal_calc.py -A flood_l1_shinsuishin_gray.tif \
  --outfile=depth_i_u32.tif --type=UInt32 --NoDataValue=0 \
  --calc="where(A==-9999, 0, clip(round(A/0.001), 0, 16777215))"

  • 欠損値処理

    • A == -9999 の画素は欠損扱いとし、出力値0を設定する。
    • 欠損と「0mm」を区別しない方針（後続で α バンドで区別する）

  • 単位変換（m → mm）

    • 実数値（m）を mm に変換するため A / 0.001 を計算する。
    • 例：0.123[m] → 0.123 / 0.001 = 123 → 123mm

  • 四捨五入

    • round(...) によって小数点以下を四捨五入する。
    • mm 単位の整数値として扱うための丸め処理を行う。

  • 値域制限

    • clip(..., 0, 16777215) で、0〜16,777,215 の範囲にクリップする。
      • 16,777,215 = 2^24 - 1
      • 後続で 24bit（R, G, B 各8bit）の 3 バンドに分解するため、24bit に収まる最大値で制限する。

  • NoData 設定

    • 出力ラスタの NoDataValue として 0 を設定する。
    • ただし、この後の RGB 分解では NoData を使わず、α バンド＋メタ情報クリアで制御する。

32bit 整数値を RGB 各バンドに分解（24bit 化）

gdal_calc.py -A depth_i_u32.tif --outfile=R.tif --type=Byte --calc="floor(A/65536)%256" --NoDataValue=None
gdal_calc.py -A depth_i_u32.tif --outfile=G.tif --type=Byte --calc="floor(A/256)%256"   --NoDataValue=None
gdal_calc.py -A depth_i_u32.tif --outfile=B.tif --type=Byte --calc="A%256"              --NoDataValue=None

  • 整数値 N を 24bit の RGB に分解するイメージ：

    • N = R * 65536 + G * 256 + B
    • 各バンドは 0〜255 の範囲（8bit）

  • 具体的な式：

R（上位 8bit）
R = floor(N / 65536) % 256
  • N / 65536 で上位 16bit を取り出し、その下位 8bit 分を % 256 で切り出し
G（中位 8bit）
G = floor(N / 256) % 256
  • N / 256 で下位 8bit を除いた 24bit を取り出し、その下位 8bit 分を % 256 で切り出し
B（下位 8bit）
B = N % 256

  • 下位 8bit をそのまま取得

  • NoData の扱い

    • ここでは NoDataValue=None としており、R/G/B の各バンドに NoData メタは付与しない。
    • 元の欠損情報は、このあと作成する α バンドで表現する

α バンド（透明度）の作成

gdal_calc.py -A flood_l1_shinsuishin_gray.tif \
  --outfile=A.tif --type=Byte --calc="(A==-9999)*0 + (A!=-9999)*255" --NoDataValue=None
  • 元のラスターの値に応じて α 値を決定
A == -9999（欠損値）の画素 → α = 0（完全透明）
A != -9999 の画素 → α = 255（完全不透明）
  • 式の意味
(A == -9999) * 0 → 欠損画素は 0
(A != -9999) * 255 → 有効画素は 255
  • 上記二つを足し合わせることで、画素ごとに 0 または 255 のいずれかの α を生成

各バンドの NoData メタ情報をクリア

gdal_edit.py -unsetnodata R.tif
gdal_edit.py -unsetnodata G.tif
gdal_edit.py -unsetnodata B.tif
gdal_edit.py -unsetnodata A.tif
  • 理由
    • 表示や後続処理で、NoData メタ情報と α バンドの両方があると挙動が複雑・不整合になる可能性がある。
    • 本処理では 欠損＝α=0 で表現する方針のため、NoData メタ情報は不要

4 バンド（R, G, B, A）を結合して GeoTIFF を作成

  • VRT で 4 バンドを仮想的に結合
gdalbuildvrt -separate rgba.vrt R.tif G.tif B.tif A.tif
  • VRT から実体の GeoTIFF へ変換
gdal_translate -of GTiff -co COMPRESS=DEFLATE -co TILED=YES rgba.vrt depth_rgba_numeric.tif

バンドごとの色解釈（Color Interpretation）の明示

gdal_edit.py \
  -colorinterp_1 red \
  -colorinterp_2 green \
  -colorinterp_3 blue \
  -colorinterp_4 alpha \
  depth_rgba_numeric.tif
  • depth_rgba_numeric.tif の各バンドに対し、以下の色解釈を設定：
バンド1：red
バンド2：green
バンド3：blue
バンド4：alpha

GeoTIFFにCRSを付与

gdal_edit.py -a_srs EPSG:6668 depth_rgba_numeric.tif

高解像度化

  • gdal2tilesのリサンプリングはあまり厳密ではないので、処理するtifを高解像度にして誤差を減らす。
  • ここでは、20倍に高解像度化している。
gdalwarp -tr 0.000003125 0.0000020833333333333335 -r near depth_rgba_numeric.tif depth_rgba_numeric_high_resolution.vrt

ラスタータイルの作成

  • ラスタータイルの作成には、GDAL/OGRのgdal2tiles.pyを使用する。
  • ズームレベル4～18で作成する。
gdal2tiles.py --xyz -r near --s_srs epsg:6668 -z 4-18 --processes=8 depth_rgba_numeric_high_resolution.vrt flood_l1_shinsuishin

下記の表は、緯度35.5度付近のラスタータイルのズームレベルと解像度の対応表である。こちらの表をもとにすると、洪水浸水想定区域図データは5mメッシュであることから、最大ズームレベルは15（解像度3.945m）程度だが、Web地図での見やすさ・操作性に配慮して最大ズームレベルは18としている。

＜緯度35.5度付近のラスタータイルのズームレベルと解像度の対応表＞

ズームレベル ①地球の周長 (m) ②タイルサイズ ③1辺のタイル数 (2^ズームレベル) （参考）総タイル数 (=③^2) ④解像度 (m/pixel) (① / (② × ③))
14 33,090,083 256ピクセル 16,384 268,435,456 7.889
15 33,090,083 256ピクセル 32,768 1,073,741,824 3.945
16 33,090,083 256ピクセル 65,536 4,294,967,296 1.972
17 33,090,083 256ピクセル 131,072 17,179,869,184 0.986
18 33,090,083 256ピクセル 262,144 68,719,476,736 0.493
19 33,090,083 256ピクセル 524,288 274,877,906,944 0.247
20 33,090,083 256ピクセル 1,048,576 1,099,511,627,776 0.123
21 33,090,083 256ピクセル 2,097,152 4,398,046,511,104 0.062
22 33,090,083 256ピクセル 4,194,304 17,592,186,044,416 0.031

作成した数値PNGタイルURL

https://public-data.geolonia.com/pref-kagawa/flood_l2_shinsuishin/{z}/{x}/{y}.png

数値PNGタイルと重ねるハザードマップの比較

  • 数値PNGタイルは、MapLibre GL JSのaddProtocolでRGB→浸水深→パレット変換して描画した。
  • 数値PNGタイルをパレット変換した結果と重ねるハザードマップを比較。
  • 本記事では、洪水浸水想定区域図データのような 高解像度なメッシュデータ（5m）に対して数値PNGを適用できることを確認した。
  • MapLibre GL JS の addProtocol を利用して、RGB → 浸水深 → パレット変換 で描画した結果は、重ねるハザードマップと 視覚的に概ね一致 することも確認できた。

デモサイト:
https://shiwaku.github.io/kagawa-numerical-png-tiles-on-maplibre/index_2#15/34.29558/133.80672

GitHubレポジトリ:
https://github.com/shiwaku/kagawa-numerical-png-tiles-on-maplibre

重ねるハザードマップ（パレットPNGタイル)

image.png (1.6 MB)

数値PNGタイルによる色表現

  • MapLibre GL JS の addProtocol を利用して、数値PNGタイルをRGB → 浸水深 → パレット変換 で動的に描画が可能である。
  • 洪水浸水想定区域図データを数値PNGタイルで整備しておけば、例えば、クライアント側で、浸水深（実数値）をもとに、国交省ガイドライン・バリアフリー配色や東京都「詳細配色版」準拠に動的に配色の変更が可能である。

20251206_数値PNGタイル⇒パレット変換イメージ.png

課題と今後の展望

課題

  • パレット変換をすべてフロントエンドで行う場合、RGB → 浸水深 → パレット変換がクライアント側の負荷となり、5mメッシュでも描画がやや重くなる。
  • このため今後の対応策としては、
    • WebGL シェーダーによるクライアント側の高速化
    • プロキシサーバやタイルサーバ側でパレット変換を行い、パレットPNGを返す方式への切り替え

などが有力と考えられる。

今後の展望

  • 筆者が把握している範囲では、数値PNGの活用事例はまだ多くなく、主に標高タイルなど 連続値の格子データ が中心である。
  • 一方で、数値PNGは、浸水深を実数値として保持できる点 や、フロント側でレンジや色分けなどを柔軟に変更できる点が大きな利点である。
  • そのため、高解像度のメッシュで実数値を持つデータ（浸水深・人口(人流)・降水量など）に対して、数値PNGは今後 有力なデータ形式の一つになり得る。

＜数値PNGの適用事例＞

データPNGの種類 特徴 主な適用事例 対応するWeb地図ライブラリ（※OSSのみ）
グリッドPNG（数値PNG） ・数値PNGは評価値の格子データを扱うためのデータPNGで、値が連続的に変化するスカラー型のデータ（標高値など）を表現する。 ・国土地理院 標高タイル
・産総研 シームレス標高タイル		 ・Maplibre GL JS
・deck.gl（国土地理院 標高タイル）
グリッドPNG（パレットPNG） ・パレットPNGはグリッドPNGの一種で、いくつかの限定された色のみが使用される（ハザードマップの浸水深ランクなど）。 ・国土交通省「重ねるハザードマップ」データ配信
・産総研 20 万分の 1 日本シームレス地質図		 ・Maplibre GL JS
・Mapbox GL JS（ver1.3）
・OpenLayers
・Leaflet
・Cesium
リストPNG（点群PNG） ・点群PNGはポイントデータを表現するためのデータPNGで、ポイントデータは 2 次元または 3 次元の座標値を持ち、それに紐づけられる何らかのデータを持つ。 ・3 次元点群データ ・Leaflet ※ただし 2 次元表示のみ
・Web地図ライブラリではないが Three.js で 3 次元表示が可能
5
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
5
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?