補完法の中でも、簡単なバイリニア補完について説明する。
バイリニア補間(Bilinear interpolation)とは、簡単に言うと
一般的なグリッドにおける2変数の補間関数による線形補間のことである。
つまり、グリッドデータを線形内挿する方法の1種である。
補間をするにあたっての前提として
・求めたいグリッドの範囲よりも広いグリッドの範囲でのデータがあること
例えば 30~35°N, 130~135°E の範囲でデータがあるとすれば、それより小さい範囲例えば 32~34°N, 132~134°Nでのグリッドを求めることができる
・元データが等間隔で並んでいること。等間隔でないと補間するグリッドデータの計算が大変(ここでは等間隔のみ扱う)
1) 求めたいグリッドの位置を囲む4角形を作る元データの4点を求める。ここで求めたいグリッドの位置に近い元データのグリッドの位置を探すのは計算が膨大となり無駄である。以下のようにして探すと容易である。
1.1) まずx方向(経度方向)について、求めたいグリッドのx座標(緯度)の値を元データのx座標(経度)の値で除し、それを元データの格子間隔で割る。
この結果、求めたいグリッドが元データのグリッドでx方向に何個目にあたるかがわかる。
最も、ぴったり一致するとは限らないので12.4個目というように小数になる。ということは、求めたいグリッドの位置を囲む点の左(西)側の点は求めた値の実部で右(東)側の点は求めた値の実部+1である。
1.2) 同様のことをy方向(緯度方向)に対しても行う。
1.3) 以上より求めたいグリッド位置を囲む4角形を作る元データの4点が求まる。
元データの経度をforlon、緯度をforlat、求めたいグリッドの経度をaftlon、緯度をaftlat
求めたいグリッドの数をx方向にP、y方向にQ
左側の個数をnumx_l、右側の個数をnumx_r、下側の個数をnumy_d、上側の個数をnumy_uとしてプログラムは、以下のように書かれる。(抜粋)
!バイリニア補間用のグリッド探索(X方向) DO i=1,P delx=aftlon(i)-forlon(1) !経度差 numx_l(i)=INT(delx/0.125d0)+1 !左側 numx_r(i)=numx_l(i)+1 !右側 a=(delx-forlon(numx_l(i)))/(forlon(numx_r(i))-forlon(numx_l(i))) IF (numx_l(i)>P.OR.numx_r(i)>P) STOP "Error at x" ENDDO !バイリニア補間用のグリッド探索(y方向) DO j=1,Q dely=aftlat(j)-forlat(1) !緯度差 numy_d(j)=INT(dely/0.1d0)+1 !下側 numy_u(j)=numy_d(j)+1 !上側 b=(dely-forlat(numy_d(i)))/(forlat(numy_u(i))-forlon(numy_d(i))) IF (numy_d(i)>Q.OR.numy_u>Q) STOP "Error at y" ENDDO
2)データの準備ができたのでバイリニア補間を行う。
具体的には、f(x,y)=f(m,n)×(1-a)×(1-b)+f(m+1,n)×a×(1-b)+f(m,n+1)×(1-a)×b+f(m+1,n+1)×a×b
ここで、f(x,y)は求めたいグリッドの値、f(m,n)は元データの左下のグリッド。
a,bは元データのグリッド間隔を求めたいグリッド位置から下ろし分割した時の比で0~1の値。あとは式に値を代入して計算していくだけ
求めたいグリッドでの値をaftval、元データのグリッドをforvalとしてプログラムは (抜粋)
DO j=1,Q
DO i=1,P
aftval(i,j)=forval(numx_l(i),numy_d(j))*(1-a)*(1-b)+ &
&forval(numx_r(i),numy_d(j))*a*(1-b)+ &
&forval(numx_l(i),numy_u(j))*(1-a)*b+ &
&forval(numx_r(i),numx_u(j))*a*b
ENDDO
ENDDO
ごく簡単な10*10のデータを用いてバイリニア補間をしてみた。
補間するグリッドは1.5~9.5まで0.5間隔で17*17のグリッドである。元データはこの項の最後にある。
PROGRAM bilinear
IMPLICIT NONE
INTEGER,PARAMETER :: M=10,N=10,P=17,Q=17
REAL(8),PARAMETER :: gridx=1.0d0,gridy=1.0d0 !元データグリッド間隔
INTEGER x,y,h,i,j,readstatus
INTEGER,DIMENSION(P) :: numx_l,numx_r
INTEGER,DIMENSION(Q) :: numy_d,numy_u
REAL(8) delx,dely,hoge
REAl(8) forlon(M),forlat(N) !元データx,y座標
REAL(8) aftlon(P),a(P),aftlat(Q),b(Q) !作成データx,y座標
REAL(8) forval(M,N),aftval(P,Q)
!元データ読み込み
OPEN(30,FILE="moto.txt",STATUS="OLD")
DO h=1,200
READ(30,FMT=*,IOSTAT=readstatus) x,y,hoge
IF (readstatus>0) STOP "Erorr"
IF (readstatus<0) EXIT
forval(x,y)=hoge
ENDDO
!元データのx,y座標割り当て
DO h=1,10
forlon(h)=DBLE(h)
forlat(h)=DBLE(h)
ENDDO
!バイリニア補間のグリッド探索(x方向)
DO i=1,P
aftlon(i)=0.5d0*DBLE(i)+1.0d0
delx=aftlon(i)-forlon(1) !経度差
numx_l(i)=INT(delx/gridx)+1 !左側
numx_r(i)=numx_l(i)+1 !右側
a(i)=(aftlon(i)-forlon(numx_l(i)))/gridx
IF (numx_l(i)>P.OR.numx_r(i)>P) STOP "Error at x"
ENDDO
!バイリニア補間のグリッド探索(y方向)
DO j=1,Q
aftlat(j)=0.5d0*DBLE(j)+1.0d0
dely=aftlat(j)-forlat(1) !緯度差
numy_d(j)=INT(dely/gridy)+1 !下側
numy_u(j)=numy_d(j)+1 !上側
b(j)=(aftlat(j)-forlon(numy_d(j)))/gridy
IF (numy_d(j)>Q.OR.numy_u(j)>Q) STOP "Error at y"
ENDDO
!求めたいグリッドの値の計算・出力
OPEN(40,FILE="aft.txt",STATUS="REPLACE")
DO j=1,Q
DO i=1,P
aftval(i,j)=forval(numx_l(i),numy_d(j))*(1-a(i))*(1-b(j))+ &
&forval(numx_r(i),numy_d(j))*a(i)*(1-b(j))+ &
&forval(numx_l(i),numy_u(j))*(1-a(i))*b(j)+ &
&forval(numx_r(i),numy_u(j))*a(i)*b(j)
WRITE(40,*) aftlon(i), aftlat(j), aftval(i,j)
ENDDO
ENDDO
STOP
END PROGRAM bilinear
