ドイツの白地図操作ができるようになるまで

hfu2007-11-02

ベクトルデータを使って、ドイツについての主題図が作成できるようになるまでのステップをまとめます。geotools.rb を使ったスクリプトを含む、ちょっと長いエントリになっています。

ソースデータの入手

WHO の SALB プロジェクトからドイツの SALB を入手しました。
http://www.who.int/whosis/database/gis/salb/salb_home.htm
通常、海外のベクトル地図データが必要な場合にはまず地球地図 ( http://www.iscgm.org, http://www.globalmap.org/ )にあたるのですが、地球地図プロジェクトにおいては現在のところヨーロッパ地域のデータは公開されていません。地球全体のベクトル地図ですでに完成しているものとしては、他に VMAP0 がありますが、これは古い上に、フォーマットが VPF というサポートの薄い MIL 規格フォーマットなので、できれば使いたくありませんでした。
そこで、国の中の第二レベル（日本で言えば、第一レベルが都道府県、第二レベルが市区町村になるでしょうか）までの行政区画のデータを整備しようとしている、SALB プロジェクトからデータを入手しました。行政区画のデータだけが必要な場合には、SALB プロジェクトは助けになるでしょう。*1
SALB プロジェクトと、SALB のデータについては、あとでもう少し書くかもしれません。

投影法の選定

SALB プロジェクトで得られるデータは、投影をしていない経緯度座標系のデータです。
ドイツはかなり高い緯度にある国ですので、投影をしないと見栄えがよくありません。そこで、地図に投影を加えることにしました。
手近にある小さい縮尺の紙の地図を何枚か見たところ、ドイツ全体を一覧するための地図には、標準緯線2本のランベルト等角円錐図法がよく使われていることが分かりました。
そこで、標準緯線2本のランベルト等角円錐図法を使うことにしました。標準緯線は、紙の地図に記載されていた、48度40分, 50度40分の緯線を使うことにしました。

投影法の再現

この投影法を、GeoTools で使えるようにします。GeoTools には、org.geotools.referencing.wkt.Parser というクラスがあり、これが WKT (Well-Known Text) 記述の投影法定義をパーズすることができます。このため、欲しい投影法の WKT 記述を作ることが当面の目標になりました。

標準緯線2本のランベルト等角円錐図法を Well-Known Text 形式で記述する方法

geotools.rb に定義しておいた Geo::find_epsg_crs を使い、標準緯線2本のランベルト等角円錐図法の既存の定義を調べることができます。いくつかある既存の定義の中で、最も癖の少ないものを探したところ、EPSG:3033 の WGS 84 / Australian Antarctic Lambert が最も癖が少ないものであるように思われました。この定義は、以下のようになります：

PROJCS["WGS 84 / Australian Antarctic Lambert", 
  GEOGCS["WGS 84", 
    DATUM["WGS_1984", 
      SPHEROID["WGS 84", 6378137.0, 298.257223563, AUTHORITY["EPSG","7030"]], 
      AUTHORITY["EPSG","6326"]], 
    PRIMEM["Greenwich", 0.0, AUTHORITY["EPSG","8901"]], 
    UNIT["degree", 0.017453292519943295], 
    AXIS["Longitude", EAST], 
    AXIS["Latitude", NORTH], 
    AUTHORITY["EPSG","4326"]], 
  PROJECTION["Lambert_Conformal_Conic_2SP"], 
  PARAMETER["central_meridian", 70.0], 
  PARAMETER["latitude_of_origin", -50.0], 
  PARAMETER["standard_parallel_1", -68.5], 
  PARAMETER["false_easting", 6000000.0], 
  PARAMETER["false_northing", 6000000.0], 
  PARAMETER["standard_parallel_2", -74.5], 
  UNIT["m", 1.0], 
  AXIS["x", EAST], 
  AXIS["y", NORTH], 
  AUTHORITY["EPSG","3033"]]

この中で手を加えるべき所は、PARAMETER の部分です。地図の形状に直接的な影響を持つのは、standard_parallel_1, standard_parallel_2, central_meridian の3つです。2つの standard_parallel には、紙の地図に記載されていた数値を、central_meridian は地図の中央付近できりの良い数字である 10 を入れることにしました。
その他の、latitude_of_origin, false_easting, false_northing は、どちらかというと数値の調整という程度のものですので、適当な数値を入れます。今回はすべて 0 を入れることにしました。
こうしてできあがった WKT 定義は、後述の Ruby スクリプトの中に記述します。

Geo::import_wkt_crs の定義

org.geotools.referencing.wkt.Parser を使う部分は、geotools.rb の中に Geo::import_wkt_crs として整備しました。
ソースは現在のところ、http://svgmapdata.sakura.ne.jp/geotools/classes/Geo.src/M000006.htmlで閲覧できます*2。

スクリプトを作成

これでソースデータも投影法も分かったので、ドイツの白地図が手に入れられるようになったことになります。あとは前回のエントリの6つの州についてフラグを立てるための属性 dgk5 を書き込む部分を加え、以下のようなスクリプトを作成することができました：

require 'geotools'

states = %w{Bremen Niedersachsen Nordrhein-Westfalen Rheinland-Pfalz Schleswig-Holstein Saarland}
wgs84 = Geo::import_epsg_crs(4326)
lambert = Geo::import_wkt_crs(<<-EOS
PROJCS["a Lambert conformal conic projection for Germany", 
  GEOGCS["WGS 84", 
    DATUM["WGS_1984", 
      SPHEROID["WGS 84", 6378137.0, 298.257223563, AUTHORITY["EPSG","7030"]], 
      AUTHORITY["EPSG","6326"]], 
    PRIMEM["Greenwich", 0.0, AUTHORITY["EPSG","8901"]], 
    UNIT["degree", 0.017453292519943295], 
    AXIS["Longitude", EAST], 
    AXIS["Latitude", NORTH], 
    AUTHORITY["EPSG","4326"]], 
  PROJECTION["Lambert_Conformal_Conic_2SP"], 
  PARAMETER["central_meridian", 10], 
  PARAMETER["latitude_of_origin", 50.6666666666667], 
  PARAMETER["standard_parallel_1", 52.6666666666667], 
  PARAMETER["false_easting", 0], 
  PARAMETER["false_northing", 0], 
  PARAMETER["standard_parallel_2", 48.6666666666667], 
  UNIT["m", 1.0], 
  AXIS["x", EAST], 
  AXIS["y", NORTH]]
EOS
)
t = Geo::Transform.new(wgs84, lambert)

Geo::Writer::open('edit.shp') do |w|
  Geo::Reader::foreach('deu_SALB_jan05_may05.shp') do |geom, attrs|
    geom = t.transform(geom)
    w.write(geom, {'dgk5' => states.member?(attrs['ADM1_NAME']).to_s,
      'name' => attrs['ADM1_NAME']})
  end
end

QGIS で画像化

前出のスクリプトによって作成した edit.shp を QGIS で読み込み、属性 dgk5 によって色分けし、その表示を画像として書き出すことにより、以下の画像を手に入れることができました：

このエントリの品質について

このエントリには誤りを含むかもしれません。

geotools.rb について

このエントリで使った geotools.rb は、http://svgmapdata.sakura.ne.jp/geotools から入手することができます。今回の作業のために、 Geo::import_wkt_crs を追加しました。

*1:但し、日本を始め、多くの国が現時点では SALB プロジェクトからデータを提供していません。

*2:このリンクは、geotools.rb の今後の改良により変化してしまう可能性があります。その場合は、http://svgmapdata.sakura.ne.jp/geotools/ から必要な箇所をたどっていただく必要があります。