これまで扱ってきたのデータはアスキー形式データ。
ここからはバイナリデータを取り扱う。
バイナリについてはここをクリック。
全球海面温度（※）データを使って上昇トレンドを調べる
※データサイエンスには関係ないけど、どうでもいいようで、案外どうでもよくない話。
本講義で使うデータは、衛星観測で得られた「Sea Surface Temperature（略してSST ）」であり、学会や気象庁のページでも「海面水温」と呼ばれるものです。
しかし、同じ衛星プロダクトのLand Surface Temperatureは、「地表面温度」と呼ばれ、英語の綴りから考えてもSST は本来「海面温度」と呼称するべきであろう、と衛星観測関連シンポジウムの懇親会でJAXA/EORCの研究者Kさんから指摘され、今に至ります。
ちなみに前東北大教授の花輪先生は、塩分という言葉にすでに濃度の意味が含まれているので「塩分濃度」はおかしいと指摘されています。頼まれてもいないのに本学で私は「塩分濃度」撲滅運動を行ってきましたが、海面水温という用語に疑問を持たなかった自分を恥じています。
というわけで、今後はSSTを海面温度と呼称します。
定義を正しく理解すること、用語を正しく使うことの大切さを考えるきっかけになればと思い、自戒の念をこめてここに書いておきます。
9回目以降ではNOAAが配布するOISST ver2を用いる。
OISSTとはOptimum Interpolation Sea Surface Temperatureの略。
つまり最適内挿法で構築された海面温度のデータセットである。
1982年の図は以下のようになる。
～GMTスクリプト（上の図を書くには）～
gmt begin OISST1982
gmt xyz2grd 1982-sst.bin -R0.5/359.5/-89.5/89.5 -I1 -ZTLf -Gsst -di-999
gmt makecpt -Chaxby -T-2/30/1 -N
gmt grdimage sst -R0.5/359.5/-89.5/89.5 -JR16c
gmt coast -B+t"1982 OISST" -B -Gwhite
gmt colorbar -DjTR+o0c+w3c/0.4c -Bx -By+loC -I -F+gwhite+p1p
gmt end show
※このスクリプトは以降の課題提出で使うので理解しておくこと。
海面温度の上昇トレンド
From : Sea surface temperature from OISST : ICDC : Universität Hamburg (uni-hamburg.de)
上の図は全球の海面温度（高解像度）のトレンドを示した例（ハンブルグ大学）である。
OISSTデータで同様の計算結果が得られる。
ただし、ここで配布するデータは1度格子（ファイルサイズを小さくしたもの）なので、
上の図と全く同じ値にはならない。
全球海面温度の上昇トレンドの計算
9回目は各月の全球平均海面温度を計算し、最小二乗法を用いて上昇トレンドを計算する。
なお、そのままでは値が小さいので「10年で何度上がるか」という値にすること。
（１）データの準備
今回使用するデータはOISSTの月平均海面温度データである。
ファイル数は494個なので、サブフォルダに解凍するほうがよい。
OISSTデータ（DLすること）：Zip形式なので各自、DL後に解凍すること。
陸マスクデータ
ファイルリスト：ファイル名をフルパスで記述したもの。
（２）データの読み込み
このデータセットは1度×1度の格子データである。
1982年以降のデータが存在するので、ＤＯループとＯＰＥＮ文を使ってデータを読む。
今回は3次元データ（緯度、経度、時間）を扱う。
複数の2次元データを一気に読むことはできないので、
各ファイルを読み込み、変数dummyに一度格納してから、
3次元配列の変数sst3dに再格納する。
xvarは時間軸のための変数。
c234567---------------------------------
character*120 ifl(494)
integer num
integer*1 lmask(360,180) ! land mask
real dummy(360,180),msst
real sst3d(360,180,494),xvar(494)
real sst1d(494),sst1d2(482)
real slope,offset
C---------------------------------------
C read the list of FILENAME
C and make time axis
C---------------------------------------
open(10,file="list.txt")
read(10,*) num
do k=1,num
read(10,'(a120)') ifl(k)
write(*,'(a120)') ifl(k)
xvar(k)=k
enddo
close(10)
C---------------------------------------
C read SST files and make 3-D dataset
C---------------------------------------
do k=1,num
open(15,file=ifl(k),form='unformatted',
& access="direct",status="old",recl=360*180*4)
read(15,rec=1) dummy
do j=1,180
do i=1,360
sst3d(i,j,k)=dummy(i,j)
enddo
enddo
close(15)
enddo
C---------------------
C read LAND MASK
C---------------------
open(20,file='lmask.bin',form='unformatted',
& access="direct",status="old",recl=360*180)
read(20,rec=1) lmask
close(20)
（３）全球の平均海面温度の計算
陸マスクデータ(0は陸、１は海)を使って、陸域を無視しつつ、各年の平均海面温度を計算する。
平均海面温度の時系列sst1dを作成。
（４）上昇トレンドの計算（最小二乗法）
説明変数をxvar、従属変数をsst1dとして、最小二乗法により線形回帰式を求める。
得られた傾きslopeは月ごとの増減トレンドを示すので、120をかけて10年での増加率を得る。
以下はサンプルコード
c234567---------------------------------
! estimate LONGTERM-TREND of GLOBAL MEAN SST
character*120 ifl(494)
integer num
integer*1 lmask(360,180) ! land mask
real dummy(360,180),msst
real sst3d(360,180,494),xvar(494)
real sst1d(494),sst1d2(482)
real sst1df(494),sst1d2f(482)
real slope,offset
C---------------------------------------
C read the list of FILENAME
C and make time axis
C---------------------------------------
open(10,file="list.txt")
read(10,*) num
do k=1,num
read(10,'(a120)') ifl(k)
write(*,'(a120)') ifl(k)
xvar(k)=k
enddo
close(10)
C---------------------------------------
C read SST files and make 3-D dataset
C---------------------------------------
do k=1,num
open(15,file=ifl(k),form='unformatted',
& access="direct",status="old",recl=360*180*4)
read(15,rec=1) dummy
do j=1,180
do i=1,360
sst3d(i,j,k)=dummy(i,j)
enddo
enddo
close(15)
enddo
C---------------------
C read LAND MASK
C---------------------
open(20,file='lmask.bin',form='unformatted',
& access="direct",status="old",recl=360*180)
read(20,rec=1) lmask
close(20)
C----------------------------
C cal mean SST from 1982 to 2022
! sst1d: timeseries of Global-mean SST
C----------------------------
do k=1,num
sst1d(k)=0
icount=0
!do j=1,180
do j=30,150
do i=1,360
if(変数lmaskが１のとき)then
sst1d(k)=sst1d(k)+sst3d(i,j,k)
icount=icount+1
endif
enddo
enddo
sst1d(k)=sst1d(k)/real(icount)
write(*,*) k,sst1d(k)
enddo
msst=0; icount2=0
do k=1,num
msst=msst+sst1d(k)
icount2=icount2+1
enddo
msst=msst/real(icount2)
write(*,*)"msst=",msst
C----------------------------
C Least-square fitting (1-D)
C X: xvar , Y:sst1d
C----------------------------
B=0; C=0; D=0; E=0;
do k=1,num
B=B+
C=C+
D=D+
E=E+
enddo
slope=
offset=
! Global SST trend degreeC/yr
write(*,200) slope*120
200 format('Global-mean SST trend =',f12.4," degC/decade")
write(*,201) offset
201 format('OFFSET =',f12.4," degC")
C------------------------------
C Build SST timeseries based on results above
C------------------------------
do k=1,num
sst1df(k)=xvar(k)*slope+offset
enddo
C----------------------------
C OUTPUT-->write ascii file
C----------------------------
decim=1./24.!(365/12/2/365)
decim2=decim*2
yr=1981
mo=12
open(80,file='GlobalmeanSST.txt')
open(81,file='GlobalmeanSST2.txt')
do k=1,num
date=yr+decim+decim2*(mo-1)
write(80,100) date,sst1d(k)
write(81,100) date,sst1df(k)
mo=mo+1
if(mo.eq.13)then
mo=1
yr=yr+1
endif
enddo
close(80)
close(81)
100 format(2f10.4)
C----------------------------
C 13-month running mean
C----------------------------
do k=7,488
stmp=0
do i=1,13
stmp=stmp+sst1d(k+i-7)
enddo
sst1d2(k-6)=stmp/13.0
enddo
msst2=0
do k=1,482
msst2=msst2+sst1d2(k)
enddo
msst2=msst2/482.0
C----------------------------
!Least-square applied to 13-mon running mean
C----------------------------
B=0; C=0; D=0; E=0;
do k=1,482
B=B+
C=C+
D=D+
E=E+
enddo
slope2=
offset2=
! Global SST trend degreeC/yr
write(*,202) slope2*120
202 format('Global-mean SST trend(13mon) =',f12.4," degC/decade")
C------------------------------
C Build SST timeseries based on results above
C------------------------------
do k=1,482
sst1d2f(k)=xvar(k)*slope2+offset2
enddo
C----------------------------
C OUTPUT-->write ascii file
C----------------------------
yr=1982
mo=6
open(82,file='GlobalmeanSSTlp13.txt')
open(83,file='GlobalmeanSSTlp132.txt')
do k=1,482
date=yr+decim+decim2*(mo-1)
write(82,100) date,sst1d2(k)
write(83,100) date,sst1d2f(k)
mo=mo+1
if(mo.eq.13)then
mo=1
yr=yr+1
endif
enddo
close(82)
close(83)
stop !<-- we nned this, as always
end
GMTでプロットする。
gmt begin SSTtimeseries
gmt set FONT_HEADING 18p
gmt basemap -R1981.9/2023.1/18.2/19.5 -JX15c/10c -Bxag5+l"Year" -Byag0.5+l"Global Mean SST(degC)"
gmt plot GlobalmeanSST.txt -R1981.9/2023.1/18.2/19.5 -W0.2p -Bxag5 -Byag0.5 -BWSne -lGlobal_Mean_SST
gmt plot GlobalmeanSST.txt -Sc0.1c -Gwhite -W0.2p,black
gmt plot GlobalmeanSSTlp13.txt -W1.6p,green -Bxag5 -Byag0.5 -BWSne -lGlobal_Mean_SST_lp13
gmt plot GlobalmeanSST2.txt -W1.6p,red -lSST_Trend
gmt plot GlobalmeanSSTlp132.txt -W0.8p,blue -lSST_Trend_lp13
gmt legend -DjTL+o0.5c -F+gwhite+pthicker --FONT_ANNOT_PRIMARY=11p,Helvetica-Bold
gmt end show
