ラインやポリゴンのデータもストリーム処理でベクトルタイルに変換する。Rubyで。

ラインやポリゴンのデータもストリーム処理でベクトルタイルに変換するプログラムを作成したので共有する。Rubyでgeoruby-ext（特にRGeo経由でGEOS）を使う。

使い方

$ ruby convert.rb [theme] [z] [shapefiles...]

• theme は、{z}/{x}/{y}.geojson が置かれるフォルダ名を想定している。例えば std という文字列を与えると、ベクトルタイルは std/{z}/{x}/{y}.geojson というフォルダに置かれることになる。
• z は、ズームレベルを指定する。例えば 14 を指定する。
• shapefile は、データソースとなる Shapefile を指定する。複数のShapefileを指定すると、立て続けに map して、まとめて reduce する。結果的にデータがひとつのタイルセットにコンポジットされる。Shapefileの文字コードはShift_JISになっていることを前提としており、処理の中で UTF-8 に変換されて出力される*1。なお、出力されるベクトルタイルには UTF-8 文字列がそのままに書きだされるようにしてあり、Unicodeエスケープはされない。

convert.rbは本エントリでこれから紹介する。これまで、ストリーム処理でMapReduce的にベクトルタイル変換する方法を実現するために、map.rbとreduce.rbを別々に作成してきたが、今回は小手先の工夫を加えて、一つのスクリプトでmapとreduceの両方が処理されるようにした。また、UNIXパイプを作って sort や reduce.rb につなぐ部分も、すべて convert.rb の中におさめてしまっている。結果として、上記のようなコマンドラインで処理できるようにした。
コマンドラインパラメータが適切でないときの処理は、現在のところほとんどまともでない。

利用例は、例えば次の通りである。

$ ruby convert.rb bus_stop 14 bus_stop.shp

ソースコード

上記の convert.rb のソースコードはこの節で示す通りである。georuby-ext が導入されていることが前提になるので、必要に応じて、

$ gem install georuby-ext

する。それ以外には特段、特別なライブラリーは必要としない。（georuby-ext が RGeo に依存し、RGeo が GEOS（GEOS = Geometry Engine, Open Source）に依存するので、結果的に GEOS が導入される。）あと、内部で sort コマンドを呼ぶので、UNIX で普通に存在するような sort が存在しない環境では動かないかもしれない。なお、将来、タイル化に並列分散処理を導入する場合には、sort コマンドのかわりに MapReduce 処理系（Hadoop等）を入れることになるのではないだろうか。

# coding: utf-8
require 'georuby-ext' # gem install georuby-ext
require 'geo_ruby/shp'
require 'fileutils'
require 'json'
include GeoRuby::Shp4r

module Math
  def self.sec(x)
    1.0 / cos(x)
  end
end

module XYZ
  def self.pt2xy(pt, z)
    lnglat2xy(pt.x, pt.y, z)
  end

  def self.lnglat2xy(lng, lat, z)
    n = 2 ** z
    rad = lat * 2 * Math::PI / 360
    [n * ((lng + 180) / 360),
      n * (1 - (Math::log(Math::tan(rad) +
        Math::sec(rad)) / Math::PI)) / 2]
  end

  def self.xyz2lnglat(x, y, z)
    n = 2 ** z
    rad = Math::atan(Math.sinh(Math::PI * (1 - 2.0 * y / n)))
    [360.0 * x / n - 180.0, rad * 180.0 / Math::PI]
  end

  def self.xyz2envelope(x, y, z)
    GeoRuby::SimpleFeatures::Envelope.from_points([
      GeoRuby::SimpleFeatures::Point.from_coordinates(
        xyz2lnglat(x, y, z)),
      GeoRuby::SimpleFeatures::Point.from_coordinates(
        xyz2lnglat(x + 1, y + 1, z))])
  end
end

class GeoRuby::SimpleFeatures::Geometry
  def tile(z)
    lower = XYZ::pt2xy(self.bounding_box[0], z)
    upper = XYZ::pt2xy(self.bounding_box[1], z)
    rg = self.to_rgeo
    lower[0].truncate.upto(upper[0].truncate) {|x|
      upper[1].truncate.upto(lower[1].truncate) {|y|
        env = XYZ::xyz2envelope(x, y, z).to_rgeo
        intersection = rg.intersection(env)
        if intersection.is_empty?
          $stderr.print "e"
          next
        end
        if intersection.respond_to?(:each)
          intersection.each {|g|
            $stderr.print "m"
            yield x, y, JSON.parse(g.to_georuby.as_geojson)
          }
        else
          $stderr.print "s"
          yield x, y, JSON.parse(intersection.to_georuby.as_geojson)
        end
      }
    }
  end
end

def write(geojson, tzxy)
  path = tzxy.join('/') + '.geojson'
  print "writing #{geojson[:features].size} features to #{path}\n"
  [File.dirname(path)].each {|d| FileUtils.mkdir_p(d) unless File.exist?(d)}
  File.open(path, 'w') {|w| w.print(JSON.dump(geojson))}
end

def map(t, z, paths)
  IO.popen("sort | ruby #{__FILE__}", 'w') {|io|
    fid = 0
    paths.each {|path|
      ShpFile.open(path) {|shp|
        shp.each{|r|
          fid += 1
          prop = r.data.attributes
          prop[:fid] = fid
          $stderr.print "[#{fid}]"
          prop.each {|k, v|
            if(v.class == String)
              prop[k] = v.encode('UTF-8', 'CP932')
            end
          }
          r.geometry.tile(z) {|x, y, g|
            f = {:type => 'Feature', :geometry => g, :properties => prop}
            io.puts([t, z, x, y, JSON.dump(f)].join("\t") + "\n")
          }
        }
      }
    }
    $stderr.print "\n"
  }
end

def reduce
  last = nil # to extend the scope of the variable
  geojson = nil # to extend the scope of the variable
  while gets
    r = $_.strip.split("\t")
    current = r[0..3]
    if current != last
      write(geojson, last) unless last.nil?
      geojson = {:type => 'FeatureCollection', :features => []}
    end
    geojson[:features] << JSON.parse(r[4])
    last = current
  end
  write(geojson, last)
end

def help
  print <<-EOS
  ruby convert.rb {theme} {z} {shapefiles...}
  EOS
end

if __FILE__ == $0
  if ARGV.size == 0
    reduce
  else
    begin
      # ruby convert.rb test 10 /Users/hfu/Downloads/A13-11_42_GML/a001420020120310.shp
      map(ARGV[0], ARGV[1].to_i, ARGV[2..-1])
    rescue
      p $!
      help
    end
  end
end

上記とまったく同じものを、gistにも置いておく。

確認用 HTML ファイル

今回のツールで、[theme] を test に設定*2して作成したベクトルタイルを確認するためのHTMLファイル index.html のサンプルを共有する。

<!doctype html>
<html>
  <head>
    <meta charset="UTF-8">
    <title></title>
    <link rel="stylesheet" href="http://cdn.leafletjs.com/leaflet-0.7.2/leaflet.css"/>
    <script src="http://cdn.leafletjs.com/leaflet-0.7.2/leaflet.js"></script>
    <script src="http://handygeospatial.github.io/mapsites/js/leaflet-hash.js"></script>
    <script src="http://handygeospatial.github.io/mapsites/js/TileLayer.GeoJSON.js"></script>
    <style>
      body {padding: 0; margin: 0}
      html, body, #mapdiv {height: 100%; width: 100%;}
      .leaflet-container {background: #fff;}
    </style>
  </head>
  <body>
    <div id="mapdiv">
    <script>
      var std = L.tileLayer(
        'http://cyberjapandata.gsi.go.jp/xyz/std/{z}/{x}/{y}.png',
        {attribution: "地理院タイル（標準地図）" });
      var vec = new L.TileLayer.GeoJSON(
        './test/{z}/{x}/{y}.geojson',
        {attribution: 'vector tile',
         minZoom: 7, maxZoom: 18,
         unique: function(feat) {return feat.properties['fid']}},
        {style: {"clickable": true, "stroke": false, "fillColor": "#3f3"},
         onEachFeature: function(feature, layer) {
           layer.bindPopup(JSON.stringify(feature.properties));
        }});

      var map = L.map('mapdiv', {
        center: [32.9263, 130.0378], zoom: 12,
        layers: [std, vec]});

      var hash = L.hash(map);
      L.control.layers({'地理院タイル（標準地図）': std},
        {'vector tile': vec}).addTo(map);
    </script>
  </body>
</html>

これを、データを作った場所と同じフォルダにおいて、例えば python の SimpleHTTPServer を用いてサーブしてブラウザで確認すれば良い。

$ python -m SimpleHTTPServer

通常 http://localhost:8000/ で表示されることになる。なお、ベクトルタイルデータすべてについてそうであるが、ベクトルタイルデータは CORS（cross-origin resource sharing）されるか又は同一オリジンから参照されることが必要なので、file: スキームでは通常表示されないことから、SimpleHTTPServer を使う必要がある。

検証状況

色塗り系ポリゴンファイルを変換してみて、上記の確認用HTMLを用い、問題ない位置に表示されることを確認している。ループの外に出せる to_rgeo 呼び出しを適切にループの外に置くことにより、ナイーブな実装にくらべてそれなりの高速化を実現したりしている。

所感

やや未熟なところはあるが、読んでくださる方が読んでくだされば、基本的なアイディアはすべて伝わると考える。

この先は、ocra など使って「普通のパソコン」でも使って頂ける形を追求する線と、速度を求めて並列分散処理を追求する線に分化する必要が出てくる。但し、それより前に、処理速度の相場観を得るためにある程度の実戦経験を積む段階が必要になる。

*1:状況が違う場合には、convert.rbの該当場所を書き換えること。

*2:状況が違う場合には、index.htmlの該当場所を書き換えること。