この記事は2024年1月13日に書かれました。
TLDR
数日前、資料を調べていると MapTiler が作成した軌跡標高データ可視化ページを見つけ、非常に良い出来でした。
公式のブログで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 のみを使用しました。
上記リポジトリの手順に従い、複数の地理データ処理ソフトウェアを含む環境を用意する必要があります：
brew install gdal geoip libspatialite librasterlite spatialite-gui spatialite-tools
以降のデータ処理では多くの Python 用地理ツールを使用します。簡単のため、Conda で Python 3.11 の仮想環境を初期化しました:
conda create -n play-terrain-rgb python=3.11
依存関係をインストールしますが、オリジナルのリポジトリとは少し異なり、rasterioはソースからビルドしてインストールする必要があります:
pip install rasterio --no-binary rasterio
そうしないと、rioがrgbifyコマンドラインプラグインを見つけられません。
その後、他の依存関係をインストールできます:
pip install rasterio rio-rgbify rio-mbtiles mbutil
rio infoコマンドで元の DEM データの情報を確認できます:
rio info --indent 2 ASTGTMV003_N35E138_dem.tif
{
"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
}
元データのカバー範囲が広すぎるため、富士山周辺だけを切り出す必要があります：
gdal_translate -projwin 138.626942 35.439672 138.899882 35.255449 ASTGTMV003_N35E138_dem.tif ASTGTMV003_N35E138_dem_clip.tif
その後、ファイルを Web Mercator 投影に変換します。これによりブラウザ上で表示できるようになります:
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
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.
変換後のファイル情報を確認します:
rio info --indent 2 ASTGTMV003_N35E138_dem_clip_EPSG3857.tif
{
"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 画像に変換します:
rio rgbify -b -10000 -i 0.001 ASTGTMV003_N35E138_dem_clip_EPSG3857.tif ASTGTMV003_N35E138_dem_clip_EPSG3857.RGB.tif
出力された画像は大まかに以下のようになります:
元のファイル内の特定の経緯度の標高を簡単に確認できます:
gdallocationinfo -wgs84 ASTGTMV003_N35E138_dem.tif 138.72739076614383 35.36067113569001
Report:
Location: (2619P,2302L)
Band 1:
Value: 3756
変換後の RGB ファイルでのデータを確認します:
gdallocationinfo -wgs84 ASTGTMV003_N35E138_dem_clip_EPSG3857.RGB.tif 138.72739076614383 35.36067113569001
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パラメータ値です）。これにより標高値を計算できます:
height = -10000 + ((209 * 256 * 256 + 234 * 256 + 72) * 0.001)
# 3757
誤差は許容範囲です。
次に、地図上に表示するために RGB 画像をタイルに変換します。ここではgdal2tiles.pyを使用しました:
time gdal2tiles.py --zoom=5-19 --processes=16 ASTGTMV003_N35E138_dem_clip_EPSG3857.RGB.tif ./tiles
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 形式に変換できます:
mb-util --image_format=png --scheme=tms ./tiles/ ./ASTGTMV003_N35E138_dem_clip_EPSG3857.RGB.mbtiles
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