紙地図からの地理情報の読みとり

紙に書かれた図形情報を読みとる装置として、 カーソルあるいはペンの位置を記録する「ディジタイザ」と 画像をラスターデータとして読みとる「スキャナ」があるが、 いずれにしても、直接得られる座標値は読みとり装置上のものである。 地理データの場合、緯度・経度などの地理座標がほしい。 そうでなくても、紙ではなく図の軸に対する座標値を得る必要がある。 このような機能は、地理情報システム(GIS)と呼ばれるパッケージ ソフトウェアの入力部分ではたいてい用意されている。 しかし、データの内容によっては、パッケージをブラックボックスで 使うだけではすまないことがある。たとえば、GISパッケージが 知らない投影法の地図を読みたい場合、地図投影変換の逆変換は 別に計算してやらなければならない。

(1) 地図(x, y)座標を得るための幾何変換

ディジタイザで取得できる位置の値は、 ディジタイザの板の上での「カーソル」(ペンあるいは虫めがねのようなもの) の横・縦の位置である。 紙の伸び縮みによる誤差を無視すれば、 これは、この地図が採用している地図投影法の定義式に 出てくる(x, y)座標と1次関数の関係にあると期待できる。 しかし、紙の置きかたをあまり精密に決められないのがふつうなので、 地図のx方向がディジタイザの横方向と厳密には一致しない。

そこで、まず、ディジタイザの座標(横方向を u、縦方向を vとする)から 地図上の直交直線座標(x, y)への変換は、 u, vそれぞれがxとyの両方の1次の項のたし合わせからなっているような 次の式(1)を採用するのがよいだろう。

      x = a u + b v + c
                                                   (1)
      y = d u + e v + f
このような変換をアフィン変換という。 これは、平行移動、回転、伸縮を特殊な場合として含んでいる。 リモートセンシングで得られた画像データの幾何変換にもよく使われる。

係数 a, b, c, d, e, fは、ここでは、対象とする図の範囲全体で定数と仮定する (図のゆがみが部分によって違うことは、とりあえず考えないことにする)。 係数を決めるには、 ディジタイザ座標と図上の座標の両方がわかっている 点がいくつか(少なくとも1直線上にない3点)必要である。 このような点を、コントロールポイントと呼ぶ。

この座標変換は、経験式の決定と適用の典型的な例と いうことができる。

  u, v ----->       [一般の点]       -----> x, y
                            ^
                            |
                    経験式の係数の決定
                            |
  u, v -----> [コントロールポイント] <----- x, y

プログラムxdigit (Fortranソースプログラムxdigit.f) では、係数を最小2乗法で求める。 つまり、コントロールポイントでの変換結果と与えた真の(x, y)座標値との 差の2乗の合計が最小になるように係数を決める。 ただし、コントロールポイントのx座標、y座標はそれぞれ (コントロールポイントだけに関する)平均値と標準偏差によって 規格化しておくことにした。したがって自動的に c = 0, f = 0となる。 (規格化はx軸とy軸の単位の違うグラフへの対応のために導入したものである。 地図では規格化しないほうがよいかもしれない。)

プログラムを実行すると、出力されるメッセージに対応して 標準入力から、読ませるデータファイル2つと、結果を書き出す ファイルの名前を入れるようにしてある。 データファイルの1つめはコントロールポイントの座標で、 各行に、u, v, dummy, x, yの5列を与える。 2つめは結果を求める点の座標で、 各行に、u, v, dummyの3列を与える。 この入力形式は、 グラフテック社のディジタイザで、 ASCII出力を設定した場合にコンピュータに伝わる情報形式に対応している。 このうちdummyは、 ディジタイザのカーソルのどのボタンを押したかによって決まる数字で、 変換には関与しないが、 あとの変換で点の性質(たとえば一連の折れ線の始め、終わりなど) の情報として使いたいことがあるので、 出力にも、x, y とともに第3列として書き出すようにしている。 なお、入力ファイルの 数値の間の区切りは、コンマまたは空白を使う。

ディジタイザは、スイッチによって送信するデータの形式を 変えられるものが多いが、ここではASCIIの数字として送信する モードに設定しておくことを想定している。 パソコンで通信ソフトウェアを起動し、速度やビット数を ディジタイザ側に合わせる。 まず、コントロールポイントの位置を読み、 続いて(紙を動かさずに)結果を求めたい点の位置を読む。 通信ソフトウェア終了後、交信記録ファイルを編集し、 コントロールポイントと結果を求めたい点を別ファイルにした うえで、コントロールポイントのファイルの各行に「真の(x, y)座標」 を追加してから、xdigitを実行する。

なお、xdigit.fの出力用のWRITE文のフォーマット(特に小数部の桁数)は、 「真の座標」の数値の範囲によって変更して使ってほしい。

(2) 緯度経度座標を得るための投影変換

データを地理情報として使うためには、位置座標を地球上の緯度(φ)・経度(λ)に 変換したほうがよいことが多い。 地図投影法のマニュアル (たとえば、Snyder and Voxland, 1989: An Album of Map Projections)には、 地図投影変換の式として、(λ, φ)を知って地図上の直交直線座標(x, y)を求める 式が書かれている。また、(x, y)から(λ, φ)への逆変換も書かれていることがある。 この式に使われている(x, y)を、前の(1)の節で述べた(x, y)として使えば、 図の縮尺の具体的数値についてはあまり気にしないで変換を行なうことができる。

なお、投影変換式は

      x = f( λ, φ )
      y = g( λ, φ )
のような、λとφが与えられたときにxとyがそれぞれ独立に求まる形(explicitな形) になっているのが望ましいが、
      F( x, y ) = f(λ, φ )
      G( x, y ) = g(λ, φ )
のようにxとyが連立された形(implicitな形)になっている場合もある。 この場合は、適当な値を仮定して逐次修正により、有効数字の範囲で正しい値が 得られるようなプログラムを作れることが多い。 逆変換の場合も, (x, y)と(λ, φ)の役割りが逆になるだけで同様である。

投影変換公式のプログラムとしては、

などがある。

ディジタイザを使ってベクター型データの緯度経度を求める手順は 次のようになる。

  1. コントロールポイント (緯度・経度がわかっていて、地図上でも同定できる地点)について、 地図投影の順変換を適用し「真の」(x, y)を求めておく。
  2. (1)節の方法により、一般の地点について(x, y)を求める。
  3. 得られた(x, y)に地図投影の逆変換を適用して(λ, φ)を求める。

実習の概略

実際の作業を行なうことができる場所は次のような形を用意した。

  1. OSとしてWindows 95が動いているPCで、 通信ソフトウェアを動かして、 このPCのシリアル(RS232)ポートに接続した ディジタイザ「グラフテックKD3800」からデータを取得し、 結果をファイルに入れる。
  2. Unix系のOSが動いているPCで、 上記のxdigit (の細部を変更したもの)と、 Evenden氏のprojを使い、緯度・経度を求める。
この手順では、1と2の間に WindowsからUnixへのファイル転送(ftp)が必要である。 これは不便だが、わたしがWindows上での数値処理にもUnixによる機器操作にも 慣れていないため、この形にしておいた。 なお、コントロールポイントの(x, y)座標値を追加するなどの ファイルの編集も必要だが、これは転送前後のどちら側で行なってもよい。

ラスター型の場合

ラスター型のデータを 緯度経度の等間隔にresamplingする場合 (リモートセンシングデータの幾何補正 として行なわれる) は、緯度経度の等間隔の位置の地図上の座標を求めて その位置のラスターデータの情報を得るため、主として使う変換は 地図投影の順変換である。


1999-06-20, 2000-07-05, 2001-04-20
増田 耕一

[地理情報学の教材の目次へ]