はじめに
こんな感じのグローバルRGB標高タイルを作ったのでメモします。
今回の作業の目的
先日、国土地理院の標高タイル（小縮尺のもの）を使って MapLibre GL JS の Terrain 3D表示を試しました（記事はこちら）。地理院地図の標高タイルはMapLibreやMapBoxのRGB標高タイルと仕様が違うので、今回はMapLibreやMapboxの仕様に準拠した小縮尺の標高RGBタイルを自分で作ります。
この記事ではMapLibreでもMapBoxでも使えるように自分でRGB標高タイルを作ることを目的にします。具体的には、地理院地図の標高タイル（小縮尺ZL 0-8)でも使っている地球地図第二版標高データをダウンロードしてきて、オープンソースのツールを使って自分でRGB標高タイルを作ってみることにします。上手く作れれば、自分で使っているグローバルなベクタータイル地図のベースとして今後使える可能性があると考えています。
なお、最近、@Kanahiro(MIERUNE)さんの記事（こちら）で、MapLibre GL JSの addProtocol() を使うことで地理院の標高タイルをそのままMapLibreで読めるという素晴らしい技があることを知りました。MapLibreを使う人はその方法が一番良いかもしれませんのでメモしておきます。
利用したデータ（地球地図第二版標高データ）について
ソースのURL （Global Map data archive）
地球地図国際運営委員会（ISCGM）は解散しましたが、以下に標高のデータアーカイブが残っています。
https://github.com/globalmaps/gm_el_v2_west
https://github.com/globalmaps/gm_el_v2_east
地球地図第二版標高データ Global Map version 2 Elevation の仕様
- 16 bit signed integer (Geotiffになっている)
- 水部は -9999
- No Dataは 9998
- 空間解像度： 15 arc seconds
- 86,400 by 43,200 pix (whole globe)
- ファイルは 48 のファイルに分かれて配布されている。
- 4 by 12 = 48 area (1-1 to 4-12)
（QGISで開いてみたところ）
上のデータは以下の様な感じで48個のファイルに分かれています。
作業で使うツール
紆余曲折があったのですが、unvt/rgbifyのDockerfileを使って作業することにします。これは、osgeo/gdal:ubuntu-small-3.4.0ベースで、python3やmapboxのrio-rgbifyを追加したイメージです。そして最後のところで、unvt/rgbifyにmapubox/mbutilを加えて作業しました。
これまでの経緯１
最初に、Water GISでやられている方法「DEMデータをTerrain RGBのラスタタイルに変換し、Mapbox GL JSで利用する」を試してみようと思いました。@JinIgarashi さんが鋭意進められていて、コードなども勉強になるからです。しかし、少し挑戦したところ、私の職場のパソコンに入っているPython3が少し古くて先に進めません。PCの管理者権限がなくてソフトのアップグレードがすぐにできないこともあり断念。
コードを見習って、自分のpythonとgdalの環境で、自力でやろうとするもpip install gdalが通らずに出来ませんでした。（gdalinfoで確認したバージョンと揃えてインストールしたのですがダメでした。）
これまでの経緯2
私のWindows環境だとPythonでgdalを使う環境を安定して構築できないので、DockerでUbuntuをベースにpythonとgdalをインストールすればいいのではないかと思い環境構築に取り組みました。しかし、Anacondaを入れてからgdalを入れようとしたためDockerfileがうまく作れない。（Anacondaのインストールの途中ではyesを入れないといけないのでDockerfileの記載だけではインストールできなかった。）
Dockerコンテナに手動でAnacondaをインストールしてもpip install gdalが通らなくて断念。
ダメだったレポジトリ https://github.com/ubukawa/satsuki
これまでの経緯3
そこで、そもそもoptgeoで標高タイルを作っていたよなぁ、と思い出し、見てみたら南極の標高タイル作成で unvt/rgbify を使っているのを思い出す（ https://github.com/optgeo/showa-terrain ）。@hfuさんが記事（標高タイルを作る rio-rgbify を Docker で動かす）まで書いてくれているので助かりました。（初めから思い出せばよかったですが、いい経験になりました）
作業環境
- Windows 10
- Powershell
- Docker （version 20.10.8）
作業手順
Step 0. データ作成範囲（ズームレベル）を考える
作業というか、作業の前の検討ですが、作成するズームレベルを考えました。
標高RGBタイルはラスタタイルになります。512×512ではなくて、256×256でタイルを準備しようと思います。ソースデータの解像度と各ズームレベルで必要なピクセル数を考えて、今回の作業対象とするズームレベルは８までとしました（理由は以下）。
- 今回使うソースデータ（地球地図第二版標高）の解像度は15秒です。
- 1分に4ピクセル、1度に240ピクセル、360度に86,400ピクセル並ぶということです。
- 緯度方向（-90度～90度）、経度方向（-180～180度）を考えると、43,200 × 86,400ピクセルになります。
- 上記のソースファイルの解像度を考えると、妥当なズームレベルは8くらいかなと思います。（南北方向では、ソースデータのピクセル数の方が少なくなるので、出来上がりをみて難しければズームレベル7までかなと思います。）
- 国土地理院でも地球地図から作ったタイルのカバー範囲はZL8までです。
- 将来、例えばALOS全球数値地表モデル (DSM) "ALOS World 3D - 30m (AW3D30)"などを使えば、ズームレベル11くらいまでいけるのかなと思います。
- 赤道が40,000kmとすると、赤道上では、30メートルのピクセルが1,333,333個並ぶ。（緯度60度ではその1/2で666,667個並ぶ）
- 北極から南極までの子午線が20,000kmとすると、30メートルのピクセルが666,667個並ぶ
（参考）表. ズームレベルと全世界のピクセル数（256ピクセルタイルの場合）
| ズームレベル | タイル数 | 全世界のピクセル数 | 備考 |
|---|---|---|---|
| 0 | 1 by 1 | 256 by 256 | |
| 1 | 2 by 2 | 512 by 512 | |
| 2 | 4 by 4 | 1,024 by 1,024 | |
| 3 | 8 by 8 | 2,048 by 2,048 | |
| 4 | 16 by 16 | 4,096 by 4,096 | |
| 5 | 32 by 32 | 8,192 by 8,192 | |
| 6 | 64 by 64 | 16,384 by 16,384 | |
| 7 | 128 by 128 | 32,768 by 32,768 | |
| 8 | 256 by 256 | 65,536 by 65,536 | 地球地図第二版標高での最大ズーム |
| 9 | 512 by 512 | 131,072 by 131,072 | |
| 10 | 1,024 by 1,024 | 262,144 by 262,144 | |
| 11 | 2,048 by 2,048 | 524,288 by 524,288 | 30メートル解像度のDEMならこのくらいはいけるのでは |
| 12 | 4,096 by 4,096 | 1,048,576 by 1,048,576 | |
| 13 | 8,192 by 8,192 | 2,097,152 by 2,097,152 |
Step 1. Dockerコンテナの起動（unvt/rgbify）
unvt/rgbifyはDocker Hubにはアップされていないようなので、GitHubレポジトリをクローンしてきて、そこのDockerfileからビルドします。
git clone https://github.com/unvt/rgbify
cd rgbify
docker build -t unvt/rgbify .
docker run -it --rm -v ${PWD}:/data unvt/rgbify
Step 2. データのダウンロードとZIP解凍
作業フォルダを作ります。そして、Powershellを立ち上げて、ソースデータをダウンロードしてくる。以下の二つのレポジトリに保存されているので、それぞれgit cloneしてくる。
git clone https://github.com/globalmaps/gm_el_v2_east
git clone https://github.com/globalmaps/gm_el_v2_west
ダウンロードしてきたzipファイルを適当なフォルダに集めて解凍する。コマンドを使って一度にunzipできるから助かりますね（以下の通り）。48個のファイルができます。srcというフォルダに入れておきました。
for f in *.zip; do unzip ${f}; done
Step 3. データをみて、作業順を考える
ダウンロードしてきたデータは10進法経緯度座標で、48ファイルに分かれているデータです。また、水部やNoData値について、少し処理が必要になりそうです。
今回は、
ラスタの計算(gdal_calc.py) → データのマージ(gdal_merge.py) → 投影変換(gdalwarp) → RGB画像タイル化(rasterio rgbify)
という手順でやってみようと思います。考えたことは以下の様なことでした。
- 今回のDEMからRGB画像タイル化するために必要な作業は、値の計算（NoDataや水部の値の処理）、ファイルの結合、投影法変換（EPSG4326 → EPSG3857）、タイル化の各ステップになります。（DEM値からRGB標高タイルへの値の変換についてはrasterio-rgbifyでやってしまうので自分で計算しなくてよさそうです。）
- 値の計算(gdal_calc.py)は、全球一括でなくて、元データの区画で複数回実施（小さいファイルサイズで計算）したほうが安定するのではないかと推測。
- RGB画像タイル化（rasterio rgbify）の前の画像は10進法経緯度ではダメなのか実験してみたが、どうもウェブメルカトルにしないとダメらしい（あるいはデータが90度に合ったからエラーだったのかもしれません）。
今回の作業ではRGB画像タイル作成の前に投影変換をすることにします。投影変換(gdal_warp)では、画像のピクセル数も変わるのでそれによる劣化が気になる。解像度もあわせないといけないので、このプロセスは変換の直前にやりたい。
Step 4. ラスタ計算・下処理（gdal_calc.py）
ソースデータは元々の仕様によるのですが、水部とNo dataの値が-9999と9998なのでここを処理します。No dataを与えようかなとも思いましたが、今回の作業ではとりあえず0にしておくことにしました。
ラスタ計算はgrassなどを使おうかなぁ、、とも少し考えましたが、今回の作業プラットフォーム unvt/rgbify にはgdal_calcも入っているのでこれを使ってやろうと思います。下記の2つの方法を考えましたが、一回の計算でできるので、option 1の計算で行うことにしました。
gdal_calc.py -A src/gm_el_v2_01_01.tif --outfile calc-out.tif --calc="A*(A!=-9999)*(A!=9998)"
gdal_calc.py -A src/gm_el_v2_01_01.tif --outfile calc-int.tif --calc="where(A==-9999,0,A)"
gdal_calc.py -A calc-int.tif --outfile calc-out.tif --calc="where(A==9998,0,A)"
48ファイルを一度にやるとこんな（↓）感じです
for f in src/*.tif; do gdal_calc.py -A src/`basename ${f}` --outfile 01_calc/`basename ${f}` --calc="A*(A!=-9999)*(A!=9998)"; done
これで水部とNo dataをゼロにしたファイルができました。本当はNodataの値を上手く扱いたいので、そこは今後検討を深める部分です。
Step 5. 画像ファイルをマージ（gdal_merge.py）
前のステップで計算した画像ファイルのサイズが一つあたりだいたい152MBになっていました。48ファイルを一つにまとめると7GBを超えるので、投影変換やタイル化が少し厳しいかなぁと想像しました。
そこで、6ファイルずつマージして、以下の様な8つの区画ごとにマージすることにしました。各画像からタイルを別々に作ってもズームレベル2以上であればタイルの重複はありません。
マージ作業はgdal_merge.pyでできました。これも今回の作業プラットフォーム unvt/rgbify に入っています。以下のスクリプトで実施しましたが、一つのファイル作成に数分かかりました。
gdal_merge.py -o 02_merge/merge01.tif 01_calc/gm_el_v2_01_01.tif 01_calc/gm_el_v2_01_02.tif 01_calc/gm_el_v2_01_03.tif 01_calc/gm_el_v2_02_01.tif 01_calc/gm_el_v2_02_02.tif 01_calc/gm_el_v2_02_03.tif;
gdal_merge.py -o 02_merge/merge02.tif 01_calc/gm_el_v2_03_01.tif 01_calc/gm_el_v2_03_02.tif 01_calc/gm_el_v2_03_03.tif 01_calc/gm_el_v2_04_01.tif 01_calc/gm_el_v2_04_02.tif 01_calc/gm_el_v2_04_03.tif;
gdal_merge.py -o 02_merge/merge03.tif 01_calc/gm_el_v2_01_04.tif 01_calc/gm_el_v2_01_05.tif 01_calc/gm_el_v2_01_06.tif 01_calc/gm_el_v2_02_04.tif 01_calc/gm_el_v2_02_05.tif 01_calc/gm_el_v2_02_06.tif;
gdal_merge.py -o 02_merge/merge04.tif 01_calc/gm_el_v2_03_04.tif 01_calc/gm_el_v2_03_05.tif 01_calc/gm_el_v2_03_06.tif 01_calc/gm_el_v2_04_04.tif 01_calc/gm_el_v2_04_05.tif 01_calc/gm_el_v2_04_06.tif;
gdal_merge.py -o 02_merge/merge05.tif 01_calc/gm_el_v2_01_07.tif 01_calc/gm_el_v2_01_08.tif 01_calc/gm_el_v2_01_09.tif 01_calc/gm_el_v2_02_07.tif 01_calc/gm_el_v2_02_08.tif 01_calc/gm_el_v2_02_09.tif;
gdal_merge.py -o 02_merge/merge06.tif 01_calc/gm_el_v2_03_07.tif 01_calc/gm_el_v2_03_08.tif 01_calc/gm_el_v2_03_09.tif 01_calc/gm_el_v2_04_07.tif 01_calc/gm_el_v2_04_08.tif 01_calc/gm_el_v2_04_09.tif;
gdal_merge.py -o 02_merge/merge07.tif 01_calc/gm_el_v2_01_10.tif 01_calc/gm_el_v2_01_11.tif 01_calc/gm_el_v2_01_12.tif 01_calc/gm_el_v2_02_10.tif 01_calc/gm_el_v2_02_11.tif 01_calc/gm_el_v2_02_12.tif;
gdal_merge.py -o 02_merge/merge08.tif 01_calc/gm_el_v2_03_10.tif 01_calc/gm_el_v2_03_11.tif 01_calc/gm_el_v2_03_12.tif 01_calc/gm_el_v2_04_10.tif 01_calc/gm_el_v2_04_11.tif 01_calc/gm_el_v2_04_12.tif;
Step 6. 投影法変換（gdalwarp）
WGS84の10進法経緯度（EPSG:4326）のデータなので、gdalwarpを使ってウェブメルカトル投影（EPSG:3857）に変換します。（もしかしてDD座標からでもそもまま変換出来るのかもしれませんが。）
gdalwarpでインプットとアプトプットの投影法を指定すればいいのですが、十分なピクセル数を確保するため解像度だけは自分で指定しておきます。今扱っている全球8分割画像はそれぞれ21,600×21,600ピクセルです。解像度の指定をしないと9,418×29,059のサイズになってしまうので、若干情報が失われるかなと思いました。解像度300mで 33,396×10,3043、500mで20,038×61,826、400mで25,047×77,283位になります。そこで、今回は解像度350mで 28,625 × 88,323 サイズの画像を得ることにしました。（ウェブメルカトルの投影方向では南北方向には一様に引き延ばされるわけではないので、元画像よりピクセル数が多ければいいという単純な問題ではない点には注意が必要）。また、リサンプルの方法はnearest neighbourにして実行しました。
for f in 02_merge/*.tif; do gdalwarp -s_srs EPSG:4326 -t_srs EPSG:3857 -r near -tr 350 350 02_merge/`basename ${f}` 03_projected/`basename ${f}`; done
実際に作業中の画像は以下の様な感じです。解像度を細かくしすぎたのか、一つのファイルが4 GBを超えてしまい、一つの投影変換に結構時間がかかりました。無事にファイルを得られたので次のステップに進みます。
Step 7. RGB画像タイル化
このステップでは、mapbox/rio-rgbify というツールを使います。ライセンスはMITライセンスです。今回の作業プラットフォーム unvt/rgbify は特にこのrio-rgbifyを使うために作られたものだと思いますので、私の作業では特に気にすることなくRGB画像タイル化を進めます。
7-1 画像タイル化（mbtilesフォーマット）
MapLibre/Mapbox様の標高RGBタイルはbaseが-10000、ステップは0.1なので、それを踏まえて以下のとおり処理しました。フォーマットをpngで出力しようとしたのですが、mbtilesだけのようだったので、まずはmbtilesで出力するようにしています。
for f in 03_projected/*.tif; do rasterio rgbify -b -10000 -i 0.1 --max-z 8 --min-z 2 --format webp 03_projected/`basename ${f}` tile/`basename ${f} .tif`.mbtiles; done
実際の画面は下の通りです。何やら警告がでているのですが、mbtilesは作成されてきています。
7-2 mbtilesからpngタイルへの変換
昭和基地のRGB画像作成では、sqlite3を使ってmbtilesを画像ファイルに展開していますが、私は既存のツールを使うことにしました。mapboxさんのmbutilをgit cloneして、インストールして使います。
git clone からインストールまで：
インストールした後に、1つめのファイルをzxy画像タイルに展開：
上のものと同じ方法で全てのファイルを展開しました。フォルダを覗くと画像タイルは出来ているようです。
全部のタイル、ZL2-8でだいたい2GBくらいのサイズでした。
7-3 画像タイルをまとめる
各地域のタイルを一つのフォルダに移動しているとき、一部の地域で重複がありました。区画の切り方は重複がないように作ったのですが、画像の端で多少の誤差があるようです。解像度の影響か、元データの影響か、微小ピクセルが隣の範囲にはみ出すことで、本来の領域からはみ出て、ほとんど空の画像タイルがタイル一区画余分に作られているようです。
例えば、南北の方向では、地区の境目（赤道）のすぐ南のタイルで重複があるようです。例えば、2/0/2, 3/x/4, 4/x/8, 5/x/16, 6/x/32, 7/x/64のタイルが重複していました（北半球の側の地区のタイルが一つ南半球の領域にはみ出している）。北半球側はほぼ空タイルなので、南側を優先すれば良いことがわかります。
-
南北方向（例：1と2、3と4、5と6、7と8。）。
- 赤道のすぐ南のタイルで重複があるようです。例えば、2/0/2, 3/x/4, 4/x/8, 5/x/16, 6/x/32, 7/x/64のタイルが重複していました（北半球の地区でタイルが余分につくられて南半球にはみ出ている。）
- これについては、南半球側の地区（偶数の地域）を優先すればよいです。
-
東西方向
- 地区1と地区7、地区2と地区8、地区5と地区7、地区6と地区8では重複がありました。1と3、2と4、3と5、4と6では重複がありませんでした。
- 地区1と地区7、地区2と地区8では、境界のすぐ東側で重複があります。これについては東側の地区（地区1、地区2）を優先すればよいです。
- 不思議なことに、地区5と地区7、地区6と地区8では境界のすぐ西側で重複があります。これについては、西側の地区（地区5と地区6）を優先するとよいです。
これらの観察結果を踏まえて、順番に気をつけて注意深く各地域のタイルをあわせていきます。
Step 8. 画像タイルのアップロード
これまでの作業で得られた画像タイルを整理しました。そして、GitHubレポジトリにアップロード出来るサイズZL2-7をアップロードしました。（ https://ubukawa.github.io/globalmap-el/rgb-el_2-7/{z}/{x}/{y}.png ）
Step 9. 見てみる
簡単な地図サイトを作ってみてみました。タイルはなんとか出来ていそうです。（ベクトル地図はNaturalEarthのもの）
※ 小縮尺だと、高さ方向は強調しないとわかりにくいですね。
https://ubukawa.github.io/globalmap-el/map.html
謝辞
今回の取組みを進めるにあたり、@hfuさん、@JinIgarashiさん、そして@Kanahiro(MIERUNE)さんの記事などに大変助けられました。いつもありがとうございます。
まとめ
今回、地理院地図の標高タイルでも使われている（ZL0-8）地球地図第二版標高データを使って、Mapbox GL JS/MapLibre GL JSでそのまま使えるRGB標高タイルを作成しました。
私のローカルな環境にはズームレベルが2～8までおいてありますが、GitHubレポジトリにはズームレベル2～7までを載せておきました。
今後の課題
- No dataの値については今後確認します。
- もう少し解像度の高いDEMを使って、ZL9～11あたりをどうするか今後検討出来たらしていきたい。
追記
2022.5.27 rgbifyの入力データの投影法について
ALOSのDEMを変換していたら、rasterio rgbifyは経緯度座標のラスタファイルからも変換ができました。地球地図データで失敗したのは85度以北・以南が入っていたからかもしれません。gdalwarpをやる代わりに、ラスタデータをクリッピングすればうまくいくかもしれません。