Qiita Teams that are logged in
You are not logged in to any team

Log in to Qiita Team
Community
OrganizationAdvent CalendarQiitadon (β)
Service
Qiita JobsQiita ZineQiita Blog
Help us understand the problem. What is going on with this article?

DEM5A2GeoJSONTile

More than 5 years have passed since last update.

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

基本方針

基本方針はこれまでと変わらず、独自形式データを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%の大きさになることになる。

handygeospatial
どのエントリにも geo や web のタグはつけません。 すべてのエントリが geo で web なことに関係することとなるであろうためです。 ※自らの所属する組織の見解を示すものでない
http://hfu.hatenablog.com
Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away