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

> **この記事は2024年1月13日に書かれました。**

## TLDR

数日前、資料を調べていると MapTiler が作成した[軌跡標高データ可視化ページ](https://labs.maptiler.com/samples/cloud/elevation-profile/)を見つけ、非常に良い出来でした。

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

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

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

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

![terrain-rgb-profile-fuji.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4073937/7b98c965-e858-4fe1-96a5-7c427a89efcc.jpeg)

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

## 作成手順

ASTER GDEM v3 のデータをダウンロードしました。これは全世界 30m 解像度の標高データで、ダウンロード先は[USGS のサイト](https://asterweb.jpl.nasa.gov/gdem.asp)です。

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

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

```bash
brew install gdal geoip libspatialite librasterlite spatialite-gui spatialite-tools
```

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

```bash
conda create -n play-terrain-rgb python=3.11
```

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

```bash
pip install rasterio --no-binary rasterio
```

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

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

```bash
pip install rasterio rio-rgbify rio-mbtiles mbutil
```

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

```bash
rio info --indent 2 ASTGTMV003_N35E138_dem.tif
```

```json
{
  "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
}
```

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

```bash
gdal_translate -projwin 138.626942 35.439672 138.899882 35.255449 ASTGTMV003_N35E138_dem.tif ASTGTMV003_N35E138_dem_clip.tif
```

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

```bash
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
```

```txt
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.
```

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

```bash
rio info --indent 2 ASTGTMV003_N35E138_dem_clip_EPSG3857.tif
```

```json
{
  "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 画像に変換します:

```bash
rio rgbify -b -10000  -i 0.001 ASTGTMV003_N35E138_dem_clip_EPSG3857.tif ASTGTMV003_N35E138_dem_clip_EPSG3857.RGB.tif
```

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

![terrain-rgb-profile-fuji-rgb.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4073937/23ec574b-e828-419d-9354-fef6e60cc4b5.jpeg)


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

```bash
gdallocationinfo -wgs84 ASTGTMV003_N35E138_dem.tif 138.72739076614383 35.36067113569001
```

```txt
Report:
  Location: (2619P,2302L)
  Band 1:
    Value: 3756
```

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

```bash
gdallocationinfo -wgs84 ASTGTMV003_N35E138_dem_clip_EPSG3857.RGB.tif 138.72739076614383 35.36067113569001
```

```txt
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`パラメータ値です）。これにより標高値を計算できます:

```python
height = -10000 + ((209 * 256 * 256 + 234 * 256 + 72) * 0.001)
# 3757
```

誤差は許容範囲です。

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

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

```txt
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 形式に変換できます:

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

```txt
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
```