サンプルプログラム (Fortran。コメントが日本語コードのことがある。)
とりあえずのテストデータとして、『理科年表』1990年版から、 アフリカの気象観測地点の年降水量を afrrain [link] というファイル名で用意する。 1行めは人のためのメモ、2行めから空白で区切られた各列が 経度(度)、緯度(度)、年降水量(mm)である。
Voronoi: ボロノイ、20世紀初めのロシアの数学者
平面内の領域を、「もよりの」(距離がいちばん近い)データ点ごとに 区分する。つまり、データ点の垂直二等分線で区分することになるが、 全部の組み合わせの垂直二等分線が使われるわけではない。
立地・配分問題でよく使われる。たとえば岡部・鈴木(1992)参照。
また、水文学で降水量の集計法として使われるThiessen (ティーセン)法 (Thiessen, 1911) は、Voronoi多角形の面積の重みをつけて集計することにあたる。
Delaunay: ドロネー [この読みは杉原厚吉先生の本に従った]、 ソ連の数学者。
データ点でかこまれた領域を、各データ点を頂点とする三角形網で 区分する方法のうち、三角形の最小の内角を最大にする (三角形をなるべく細長くしない)という基準に従って最適なもの。
Voronoi多角形が相互に接する辺があるようなデータ点どうしを結ぶ ことに相当している。
Voronoi分割を計算して作図するプログラムの例文(Pascal)が、 永野ほか(1987)の本にのっていたので、 このアルゴリズムにしたがって、Fortranになおした 例文 voronoic.f を作った。 ただし、 データ点(Voronoi分割の理論では「母点」という)の座標を、 乱数ではなく、 ファイルから読みとるようにした。 また、作図はせずに、Voronoi多角形の座標値を、辺がつながる ように並べかえてテキストファイルに書き出すようにした。
このプログラムのアルゴリズムでは、 それぞれの母点について、他のすべての母点との間の垂直二等分線 を調べている。したがって、母点の数がふえていくと、少なくとも 点の数の2乗に比例して計算量がふえてしまう。 Voronoi図を構成するアルゴリズムについては、「計算幾何学」と呼ばれる 研究分野で多くの研究があり、
プログラム voronoic は、何もプロンプトメッセージは出さないが、 標準入力から1行めに入力ファイル名、2行めに出力ファイル名を入れてやると 動き出す。
テストデータ afrrain をこのプログラムに入れると、 経度・緯度をそのままx, yの座標値とみなして、 平面の距離を計算してVoronoi分割してしまう。 これは、正しくない。 地球上の大地域をVoronoi分割するには、 球面上の距離を計算し、球面多角形(各辺は大円の弧)で分割すべきである。 あるいは、完全ではないが、距離のゆがみの少ない地図投影をして、 その地図上で分割する。 そのような道具を作ろうと考えているが、まだできていないので、 正しくないことを承知のうえで例文としてある。 杉原(1998)の本には、 球面上の領域のVoronoi分割のプログラム例(Fortran)があるので、 それを参考に作るとよさそうである。
プログラムを修正し、 分布を計算する領域は母点の分布領域よりじゅうぶん大きく (たとえばプラスマイナス50倍の幅に)とっておき、 あとで、分布領域とその付近(たとえばプラスマイナス0.1倍)だけ を切断した多角形として出力できるようにした。 また、多角形を構成する線分を別の直線で打ち切る操作、 およびその中で線分と直線の交点を求める操作を サブルーチンとしてプログラムを分割した。
GKSを使って、 voronoic.f で計算されたVoronoi多角形を、 相互に区別できるように塗り分けるプログラムが voronoig.f である。色は順番に変えるようにしたが、 6色だけを使っているので、ときにはとなりどうしが同じ色に なる場合もある。 また、Delaunay三角形網も重ねてかいてみた。
voronoig.fの結果
また、同じVoronoi分割の各多角形を、観測地点(母点)の降水量の階級に よって色を変えて塗るプログラムが voronoip.f である。 (地理情報を見る人に伝えるためには、これでは未完成であり、 色の凡例や緯度・経度のめやすを入れなければならない。)
voronoip.fの結果
不規則に分布する観測地点でのデータをもとに等値線をひく方法の ひとつは、対象領域を観測地点を頂点とする三角形で分割することである。 対象となる量が、三角形の辺上で、距離の1次関数で変化すると 仮定すれば、 格子点データの1つのマス目の場合 (塩野, 1988a)と同様に、その三角形の内部の等値線をひくことができる。 三角形では、(観測値が等値線レベルに一致する場合、数値をずらして ごまかすことにすれば)等値線と辺との交点数は2か0のいずれかに 限られるので、四角形の場合より簡単になる。 サブルーチン drawc3.f [link]は、格子点データの場合の drawc4.f を なおして作ったものである。
このための三角形分割は、必ずしもDelaunay分割である必要はない。 プログラム conttri2.f [link] では、Akima のプログラム BIVAR (別記事「格子点に内挿して等値線をひく」参照) の三角形分割の部分 (ファイル idtang.f )を利用した。 この三角形分割は最小の内角を最大にするものなので、Delaunay分割と 一致すると思われるが確かめていない。 この例文の実行では、データファイル名 afrrain を標準入力から 与えてやればよい。
一方、プログラム conttri.f [残っていない] では、 voronoic.f で得られたVoronoi分割に対応するDelaunay分割 によって、三角形を決めることにした。 そのため、例文の実行では、データファイル名 afrrain と afrrain.vorの 両方を読ませる必要がある。 voronoic.fの出力には、どの母点(観測地点)とどの母点を結ぶ辺があるかは 書かれているが、どの辺の組み合わせの三角形があるかを直接には示していない。 そこで、森(1991)のサブルーチン mkntw と lstotr をこの順に使って、 辺の羅列から三角形のリストに変換することにした。 (この部分はもっと簡単な操作になおせると思う。)
三角形網が得られたところで、作図にはいり、 それぞれの三角形について drawc3 を呼んで、等値線をひいている。 例文では、降水量100 mm ごとに線をひき(これだけの地点数にしてはやや 細かすぎるようだが)、500 mm、1000 mm、2000 mmの倍数のところでは それぞれ色を変えるようにした。 各辺上では直線で空間内挿しているので、なめらかではない。
等値線をかき終わったあと、観測地点の位置をマーカーで、 また、補助情報として地図(海岸線と国境線)を線で重ねた。 海岸線などをかく部分は drwmap.f [link] で、 実質的には格子点気候データの等値線図 の例文のプログラムの地図をかいている部分と同じである。 地図としてはこのほか、経緯度のワク(あるいは格子)も重ねるべきだろうが、 まだ入れていない。
conttri2.fの結果
