GMTで海面水温図＆動画を作る方法 - 自然環境保全のための周辺技術

気象庁の日別海面水温図のようなものをGMTで作りたいと思い、その方法をハッキング（？）試行錯誤したのでメモっときます。

↑できあがり　その1

[:W360]
↑できあがり　その2

↑できあがり　その3

1.水温データの取得

NEAR-GOOS 地域遅延モードデータベース (RDMDB)からdailysst(JMA)を取ってきます。
データを取得するには、簡単なユーザー登録が必要です。

2.水温データの変換

以下のgawk用スクリプとをgoos.awkという名前で保存してください。

BEGIN {
  FIELDWIDTHS = "3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3";
  j=0;
}
{
  if(2+(121*(day-1)) <= NR && NR <=121+(121*(day-1))){
   for (i = 0; i < NF; i++) {
    val = $(i+1);
     if(val==999 || val==888 || val==777){
      printf("%s,%s,NaN",120.125+i*0.25,49.875-j*0.25);
     }else{
      printf("%s,%s,%.1f",120.125+i*0.25,49.875-j*0.25, $(i+1)/10);
     }
    printf("\n");
   }
   j++;
  }
}

以下のコマンドでgawk スクリプトを実行します。引数dayで変換したい日を指定します。以下の例は、2011年8月1日のデータをout.txtに変換しています。（awkじゃなくて、gawkです。(FIELDWIDTHSを利用するため)。cygwinならawk==gawkかも）

gawk -v day=1 -f goos.awk 2011_08.dailysst > out.txt

3.GMTで作図

以下のGMTのコマンドで作図します。

gmtset BASEMAP_TYPE plain PLOT_DEGREE_FORMAT ddd:mm:ssF ANNOT_FONT_PRIMARY Times-Bold ANNOT_FONT_SIZE_PRIMARY 10p
xyz2grd out.txt -Gout.grd -I15m -R120.125/159.875/20.125/49.875
makecpt -T-2/32/1 -Z > test.cpt
grdimage out.grd -R120/160/20/50 -Jm0.3 -Ctest.cpt -K -P > test.ps
grdcontour out.grd -Jm0.3 -Gd4c -C1 -A5+gwhite+au -Wc1,50/50/50 -Wa4,20/20/20 -K -O  >> test.ps
pscoast -R120/160/20/50 -Bg5a10f1WSne -Jm0.3 -W -Cwhite -G100/100/100 -K -O >> test.ps
psscale -Bg1a5/:"[C]": -D13/5.5/10/0.5 -Ctest.cpt  -O >> test.ps
ps2raster  test.ps -A -Qt -Qg

test.jpgが出来上がっています。

その２

陸域を標高データにしてみました。
1分間隔の標高データはこちらから取得。http://topex.ucsd.edu/cgi-bin/get_data.cgi

gawk -v day=1 -f goos.awk 2011_08.dailysst > out.txt
gmtset BASEMAP_TYPE plain PLOT_DEGREE_FORMAT ddd:mm:ssF ANNOT_FONT_PRIMARY Times-Bold ANNOT_FONT_SIZE_PRIMARY 10p
surface out.txt -Gout.grd -I1m -R120/160/20/50 -T0.5
makecpt -T-2/32/1 -Z > test.cpt
surface srtm60.txt -I1m -Gsrtm60.grd -R120/160/20/50
grdgradient srtm60.grd -A270 -Gsrtm60.int -Ne0.9 -M
grdimage out.grd -R120/160/20/50 -Jm0.3 -Ctest.cpt -K -P > test.ps
pscoast -R120/160/20/50  -Cwhite -G100/100/100 -Jm0.3 -W -K -O >> test.ps
grdcontour out.grd -Jm0.3 -Gd4c -C1 -A5+gwhite+au -Wc1,50/50/50 -Wa4,20/20/20 -K -O  >> test.ps
echo "-10000 150 10000 150" > gray.cpt
pscoast -R120/160/20/50 -J -Di -Gc -O -K >> test.ps
grdimage srtm60.grd -Isrtm60.int -R120/160/20/50 -Jm0.3 -CGMT_relief.cpt -K -O >> test.ps
pscoast -R -J -Q -Bg5a10f1WSne:."2011/08/01": -K -O  >> test.ps
psscale -Bg1a5/:"[C]": -D13/5.5/10/0.5 -Ctest.cpt  -O >> test.ps
ps2raster test.ps -A -Qt -Qg

その３

2011年の海面水温を動画にしてみました。

シェルスクリプトで2011年の海水温データを1日ずつのjpg画像に変換して、それをアニメーションgifや動画に変換します。
動画への変換は、これを使用するといいかもしれません。http://www.tmpgenc.net/ja/j_download.html　
※動画にしたいjpgのファイル名は連番にしておく必要があります。
アニメーションgifへの変換は、ImageMagicのconvertコマンドでできます。

あらかじめ、海水温データを1年分とってきておいてください。標高データも。

surface srtm60.txt -I1m -Gsrtm60.grd -R120/160/20/50
grdgradient srtm60.grd -A270 -Gsrtm60.int -Ne0.9 -M

for j in `seq -w 1 12`
do
	for i in `seq -w 1 31`
	do
	gawk -v day=$i -f goos.awk 2011_$j.dailysst > out.txt
	fsize=`ls -s out.txt | cut -c 1-1`
	if [ $fsize -ne 0 ]; then
		gmtset BASEMAP_TYPE plain PLOT_DEGREE_FORMAT ddd:mm:ssF ANNOT_FONT_PRIMARY Times-Bold ANNOT_FONT_SIZE_PRIMARY 10p
		surface out.txt -Gout.grd -I1m -R120/160/20/50 -T0.5
		makecpt -T-2/32/1 -Z > test.cpt		
		grdimage out.grd -R120/160/20/50 -Jm0.3 -Ctest.cpt -K -P > sst$j$i.ps
		####for drawing lake with white color
		pscoast -R120/160/20/50  -Cwhite -G100/100/100 -Jm0.3 -W -K -O >> sst$j$i.ps
		grdcontour out.grd -Jm0.3 -Gd4c -C1 -Wc1,50/50/50 -Wa4,20/20/20 -K -O  >> sst$j$i.ps
		echo "-10000 150 10000 150" > gray.cpt
		pscoast -R120/160/20/50 -J -Di -Gc -O -K >> sst$j$i.ps
		grdimage srtm60.grd -Isrtm60.int -R120/160/20/50 -Jm0.3 -CGMT_relief.cpt -K -O >> sst$j$i.ps
		pscoast -R -J -Q -Bg5a10f1WSne:."2011/$j/$i": -K -O >> sst$j$i.ps
		psscale -Bg1a5/:"[C]": -D13/5.5/10/0.5 -Ctest.cpt -O >> sst$j$i.ps
		ps2raster sst$j$i.ps -A -Qt -Qg
		# for animation gif
#		convert sst$j$i.jpg sst$j$i.gif 
	fi
	done
done

#for animation gif
#convert sst*.gif sst_anime.gif

【参考】海面水温データについて

海面水温（Sea surface temperature(SST)のデータには、人工衛星NOAA、船舶やブイで計測、人工衛星AQUAに搭載されたAMSR-Eセンサからのデータがあるようです。先日打ち上げられた「しずく」にはAMSR-Eの後継となるAMSR2センサが搭載されています。
これらのデータをすべて考慮して日本周辺だけ取り出したのがNEAR-GOOS 地域遅延モードデータベース (RDMDB)ということのようです。ただし、このデータはリアルタイムでは取得できず、１ヶ月前までとなっています。

NOAAのデータ
http://www.ncdc.noaa.gov/cdr/operationalcdrs.html

GCOM-W1データ提供サービス（AMSR-E,AMSR2のデータ）
https://gcom-w1.jaxa.jp/auth.html

漁業情報サービスセンター（リアルタイムの海面水温データ(有料)）
http://www.jafic.or.jp/

海洋政策情報支援ツール（最近公開されたみたいです）
http://www5.kaiho.mlit.go.jp/kaiyo/