世界の測量

Sibling of "Relevant, Timely, and Accurate, " but much lighter and shorter ※自らの所属する組織の見解を示すものでない

MapReduce的ベクトルタイル作成(地価公示)Pt. 2

ついでreduce.rbを紹介し、実際にタイルを作成する。

reduce.rb

ソートされ、タイルごとに並んでいるデータをTileJSONに書き出すだけの簡単なプログラム。

require 'fileutils'
require 'json'

last = nil
geojson = nil
$n_tiles = 0

def write(geojson, zxy)
  path = zxy.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))}
  $n_tiles += 1
end

while gets
  r = $_.strip.split("\t")
  current = r[0..2].map{|v| v.to_i}
  if current != last
    write(geojson, last) unless last.nil?
    geojson = {:type => 'FeatureCollection', :features => []}
  end
  geojson[:features] << JSON.parse(r[3])
  last = current
end
write(geojson, last)
print "finished. #{$n_tiles} tiles written.\n"

map.rb(多少改造)

前回紹介したmap.rbについて、ズームレベルをコマンドライン指定できるよう多少改造。

require 'find'
require 'json'
require 'rubygems'
require 'geo_ruby'
require 'geo_ruby/shp'
include GeoRuby::Shp4r

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

def process(r, z)
  n = 2 ** z
  lon_deg = r.geometry.x
  lat_rad = r.geometry.y * 2 * Math::PI / 360
  x = n * ((lon_deg + 180) / 360)
  y = n *
    (1 - (Math::log(Math::tan(lat_rad) + Math::sec(lat_rad)) / Math::PI)) / 2
  prop = {}
  r.data.attributes.each {|k, v|
    prop[k] = case v
    when 'true'
      true
    when 'false'
      false
    else
      s = v.encode('UTF-8', 'CP932')
      /^\d*$/.match(s) ? s.to_i : s
    end
  }
  f = {:type => 'Feature', :geometry => {:type => 'Point', 
    :coordinates => [r.geometry.x, r.geometry.y]}, 
    :properties => prop}
  print [z, x.to_i, y.to_i, JSON.dump(f)].join("\t"), "\n"
end

ShpFile.open('L01-13.shp') {|shp|
  shp.each {|r|
    ARGV.each {|z| process(r, z.to_i)}
  }
}

補足であるが、x,yの算出式については、 Slippy map tilenames - OpenStreetMap Wiki が出典である。

実行

ズームレベル14のタイルを作成するには、このようにする。

$ ruby map.rb 14 | sort | ruby reduce.rb

手元のMacBook Airの場合、21秒程度で9373タイルが生成された。