富士山周辺の標高RGBタイルの作成

TLDR

数日前、資料を調べていると MapTiler が作成した軌跡標高データ可視化ページを見つけ、非常に良い出来でした：

这个做的是真好的啊，而且完全在浏览器实现的https://t.co/R1vQMbdYAo pic.twitter.com/IyTYy5Fxo1
— ringsaturn.me (@ringsaturn_me) January 10, 2024

公式のブログで2019 年のブログ記事を見つけ、使用している技術は標高データを画像の RGB チャンネルの値にエンコードし、ブラウザでデコードするというものでした。この方法によりブラウザ上で標高データの可視化が可能になり、とくに軌跡データのような連続的なデータを、数枚の画像をダウンロードするだけで連続した標高プロファイルを描画できます。

GitHub でsyncpoint/terrain-rgbプロジェクトを見つけ、標高データを RGB 画像に変換するための完全な手順が提供されていました。

整理したコードをリポジトリringsaturn/play-terrain-rgbに置きましたが、処理の流れは@syncpointのものとまったく同じで、いくつかの詳細なパラメータと元データだけが異なります。

また、https://github.com/maptiler/samples/tree/main/cloud/elevation-profileのコードを修正して、富士山周辺の標高データ可視化ページを公開しました。結果は以下のようになります:

デモはこちらで確認できます：https://ringsaturn.github.io/play-terrain-rgb

作成手順

ASTER GDEM v3 のデータをダウンロードしました。これは全世界 30m 解像度の標高データで、ダウンロード先はUSGS のサイトです。

デモなので富士山を覆うファイル ASTGTMV003_N35E138_dem.tif のみを使用しました。

上記リポジトリの手順に従い、複数の地理データ処理ソフトウェアを含む環境を用意する必要があります：

1
brew install gdal geoip libspatialite librasterlite spatialite-gui spatialite-tools

以降のデータ処理では多くの Python 用地理ツールを使用します。簡単のため、Conda で Python 3.11 の仮想環境を初期化しました:

1
conda create -n play-terrain-rgb python=3.11

依存関係をインストールしますが、オリジナルのリポジトリとは少し異なり、rasterioはソースからビルドしてインストールする必要があります:

1
pip install rasterio --no-binary rasterio

そうしないと、rioがrgbifyコマンドラインプラグインを見つけられません。

その後、他の依存関係をインストールできます:

1
pip install rasterio rio-rgbify rio-mbtiles mbutil

rio infoコマンドで元の DEM データの情報を確認できます:

1
rio info --indent 2 ASTGTMV003_N35E138_dem.tif

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
{
  "blockxsize": 256,
  "blockysize": 256,
  "bounds": [
    137.999861111111, 34.99986111111112, 139.00013888888876, 36.0001388888889
  ],
  "colorinterp": ["gray"],
  "compress": "lzw",
  "count": 1,
  "crs": "EPSG:4326",
  "descriptions": ["Band 1"],
  "driver": "GTiff",
  "dtype": "int16",
  "height": 3601,
  "indexes": [1],
  "interleave": "band",
  "lnglat": [138.4999999999999, 35.500000000000014],
  "mask_flags": [["all_valid"]],
  "nodata": null,
  "res": [0.000277777777777778, 0.000277777777777778],
  "shape": [3601, 3601],
  "tiled": true,
  "transform": [
    0.000277777777777778, 0.0, 137.999861111111, 0.0, -0.000277777777777778,
    36.0001388888889, 0.0, 0.0, 1.0
  ],
  "units": [null],
  "width": 3601
}

元データのカバー範囲が広すぎるため、富士山周辺だけを切り出す必要があります：

1
gdal_translate -projwin 138.626942 35.439672 138.899882 35.255449 ASTGTMV003_N35E138_dem.tif ASTGTMV003_N35E138_dem_clip.tif

その後、ファイルを Web Mercator 投影に変換します。これによりブラウザ上で表示できるようになります:

1
gdalwarp -t_srs EPSG:3857  -dstnodata None  -novshiftgrid -co TILED=YES  -co COMPRESS=DEFLATE  -co BIGTIFF=IF_NEEDED -r lanczos ASTGTMV003_N35E138_dem_clip.tif  ASTGTMV003_N35E138_dem_clip_EPSG3857.tif

1
2
Creating output file that is 914P x 756L.
Processing ASTGTMV003_N35E138_dem_clip.tif [1/1] : 0...10...20...30...40...50...60...70...80...90...100 - done.

変換後のファイル情報を確認します:

1
rio info --indent 2 ASTGTMV003_N35E138_dem_clip_EPSG3857.tif

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
{
  "blockxsize": 256,
  "blockysize": 256,
  "bounds": [
    15431865.404742576, 4198669.725048377, 15462269.950626153, 4223818.342868929
  ],
  "colorinterp": ["gray"],
  "compress": "deflate",
  "count": 1,
  "crs": "EPSG:3857",
  "descriptions": ["Band 1"],
  "driver": "GTiff",
  "dtype": "int16",
  "height": 756,
  "indexes": [1],
  "interleave": "band",
  "lnglat": [138.76336989692504, 35.347779733715655],
  "mask_flags": [["all_valid"]],
  "nodata": null,
  "res": [33.265367487502644, 33.265367487502644],
  "shape": [756, 914],
  "tiled": true,
  "transform": [
    33.265367487502644, 0.0, 15431865.404742576, 0.0, -33.265367487502644,
    4223818.342868929, 0.0, 0.0, 1.0
  ],
  "units": [null],
  "width": 914
}

その後、rio rgbifyコマンドで標高データを RGB 画像に変換します:

1
rio rgbify -b -10000  -i 0.001 ASTGTMV003_N35E138_dem_clip_EPSG3857.tif ASTGTMV003_N35E138_dem_clip_EPSG3857.RGB.tif

出力された画像は大まかに以下のようになります:

元のファイル内の特定の経緯度の標高を簡単に確認できます:

1
gdallocationinfo -wgs84 ASTGTMV003_N35E138_dem.tif 138.72739076614383 35.36067113569001

1
2
3
4
Report:
  Location: (2619P,2302L)
  Band 1:
    Value: 3756

変換後の RGB ファイルでのデータを確認します:

1
gdallocationinfo -wgs84 ASTGTMV003_N35E138_dem_clip_EPSG3857.RGB.tif 138.72739076614383 35.36067113569001

1
2
3
4
5
6
7
8
Report:
  Location: (336P,325L)
  Band 1:
    Value: 209
  Band 2:
    Value: 234
  Band 3:
    Value: 72

読み出された 3 つの値はそれぞれ 209、234、72 で、変換式は次の通りです：height = -10000 + ((R × 256 × 256 + G × 256 + B) × 0.001)（注: 0.0001 という係数はrio rgbifyの-iパラメータ値です）。これにより標高値を計算できます:

1
2
height = -10000 + ((209 * 256 * 256 + 234 * 256 + 72) * 0.001)
# 3757

誤差は許容範囲です。

次に、地図上に表示するために RGB 画像をタイルに変換します。ここではgdal2tiles.pyを使用しました:

1
time gdal2tiles.py --zoom=5-19 --processes=16 ASTGTMV003_N35E138_dem_clip_EPSG3857.RGB.tif ./tiles

1
2
3
4
5
Generating Base Tiles:
...10...20...30...40...50...60...70...80...90...100 - done.
Generating Overview Tiles:
0...10...20...30...40...50...60...70...80...90...100 - done.
gdal2tiles.py --zoom=5-19 --processes=16  ./tiles  1332.07s user 360.34s system 1297% cpu 2:10.39 total

PS: 複数回タイルを生成すると、gdal2tiles.pyが生成される HTML ファイル内のズーム範囲を変更するため、実際の状況に応じて手動で修正する必要があります。

最後に、mb-utilを使用してタイルを MBTiles 形式に変換できます:

1
mb-util --image_format=png --scheme=tms ./tiles/ ./ASTGTMV003_N35E138_dem_clip_EPSG3857.RGB.mbtiles

1
176000 tiles inserted (24076 tiles/sec)mb-util --image_format=png --scheme=tms ./tiles/   2.12s user 3.92s system 75% cpu 7.995 total

TLDR#

作成手順#

TLDR

作成手順