等値線(地形図の等高線や、天気図の等圧線など)は、 2次元の座標 x, y に対して与えられる量 Z が、 x, yに対して連続に分布しているときに定義されます。 ただし、有限の(無限小でない)長さにわたって Z が同じ値をとる場合は 定義できないことがあります。
等値線とは、Z が ある値 Z1をとる位置を結んだ線です。 線の(たとえば)右側で Z が Z1よりも大きいとすれば、 左側では Z は Z1よりも小さい値をとります。 したがって、1つの等値線は、(x, y)平面を、ふたつの領域に分離する 閉じた曲線であるはずです。 ただし、対象となる領域全体が有限のときには、等値線は領域の端で 止まることがあります。この場合、止まるところは偶数個あるはずで、 仮に端に沿って結んで考えれば、閉じた領域を定義することができます。
ふつう、いくつかの Z の値についての等値線を同じ紙面に重ねてかきます。 このようにすると、等値線の間隔から、斜面の傾斜(一般には Z の空間座標による 微分)を視覚的に読みとることができます。
計算機でデータを処理する場合、 データが等間隔の格子で与えられているというところから出発することが多いです。 数値実験(シミュレーション)では実際に格子上に値があることが多いのですが、 現実の観測(測量)では、必ずしもそうではありません。しかし、じゅうぶん密に 観測値が得られて、しかも、値の空間分布が連続であると仮定することができれば、 統計的方法で、格子点上の値を推定することができます。ここでは、その方法には ふれないことにします。
Zの値は格子点でだけ与えられているのですが、等値線は格子点の中間を 通ることができます。どこを通るかを推定しなければなりません。 いろいろな考え方があるのですが、 いちばん簡単なのは(必ずしもいちばんきれいな結果にはなりませんが)、
このように問題を決めても、計算する手順は何通りか考えられます。 ここでは、塩野ほか(1988a)にならって例文を作りました。 例文プログラムは次のものです。 (/pub/sfc/earthobs/example/ncargks/{f,c}/cont/ の下に置きます。 (waku2, wc2ndc, senbun は、すでに[GKSを使ったプログラム例]使ったものです。)
contmain.f [link] ... メインプログラム
drawc4.f [link] ... 主要部分のサブルーチン
rdelev.f [link] ... データを読むサブルーチン
waku2.f [link] ... 目盛りつきのワクをかくサブルーチン
wc2ndc.f [link] ... 座標変換(waku2の補助)サブルーチン
senbun.f [link] ... 線分をかく(GKSのGPLの応用)サブルーチン
leglin.f [link] ... 線の凡例をかくサブルーチン
なお、サンプルデータとして、
elevtibet.dat [link] ... チベット高原を中心とする地域の標高
を置きます。これの由来はこのページの下のほうで説明します。
実行のときは、このデータは、datapath.fiの dapath という文字列で
指定されたファイルパスの下の elev というディレクトリの置かれることを想定しています。
実行結果は次の図のようになります。
ここでの説明は、Fortranの例文にしたがって進めますが、 Cの例文も対応するように作ってあります。
関係するファイルを自分の適当なディレクトリにコピーし、
% ncargf77 contmain.f % ./a.out % ctrans gmetaとすると作図されるはずです。
計算機で何かをするプログラムを作る場合に考えることは、 同じ形式の仕事を繰り返すように分割することです。 ここでは、
Zの値の各レベル Z1, Z2, ... Zk (たとえば 0 m, 1000 m, ... ) x座標の格子 x1, x2, ... xi (たとえば 東経70度, 71度, ...) y座標の格子 y1, y2, ... yj (たとえば 北緯20度, 21度, ...)という3つのループ(いいかえれば、数え上げ)があります。 どちらを先にしてもよいのですが、例文のプログラムでは、 まず、「1つのレベルの等値線をかく」サブルーチン「drawc4」を作り、 それを、各レベルごとに CALL するようにしています。 drawc4の中に、xの格子番号、yの格子番号についての2重ループがあります。
プログラムの構造を理解するのに、「トップダウン」つまり、大きな構造から だんだん細かく見て行くのと、「ボトムアップ」つまり、小さな部品から 組み立てて行くのとあります。ここでは、等値線をかく仕事に本質的な 部分に限って、ボトムアップで見て行くことにします。
まず、この3つのループのいちばん内側でやるべきことを考えましょう。
説明用の図
(この図は、
onebox.f [link]
を自分のディレクトリにコピーし、
% ncargf77 onebox.f % ./a.out % ctrans gmetaで作ることができます。)
今、x方向はi番めと(i+1)番め、y方向はj番め(yj)と(j+1)番めの格子線に 注目しているとします。ここには、(1)(2)(3)(4)の4つの格子点があります。 (C言語の場合、配列添字は0から始まるので、ここでの(1)..(4)が 配列添字 0..3に対応しています。)
ここで、話を簡単にするために、どの格子点の Z の値も、等値線を引きたい 値(drawc4.f の 変数 cval) に等しくないとします。
1-2, 2-3, 3-4, 4-1 の4つの辺上を等値線が通るかどうか考えます。 すると、各辺上で、片側の格子点の Z の値が cval よりも大きく、反対側の Z の値が cval よりも小さい場合に限って、等値線が1回通ることになります。 (2回以上通ったり、この条件を満たさないのに通ったりすることがないのは、 辺上では Z の値は座標の単調な関数(この場合1次関数)と仮定したためです。)
4つの頂点の値が cval よりも大きい[+]か小さいか[-]の組み合わせを考えると、 等値線が通る辺の数は、0, 2, 4 のいずれかになります。
+--C--+ +--C--+
+/ + + \+
D B D B
+ / + +\ +
+--A--+ +--A--+
どちらを結んでも等値線の理屈には合うのですが、全体としてなめらかな面
になったほうが望ましいことが多いので、プログラム drawc4 では、線分の
長さの合計を比較して、短いほうをとることにしています。
なお、「一方が cval より大きく、他方が cval より小さい」かどうかを判定 する方法として、それぞれの Z の値から cval を引いたものどうしの積をとり その符号で判断しています。(数値のかけ算をする必要はなく、符号だけ演算す ればいいのですが、Fortranでは符号の演算があまり便利でないのでこのように書 きました。)プログラム中「IF(fc(1)*fc(2) .LE. 0) THEN ...」のところです。
現実には、格子点の Z の値が cval に一致することがありえます。 ここでは、そのような場合にはZの値に小さな値 deltaf をたすことによって 逃げています。
サブルーチン drawc4 は、 それを CALL するプログラムから、次の情報を引き渡されます。
nxdim (整数) 配列zの第1次元の寸法、nxdim >= nx
nydim (整数) 配列zの第2次元の寸法、nydim >= ny
nx (整数) x 方向の格子線の数
ny (整数) y 方向の格子線の数
x (実数1次元配列、寸法 nx以上 ) x 方向の格子線のx座標値
y (実数1次元配列、寸法 ny以上 ) y 方向の格子線のy座標値
z (実数2次元配列、寸法 nxdim, nydim ) 格子点での Z の値
X=x(ix), Y=y(jy) で Z=z(ix,jy)
cval (実数) 等値線をひく値
itype (整数) 線の種類 GKSのGSLNへ
icol (整数) 線の色番号 GKSのGSPLCIへ
deltaf (実数) Zがcvalに一致したとき便宜上ずらす値
vnull (実数) Zの値がvnull以上ならばデータ欠損とみ
なす
サブルーチンの大部分は、
DO jy=1, ny-1 DO ix=1, nx-1の2重DOループになっていて、 その中で、この記事の3節で書いたように、4点に注目します。 プログラムを見やすくするため、この4点の x, y, z の値を xp, yp, fcという 1次元、寸法4の配列につめかえて作業しています。 実際に fc に入れるのは、zの値から cval をひいたものです。 以下、ループ内の作業内容は、3節を見てください。
実際に線を引くところでは、senbun というサブルーチンを呼んでいます。 これは、
CALL senbun (x1, y1, x2, y2)
の形で、点 (x1, y1) と 点 (x2, y2) を結ぶ直線の線分をひくものです。
内容は、senbun.f にあり、GKSのGPLを点の数2で呼んでいるだけです。
まず宣言部がきますが、実行部で使う変数を宣言しているだけなので 説明は省略します。 ただし、PARAMETER文は定数の宣言、 DATA文は変数の値(初期値)を与える文です。 機能は似ていますが、配列の寸法は、定数で宣言しなければなりません。 (ただし、サブルーチンの引き数となった配列の寸法の、サブルーチン側での 宣言は変数でかまいません。drawc4のx, y, zの場合がそうです。)
実行部の最初で、rdelevというサブルーチンを呼んで、データを読み込みます。 rdelevの主要部分はFortranのREAD文です。1つめの引き数でファイル名を指定 しています。rdelevは、ファイルを読んで、nx, ny, x, y, z の値をメイン プログラムに返してきます。 (rdelevの中のREAD文とデータファイルの構造は対応していなければ なりませんが、ここでは詳しい説明を省略します。)
次に、GKSを開始します。opngksはNCAR Graphicsのルーチンで、 GKS本来の開始に必要な3つのルーチンを中で呼んでいます。 続いて、座標変換(変換番号 1)を定義します。 この例では、次のような window と viewport を同一視します。 (次の図は、<pre> の中ではASCII文字(空白を含む)の幅が一定である ことを前提としてかいてあります。)
WC (世界座標) NDC (正規化装置座標)
42.0 +--------+ 0.7 +--------+
| | | |
| window | |viewport|
| | | |
18.0 +--------+ 0.2 +--------+
68.0 112.0 0.2 0.9
次に、目盛りつきのワクをかきます。
それからいよいよ、等値線をかくのですが、 drawc4を1回呼ぶと、1つの cval についてかくので、 変数 k の DOループで、cvalを変えながら繰り返しdrawc4を呼んでいます。
そのあと、図の上の表題と、左側、下側の軸の見出しを、 GKSの文字列作図ルーチン GTX を呼んでかいています。 ただし、文字の位置は NDC で指定したほうがわかりやすいので、 GTXを呼ぶ前に、CALL GSELNT(0) で、変換番号 0、つまりWCをNDCと 等しくする変換を指定しています。 なお、縦軸の目盛りを横倒しにするため GSCHUP で「文字の上向きの方向」 を指定しています。
最後に clsgks (NCAR Graphicsのルーチン) を利用して GKSを終了します。
ここで必要なサブルーチンは、GKSやNCAR Graphicsのものを除いて、 Fortranプログラムで用意してあり、メインプログラムの最後にかかれた INCLUDE文で取り込まれていっしょにコンパイルされるようにしてあります。 INCLUDE文はFortranの標準にはありませんが、かなり多くのコンパイラで 利用可能です。
なお、ここで例として使った標高データは、 アメリカのスクリプス海洋研究所で作られ、ランド(Rand)研究所で改訂された 世界の緯度・経度1度方眼の標高・水深データの一部です。
data : 標高データ [link]MASTER : データの説明 (英語、NCARによる) [link]format : データの記録形式の説明 (英語、NCARによる) [link]read.f : データを読むプログラム例 (Fortran、NCARによる) [link]これを使って全地球の標高の等高線図を作るプログラム例を
contglob.f [link] として置きます。
別ファイルにしたサブルーチンとしては、すでにのべたもののほかに、
Randの形式の標高データを読むための
relevg.f [link] を使います。
2000年現在、世界規模では、緯度・経度30秒つまり(1/120)度、約1 kmの標高データがあり、 自由に使えます。 ただし、実質的にはもっとあらい分解能しかないところもあります。 【2020年現在、緯度・経度 3秒、約 90 m のデータもあります。 】