電子国土「背景地図の更新情報」を geotools.rb でハンドリング

電子国土「背景地図の更新情報」で公開されている情報を geotools.rb で受け取ってみました。

電子国土プロファイル形式のデータを読み取る

電子国土の背景地図更新状況が、下記のように公開されています。

電子国土の背景地図更新状況は、このページでご確認いただけます。（中略）背景地図の更新は随時行っていますが、更新状況は月に一度、このページでお知らせします。

http://portal.cyberjapan.jp/maprireki.htm

このページのソースを見ることで、背景地図更新状況は

top.map.openJSGIXML('http://cyberjapan.jp/koushin_rireki/200711.xml',0);

といったように、cyberjapan.jp/koushin_rireki/ の下の URL から公開されていることが分かります。地理情報標準電子国土プロファイル形式でデータを公開されているということは、多様なシステムでこのデータを利用して欲しいということでしょうし、また、データに内蔵されているメタデータを見ても、特に別システムでの利用を制限する記述はなかったので、geotools.rb でこのデータを読み取ってみることにしました。

作成したスクリプト

作成したスクリプトは、以下のようになりました。このスクリプトは、別途 http://cyberjapan.jp/koushin_rireki/200708.xml からダウンロードした 200708.xml からデータを取り出し、200708_polygon.shp などの Shapefile を作成するものです。

主スクリプト：cyberjapan.rb

require 'geotools'
require 'rexml/document'
require 'open-uri'
require 'combinadic'

polygon_w = Geo::FeatureList.open('200708_polygon.shp') 
line_string_w = Geo::FeatureList.open('200708_line_string.shp')
point_w = Geo::FeatureList.open('200708_point.shp')

def get_polygon(line_string, dups)
  polygonizer = Geo::Tools::Polygonizer.new

  # つなぎの線を消す操作
  dups.c(2) do |c|
    seg = Geo::import_wkt_geometry("LINESTRING (#{c[0]}, #{c[1]})")
    if line_string.intersection(seg).getLength == seg.getLength
      line_string = line_string.difference(seg)
    end
  end if dups.size > 2

  # 一点で接する穴を invalid ring lines と判断してしまうことへの対応 (Noding)
  line_string = line_string.union(line_string.getGeometryN(0).getStartPoint)

  polygonizer.add(line_string)
  polygons = polygonizer.getPolygons

  # 包含矩形の面積が最大のものが必要なポリゴンである。
  iter = polygons.iterator
  biggest_area = 0
  biggest_geometry = nil
  while(iter.hasNext)
    pg = iter.next
    area = pg.getEnvelope.getArea
    if(area > biggest_area)
      biggest_area = area
      biggest_geometry = pg
    end
  end
  biggest_geometry
end

open('200708.xml') do |p|
  doc = REXML::Document.new(p.read).root
  REXML::XPath.each(doc, '/GI/dataset/layer/Surface/GM_Surface/patch/GM_Polygon/boundary/exterior/generator/GM_Curve/segment/GM_LineString/controlPoint') do |cp|
    coordinates = []
    REXML::XPath.each(cp, 'column/direct/coordinate') do |coordinate|
      coordinates << coordinate.text
    end

    line_string_wkt = 'LINESTRING ('
    coordinates.each_index do |i|
      point_w.write({:the_geom => 
                      Geo::import_wkt_geometry("POINT (#{coordinates[i]})"),
                      :order => i})
      line_string_wkt += coordinates[i] + ','
    end
    line_string_wkt = line_string_wkt.chop + ')'
    line_string = Geo::import_wkt_geometry(line_string_wkt)
    line_string_w.write({:the_geom => line_string})

    dups = coordinates.clone
    coordinates.each do |e|
      dups.delete(e) if dups.index(e) == dups.rindex(e)
    end
    dups.uniq!

    polygon_w.write({:the_geom => get_polygon(line_string, dups)})
  end
end

polygon_w.close
line_string_w.close
point_w.close

本質的ではない XML まわりの処理や Shapefile まわりの処理がうまく隠蔽されていて、うれしいと思っています。

combinadic.rb 組合せ作成スクリプト

上記 cyberjapan.rb で require する combinadic は、d:id:hfu:20071109 で作成したものです。再掲します：

class Array
  # get combinations.
  # when called without a block, returns the number of combinations.
  # when called with a block, pass each combination to the block.
  def c(n_token, &block)
    if n_token > size
      raise "combination size (#{n_token}) larger than the array size (#{size})."
    end
    if n_token < 0
      raise "combination size (#{n_token}) should not be negative."
    end
    if block == nil
      return n_combination(size, n_token)
    else
      n_combination(size, n_token).times do |i|
        yield combination(n_token, i)
      end
    end
  end

  # get the number of combinations.
  def n_combination(n_slot, n_token)
    if(n_token == 0 || n_token == n_slot)
      1
    else
      n_combination(n_slot - 1, n_token - 1) + 
        n_combination(n_slot - 1, n_token)
    end
  end
  private :n_combination

  # return the combination at index i of combinadic
  def combination(n_token, i)
    c = []
    size.times do |sl|
      break if n_token == 0
      threshold = n_combination(size - sl - 1, n_token - 1)
      if i < threshold
        c << self[sl]
        n_token -= 1
      else
        i = i - threshold
      end
    end
    c
  end
  private :combination
end

作成したデータ

作成したデータは、以下のようになりました。

電子国土 Web システムで表示されるものが抜けなく取り出されていると思います。

苦労したところ

この電子国土プロファイルのポリゴンデータ、穴の処理がかなり特殊でした。

内周が外周に接する場合、一筆書き的に書き進む

しかも、ポリゴンの始点が内周から始まる場合があります。

接しない内周は、外周を書き終わった後、おもむろに書かれる

本来、内周を表現するものとして、exterior に対比した interior が用意されているのですが、その区別は使われず、exterior の中に、内周の情報が続けて書き込まれています。そのため、外周の終点と内周の始点を結ぶ線、あるいは内周の終点と別の内周の始点を結ぶ線が現れます。

また、線の向きは、内周・外周に関わらず、まちまちのようで、線の向きには特に意味はないようでした。

Polygonizer を使ってポリゴンを構築し直す

結局、点の順番を考慮していろいろ悩むことはやめ、JTS の Polygonizer に頼ることにしました。つまり、以下のようなアルゴリズムでポリゴンを得るようにしています。

• 電子国土プロファイルのデータで提供されている座標列から LineString を作る
• そのうち、外周の終点から内周の始点への線、内周の終点から内周の終点への線は、「2回訪れた点から2回訪れた点へ直接接続されている線」として特徴付けられると見なし、これらの線を消去
• Polygonizer で LineString からポリゴンを作成
• 目的とするポリゴンは、そのエンベロープの面積が最大になるもの*1として取得

*1:内周の中に作成されてしまうポリゴンは、エンベロープの面積が目的のポリゴンのエンベロープの面積よりも小さくなる。