基盤地図情報（数値標高モデル）のベクトルタイルへの変換の続き。

基本方針

基本方針はこれまでと変わらず、独自形式データをzxyつきのGeoJSONLストリームにし、MapReduce的に sort して reduce する。

mapper

DEM5AのGMLデータをGeoJSONLに変換するプログラムについて、少し実装を改善しつつ zxy を付与するようにした。

# coding: utf-8
# convert.rb
require 'json'

# Meshcode is probably Japanese English.
module Meshcode
  def self.width(code)
    case code.size
    when 8
      45.0 / 60 / 60
    else
      raise 'not implemented.'
    end
  end

  def self.height(code)
    case code.size
    when 8
      30.0 / 60 / 60
    else
      raise 'not implemented.'
    end
  end

  def self.lefttop(code)
    case code.size
    when 8
      [(code[2..3].to_f + code[5].to_f / 8 + code[7].to_f / 80) + 100, 
       (code[0..1].to_f + code[4].to_f / 8 + (code[6].to_f + 1) / 80) * 2 / 3]
    else
      raise 'not implemented.'
    end
  end
end

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

class DEM
  def parse(params)
    (left, top) = Meshcode::lefttop(params[:meshcode])
    skip = true
    count = 0
    File.foreach(params[:path], encoding: 'cp932') {|l|
      if l.include?('<gml:tupleList>')
        skip = false
        next
      elsif l.include?('</gml:tupleList>')
        skip = true
        next
      elsif !skip
        (i, j) = [count % @n_lng, count / @n_lng]
        lng = left + @d_lng * (i + 0.5)
        lat = top - @d_lat * (j + 0.5)
        lat_rad = lat * 2 * Math::PI / 360
        n = 2 ** params[:z]
        x = n * ((lng + 180) / 360)
        y = n *
          (1 - (Math::log(Math::tan(lat_rad) + Math::sec(lat_rad)) /
          Math::PI)) / 2
        x = x.to_i
        y = y.to_i
        (type, height) = l.encode('UTF-8').strip.split(',')
        f = {:type => 'Feature', 
          :geometry => {:type => 'Point', :coordinates => [lng, lat]},
          :properties => {:type => type, :height => height.to_f}}
        print [params[:z], x, y, JSON::dump(f)].join("\t"), "\n"
        count += 1
      end
    }
  end
end

class DEM5A < DEM
  def initialize
    @n_lng = 225
    @n_lat = 150
    @d_lng = 1.0 / 80 / 225
    @d_lat = 2.0 / 3 / 80 / 150
  end
end

ARGV.each {|path|
  next unless /xml$/.match path
  r = File.basename(path, '.xml').split('-')
  r.pop if r[-1].size == 4
  next unless r.shift == 'FG'
  next unless r.shift == 'GML'
  date = r.pop
  type = r.pop
  meshcode = r.join
  params = {:path => path, :type => type, :meshcode => meshcode, :z => 18}
  case type
  when 'DEM5A'
    Kernel.const_get(type).new.parse(params)
  else
    # print "converter for #{type} not implemented.\n"
  end
}

reducer

これまでつかってきた reducer の使い回し。

#reduce.rb
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"

変換実行例

$ ruby convert.rb FG-GML-5339-45-DEM5A/FG-GML-5339-45-22-DEM5A-20130702.xml | sort | ruby reduce.rb

これで、ズームレベル18のタイルとして作成される。

生成されたタイルの大きさについて分析

標準出力を見ると、タイルあたりポイント数が400〜525あたりになるようだ。525ポイント入ったものについて、サイズは76KB程度になっている。
ズームレベルを一つあげると、サイズは1/4程度になるであろうことは言うまでもない。

また、基盤地図情報DEM5Aの1zipファイルあたり生成されるz=18のタイル数は6992程度である。

TopoJSON に変換した場合

$ /usr/local/share/npm/bin/topojson 18/232765/103220.geojson > a.topojson

としてa.topojson のサイズを計測したところ21KBとなった。もとのGeoJSONデータが72KBということだから、約30%の大きさになることになる。

DEM5A2GeoJSONTile