無料のSRTM（1秒）標高データからRGB標高タイルを作る(docker, rasterio, gdal_merge)

先日、以下の記事を書いたのですが、あんがい便利ですので、日本語でも書いておこうと思います。

https://qiita.com/T-ubu/items/09720f25d5278d85222c

:::note 
（2022.10.31 追記）
このメモはSRTMデータ用ですが、SRTMは西アルプスなど山岳地帯でデータ欠損が目立つ場合があります。その場合、SRTMをASTER G-DEMやALOS/PRISIM DEMで補完した**NASADEM**の方がよいことがあります。（EarthExplorerでもダウンどーどできます。）
　　　参考　https://qiita.com/T-ubu/items/c7bb5697b51f245ef163
:::

# はじめに
3D Terrain Mapを作るために、MapLibre GL JS にしろ Mapbox GL JS にしろ、RGB標高タイルを使っています。世の中には、既存のRGB標高タイルサービスもありますが、ここでは自分でRGB標高タイルを作る方法を紹介します。

USGSからパブリックドメインとして公開されているSRTMデータや、JAXAさんが公開しているALOSの30メートルデータを使えば、**無料で自由に使えるRGB標高タイル**を作ることができます。

# 私の環境
* Windows 10 Enterprise
* Docker version 20.10.8
* PowerShell version 5.1.19041.1237

# 手順
## Step 1: USGSからSRTM データのダウンロード

[USGS EarthExplorer](https://earthexplorer.usgs.gov/) からデータをダウンロードします。今回は SRTM の 1 arc second DEM を tiff 形式でダウンロードしました。

今回の対象エリアは、以下の図のように10シーンでカバーされます。西経80－75度、北緯42－44度の地域です。ナイアガラの滝やフィンガー・レイクがあるところです。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/679143/f9d80e54-cbb2-f8d4-d60b-7db41f591879.png)

そして、ダウンロードしてきたTIFFファイルを、作業フォルダにsrcというフォルダを作って保存しておきます。（なお、パソコンの処理能力にもよりますが、私の環境では40ファイルを超えるくらい（マージして1GBを超えるくらい）になると一括での処理が難しくなるので、ほかの方法が必要になりました。）
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/679143/6fe973ca-af70-2be9-382b-59baab11a179.png)


参考ですが、もっとたくさん、一括でダウンロードしたい場合にはこちらの記事が参考になるかもしれません。 https://qiita.com/T-ubu/items/bb353935085eec2ab7a0

## Step 2: Dockerコンテナの起動
この作業中、いくつかの地理空間情報ツール（例えばgdal_merge.py、rasterio、mb-utilなど）を使います。個別ツールをインストールする手間を避ける意味もあり、ここでは既存のDockerイメージからコンテナを作って作業します。**unvt/rgbify:ubukawa**があるので、プルしてきて使えばよいと思います。以下のようにDockerを始めます。作業フォルダをDocker コンテナ中の/dataとしてつなげておくことも忘れないようにします。


```bash
docker pull unvt/rgbify:ubukawa
docker run -it --rm -v ${PWD}:/data unvt/rgbify:ubukawa
cd /data
```

## Step 3: DEMファイルをマージしてインプット用のファイルにする
もしファイルサイズが大きすぎると分割した方がよいですが、ファイルサイズが適切なら（後続の処理でトラブルがでないなら）、全部マージしてしまって一つのファイルにします。gdal_merge.py を使って以下のようにできます。（gdal_pyはdockerコンテナに入っています。）

```bash
gdal_merge.py -o input.tif src/*.tif
```
マージのプログレスが見られますね。↓
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/679143/fdf935ac-9276-2ecb-ecb4-ddf31bbd536f.png)

(Note for Advanced users)
From my experience, if the total size of your SRTM data is greater than a few GB, it would be a good idea to prepare separate input files. When I worked for the whole SRTM data, I made 932 inputs files from 14,277 SRTM files (I used the ZL6 tile extent to devide the input). If you are interested in, use translator and check this text: https://qiita.com/T-ubu/items/a90c3d42f809725f50ab.


## Step 4: RGB標高タイル作成
この次に、rasterioというツールを使ってRGB標高タイルを作ります。以下のコマンドを実行すればよいです。ズームレベル（max zoom / min zoom）は適宜調整して大丈夫ですが、オリジナルデータの解像度が1秒（≒30メートル）と考えると最大ズームは12とか13でよいと思います。（当たり前ですが、最大ズームを大きくするとデータが大きくなります。）

```bash
rasterio rgbify -b -10000 -i 0.1 --max-z 12 --min-z 9 --format webp input.tif out.mbtiles
mb-util out.mbtiles zxy
```

上のコマンドでは、まずはmbtilesを作って、その次にそれをzxy構造をもったpngファイルにmb-utilを使って変換しています。これが終わると、下の図のようなフォルダ構造をもったデータが出てきます。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/679143/812d6dd3-04df-58e4-4bbb-1231a9e2c3df.png)

これで、RGB標高タイルができました。おめでとうございます。
（あとはページにホストして使うだけ。）

# まとめ
ここでは、どのようにRGB標高タイルをSRTMデータから作るか書きました。

なお、この方法はALOSのdem（30メートルは無料です）にも使うことができると思います。

# 謝辞
素晴らしいSRTMデータを公開しているUSGSに感謝します。
また、ここで使っているツールの開発者にも感謝します。

# 参考
https://earthexplorer.usgs.gov/

https://github.com/mapbox/rio-rgbify

https://gdal.org/programs/gdal_merge.html

https://github.com/ubukawa/rgbify

https://github.com/unvt/rgbify