GDALの外部実行ファイルとPerlでパノラマ写真をタイル画像にしてやる

  • 3
    いいね
  • 1
    コメント
この記事は最終更新日から1年以上が経過しています。

はじめに

FOSS4G Advent Calendar

この記事はFOSS4G Advent Calendar 2014 11日目の記事です。

さあハードル下げるよー

Advent Calendarでは、前日にすごいことをされると次の日の人は、複雑な気分になるものです。少なくとも、FOSS4G Advent Calendarでは、数年前から「ハードル上がった」等の言葉がときどき出ています。

現に前日の記事で、ハードルが上がってますね。
しかしここでハードルが下がります。

若干「ジオ」なものとは違うかも知れませんし、なにより、ネタかよと思われるようなプログラムを書いています(なのでGitHubに上げる気もない)。

この記事について

某所で巨大なパノラマ写真をもらいました。そういうものをもらったら、とりあえずスクロール地図のライブラリを使ってブラウザで負担なく見れるようにしようと思うところですね。OpenLayers 2で見れるようにしてみようと思いました。

しかしここで問題が。当然パノラマ写真画像に対してジオリファレンスが取れるようなものではありません。gdal2tiles.pyは、ジオリファレンスが取れているラスタ画像について動かすようにできています。MapTilerを使うと、ジオリファレンスが取れてないラスタ画像でもタイルを作ってくれるのですが、リサンプリングメソッドが決められない。そしてなにより、GUIなアプリケーションなので、複数ファイルがある場合には向きません。やはりバッチ処理的にCUIなツールを実行していくように仕込んで「後はコンピュータさんよろしく」で済むようにしたいところです。

ということで、GDALの外部実行ファイルとPerlを無理やり使って、タイル作成とHTML生成とを行う、いい加減なスクリプトを書きました。

ImageMagickが使えなかった

80000x20000ピクセルぐらいのTIFFをconvertで若干の拡大をしながらPNGに変換しようとしたところ、カーネルパニックで落ちました。

FreeBSD機では、カーネルパニックは、ほとんど遭遇しません。

ちゃんと調べていないので、環境変数をちょこっとだけ変えたり、いったんPNGにした後拡大するとかすると済んだのかも知れませんが、調査を切り上げ、大きなサイズの画像が得意なGDALを使うこととしました。

ご注意

GDALを使いますが、Perlから外部プログラムをsystem()で実行したり、popen()で実行した結果をパイプで読んだりで動かしたりしてます。GDALパッケージ不使用のピュアPerlという、何したいんだこいつやる気あるのか、というのをやってしまっているのがポイントです。

吐き気を催す等したらすぐさまスクリプトを追うのを止めてください。

GDALの注意点

  • GDAL_CACHEMAXはデフォルト設定から上に引き上げておいた方が幸せになれます(http://d.hatena.ne.jp/yellow_73/20110427)。
  • gdalwarp と gdal2tiles.py につっこむ画像には空間参照系を指定していないといけない。
  • gdalwarpでPNGを出力フォーマットに指定できないがgdal_translateならOK

方針

  • ズームレベル0は256*256ピクセルとなるようにする
    • 画像は左右方向、上下方向のいずれについても中央寄せにする
  • 最大ズームレベルは元画像より大きなズーム値の最小とする
    • 元画像の解像度は最大ズームレベルの解像度と最大ズームレベル-1の解像度の間になる
  • 最大ズームレベルはこだわりのlanczos、最大ズームレベル-1もこだわりのlanczos、それより小さいのはcubicとする
  • タイル付番は左下原点とする
    • gdal2tiles.pyはこれがデフォルト
    • Google Maps等は左上原点

空間参照系

今回の場合、元画像に空間参照系は付くはずがないものです。が、GDALの注意点の二つ目から、無理やりでも空間参照系は付けたいところ。

じゃあどういう空間参照系にすればよいか?

方針の一つ目から、Webメルカトルと同じにできそうですね。じゃあWebメルカトルとしましょう。

地球全体のサイズの半分は、おおむね 20037508.34 です。となるとですね、元画像のサイズは 40075016.68m もある巨大なものということになります。

どれだけ無茶か

ところで巨大パノラマ写真は自分で撮ったものでもないので、今回は、オリジナルの写真を使います。

ウサギさんの写真

このウサギさん、横幅40075016.68mの地図画像であると考えると、かるく2万キロメートルを超えることになります。

ためしに重ね合わせるとこんなかんじです。

p1.jpg

(C) OpenStreetMap contributors, CC-BY-SA

このウサギさんによって人類が滅んでも仕方ないのではないでしょうか。

スクリプトを作っていく

画像サイズを読む

画像サイズはgdalinfoから読めます。

my($w) = -1;
my($h) = -1;
my($file) = ''; # どこかでファイル名を取得してここに納める
...
open(P,"gdalinfo ${file} | grep 'Size is' | sed 's/Size is //' | sed 's/,//'|");
while(<P>) {
  chomp;
  ($w, $h) = split;
}
if( !($w > 0 && $h > 0) ) {
  die "Cannot get size of ${file}\n";
}

gdalinfoを呼び出して、"Size is "となっている行のみ抜き出し、"Size is "を切り捨て、","を切り捨てると、"(数字) (数字)"の文字列になるので、これをsplitで切ってます。

最大ズームレベルを得る

幅を256ピクセルにするか高さを256ピクセルにするかは、幅または高さのうちピクセル数の多いほうを256ピクセルにします。

次に示すのは、幅優先とした場合です。

$sh_world = 20037508.34;
...
  $maxz = ceil(log($w)/log(2)) - 8;
  $w_maxz = 256 * (1<<$maxz);
  $h_maxz = floor($w_maxz / $w * $h);
  $wh_world = $sh_world;
  $hh_world = $sh_world / $w * $h;
...

log($w)/log(2)は、幅ピクセル数が2の何乗(整数とは限らない)になるかを見ています。さらにceil(log($w)/log(2))''で、log($w)/log(2)``以上の最小値を得ています。

後は、ズームレベルが$maxzの時のピクセル数($w_maxs'', ''$h_maxz'')や、画像の「メートル単位」のサイズの半分の数($wh_world,$hh_world``)を計算しています。

空間参照系を付与する

my($file); # 上で説明済み
my($f_tmps); # テンポラリの出力ファイル名を入れておく
...
  my($cmd) = "gdal_translate -of VRT ".
    "-a_nodata \"255 255 255 0\" ".
    "-a_srs epsg:3857 ".
    "-a_ullr -${wh_world} ${hh_world} ${wh_world} -${hh_world} ".
    "\"${file}\" ".
    "\"${f_tmps}\"";

  system($cmd);
...

これで、$f_tmpsで示されたファイルに、「適切な」空間参照系を付与されたVRT形式のファイル(実際には何も処理していない)として出力しています。

また、NODATAの組み合わせを (R,G,B,A)=(255,255,255,0)にしています。

指定したズームレベルに対応した画像を作成する

  # $f_tmpsはジオリファレンスを与えたVRTファイル(上述、入力)
  # $f_tmpzは拡大したVRTファイル(出力)
  # $zは、ズームレベル
  my($w_z) = int($w_maxz / (1<<($maxz-$z)));
  my($h_z) = int($h_maxz / (1<<($maxz-$z)));

  $cmd = "gdalwarp -of VRT ".
    "-dstnodata \"255 255 255 0\" ".
    "-r lanczos ".
    "-ts $w_z $h_z ".
    "\"${f_tmps}\" ".
    "\"${f_tmpz}\" ";
  system($cmd);

実際にはVRTなので作成していません。

VRTファイル?

VRTファイルの詳細については http://www.gdal.org/gdal_vrttut.html 参照。

元ファイルがこれで、元ファイルに対してこういうオフセットを取って、これだけ拡大して、空間参照系はこうやって、…といった計算をgdal_translateとかで行いますが、実際のリサンプリングは行わずに、操作記録だけを持たせているファイルです。

もちろん、最終成果物を出すまでに一度以上は実際のラスタ画像を吐き出す必要があり、その際はそれなりの時間がかかります。

ただ、VRTを中間ファイルで使うと、うれしいこともあります。

  • GeoTIFF生成の際にNODATAの組み合わせを(R,G,B,A)=(255,255,255,0)に設定しても、生成されるGeoTIFF画像は(R,G,B,A)=(0,0,0,0)になってたりしますが、そういうことがありません。
  • とにかく速いです。元ファイルへのリンクと操作の記録を書き込んでるだけですから。
  • どのGDALツールでも入出力を受けてくれます。PNGを吐けないgdalwarpの結果をPNGで得たい場合、いったんVRTで吐かせて、gdal_translateでPNGに変換する(ここで実際のwarpが実行される)ことで目的を達成できます。

PNGで吐く

VRTは実際には何もリサンプリング処理等をしていないので、ここでやってしまいましょう。

ただし、上述のNODATA問題がGeoTIFFにはあるので、あえてPNGにしちゃいましょう。

  # $f_tmpzは拡大したVRTファイル(上述、入力)
  # $f_tmpzは拡大したPNGファイル(出力)
  $cmd = "gdal_translate -of PNG \"${f_tmpz}\" \"${f_out}\"";
  system($cmd);

タイル化する

当然ながらgdal2tiles.pyを使いますが、最大ズームレベルだけをnearestで分けます。これは既に最大ズームレベルとほぼ同じ解像度でリサンプリングしているため、無駄な補完計算を省略するためです。

最大ズームレベル未満は、cubic補完で縮小しつつタイル化します。

最大ズームレベル

    $cmd = "gdal2tiles.py -a '255,255,255,0' -z $maxz \"${f_out}\"";
    ...
  system($cmd);

最大ズームレベル未満

    my($zarg);
    if( $z == 0 ) {
      $zarg = '0';
    }
    else {
      $zarg = "0-${z}";
    }
    $cmd = "gdal2tiles.py -r cubic -z ${zarg} \"${f_out}\"";
  ...
  system($cmd);

HTML,JavaScriptを吐く

HTML,JavaScriptでは、Webメルカトルが完全に消えています。このスクリプトで前提としているのは、256x256ピクセルのタイルをズームレベル0とし、ズームレベルが1増えるごとに2倍に拡大され、それにあわせてタイル番号が付されている、ということだけです。それがWebメルカトルとか、UTMとか、まともな空間参照系のことは考えていません。

たいたい、ふつうのスナップ写真ひとつが全世界を覆うだなんて、まともではありません。上で「Webメルカトルを使います」とか言ってたのを自己全否定しているようですが、上の記述はあくまで方便であって、タイルができればそれでいいのです。

HTML, JavaScriptを出力する部分は次のようになります。

  # $maxzは最大ズームレベル
  # $w_maxz, $h_maxzは、最大ズームレベルの画像ピクセルサイズ
  # $fbodyは、ファイル名のプリフィックス部(ディレクトリ名に使用)
  my($maxz, $w_maxz, $h_maxz, $fbody);
  # Wはファイルハンドル
  # 
  print W << "_EOL_";
<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Strict//EN" "http://www.w3.org/TR/xhtml1/DTD/xhtml1-strict.dtd">
<html xmlns="http://www.w3.org/1999/xhtml>"
<head>
<title>$fbody</title>
<style type="text/css">
html, body, #MAP {
  margin: 0;
  padding: 0;
  position: absolute;
  width: 100%;
  height: 100%;
  overflow: hidden;
}

.olImageLoadError {
  display: none;
}

</style>
<script src="http://www.openlayers.org/api/OpenLayers.js" type="text/javascript"></script>
<script type="text/javascript">
// specific
var w_maxz = $w_maxz;
var h_maxz = $h_maxz;
var maxz = $maxz;
var s_world = 256 * (1<<maxz);

var w_maxz_h = w_maxz / 2;
var h_maxz_h = h_maxz / 2;
var s_world_h = s_world / 2;

var map;
var mapBounds = new OpenLayers.Bounds(-w_maxz_h,-h_maxz_h,w_maxz_h,h_maxz_h);

OpenLayers.Util.onImageLoadErrorColor = "transparent";

function init(){
  var options = {
    controls: [],
    maxExtent: new OpenLayers.Bounds(-s_world_h, -s_world_h, s_world_h, s_world_h),
    maxResolution: s_world/256,
//    restrictedExtent: mapBounds,
    numZoomLevels: maxz + 1
  };
  map = new OpenLayers.Map('MAP', options);
  var layer = new OpenLayers.Layer.TMS(
    "TMS Layer", "", {
      url: '', serviceVersion: '.', layername: '.', alpha: true,
      type: 'jpg'
    }
  );
  map.addLayer(layer);
  map.zoomToExtent( mapBounds );  
  map.addControl(new OpenLayers.Control.PanZoom());
  map.addControl(new OpenLayers.Control.MousePosition());
  map.addControl(new OpenLayers.Control.Navigation());
}


function resize() {  
  var map = document.getElementById("MAP");
  map.style.width = "100%";
  map.style.height = "100%";
}
</script>
</head>
<body onload="init()">
<div id="MAP"></div>
<script type="text/javascript">resize();</script>
</body>
</html>
_EOL_

完全版ソース

好きにお使い下さい、使えるものならなー。

#!/usr/bin/perl

use POSIX qw(ceil floor);

my($w) = -1;
my($h) = -1;
# $w = 73149;
# $h = 9896;
my($file) = '';

foreach my $arg (@ARGV) {
  if( $arg eq '-h' || $arg eq '--help' ) {
    print_help();
    exit 0;
  }
  if( $file eq '' ) {
    $file = $arg;
  }
}

if( $file eq '' ) {
  print_help();
  exit 64;
}

sub print_help {
  print STDERR "$0: <inputfilename>\n";
}

open(P,"gdalinfo ${file} | grep 'Size is' | sed 's/Size is //' | sed 's/,//'|");
while(<P>) {
  chomp;
  ($w, $h) = split;
}

if( !($w > 0 && $h > 0) ) {
  die "Cannot get size of ${file}\n";
}

print STDERR "Size is $w $h\n";

my(@fsegs) = split(/\//, $file);
my($fbody) = $fsegs[$#fsegs];
$fbody =~ s/\.[^\.]+$//;


$sh_world = 20037508.34;

if( $w >= $h ) {
  $maxz = ceil(log($w)/log(2)) - 8;
  $w_maxz = 256 * (1<<$maxz);
  $h_maxz = floor($w_maxz / $w * $h);
  $wh_world = $sh_world;
  $hh_world = $sh_world / $w * $h;
}
else {
  $maxz = ceil(log($h)/log(2)) - 8;
  $h_maxz = 256 * (1<<$h_maxz);
  $w_maxz = floor($h_maxz / $h * $w);
  $hh_world = $sh_world;
  $wh_world = $sh_world / $h * $w;
}

# setenv GDAL_CACHEMAX 8192

if( !-d "tmp" ) {
  mkdir("tmp") || die "tmp: $!\n";
}

my($f_tmps) = "tmp/s-${fbody}.vrt";

# --------

mksrc($file, $f_tmps, $wh_world, $hh_world);
warp($maxz, $maxz, $w_maxz, $h_maxz, $fbody, $f_tmps);
if( $maxz >= 1 ) {
  warp($maxz-1, $maxz, $w_maxz, $h_maxz, $fbody, $f_tmps);
}
png2jpg($fbody);
write_html($maxz, $w_maxz, $h_maxz, $fbody);

# ----------------

# georeferences and converts into VRT
sub mksrc {
  my($file, $f_tmps, $wh_world, $hh_world) = @_;
  my($cmd) = "gdal_translate -of VRT ".
    "-a_nodata \"255 255 255 0\" ".
    "-a_srs epsg:3857 ".
    "-a_ullr -${wh_world} ${hh_world} ${wh_world} -${hh_world} ".
    "\"${file}\" ".
    "\"${f_tmps}\"";
  &exec($cmd);
}

sub png2jpg {
  my($fbody) = @_;
  system("find \"$fbody\" -name '*.png' -exec ./convert.sh {} \\; -print");
}

sub write_html {
  my($maxz, $w_maxz, $h_maxz, $fbody) = @_;
  my($htmlpath) = $fbody.'/index.html';
  open(W, ">$htmlpath");
  print W << "_EOL_";
<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Strict//EN" "http://www.w3.org/TR/xhtml1/DTD/xhtml1-strict.dtd">
<html xmlns="http://www.w3.org/1999/xhtml>"
<head>
<title>$fbody</title>
<style type="text/css">
html, body, #MAP {
  margin: 0;
  padding: 0;
  position: absolute;
  width: 100%;
  height: 100%;
  overflow: hidden;
}

.olImageLoadError {
  display: none;
}

</style>
<script src="http://www.openlayers.org/api/OpenLayers.js" type="text/javascript"></script>
<script type="text/javascript">
// specific
var w_maxz = $w_maxz;
var h_maxz = $h_maxz;
var maxz = $maxz;
var s_world = 256 * (1<<maxz);

var w_maxz_h = w_maxz / 2;
var h_maxz_h = h_maxz / 2;
var s_world_h = s_world / 2;

var map;
var mapBounds = new OpenLayers.Bounds(-w_maxz_h,-h_maxz_h,w_maxz_h,h_maxz_h);

OpenLayers.Util.onImageLoadErrorColor = "transparent";

function init(){
  var options = {
    controls: [],
    maxExtent: new OpenLayers.Bounds(-s_world_h, -s_world_h, s_world_h, s_world_h),
    maxResolution: s_world/256,
//    restrictedExtent: mapBounds,
    numZoomLevels: maxz + 1
  };
  map = new OpenLayers.Map('MAP', options);
  var layer = new OpenLayers.Layer.TMS(
    "TMS Layer", "", {
      url: '', serviceVersion: '.', layername: '.', alpha: true,
      type: 'jpg'
    }
  );
  map.addLayer(layer);
  map.zoomToExtent( mapBounds );  
  map.addControl(new OpenLayers.Control.PanZoom());
  map.addControl(new OpenLayers.Control.MousePosition());
  map.addControl(new OpenLayers.Control.Navigation());
}


function resize() {  
  var map = document.getElementById("MAP");
  map.style.width = "100%";
  map.style.height = "100%";
}
</script>
</head>
<body onload="init()">
<div id="MAP"></div>
<script type="text/javascript">resize();</script>
</body>
</html>
_EOL_
  close(W);
}

sub warp {
  my($z,$maxz,$w_maxz,$h_maxz,$fbody,$f_tmps) = @_;

  my($w_z) = int($w_maxz / (1<<($maxz-$z)));
  my($h_z) = int($h_maxz / (1<<($maxz-$z)));

  my($f_outdir) = "tmp/${z}/";
  my($f_out) = $f_outdir."${fbody}.png";
  my($f_tmpz) = "tmp/${z}-${fbody}.vrt";

  my($cmd);

  if( !-d $f_outdir ) {
    mkdir($f_outdir) || die "${f_outdir}: $!\n";
  }

  # warps
  $cmd = "gdalwarp -of VRT ".
    "-dstnodata \"255 255 255 0\" ".
    "-r lanczos ".
    "-ts $w_z $h_z ".
    "\"${f_tmps}\" ".
    "\"${f_tmpz}\" ";
  &exec($cmd);

  # to PNG
  $cmd = "gdal_translate -of PNG \"${f_tmpz}\" \"${f_out}\"";
  &exec($cmd);

  # tile
  if( $z == $maxz ) {
    $cmd = "gdal2tiles.py -a '255,255,255,0' -z $maxz \"${f_out}\"";
  }
  else {
    my($zarg);
    if( $z == 0 ) {
      $zarg = '0';
    }
    else {
      $zarg = "0-${z}";
    }
    $cmd = "gdal2tiles.py -r cubic -z ${zarg} \"${f_out}\"";
  }
  &exec($cmd);
}

sub exec {
  my($cmd) = @_;
  print STDERR "--------\n$cmd\n";
  my($stat) = system($cmd);
  if( $stat != 0 ) {
    exit($stat);
  }
}

おわりに

無理にジオリファレンスを取ることで、gdalの外部プログラムを駆使してタイル画像を作ることに成功しました。

中身はネタのようなスクリプトをご紹介しましたが、実際にタイルとHTMLが作成できるから不思議に思うかもしれません。

しかし、ここまでさんざん「ネタ」と言ってますが、gdal_trasnlate, gdalwarp, VRT, GDAL_CACHEMAX を説明しているので、決して「ネタ」だけで終わってないので、まあハードルを下げた記事であった程度と認識頂ければありがたいです。

ためしにやってみた

最後の最後に、サンプルのスナップショットを出しておきます。

1.png

2.png