事象発生日：2024-12-22
記事公開日：2024-12-22
アクセス数：859
オープンデータである「ALOS 全球数値地表モデル」を利用して，全球の 30 m 解像度の地形データ (DSM) を pmtiles でタイル化して MapLibre GL JS で表示できるようにしました．
というだけだととても簡単なのですが，全球のデータを加工するのはなかなか骨の折れる作業でした．
また，苦労しただけ様々な学びが得られたので，それらをまとめておきたいと思います．
これは FOSS4G Advent Calendar 2024 の12月22日の記事です．
そこでは，事業 PoC のために Web GIS 技術を用いたアプリケーションのモックアップを開発し，プロトタイピングを重ねています．
ありがたいことに，2024 年の今日，様々な FOSS4G (Free Open Source Software for GeoSpatial) のお陰で，短期間のうちにいくつかのモックアップを作り上げることができました．
私の所属している株式会社アークエッジ・スペース (ArkEdge Space Inc.) では，OSS 活動に力を入れており，たとえば衛星のフライトソフトウェアや地上システム，様々な社内ツールなどを積極的に OSS として公開しています（一例として のような取り組みなど）．
FOSS4G の強力な恩恵を実感した以上，なにかしらこのコミュニティにも貢献したいと考え，今年はじめて FOSS4G 2024 Japan に参加，また，仕事で南米に行くついでにブラジル ベレンで開催される FOSS4G 2024 Belém Brasil にスポンサー参加することにしました．
FOSS4G 2024 Japan に参加すると，そのオープニングでこの Advent Calendar が紹介されており，「まあじゃあ新参者ですが，書くしかないな」ということで，筆を執っています．
というわけで，初心者が DSM のタイルを生成しようとして得た Lessons Learned を紹介しようと思います．
アークエッジ・スペース (ArkEdge Space Inc.) は，オープンソース地理空間情報ソフトウェアとそのコミュニティへの貢献のために，FOSS4G 2024 Belém にスポンサー参加しています！#foss4g に限らず OSS には日々お世話になっているし，ちょうど南米で仕事があったので，参加するしかないよね，と pic.twitter.com/lEzX4EnSGe
— 溶けかけてるうさぎ (@_meltingrabbit) December 4, 2024
| オープンデータである「ALOS全球数値地表モデル (DSM)」から，全球の地形タイルを pmtiles で生成しようとした | |
| GDAL や mapbox/rio-rgbify, go-pmtiles などを駆使してゴニョゴニョ生成しようとしたら，なんと 51 日も費やしてしまった（試行錯誤含む．また，ほとんどがファイル処理実行時間であり，作業時間はそこまでではないことに注意．） | |
| 様々な気づきや学びがあったので，まとめる |
下図のような感じに MapLibre GL JS で全球どこでも地形を表示できるようにするために，地形データ (DSM) の全球 pmtiles をつくりました．
なお，下の例では，ベースマップは HERE を使っています．
また，今回は 30 m 分解能の（DEM ではなく）DSM ということで，下図の新橋の結果のように，東京のビル群はちょっとイマイチですね（DEM ならよかった / DSM としてはビルに対して解像度不足）．
※ 上図のスクリーンショットを撮る際に，右下の Attribution Toggle を畳んだ状態で撮ってしまいました． Attribution は「AW3D30 (JAXA) | MapLibre | © 2024 HERE」です．後半のスクリーンショットも同様です．
ざっとこんな感じです．詳細と学びは後ろでそれぞれ記載します．
| ALOS全球数値地表モデル (DSM) "ALOS World 3D - 30m (AW3D30)" の全タイル (GeoTIFF) を取得する | |
| GDAL ですべてのタイルを結合した 1 つの全球 GeoTIFF を作成する | |
| mapbox/rio-rgbify で GeoTIFF から Terrain-RGB な mbtiles に変換する | |
| protomaps/go-pmtiles で mbtiles から pmtiles へ変換する |
各処理の詳細，気づきや学び，失敗談などを記載していきます．
なお，特に言及がない場合，処理時間はそこそこのスペックのデスクトップ PC (Windows) 上で動作する WSL2 での実行時間です．
ALOS全球数値地表モデル (DSM) "ALOS World 3D - 30m (AW3D30)" とは，陸域観測技術衛星「ALOS」搭載の光学センサ・パンクロマチック立体視センサ (PRISM) による解像度 30 m 相当（緯度経度 1 秒間隔を基本）の全球数値地表モデル (DSM) です．
より高解像度なデータセットもありますが，この 30 m 精度であれば JAXA 地球観測研究センター (EORC) から無償で利用できます．
この AW3D30 のデータダウンロードページが下図のような感じなのですが，見ての通り，緯度経度 30° ごとにサブページがあり，さらに 5° ごとに zip になっていて，一つづつポチポチするとダウンロードできます．
海しかないグリッドはデータがないので除いて，全部で 1419 グリッド（ファイル）をダウンロードする必要がある という感じのとてつもなくだるいダウンロードサイトなので，サーバー負荷を考慮しながら適当なスクリプトを書いて落とします ので，概ね 1500 回クリックしてダウンロードしましょう！
ちなみに，全部落とすとだいたい 140 GB 程度になりました．
$ ls -lh total 140G -rwxrwxrwx 0 root root 45M Jun 12 11:27 N000E005_N005E010.zip -rwxrwxrwx 0 root root 274M Jun 12 11:27 N000E010_N005E015.zip -rwxrwxrwx 0 root root 252M Jun 12 11:27 N000E015_N005E020.zip -rwxrwxrwx 0 root root 274M Jun 12 11:27 N000E020_N005E025.zip ...
取得した zip を解凍すると，以下のようにいくつかのファイルが入っているので，_DSM.tif で終わるファイルだけ抽出します．
これが，標高が Int16 で格納された GeoTIFF となっています．
ちなみに，全ファイルをかき集めると，約 431 GB の 23,993 ファイルとなりました．
$ ls N035E135_N040E140/ -lh total 843M -rwxrwxrwx 1 root root 25M Mar 10 2020 ALPSMLC30_N035E135_DSM.tif -rwxrwxrwx 1 root root 1.1K Mar 10 2020 ALPSMLC30_N035E135_HDR.txt -rwxrwxrwx 1 root root 9.9K Mar 10 2020 ALPSMLC30_N035E135_LST.txt -rwxrwxrwx 1 root root 13M Mar 10 2020 ALPSMLC30_N035E135_MSK.tif -rwxrwxrwx 1 root root 7.3K Mar 10 2020 ALPSMLC30_N035E135_QAI.txt -rwxrwxrwx 1 root root 13M Mar 10 2020 ALPSMLC30_N035E135_STK.tif -rwxrwxrwx 1 root root 25M Mar 10 2020 ALPSMLC30_N035E136_DSM.tif -rwxrwxrwx 1 root root 1.1K Mar 10 2020 ALPSMLC30_N035E136_HDR.txt -rwxrwxrwx 1 root root 14K Mar 10 2020 ALPSMLC30_N035E136_LST.txt -rwxrwxrwx 1 root root 13M Mar 10 2020 ALPSMLC30_N035E136_MSK.tif -rwxrwxrwx 1 root root 7.3K Mar 10 2020 ALPSMLC30_N035E136_QAI.txt -rwxrwxrwx 1 root root 13M Mar 10 2020 ALPSMLC30_N035E136_STK.tif ...
緯度経度 5° ごとにバラバラの GeoTIFF になっているので，これを統合して全球ファイルにします．
その時の気づきは次のようなものがありました．
| 解像度の異なるファイルを結合するときの注意 | |
| GeoTIFF を圧縮しないと処理時間もファイルサイズも大変なことになる | |
| TIFF は横方向単位で圧縮するため，巨大なファイルを扱う場合は tiled TIFF にしないと処理時間が大変なことになる | |
| 生成した巨大 TIFF は，Pyramid TIFF 化しないと QGIS などではまともに扱えない |
複数の GeoTIFF の結合には，ラスターおよびベクター地理空間情報データフォーマットのための変換用ライブラリ（FOSS4G ですね！）である GDAL の gdal_merge を使います．
最終的に用いたコマンド引数は次になりました．
$ gdal_merge.py -init 0 -co PREDICTOR=2 -co TILED=YES -co COMPRESS=DEFLATE -co BIGTIFF=YES -o ./dsm.tif --optfile ./tiff_list.txt
コマンドライン引数は，頭から「結合時の初期値（つまりデータがないピクセル）の値は 0 （欠損タイルは海なので標高 0 m）」，「水平方向の差分符号化を使う」，「Tiled TIFF を使う」，「圧縮アルゴリズムは Deflate」，「BigTIFF を有効に」，「出力ファイル名」，「結合する入力ファイルのリスト」です．
GeoTIFF に関するオプション（-co のオプション）は，GDAL の GTiff Raster Driver に詳しく書いてあります．
以下より，気づきの詳細を見ていきます．
『ALOS World 3D-30m（AW3D30） Version 4.1 プロダクト説明書 第 1 版』を読むと，次のような記載があります．
つまり，地上は球面であり，それを緯度経度グリッドでメッシュを切ると高緯度地域では解像度が高くなってしまうため，高緯度側ほどメッシュを荒くしている，ということです．
gdal_merge では，最初に入力されるファイルの解像度が出力ファイルの解像度となります．
したがって，今回の例では tiff_list.txt の先頭行には緯度が 60° 以下のファイルを記載する必要があります．
その状態で結合を実行すると，たとえば緯度 70° ～ 80°（ゾーン III）のファイルは 1 ピクセルが東西方向に（同じ値で）3 分割に高解像度化されて結合されます（下 2 図参照）．
逆に，緯度 70° ～ 80° のファイルを先頭行に記載してしまうと，赤道付近のピクセルは東西方向 3 ピクセルが 1 ピクセルに丸められてしまいます．
さて，上の図で結合前の 1 ピクセルが東西に長細く見えますが，下図に示すように実際はほぼ正方形の領域です．
これは EPSG:4326 に投影しているために極に近い高緯度ほど横に引き伸ばされているだけであって（そもそも地理座標系では高緯度地域で東西方向が引き伸ばされるので，先述の通り AW3D30 ではゾーンによって東西方向のピクセルスペーシングを段階的に変更している），もちろん EPSG:3857（いわゆる Webメルカトル）でみると（ほぼ）正方形にみえます．
今はただ純粋に TIFF ファイルを結合しているという気持ちなので，画像ファイル（TIFF ファイル）の 1 ピクセルが正方形であったほうがわかりやすいかなと思って EPSG:4326 で作業してます．
ここで，このような解像度が異なるタイルを結合することによる弊害に一つ気づいてしまったので記しておきます．
たとえば MapLibre GL JS のような描画ライブラリでは，データをいい感じに補間（やスムージング）して表示してくれます．
したがって，あるデータを上のように雑に高解像度化すると，その処理に悪影響が出ることが予想されます．
イメージとしては次のような図で考えるとわかりやすいです．
本来は青線で示すように補間されるところが，雑に高解像度化したことによって赤線のようになってしまう，ということです（なお，図では雑に折れ線で示していますが実際の補間方法は実装依存）．
実際，下図は北緯 81° 付近の生成した DSM の様子ですが，高解像度化されなかった南北方向の稜線は滑らかであるものの，東西方向の稜線はおよそ 30 m ごとに階段状になってしまっています（衛星写真の解像度が悪く見えるのは，極付近は需要が低くそもそも提供されるデータの解像度が悪いだけ．）．
この問題がぱっと解決できるか調べるために，gdal_merge.py の実装 を眺めてみました．
内部で入力ファイルリストの先頭ファイルを元に出力となるデータの解像度が計算された後，ファイルリストに記載されいてるファイルを 1 枚ずつ Python Raster API の osgeo.gdal.Band class の ReadRaster で読み出し，WriteRaster で書き込む，という処理を繰り返しています．
その近傍のコメント には
if the destination file is a different resolution, or different
image pixel type, the appropriate resampling and conversions will
be done (using normal GDAL promotion/demotion rules).
とあり，normal GDAL promotion/demotion rules とやらの方法で，複数の入力ファイルの解像度をあわせているようです．
何が normal なのか知らないので更にコードを読み進めると，ReadRaster の実行時に出力解像度に合わせて読み出しが行われており，最終的に GDAL の C++ 実装の GDALRasterBand::ReadRaster が呼ばれていることがわかりました．
このリファレンスを読むと，
eResampleAlg -- Resampling algorithm. Defaults to GRIORA_NearestNeighbour.
とあり，なるほど，nearest で入力ファイルのリサンプリングを行っている，ということがわかりました．
うん，気になったらすぐ実装を見に行けるという意味で，OSS は本当に便利ですね！
とはいえ，残念ながら gdal_merge ではこのリサンプリングアルゴリズムを変更するオプションが与えられていないのと，そのためだけに gdal_merge そのものを修正するのもちょっとめんどくさい（多分，ReadRaster の代わりに ReadRaster1 を使って適切な GDALRIOResampleAlg を指定してあげればよさそう，というだけではありそう．）なぁという気持ちです．
そうすると，本来はタイルの結合前に適当に補間して高解像度化することで解像度を揃えておくことが適切のように思えます．
最適な補間方法も，ケース・バイ・ケースだと思いますので（たとえば今回のケースは，南北方向の補間は不要なので，東西方向のピクセルだけで補間処理するのが良さそう，など．）．
今回は高緯度地域のこの問題は一旦無視することにしましたが，もし，ぱっと解決できる方法をご存じの方がいらっしゃったら教えて下さい．
ここ処理は試行錯誤等で多くの時間を浪費してしまったところです．
最初，何も考えずに以下のコマンドラインオプションで処理していました（TIFF の圧縮をしないと大変なことになることは知っていた）．
この項の冒頭で記載したものとの差分は -co PREDICTOR=2 と -co TILED=YES がないことのみです．
$ gdal_merge.py -init 0 -co COMPRESS=DEFLATE -co BIGTIFF=YES -o ./dsm.tif --optfile ./tiff_list.txt
ベンチマークとして南米を包含するエリア（北緯 15° ～南緯 60°，西経 30° ～ 85°）を処理すると，21 分ほどで 1794 ファイルの 合計 43.4 GB の GeoTIFF が 17 GB の GeoTIFF に結合されました．
これは，AW3D30 の定義域（緯度 85° 以下）のおよそ 6.8% に相当するため，全球 23,993 ファイル (431 GB) の結合は 5 時間ほどで完了し，ファイルサイズも 167 GB 程度か（AW3D30 のタイルのない海の部分もかなりあるので，多少これよりも増えそう），と予想できます．
満を持して全球の結合を開始すると，32 時間かかって 500 GB の出力 TIFF ファイルを書き込んだところで，disk full のため処理が中断．
また，ディスクアクセスログを見るとどうやら処理の進捗はたった 6% 程度であったことも判明．
一体何が起きてるんだ...？
いくつか行った実験で解決につながったのは，次のようなものでした．
| 東経 130° ～ 140° の南北に長細い領域 875 ファイル (16.0 GB) を結合すると，23 分ほどで処理が完了し，結合後のファイルサイズは 4.9 GB． | |
| 北緯 0° ～ 10° の東西に細長い領域 1144 ファイル (27.6 GB) を結合しようとすると，処理開始後 5 時間弱で進捗は 40% ほど，出力ファイルサイズは 160 GB を超過． |
結合後の領域が東西に伸びると，指数関数的に処理時間が伸びているように感じられました．
そこで，GTiff Raster Driver をより詳細に読んでいくと，次のような記載を見つけました．
TILED=[YES/NO]: Defaults to NO. By default striped TIFF files are created.（中略）
BLOCKXSIZE=256. Sets tile width.（中略）
BLOCKYSIZE=
（中略）
PREDICTOR=[1/2/3]: Defaults to 1. Set the predictor for LZW, DEFLATE and ZSTD compression. The default is 1 (no predictor), 2 is horizontal differencing and 3 is floating point prediction.（中略）
ここで思い出す，Striped TIFF と Tiled TIFF...．
TIFF を含め通常の画像ファイルは，ピクセルのデータは画像の左上から行（横）方向に向かって格納されていきます（下図左）．
そして，圧縮もこの行を数本まとめた Strip 単位で行われます．
こうすることで，シーケンシャルにアクセスできます（通信回線の細い環境でインターネットを見ると，画像が上からちょっとずつ表示されるアレ）．
一方で画像編集などの文脈では，画像の一部を部分的にアクセスして書き換えたり，大きな画像をメモリ効率よく処理したりするために，データの格納順を Strip 単位ではなくタイル化したブロック単位で行うことがあり，これが Tiled TIFF です（下図右．なお，地図配信におけるタイル化とは別の概念）．
全球では，東西方向のピクセル数は 1,296,000 ピクセルに達します．
一部分のデータを追記するたびに，この行をすべて読み込み，再圧縮（当然 Deflate の処理もデータ列が部分的にでも変わるとやり直し），再書き込みを行うとなると，ありえないほどの非効率な処理となるわけです．
次に Predictor です．
今回のデータは標高（メートル単位の整数値）なので，値が急激に変化することはほぼないはず，つまり隣接するピクセルの値の差分はある程度の小さい値で固まっている（差分データの度数分布は 0 付近に集中する）はずです．
となると，度数が上位のデータは Deflate のハフマン符号化でごっそり潰れることを考えると，「2 is horizontal differencing」が有利です（7 割の領域を占める海は 0 が連続するので，どちらでも良さそう）．
というわけで，冒頭に示した
$ gdal_merge.py -init 0 -co PREDICTOR=2 -co TILED=YES -co COMPRESS=DEFLATE -co BIGTIFF=YES -o ./dsm.tif --optfile ./tiff_list.txt
で処理すると，たったの 11 時間で処理が完了し，生成された GeoTIFF のファイルサイズはたったの 90.3 GB でした．
なお，全球で -co PREDICTOR=1（デフォルト）と -co TILED=YES を指定すると 33 時間ほどの処理で 113 GB の出力が得られました．
定性的には，-co TILED=YES によってディスクアクセスがまともになり，-co PREDICTOR=2 によって Deflate の圧縮効率が上昇する，というところです．
全球において，-co TILED=YES を指定しないとファイルサイズが爆発する理由はいまいちわかりませんでしたが，同様の問題は Geographic Information Systems Stack Exchange でも報告されていました．
そんなこんなで無事に完成した全球の GeoTIFF はこんな感じです．
ここの PREDICTOR と TILED のオプションを発見するための試行錯誤で 3 日潰しました．
データサイズが大きいと，一回の施行にかかる時間も長いので大変ですね．
生成された全球の GeoTIFF を QGIS などの GIS ソフトウェアで扱おうとすると，特にズームレベルが低いときには広大な面積（ピクセル数）のデータを読みに行くため，非常に動作が遅くなります．
TIFF には，一つのファイルに複数の画像を格納できることを利用した Pyramid TIFF という形式（下図参照）があります．
これは異なる解像度の画像を重ねて格納することで，大きなファイルでも効率的に表示・処理することができるようになるものです（つまりズームレベルが低いときにはわざわざ高解像度なデータをかき集めるのではなく，解像度の低い画像をもってきて表示する．）．
Pyramid 化には gdaladdo を使います．
以下のコマンドで，元画像を 1/2, 1/4, ... , 1/8192 にした画像を格納したファイルを生成することができます．
注意点としては，入力ファイルを Pyramid 化する形，つまり，入力ファイルを更新する形（破壊的操作）になるため，オリジナルのファイルのバックアップをとっておくことが重要です．
$ gdaladdo -r average dsm.tif 2 4 8 16 32 64 128 256 512 1024 2048 4096 8192
生成された Pyramid TIFF を QGIS で開いた様子が下図です．
元画像の解像度が 1,296,000 x 604,800 であり，それを縮小した 13 枚の画像が生成されていることがわかります（QGIS では 14 枚つくるといいね，とレコメンドされてますね）．
また，生成された結果は gdalinfo などでも確認できます（まあ，QGIS の内部で動いてるのも GDAL なので，同様ですが．）
これによって，QGIS でストレスなく操作できるようになりました．
一枚の GeoTIFF になっているので，あとは普通に処理していくだけです．
そう，処理にかかる時間以外は...．
なお，GeoTIFF から一発で pmtiles へ変換する方法を知らない（誰か教えて下さい）ので， mbtiles を経由します．
今回生成する pmtiles は，最終的には MapLibre GL JS（FOSS4G です！）で利用しようと考えており，公式ドキュメントのサンプルコードや MapLibre Style Spec の Sources などを見ると，「A raster DEM source. Only supports Mapbox Terrain RGB and Mapzen Terrarium tiles.」とあるので，Mapbox Terrain RGB 形式で作っていきます．
Mapbox Terrain RGB は標高データを RGB 画像にエンコードしたものです．
デコードの定義は，公式ドキュメントより，
elevation = -10000 + (({R} * 256 * 256 + {G} * 256 + {B}) * 0.1)
であり，
| RGB のビット列をそのまま 24 bit 整数として読む | |
| 量子化の分解能は 0.1 m | |
| -10,000 m でオフセットすることで，正数としてエンコード（海抜未満の標高に対応） |
といったものです．
ちなみに，ここは盛大につまずいたポイントの一つでした．
最初は普通に GeoTIFF をそのまま gdal_translate と protomaps/go-pmtiles で pmtiles を作っていたのですが，生成されたタイルを見ると標高の高いところのデータが潰れてしまいました（下図参照）．
データを眺めていると， mbtiles がもともとの GeoTIFF のデータ型である Int16 の表現力がないことに気づきます（こういうとき，思わず，「あ゛゛゛～～」ってなりますよね．）．
GDAL の MBTiles Raster Driver のドキュメントの「Tile formats」の項に
MBTiles can store tiles in PNG, JPEG or WEBP (since 3.8). Support for those tile formats depend if the underlying drivers are available in GDAL. By default, GDAL will produce PNG tiles.
It is possible to select the tile format by setting the creation/open option TILE_FORMAT to one of PNG, PNG8, JPEG or WEBP. When using JPEG, the alpha channel will not be stored.
（中略）Note that at that time, such an 8-bit PNG formulation is only used for fully opaque tiles, as the median-cut algorithm currently implemented to compute the optimal color table does not support alpha channel (even if PNG8 format would potentially allow color table with transparency). So when selecting PNG8, non fully opaque tiles will be stored as 32-bit PNG.
とあり，PNG (PNG32) でのグレースケールであれば 8 bit しか表現能力がないので，うまくいかないのも納得です（16 bit グレースケール PNG は使えないのね．）．
根本的に間違えていることを悟った後，MapLibre GL JS の公式ドキュメントのサンプルコードを手元で動かし取得されるタイルを眺めていると，そのタイルがグレースケールではなく RGB Truecolor であることに気づき，上記の正しい仕様にたどり着くことができました．
まあ，最初から公式ドキュメントを "きちんと" 読みましょう，ということではあります．
結果的に不必要なデータ変換などで，ここで 1 日潰しました．
さて，生成すべきファイルフォーマットがわかったので，変換していきます．
Rasterio などを使って変換スクリプトを自作しても良いのですが，Mapbox から変換ツール mapbox/rio-rgbify が公開されているので，これを使っていきます（FOSS4G に感謝）．
使用したコマンドは以下です．
$ rio rgbify -j 1 -b -10000 -i 0.1 --max-z 12 --min-z 2 --format webp dsm.tif dsm.mbtiles
ズームレベルは 2 ~ 12 を指定していますが，AW3D30 の赤道上での分解能は 30.922 m であるので，まあ若干足りてないですが Lv.12 (40,000,000 [m] / 2^12 / 256 [pixel] = 38.2 [m/pixel]) にしました．
本当は Lv.13 まで入れておくと良いかもなのですが，計算時間がさらに 5 倍になると考えると...，ね？
さて，上記コマンドの実行時間についてです．
ベンチマークとして南米を包含するエリア（北緯 15° ～南緯 60°，西経 30° ～ 85°）を，-j 4 と引数を修正（4 並列で処理）して実行した結果，実行時間は 43 時間程度でした．
GeoTIFF の結合のときと同様に，これは，AW3D30 の定義域のおよそ 6.8% に相当するため，つまり全球の処理時間は 26 日か...，海しかないところはもう少し処理は早くなるよな...，などと思いながら恐る恐る全球での実行を開始します．
．．． 最初の試みから 36 日後に，手元の Windows PC での処理を諦めました．
理由はざっとこんな感じ．
処理が半分を超えてくると，しばしば BrokenPipeError: [Errno 32] Broken pipe や numpy._core._exceptions._ArrayMemoryError: Unable to allocate 2.00 MiB for an array with shape (512, 512) and data type float64 が発生する | |
| 業務 PC のバックグラウンドで動かすには，あまりにも邪魔 | |
| しばしばミスって disk full させちゃう（出社して PC をつけたら disk full してる，など） |
一つ目のエラーは原因の特定には至らなかったのですが，マルチプロセス処理部分で発生していたので -j 1 とすることで回避できました．
まぁ，残りがどうしようもないので，要件を満たす AWS EC2 インスタンスを立ててそこで作業することにしました（札束で解決ってやつですね）．
なお，この処理はディスクアクセスではなく CPU が律速なので安い HDD をアタッチし，また，並列化もしないので CPU コア数も下げたインスタンスを活用しました．
その結果，21 日弱で変換が終わりました．
出力されたのは 288 GB の mbtiles で，含まれるタイル数は 21,008,336 でした．
以下に，その時の EC2 のメトリクスを掲載します．色々と興味深いですね．
まず，CPU は 8 コアのインスタンスを使った一方で並列処理をさせていないので，1/8 である 12.5% を上限に張り付いています．
書き込みスループットをみると，ズームレベルが上がっていくにつれ書き込み速度が上がっていることがわかります（処理は z/x/z の順にインクリメントしながら処理されることが，mapbox/rio-rgbify のソースコードの「ここ」を見るとわかる）．
これはズームレベルが上がるごとに一つのタイルの生成時間は短くなっていくことを示しており，ズームレベルが小さいと多くのデータを読んでいる事がわかります（多分，1 タイルあたりの元データ領域が広くなるため）．
途中で 2 回ほど書き込みレイテンシがバーストして CPU 使用率が下がっている部分もありますが，ここは何なんでしょうね？
ちなみに，こんなにも処理時間がかかるのにもかかわらず，mapbox/rio-rgbify はプログレスバーなどの進捗を表示するオプションがありません．
そこで手元では，元のソースコードを編集し進捗がわかるようにしていました．
自由に改変できたり内部仕様を読み取れたりするので，OSS って最高ですね！
さて，全球 DSM の pmtiles を生成する作業もこれで最後です．
この最後の処理では，絶望と感動（データベースってすごいね，という再認識）がありました．
変換作業自体は，protomaps/go-pmtiles を使って一発です．
$ pmtiles convert dsm.mbtiles dsm.pmtiles --tmpdir=./tmp/
ここでは --tmpdir を指定しています．
pmtiles convert の挙動としては，os.CreateTemp で生成される一時フォルダに生成されたタイルをためていき，最後に一時フォルダにあるタイルにヘッダ類をくっつけて出力ファイルとしてコピーする，というものだったからです．
使っていた EC2 ではシステムストレージをケチっていたので，最初何もオプションをつけずに実行したところ，タイル変換作業が始まってすぐ，disk full で死にました．
（OSS はコードを自ら読みに行けるので，本当に便利！）
問題は実行時間でした．
ベンチマークの南米（北緯 15° ～南緯 60°，西経 30° ～ 85°）は 39 時間で処理できたので，全球は 24 日で終わるだろう，と単純には思うのですが，きちんと検討すると 1350 日，つまり 3.8 年かかることが判明しました．
3.8 年...？ 絶望です．
ちなみにこの概算は，処理時間は「前処理時間 + 1 タイルの処理の固定時間 + 1 タイルの処理の可変時間」であり，可変時間は処理済みタイル数に比例する（処理が進むと線形に処理時間が増えるということ，これは後述するように後で正しい仮定であることが判明する）ということを仮定して，瞬時処理時間などの観測結果とフィッティングして求めただけです．
観測される挙動としては，初速は秒速数数 1000 ～ 数 10,000 タイルで処理されていたものが，数日後などには秒速 2 タイルなどに速度低下していました．
更にこのとき，頻繁なディスクアクセスが発生していることも観測できました．
同僚の KOBA789 に相談すると，「pmtiles convert、1 tile 変換するごとに mbtiles (SQLite3) の tiles テーブルをフルスキャンしてる！！！！ 」と．
どうやら，最初に pmtiles のメタデータを構築し，その後各タイルを生成するときに mbtiles 上の対応する各タイルを検索して読み出すのですが，その検索にものすごい時間がかかっていました．
そして，pmtiles convert で生成するタイルの順番が mbtiles を生成した時に処理した順番と（たまたま）一緒だったために，mbtiles のタイルを検索するために tiles テーブルを読むとき，序盤に処理されるタイルはそのテーブルの上位にあるためフルスキャン時にすぐにヒットするが，徐々にテーブルの下のタイルになるにつれ検索時間が線形に増加していく，というものでした．
この挙動は mbtiles にインデックスが貼られていないせいでした．
というわけで，
create index tiles_idx on tiles(zoom_level, tile_column, tile_row);
とインデックスを張ってから処理してあげると，2 時間かからずあっさりと変換が終了しました（めでたしめでたし）．
こんなところで，インデックスによる O(N) → O(logN) を体感することになるとは...．
後日落ち着いて mbtiles のドキュメントを読むと，
The database MAY contain an index for efficient access to this table:
CREATE UNIQUE INDEX tile_index on tiles (zoom_level, tile_column, tile_row);
とあり，「MAY」かぁ～，という気持ちになりました．
まあ確かに，インデックスの貼り方は use case 次第，という気もします．
さらに，タイルは重複することはありえないので，公式ドキュメントのように UNIQUE をつけるほうが良さそうですね．
なにはともあれ，これでファイルサイズが 282 GB 程度の，目的の pmtiles が完成（下図）しました．
これで，「」に示すような 3D 地形表示ができるようになりました．
というわけで，初学者が全球 DSM pmtiles を作ろうとしたら想像以上に時間がかかりました，というお話でした．
なにはともあれ，巨大なデータを扱うときはきちんと見積もりをすることが大事ですね．
あれ？ 数 100 GB って，「小さい」のか...．
はい令和版 pic.twitter.com/Aa222fEw4S
— Minero Aoki (@mineroaoki) November 19, 2021
いろいろつまずきもありましたが，学びも多かったです．
みなさんもつまずきながら頑張っていきましょう！
名前
Email (※公開されることはありません)
コメント