タイル画像を作った記録 #1 - ここのことはなかったことにするかも

はじめに

某方面からあるラスタを頂きました。
ピクセルサイズは忘れたのですが、Google Maps等のズームレベル16よりもうちょっと細かい解像度で関東平野全体のラスタで、圧縮なしTIFFで40G程度だったと記憶しています。
前からクレクレ言ってて向こうから「送った」とか連絡をもらった際、ブルーレイだろうと思ってたらUSB接続型のハードディスクでした。

このデータでタイル画像を作ろうとしました。ただタイル画像を作るだけでなく、もうちょっとこだわりのタイル画像を作るべきだろうと考えました。
そのこだわりを持つことこそ間違いだったのではないかと途中何度も思いましたが、おおむねできあがったので、今となっては良い思い出です。

このエントリは、要領悪い人が苦労して巨大ラスタ画像からタイル画像を作成した実話をあますところなく文書化したものです。

何を作るか

作成するタイルは、座標系はWebメルカトルと正距円筒、画像形式はPNGとJPEG、の組み合わせで4種類あります。また (正距円筒, PNG) の組み合わせについてはKMLもあわせて生成します。あと、リサンプリングアルゴリズムには、今回はこだわりのランチョスでやったろうと決意しました。

gdal2tiles.py を実行してそっと終了した

とりあえず gdal2tiles.py でタイル画像が作れますので実行。…1日放置して10%。単純計算したら10日はかかる。しかもWebメルカトルだけで10日。正距円筒版も作ったりと考えたら、かなりな時間が取られる。これでは話にならない。ということで、そっとCtrl-C。

先にwarpを実行しておく

元画像はUTMだったので、投影変換が必要。その際に、ベースタイル(最大ズームのタイル)と同じ解像度(この場合はUnit-Per-Pixel)で元データを変換して、その後で、near (nearest neighbour) のリサンプリングでタイル化してやろう、と考えました。

理由は、手戻りできることと、タイルごとにリサンプリングを実行すると効率が落ちるような気がしたためです。

ちなみにWebメルカトルのズームレベル17を例にとると、

2*PI*6378137/(256 * 2**17) = 1.19432856696

ですので、コマンドラインは

% gdalwarp -t_srs "epsg:3857" -tr 1.19432856696 1.19432856696 (元ファイル) (出力ファイル)

となります。

gdal2tiles.py の改造

マルチプロセスで動くようにしました。

管理プロセスは存在させず、プロセス数とプロセス番号とをコマンドライン引数で渡すものです。ので、複数のシェルを上げて順次プロセス番号をずらしたコマンドを実行します。

単純にXのループの開始をプロセス番号でずらして、1ずつ増やしているのをプロセス数ずつ増やすようにしただけです。

ベース(最大ズームのタイル)を作成する際は全く問題がありませんでしたが、それが終わってオーバビュータイル(最大ズーム以外のズームのタイル)作成段階で、各プロセスが例外を出してバタバタと落ちていき、1プロセスだけ残りましたが、全部動かさないといけないので、生き残ったプロセスも泣く泣くkill。

原因は、他プロセス依存タイルを使用しなければならない点でした。
オーバビュータイルは、ズームレベルの1段大きいタイル2*2=4個から生成しています。Xごとにプロセスが割り当てられているので、左のタイルと右のタイルは担当プロセスが異なります。プロセスの進行にばらつきがあって、不幸にも1段下のタイルがまだできていない場合には、本来あるべきタイルが無くて、例外が発生してしまったのです。

ズームレベルごとに逐次プロセスを起動させることで回避しましたが、やはり作業プロセスを管理するプロセスが必要だろうと思います。