gdal2tiles.pyのマルチプロセス改造は既に存在していた

id:yellow_73:20130809とかid:yellow_73:20130814とかの続き。

https://twitter.com/tmizu23/statuses/363253994665152512 参照。

http://trac.osgeo.org/gdal/ticket/4379 で、マルチプロセス用パッチが出てたという。

同じようなこと考えてるんだ、20ヶ月も前にな。
困ってたガワも作ってらっしゃるし…。

あと発見が1ヶ月早ければよかったんだが仕方ないorz

タイル画像を作った記録 #2

前回のあらすじ

id:yellow_73:20130809 の続き。

先にgdalwarpでリサンプリングすることでgdal2tiles.pyではニアレストですむようにし、また gdal2tiles.py の魔改造によりマルチプロセス化した。

魔改造スクリプトはガワが全く無い状態なのでガワを被せたいなと思ってますが、マルチプロセスでプログレスをキャラクタベースで出すなんてTERMCAPなんかとうのむかしに忘れたぜヒャッハー、となっているため無理です。

JPEG の生成で透過部で白と黒が出てあせる

gdal2tiles.py はPNGを生成しますが、ここはJPEGも生成しないといかんでしょう、と考えました。

JPEGはファイルサイズが圧倒的に小さいという利点があります。ただし、アルファチャネルをサポートしてない、コントラストが強い箇所の近傍ではモスキートノイズが出る、という欠点があります。

今回は、対象が背景図に使うので透過は必須でない、地図だけどにじみ等も結構あるのでコントラストが強すぎるとは言えない、むしろ3G回線といった高速とは言えない回線でも使うことを考えたらJPEGも用意するメリットはある、と判断しました。

maptiler ではJPEGもあわせて生成してくれますが、gdal2tiles.py ではできないので、ImageMagick を使いましたPNGファイルからJPEGファイルに変換するのは次の通り。

% convert a.png a.jpg

サフィックスを見て判断してくれます。ImageMagickばんざい…。

find コマンドで *.png を見つけたら jpeg に変換するようにしました。

で結果を見ると…。

pngで透過になっている領域が、黒になった箇所と白になった箇所が出ました。

たぶん gdalwarp で失敗してたんだと思います。

検討しきってないのですが、考えられる原因を。元データはアルファチャネルなし圧縮なしのTIFFで 255,255,255(白) をNoDataにしてました。gdalwarpを実行する際に、アルファチャネルを追加、NoDataを255,255,255,0と指定しました。アルファチャネルはいいとして、NoData指定がgdalwarpによる回転(元データはUTM)による現れる元データが存在していないピクセルに255,255,255,0が入るかなと思ったら入ってない。gdalinfoでwarp後のファイルを見たら全チャネルのNoDataが255になってる。

これ、盛大に勘違いしてる可能性がありますね…。

事前にチェックしてたはずがチェックしてなかった

チェックしてたんですよ、アルファチャネルを取ったら白以外が現れないかは。

GIMPで。

gdalwarpで生成した画像の一部をgdalwarpで取り出して、GIMPで開いてアルファチャネルを除いてやったんです、そうすると透過部が背景色になるのを確認してたんですよ…これ意味無いですからねorz 後で背景色を黒にした上でアルファチャネルを除いたら黒になったところを見ると、元のピクセルが持つR,G,B値に関係なくGIMPで指定した背景色が使われます。

pngから透過をサポートして無いjpegに変換するとチェックできると思います。

JPEG の生成で透過部で白だけが出るようにしたがこれはおすすめしない

手戻りせんでええように考えてたのが、gdalwarp以前まで手戻りしなければならない危機。これは避けたいので、256*256の白一色のpngファイル (ffff.png)を用意して、

convert -composit ffff.png (元png) (出力jpg)

で変換するようにしました。

これでとりあえずは透過領域は白で出るようになった。ほっとした。gdalwarpははっきりいって億ピクセルかつランチョスでりんサンプリング、となってくると1日では終わらんのですから。

で、PNGはどうせ透過サポートですから放置。透過領域の隠れた色問題はpngでも響いてくることになります。

Google Earthさんのばかぁ

その後KML生成もうまくいったので、とりあえずGoogle Earthで見てみることに。

…なんか変なのが出てる。

透過されている部分が表示されていました。前述と同じように、黒領域と白領域が表示されます。

…原因不明の怪現象に呆然とするしかないですね。

呆然としたままいても話が始まらないのですが、とりあえず放置しようと思って別作業にかかっていたところ、ふと前に作成したpngと比較してみようと思いました。

gdalinfo で新旧PNGを見ると、新しい方はインデクス化されていて、古い方はインデクス化されていないのが判明。地理情報を持たない画像ファイルでもインデクスされているかどうかが分かったりするので、gdalinfoを利用してます。

そうすると、旧版は32ビットPNG、新版は気合入れてoptipng -o7で圧縮したカラーインデックスPNG

もしかして、32ビットPNGにしたら問題解決? まさかね… など思いながら、32ビットPNGにしてみたら、ああ、うまくいくやないの…

透過領域を背景色に統一しつつインデックスPNGをG32ビットPNGに変換する

convertで32ビットPNGで出力するには、ファイル名の前に "png32:"を付けるとOK。あと、黒領域問題を解決するため、透過部分は背景色に変更するオプション(-alpha background)を付けることにしました。

convert -alpha background (pngファイル) png32:(中間ファイル)
mv (中間ファイル) (pngファイル)

を繰り返し実行すると、インデクスを使わず、透過部が背景色となった PNG に変換できました。

前述のJPEG生成において、ffff.png と重ねて使うより、先にこちらを実行したうえで JPEG を生成するともっとスマートに行ったと思います。

あと余談ですが、KMLは正距円筒でのみ対応しているので、そちらのPNGしか処理しておらず、Webメルカトルの方は現在作成しているところだったりしますので、隅っこ近くのタイルを落として convert なりでJPEGとかに変換すると「出る」と思いますorz

感動のフィナーレ

ここでは出てこなかった Windows x64 用バイナリを作ってみたが地理情報がうまく載らずにあーだこーだしてた時期を除くと、着手から2週間で、億ピクセル画像のWebメルカトル(z<=17)と正距円筒(z<=16)のタイルを作りサーバにアップできたました。まじめにgdal2tiles.pyを駆動してたら正距円筒のPNGができあがらなかったというのを考えると、かなり速度向上が図れたと思います。

問題点は

  • 魔改造gdal2tiles.pyを外に出せる程度にガワ付けたい
  • 透過領域を持たせたいなら Google Earth に渡す PNG は png32 にしとくべき
  • ImageMagickばんざい optipngばんざい

タイル画像を作った記録 #1

はじめに

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

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

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

何を作るか

作成するタイルは、座標系はWebメルカトルと正距円筒、画像形式はPNGJPEG、の組み合わせで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段下のタイルがまだできていない場合には、本来あるべきタイルが無くて、例外が発生してしまったのです。

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

センターマーカーでホイールイベントを邪魔しない別の方法

divなりの要素を作ってborderで細い横線だけ持つ要素を作る。同じように縦線だけ持つ要素を作る。これらをセンターマーカーとなるような位置に持っていく。そうすると、縦線、横線の真上だけホイールイベントを邪魔するので、以外と分かりにくくなります。

画像でマーカーを作るから、透明な部分にマウスカーソルがあっても邪魔されてしまう、と。

電子国土さんがそんなかんじです。

センターマーカーを表示するコントロール

id:yellow_73:20120909 の続き。id:yellow_73:20130717 の通り、マウスホイールイベントが効かないのが解決されていましたので、あらためてセンターマーカーのコントロールを作ってみます。
OpenLayers 2.13以上でのみ有効です。

resizeイベントが取れない云々と言ってたことについては、 CSS でどうにかしちゃうことにして、resizeイベントを拾わなくて良いようにしました。

<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=UTF-8">
<script src="http://openlayers.org/api/OpenLayers.js"></script>
<script>
function init() {
  // コントロールクラスを作成
  // -------- ここまでクラス
  var CenterMarker = OpenLayers.Class(OpenLayers.Control, {
    // 画像ソース
    url: null,
    // 画像オフセット デフォルトでは画像左上隅がマップ中心に来る
    offset: null,
    // 初期化
    initialize: function(options) {
      OpenLayers.Control.prototype.initialize.apply(this, arguments);
      this.url = options.url;
      this.offset = options.offset;
    },
    // 廃棄
    destroy: function() {
      OpenLayers.Control.prototype.destroy.apply(this, arguments);
    },
    // 表示
    draw: function() {
      OpenLayers.Control.prototype.draw.apply(this, arguments);
      var ox = 0, oy = 0;
      if( this.offset ) {
        ox = this.offset.x;
        oy = this.offset.y;
      }
      // コントロールの左上をマップ中心に持っていく
      this.div.style.position = "absolute";
      this.div.style.left = "50%";
      this.div.style.top = "50%";
      // マーカ作成
      // left, top はマップ中心からのオフセット
      var img = document.createElement("img");
      img.src = this.url;
      img.style.position = "absolute";
      img.style.left = ox + "px";
      img.style.top = oy + "px";
      img.className = "olScrollable";
      // マーカをコントロールに追加
      this.div.appendChild(img);
      return this.div;
    }
  });
  // -------- ここまでクラス

  var centerMarkerControl = new CenterMarker({url: "cmark.png", offset: new OpenLayers.Pixel(-16,-16)});
  var options = {
    maxResolution:1.40625/2,
    controls: [
      new OpenLayers.Control.PanZoomBar(),
      new OpenLayers.Control.LayerSwitcher(),
      new OpenLayers.Control.Navigation(),
      centerMarkerControl,
      new OpenLayers.Control.Attribution({displayClass: 'prmtcd'})
    ],
    numZoomLevels: 18
  };

  var map = new OpenLayers.Map('MAP',options);
  var kibanwms = new OpenLayers.Layer.TMS(
    "KIBAN 25000 TMS",
    "http://www.finds.jp/ws/tmc/",
    {
      layername: "KBN25000ANF-4326",
      type: "png",
      attribution: '<a target=\"_blank\" href="http://www.finds.jp/wsdocs/kibanwms/index.html">基盤地図情報(平20業使、第449号)</a>',
      isBaseLayer: true
    }
  );

  map.addLayer(kibanwms);
  map.setCenter(new OpenLayers.LonLat(133.39, 34.5), 15);
}
</script>
<style type="text/css">
html, body, #MAP {
  margin: 0;
  padding: 0;
  position: relative;
  width: 100%;
  height: 100%;
}

.prmtcd {
  bottom: 4px;
  right: 4px;
  background: #fff;
}
</style>
</head>
<body onload="init()">
<div id="MAP">
</div>
</body>
</html>