時空間方向のラグ相関
対象海域における「ある変数」の空間代表性、卓越する時空間変動を理解する。
「ある変数」はなんでもよいが、ここでは海洋のデータを扱っているので、
衛星観測値の中でも密度（水温と塩分の関数）の逆数の積分、つまり水柱の情報を
反映する海面高度を用いることにする。
（実践）
まずは図を示す。以下の図は北緯27度、東経127度の海面高度偏差と周囲の海域における海面高度偏差との相関係数の分布図である。混乱すると思うので、今回は南北方向のみ考える。
X軸は空間方向のラグ、Y軸は時間方向のラグを示している。
たとえば
X=0、Y=0（時間のラグなし）における相関係数は「１」
X=100km、Y=0（時間のラグなし）における相関係数は「0.7」程度
X=100km、Y=２（2か月のラグ）における相関係数は「0.１」程度
時間方向のラグで考えると同じ場所でも1か月ズレると相関係数は落ちる。
ここでは、黒潮＆黒潮続流域で切り出した海面高度データを使用する（Teamsで配布済み）
端的に言えば、自分が指定したグリッドの海面高度偏差時系列と、南北方向にズレた海域における海面高度偏差時系列との相関係数（同時相関・ラグ相関の双方）を計算する。
以下はFortranのサンプルコード。例によって穴埋め形式。
あくまで説明用なので、Teamsで配布したコードを使うこと。
c234567
! This code calculates autocorrelation in time & space domain
character*120 header(364)
character*120 ifl
real sla3d(80,80,364),lon(80,80),lat(80,80),dummy(80,80)
real COR(23,13),ilon,ilat,DIST(23,13),TIME(23,13),dif1,dif2
integer year(364),month(364),xvar(364),smo,emo
real ave1,ave2,var1,var2,std1,std2,covar,dnum
!------------ parameter! --------------
g=9.8; ! m/s^2
pi=3.141592;
omega=7.29*1e-5;
a=6378136; ! m
e2=0.006694470
!----------------------------------------------
! read SLA dataset; build SLA 3d(lon,lat,time)
! ファイル名を作成して、海面高度データを読み込む
!----------------------------------------------
inum=1
smo=1;emo=12;
do iyr=1993,2023
if(iyr.eq.1993)then
smo=2
endif
if(iyr.eq.2023)then
emo=5
endif
do imo=smo,emo
if(imo.ge.10)then
write(ifl,'("./KuroshioBIN/",i4,i2,"-80x80.bin")') iyr,imo
else
write(ifl,'("./KuroshioBIN/",i4,"0",i1,"-80x80.bin")') iyr,imo
endif
write(*,*) ifl
200 format(a120) !--make file's name
open(10,file=ifl,form='unformatted',
& access="direct",status="old",recl=80*80*4)
read(10,rec=1) dummy
do j=1,80
do i=1,80
sla3d(i,j,inum)=dummy(i,j)
enddo
enddo
year(inum)=iyr ! we will use it later!
month(inum)=imo ! we will use it later!
xvar(inum)=inum
write(*,*) inum,year(inum),month(inum),xvar(inum)
inum=inum+1
enddo
enddo
inum=inum-1
write(*,*) "inum=",inum
!----LONLAT info
open(10,file="./KuroshioBIN/lon-80x80.bin",form='unformatted',
& access="direct",status="old",recl=80*80*4)
read(10,rec=1) lon
close(10)
open(10,file="./KuroshioBIN/lat-80x80.bin",form='unformatted',
& access="direct",status="old",recl=80*80*4)
read(10,rec=1) lat
close(10)
!----------------------------------------------
! specify x,y coordinate based on my input
! 調べたい場所の経度、緯度を入力する。
! そのあと、入力した緯度経度に最も近いグリッドの座標をix,iyに格納する。
!----------------------------------------------
write(*,*) "Enter Lon/Lat info"
write(*,*) "Enter longitude"
read(*,*) ilon
write(*,*) "Enter lat"
read(*,*) ilat
write(*,*) "You entered",ilon,ilat
dif1=100000; dif2=100000
do j=1,80
do i=1,80
write(*,*) i,abs(ilon-lon(i,j)),ilon,lon(i,j)
if(abs(ilon-lon(i,j)).lt.dif1)then
dif1=abs(ilon-lon(i,j))
ix=i
endif
if(abs(ilat-lat(i,j)).lt.dif2)then
dif2=abs(ilat-lat(i,j))
iy=j
endif
enddo
enddo
write(*,*) "ix and iy are ",ix,iy
!----------------------------------------------
! cal autocorrelation in time-space domain
!----------------------------------------------
! consider in latitudinal direction only, this time
! COR=nan(23,13)-->(space and time)
!----------------------------------------------
do it=0,12 ! temporal 時間方向のラグを決める
istart=1; iend=364;
do j=-11,11 ! spatial 空間方向のラグを決める
!~~~~ cal distance ~~~~! ここは前回行った距離計算！
mlat=(lat(ix,iy)+lat(ix,iy+j))*0.5
dlat=
coli=
AA=
BB=
CC=
DD=
disty=(pi/180)*(AA/DD)*dlat !! distance (m) between two grids
if(lat(ix,iy).gt.lat(ix,iy+j))then !!!!! caution !!!!!
disty=-1*disty
endif
! write(*,*) mlat,coli,AA,BB,CC,DD,"disty=",disty
!------ Cal. Autocorrelation
! 1. get "MEAN" values
ave1=0
ave2=0
iend2=iend-it
do k=istart,iend2 !--> ex) 1-363 2-364
ave1=ave1+sla3d(ix,iy,k) 自分が入力した緯度経度の海面高度偏差を使う
ave2=ave2+ 入力した緯度経度から空間・時間方向にラグを持つ海面高度偏差を使う、変数itとｊも使う
enddo
dnum=iend2-istart+1
ave1=ave1/dnum
ave2=ave2/dnum
write(*,*) "average",ave1,ave2,iend2
! 2. get variance and std
var1=0; var2=0;covar=0
do k=istart,iend2
var1=var1+
var2=var2+
covar=covar+
enddo
var1=var1/dnum
var2=var2/dnum
covar=covar/dnum
std1=sqrt(var1)
std2=sqrt(var2)
! 3. cal correlation
write(*,*) "var1,var2,covar",var1,var2,covar
write(*,*) "std1,std2",std1,std2
COR(12+j,1+it)=covar/std1/std2 相関係数の行列
DIST(12+j,1+it)=disty; 空間方向のラグの行列
TIME(12+j,1+it)=it; 時間方向のラグの行列
write(*,202) j,it,COR(12+j,1+it)
202 format("Area and Correlation-->",i4,i4,2f12.2)
enddo
enddo
open(30,file="SLAautocor.txt")
do j=1,13
do i=1,23
write(*,300) DIST(i,j)/1000,TIME(i,j),COR(i,j) 距離の単位は「m」なので「km」に換算
write(30,300) DIST(i,j)/1000,TIME(i,j),COR(i,j)
enddo
enddo
300 format(3f14.4)
stop
end
※赤字の箇所は自分が指定した緯度・経度に書き換えること
GMTスクリプト
gmt begin TimeSpaceCor
gmt set FONT 18p,Helvetica-Bold フォントとフォントサイズ
gmt surface SLAautocor.txt -R-305/305/0/12 -I5/0.1 -GSLAautocor.nc テキストデータからNetCDFを作成
gmt makecpt -Crainbow -T-0.5/1/0.1 -N カラーパレットを作成
gmt grdimage SLAautocor.nc -R-305/305/0/12 -Jx0.05/2 -BWSne+t"27N, 127E Correlation matrix" 塗りつぶし
gmt grdcontour SLAautocor.nc -C0.1 -A0.2+f16p -Gd4c -S4 等値線の描画
gmt basemap -R-305/305/0/12 -Jx0.05/2 -Bxaf5+l"Lag(distance, km)" -Byaf0.5+l"Lag(time, month)" X軸・Y軸の描画
gmt colorbar -DJRM+o0.3c/0+e+mc -Bx0.2 -By0.2 カラーパレットの描画
gmt end show