基盤地図情報（数値標高モデル）（5mメッシュ（標高））のGeoJSONストリームへのコンバータ

基盤地図情報（数値標高モデル）（5mメッシュ（標高））のGeoJSONストリームへのコンバータをつくった。かなりひねくれたコードだ。

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

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

  def Meshcode::height(meshcode)
    case meshcode.size
    when 8
      30.0 / 60 / 60
    else
      raise 'not implemented.'
    end
  end

  def Meshcode::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 DEM5A
  N_LNG = 225
  N_LAT = 150
  D_LNG = 1.0 / 80 / 225
  D_LAT = 2.0 / 3 / 80 / 150
  def DEM5A::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
        (x, y) = [count % N_LNG, count / N_LNG]
        lng = left + D_LNG * (x + 0.5)
        lat = top - D_LAT * (y + 0.5) 
        (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 JSON::dump(f), "\n"
        count += 1
      end
    }
  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}
  case type
  when 'DEM5A'
    Kernel.const_get(type)::parse(params)
  else
    # print "converter for #{type} not implemented.\n"
  end
}

これで、JSONLっぽく、一行一レコードでデータが出てくる。QGISとかで読み込めるgeojsonファイルにするには、次のスクリプトにパイプする。

print <<-EOS
{"type": "FeatureCollection",
 "features": [
EOS
buf = ''
while gets
  print buf, ",\n" unless buf == ''
  buf = $_.strip 
end
print buf
print "]}\n"

具体的には、次のようなコマンドになる。

$ ruby convert.rb FG-GML-5339-45-DEM5A/FG-GML-5339-45-22-DEM5A-20130702.xml | ruby fc.rb > 53394522.geojson

行指向なので、convert.rb のコマンドライン引数に複数の xml ファイルを指定しても OK である。

実用的には、convert.rb と fc.rb の間にタイル化（タイル番号での Map & Reduce）が入る。これで、ポイントクラウドタイルができることになるはず。