はじめに

今回の記事の趣旨

前回、地球地図第二版標高（解像度15秒）から、ズームレベル2～8までの標高RGBタイルを作りました（記事こちら）。
今回はもう少し大きい縮尺で標高タイルを作るため、JAXAさんから無料で公開されているALOS World 3D - 30m (AW3D30)（注：解像度は1秒）をつかって標高RGBタイルを作成しました。ズームレベル9～11まで、範囲は全球ではなく東経5度～15度、北緯35度から45度の範囲を試験的に行いました）。

その作業の経緯をメモします。

2022.5.29 追記
この記事では、どうやったらAW3D30からRGB標高DEMを作ることが出来るか書いていますが、テンプレート（Docker）を使った簡単なRGB標高タイル作成方法も以下にまとめています。
https://qiita.com/T-ubu/items/dc5022e15f2b770429b3

ALOS World 3D - 30m (AW3D30)　について

JAXAさんがALOS/PRISMデータから作成したグローバルDEMで、30メートル解像度のものです。無料で公開されています。最新のものは第3.2版ですが、2016年4月の第1版公開以降、着実に公開範囲を広げ、またデータの改善を行ってきているようです。

本データセットは、陸域観測技術衛星「ALOS」搭載の光学センサ・パンクロマチック立体視センサ（PRISM）による解像度30m相当（緯度経度1秒間隔を基本）の全球数値地表モデル（DSM）です。最新版は第3.2版と第3.1版です。本データセットは、下記「5. 利用規約」の条件において、商用・非商用目的を問わずどなたでも無償でご利用いただけます。科学研究や教育分野、地理空間情報を活用した新しいサービス等に、本データセットを広くご活用いただければ幸いです。
（JAXAさんのホームページより　https://www.eorc.jaxa.jp/ALOS/jp/dataset/aw3d30/aw3d30_j.htm　）

手順

Step 1: データのダウンロード、データの確認

ユーザー登録が必要ですが、以下のページからリンクをたどってダウンロードしました。（以下のページは英語ページですが、日本語ページもあります。）
https://www.eorc.jaxa.jp/ALOS/en/dataset/aw3d30/aw3d30_e.htm

1度×1度のデータ、あるいは5度×5度のデータパッケージをHPからダウンロードできます。
今回はイタリア周辺の6つの5度×5度の区画をダウンロードしました。東経5度～20度、北緯40度から50度の範囲をカバーします。

ダウンロードしてきた5度×5度の区画のZIPファイルを見てみます。この中には、1度×1度の地域ごとにデータ（DSM）やマスク情報などが入っていました。水部など、データがない所もあります。解像度は1秒なので、1度×1度の画像のピクセル数は 3，600 × 3,600 で、サイズは25,342kbでした。

マスク情報も参考になります。水部のマスクの部分の標高は0ということを確認しました。 なお、仕様書は二進法での表記ですが、QGISだと10進法の表現になっているので、気をつけましょう。

技術情報も公開されています。
https://www.eorc.jaxa.jp/ALOS/en/dataset/aw3d30/data/aw3d30v3.2_product_e_e1.2.pdf　
（英語版ですが、日本語版もあると思います。）

Step 2: Docker スタート

今回もRGB標高タイルを作るためにunvt/rgbifyを使います。ただし、mapbox/mbutilを入れるために、私のGitHubアカウントでunvt/rgbifyをフォークして、mbutilをインストールしたイメージ（ubukawa/rgbify:ubukawa）を作成しています。

以下の様に、ubukawa/rgbifyにあるDockerfileのイメージを作って、今回の作業フォルダでDockerコンテナをスタートします。

git clone https://github.com/ubukawa/rgbify
cd rgbify
docker build -t unvt/rgbify:ubukawa .
cd ..
docker run -it --rm -v ${PWD}:/data unvt/rgbify:ubukawa
cd /data

Step 3: ZIPファイルの解凍

srcというフォルダをつくって、ダウンロードしてきた5度×5度の地域ごとのZIPファイルを入れておきます。（今回の場合は6つのZIPファイル）。

docker コンテナをスタートしており、Linuxのコマンドが使えるので、以下のように解凍用フォルダ01_unzipを準備してから、解凍します。

mkdir 01_unzip
for f in src/*.zip; do unzip ${f} -d 01_unzip; done

↓こんな感じになります

Step 4：　画像ファイルのマージ、標高タイル作成

解凍した5度×5度のフォルダには、1度×1度のファイルがたくさん入っています。DSMデータはファイル名の終わりがDSM.tifになっていますので、これを全てマージしてしまうことにします。タイル作成のインプット用のtifなので、とりあえずinput.tifという名前で作ることにします。

入力ファイルの作成（マージ）

gdal_merge.py -o input.tif 01_unzip/*/*_DSM.tif

（メモ）
今回は、5度×5度の区画6枚分なので、1度×1度のファイル（約25MB）が25×6 = 150枚相当あることになります（水部はtifファイルがないものがあるのですが、マージするとそこもサイズに含まれる）。サイズはおおまかに見積もって4GB前後だろうと想像しましたが、このくらいならなんとか処理できると考えました。
※実際にできたinput.tifのサイズは約3.8GBでしたし、それを変換したmbtilesファイルはZL9-11で879MBでした。

入力用ファイルを作ってしまえば、あとは変換するだけです。どうやら10進法の経緯度座標でもそのまま変換できるので、gdalwarpなどでの投影変換は考えなくてよいです。

RGB標高タイル作成

rasterio rgbify -b -10000 -i 0.1 --max-z 11 --min-z 9 --format webp input.tif out.mbtiles

↓実際の画面（少し時間がかかります）

上記はダウンロードしてきたファイルを全てまとめてから変換する方法ですが、例えば、まずは5度×5度の区画ごとにマージしたければ以下のような感じでいいと思います。

mkdir merge
for f in 01_unzip/*; do gdal_merge.py -o merge/`basename ${f}`.tif ${f}/*_DSM.tif; done

そのほかにも、5度×5度の区画ごとに中身のDSMファイルをコピーしてきたいときは以下の様な感じでできます。

for f in 01_unzip/*;do mkdir 02_extract/`basename ${f}`; for m in ${f}/*_DSM.tif; do cp ${m} 02_extract/`basename ${f}`/  ;done  ; done

データの数が多いとサイズが大きくなりますので、あらかじめサイズを検討してからマージをするとよいと思います。

Step 5: mbtiles から pngタイルへの展開

簡単にこのプロセスをおこなうために、mbutilを入れてあります。以下のコマンドでmbtilesをpngタイルに展開します。

mb-util out.mbtiles zxy

Step 6: タイルの確認

フォルダを覗いてみると出てきたPNGファイルに中身があります。z/x/yの位置も正しいところにきています。

出来たタイルは、テストの目的でGitHubレポジトリにアップしました。MapLibreで簡単なWeb地図にしてみてみました。
下図は、前回作った地球地図のDEMの上にALOSのDEMを重ねたものです。ALOSのDEMの作成範囲の南端なので背景にある地球地図標高と比較ができると思います。15秒解像度の地球地図と、1秒解像度のALOSでは解像度の違いがよくわかりますね。（Hillshadeはきちんとでているのですが、Terrainの感じがあまり上手く出ていないような気もするのでもう少し確認します。）
https://ubukawa.github.io/globalmap-el/map-withALOS.html#9.17/39.9657/15.7677/0/29

謝辞

AW3D30 (ALOS World 3D - 30m) を作成・公開してくれているJAXAさんに感謝しています。unvt/rgbifyの作成者のhfuさんに感謝しています。

まとめ

AW3D30 (ALOS World 3D - 30m) をRGB標高タイルにする方法をためしました。
少し時間はかかりますが、手順的には結構簡単に作成できるように思いました。

References