標高の入った GeoTIFF から標高PNGタイルをつくるまで

はじめに

標高の入った GeoTIFF から 地理院仕様の標高 PNG タイル をつくるまでの手順をまとめます。おそらくもっといい方法があるのだと思いますが、記録として。

1. 投影変換

最初に gdalwarp で GeoTIFF を EPSG:3857 に変換しておきます。

$ gdalwarp -t_srs epsg:3857 -r bilinear input.tiff output.tiff

• -r には画像の拡大縮小のためのアルゴリズムを指定します。ここでは bilinear を指定。

2. タイル生成

zxy ごとに以下の手順を必要なだけ繰り返します

2.1 クリップ

gdal_translate で必要な領域を切り抜きます

$ gdal_translate -outsize 256 256 -r bilinear -projwin [left] [top] [right] [bottom] input.tiff [zxy].tiff

• -outsize にはクリップ画像のサイズを指定します。地理院標高PNGタイルは256x256なのでこの指定です。
• -r には画像の拡大縮小のためのアルゴリズムを指定します。ここでは bilinear を指定。
• -projwin にはクリップ領域を EPSG:3857 の座標系で指定します

2.2. 着色とPNG化

gdaldem の color-relief を使って、段彩図を作る要領で GeoTIFF の各ピクセルの値を色に変換していきます

$ gdaldem color-relief -nearest_color_entry -of png [zxy].tiff color.txt [zxy].png

• color.txt は値と色の対応を格納したファイルです
• -nearest_color_entry オプションによって、color.txt で定義された値・色に寄せる形で色が選択されます。標高PNGでは中間色を作られると意味がないので、これは重要です。
• -of オプションでは出力フォーマットを指定します。標高PNGなのでこれは png 固定です。

3.補足

クリップ領域を確定したり座標・色の対応ファイルを作ったりなど、以下のような使い捨てのコードを作って対応しました。

zxy と epsg:3857 座標系の対応

とりあえず min_lat,min_lon,max_lat,max_lon,min_z,max_z を入力として、zxy と EPSG:3857の min_x,min_y,max_x,max_y を返すものを Leaflet で作ったので置いておきます。Proj4js + node.js でリライトしたほうがいいでしょうけど。

index.html

  <script>
    var map = L.map("map", {
      maxZoom: 18
    });
    L.tileLayer("https://cyberjapandata.gsi.go.jp/xyz/std/{z}/{x}/{y}.png").addTo(map);

    var points = [];
    points.push(L.latLng(35.6607278, 139.2369119));
    points.push(L.latLng(35.6607278, 139.2813654));
    points.push(L.latLng(35.6199538, 139.2813654));
    points.push(L.latLng(35.6199538, 139.2369119));

    var bounds = L.latLngBounds(points);
    map.fitBounds(bounds);
    L.polygon(points).addTo(map);

    for (var z = 14; z <= 17; z++) {
      var nw = map.project(bounds.getNorthWest(), z).divideBy(256).floor();
      var se = map.project(bounds.getSouthEast(), z).divideBy(256).floor();
      for (var y = nw.y; y <= se.y; y++) {
        for (var x = nw.x; x <= se.x; x++) {
          var box = L.latLngBounds([
            map.unproject(L.point(x, y).multiplyBy(256), z),
            map.unproject(L.point(x + 1, y).multiplyBy(256), z),
            map.unproject(L.point(x + 1, y + 1).multiplyBy(256), z),
            map.unproject(L.point(x, y + 1).multiplyBy(256), z)
          ]);
          L.rectangle(box, {
            color: "red",
            weight: 1
          }).addTo(map);
          L.marker(box.getCenter()).addTo(map);

          var p1 = L.CRS.EPSG3857.project(box.getNorthWest());
          var p2 = L.CRS.EPSG3857.project(box.getSouthEast());

          var win = [p1.x, p1.y, p2.x, p2.y].join(" ");
          var file = [z, x, y].join("_");
          console.log("gdal_translate -outsize 256 256 -r bilinear -projwin " + win + " alldem_1m_3857.tif " + file + ".tif");
          console.log("gdaldem color-relief -nearest_color_entry -of png " + file + ".tif dem.mini.txt " + file + ".png");
        }
      }
    }
  </script>

実際にはこんな雰囲気のコマンド群が出力されるので、これをシェルに貼って実行、という雑な対応です。

out.txt

gdal_translate -outsize 256 256 -r bilinear -projwin 15497760.358876059 4256013.734918613 15500206.343781184 4253567.750013488 alldem_1m_3857.tif 14_14528_6452.tif
gdaldem color-relief -nearest_color_entry -of png 14_14528_6452.tif dem.mini.txt 14_14528_6452.png
gdal_translate -outsize 256 256 -r bilinear -projwin 15500206.343781184 4256013.734918613 15502652.32868631 4253567.750013488 alldem_1m_3857.tif 14_14529_6452.tif
gdaldem color-relief -nearest_color_entry -of png 14_14529_6452.tif dem.mini.txt 14_14529_6452.png
gdal_translate -outsize 256 256 -r bilinear -projwin 15502652.32868631 4256013.734918613 15505098.313591436 4253567.750013488 alldem_1m_3857.tif 14_14530_6452.tif
gdaldem color-relief -nearest_color_entry -of png 14_14530_6452.tif dem.mini.txt 14_14530_6452.png
gdal_translate -outsize 256 256 -r bilinear -projwin 15497760.358876059 4253567.750013488 15500206.343781184 4251121.765108366 alldem_1m_3857.tif 14_14528_6453.tif
gdaldem color-relief -nearest_color_entry -of png 14_14528_6453.tif dem.mini.txt 14_14528_6453.png
gdal_translate -outsize 256 256 -r bilinear -projwin 15500206.343781184 4253567.750013488 15502652.32868631 4251121.765108366 alldem_1m_3857.tif 14_14529_6453.tif
gdaldem color-relief -nearest_color_entry -of png 14_14529_6453.tif dem.mini.txt 14_14529_6453.png
gdal_translate -outsize 256 256 -r bilinear -projwin 15502652.32868631 4253567.750013488 15505098.313591436 4251121.765108366 alldem_1m_3857.tif 14_14530_6453.tif
gdaldem color-relief -nearest_color_entry -of png 14_14530_6453.tif dem.mini.txt 14_14530_6453.png
...

color.txt

適当な既存のスクリプトがあるのかもしれませんが、こんな node.js スクリプトを作って作成しました。gdaldem は color.txt のサイズが大きいと遅くなるので、min と max の範囲の標高値に対応する色だけを出すようにしています。

main.js

var a = Math.pow(2, 23);

var min = 180;
var max = 410;

for (var r = 0; r < 0x100; r++) {
  for (var g = 0; g < 0x100; g++) {
    for (var b = 0; b < 0x100; b++) {
      var x = r * 0x10000 + g * 0x100 + b;
      if (x < a) {
        var h = x * 0.01;
        if (min < h && h < max) console.log(h, r, g, b);
      } else if (x > a) {
        var h = (x - Math.pow(2, 24)) * 0.01;
        if (min < h && h < max) console.log(h, r, g, b);
      } else {
        console.log("nv", r, g, b)
      }
    }
  }
}

実際の出力はこんなかんじです。

color.txt

180.01 0 70 81
180.02 0 70 82
180.03 0 70 83
180.04 0 70 84
180.05 0 70 85
180.06 0 70 86
180.07 0 70 87
180.08 0 70 88
180.09 0 70 89
180.1 0 70 90
180.11 0 70 91
...

4.まとめ

とりあえず上記の手順で目標は達成しました。

• 普通は gdal2tiles.py を使えばいいのでしょうが、gdal2tiles.py だと着色結果をリサンプリングするような手順になってしまい、標高PNGとしては意味のない中間色が作られてしまって望む結果が得られませんでした
• gdaldem は color.txt が大きいとオーバーヘッドが大きいので、処理するタイル数が多くなると実用的ではないかも
• 「epsg:3857 な TIFF を入力として標高PNGタイルを出力するスクリプト」を作っておくと役に立つかも（いつかやるかも）

追記 (20180731)

「epsg:3857 な TIFF を入力として標高PNGタイルを出力するスクリプト」

を Node.js で書くとこんなかんじです。参考程度に。

tiff2dempng.js

const fs = require("fs");
const TIFF = require("tiff");
const PNG = require("pngjs").PNG;

const src = TIFF.decode(fs.readFileSync(process.argv[2]))[0];
const dst = new PNG({
  width: src.width,
  height: src.height
});

src.data.forEach((h, i) => {
  var x = isNaN(h) ? 0x800000 : h * 100 + (h < 0 ? Math.pow(2, 24) : 0);
  dst.data.writeInt32BE((x << 8) + 0xff, i * 4);
});

dst.pack().pipe(process.stdout);