不規則に分布するデータをもとに等値線をひく方法として、 三角形分割して各三角形要素についてひく方法 もあるが、いったん格子点データを作り、あとは 格子点データから等値線をひく方法 を適用することが多い。
データ点から格子点への内挿の方法はいろいろ考えられる。 内挿についての概論は 別の授業用の記事にまとめておいた。
実際のデータを操作するプログラム例のうちいくつかを [link]と書いたリンクさきに用意する[この段落 2020-08-16 変更]。 ただし、Fortran言語によるものだけである。
bivar.f [link]は、秋間浩氏 (アメリカ合衆国Institute of Telecommunication Science)による 2次元の空間内挿サブルーチンである。 著作権上はパブリックドメインとなっており、 National Center for Atmospheric Researchその他から公開されている。 (NCAR Graphicsの例文でも使われているが、NCAR Graphicsの一部ではない)。 アルゴリズムの説明は Akima (1978a, 1978b)、部分的修正の説明が Akima (1984)の論文にある。概略を増田(1987)の解説中に紹介した。
これは、(数値計算の誤差を別として)データ点では観測値に一致する内挿である。 領域をデータ点を節点とした三角形要素に区分し、各三角形内で 5次式を使い、辺上で1階微分が連続になるようにしたものである。 三角形分割は、最小の内角が大きいほうを選ぶという基準で行なっており、 ほぼDelaunay分割と同じものだが、全く同じかどうかは確認していない。
秋間氏のプログラム(NCARで少し編集されているが) bivar_orig.f [link] は、1970年代(Fortran 77標準より前)に作られたものであるため、 今からみれば古いスタイルである。 bivar.fは、見通しをよくするため、 Fortran 77のIF..END IFや、Fortran 90のDO..END DOなどの構文を使って 書きなおしを試みたものである。GO TO文がいくらか残ったが、これは ループを脱出するためのもので、おそらくFortran 90のEXIT文に変える ことによってなくすことができるだろう。 計算結果は、わたしが試みた実行例に関する限り同じになるが、 どんな入力に対しても同じになるかどうかは確認していない。 また、コメント行の形でプログラムの初めについている説明は 古いままである。
bivmain.f [link] は、bivarのうちのサブルーチン idsfftを使って、 地点ごとのデータxd, yd, zdを読み、 格子点値xg, yg, zgを別のファイルに書き出すものである。 作図を含んでいないので、f77でコンパイル・リンクできる。 このプログラムには、コメントが日本語EUC文字コードで書かれている。 このプログラムでは、bivar.f と、つぎの別ファイルに置かれているサブルーチンを利用する。
実行するときは、標準入力から
サンプルデータとしては、 三角形分割の例文と共通に、 アフリカの年降水量のデータ afrrain [link] (『理科年表 1990』による)を使うようにしてある。 このファイルのx, y座標は実は経度、緯度なので、 bivarで使われる距離は地球上の距離としては正しくないという 問題点を残している。
fit2.f [link]は、塩野ほか(1988b)のN88 Basicによる例文 「FIT(2).FIL」をFortranで書きなおしたものである。 この本では「制約条件つき最適化」と呼んでいるが、わたしの分類 (別記事)では変分法に含まれる方法である。 なめらかな場を求めるので、データ点でも必ずしも観測値に一致しない。
目的関数は、変量のx, y座標による1階微分のノルム、2階微分のノルム、 観測値と場の値との差のノルムの線形結合である。 ただし、例文のプログラム中で与えた係数の値では、1階微分の項の重み (pm1)が0に なっているので、 「データに近いという制約のもとで、薄板の変形エネルギーを最小化する」 という定式化だと考えてよい。(曲率を最小化するのに似ているが、少し違う。)
場の値を求めるために連立1次方程式をとく必要がある。 例文fit2では、くり返し過緩和法(successive over-relaxation, SOR)という方法を使っている。 MS-DOS上でも実用的に動くようにするため、係数行列の特殊な形を利用して メモリーを節約できる方法を採用したのである。 ワークステーション教室の機械・OSを前提とすれば、 塩野ほか(1988b)の例文「FIT(1).FIL」のように、 連立一次方程式を消去法(具体的にはコレスキー法)で 直接といたほうがよいかもしれない。
入出力データの形式は、bivmainの場合と同じにした。 (出力ファイル名の例は、afrrain.fit)。 作図は含んでいない。
contgrd.f [link]は、 bivmainまたはfit2で作られた格子点データを 等値線表示するプログラムである。 ncargf77でコンパイル・リンクする。 基本的に格子点データの等値線の説明で 使ったcontmain.fと同じであり、 同じdrawc4.fを呼び出している。
また、三角形分割の例文conttri.f の場合と同様に、観測地点の位置に記号をかき、 さらに、(大陸規模の気候データの例を想定して) 地図の海岸線・国境線を重ねた。
実際にデータの表示に使う場合は、等値線の凡例と、座標軸の方向や尺度を 示す線(この例では緯線・経線)を入れるべきであるが、今のところ省略している。 また、いずれの内挿方法も、データ点に囲まれた領域以外への外挿には あまり適さないので、内挿のプログラムに手を加えて、外挿部分を欠損を示す値 (contgrdのvnullかそれ以上の値)に 置きかえて書き出しておくのがよいだろう。
