気象庁の「地震月報(カタログ編)」の震源データを使って震央の時空間分布図を作ってみました．これも最初はresearchmapに書いたものです．

前の記事と同じく，『歴史のなかの地震・噴火――過去がしめす未来』に載せた図の描き方を書いてみよう，というものです．この本のもとになった授業で2016年熊本地震の解説をするにあたり，気象庁の報道発表資料から図を転載していました．本にするにあたっては，余分な情報もあり，また白黒印刷対応のため，作りなおすことにしました．

対象とするデータ

気象庁「地震月報(カタログ編)」の震源データ

使ったツール・環境

gmt6 https://www.generic-mapping-tools.org/
awk
WSL2上のubuntu Ubuntu 20.04.2 LTS 「Windows 10 用 Windows Subsystem for Linux のインストール ガイド」

描き方

まずデータをダウンロードします．年ごとのファイル（新しいところは月ごとの場合もあり）になっています．

今回は熊本地震のあった2016年4月の図をつくりますので，2016年のデータ（h2016）をダウンロードしました．適宜解凍すると，震源ファイルフォーマットに従ったテキストファイルです．今後も自分でデータを使ってみようという方は，フォーマットや気象庁の地震カタログの解説もぜひ読んでみてください．

以下のスクリプト中では，気象庁の震源ファイルフォーマットをいったんWINシステムの震源データベースファイルの形式に変換しています．必ずしもこのフォーマットでなくてもいいのですが，慣れているのでそうしています．

時空間の範囲でフィルタして一時ファイルに書き出していますが，このくらいのデータ量であれば，毎回もとのデータを使ってもたいした違いはないと思います．4月14日から31日までの地震を描いています．

震源データへのパスなどは適宜変更してください．ラベルなどを日本語で書くときは，スクリプトの文字コードはEUC-JPにしましょう．震央の◯の大きさはいろんな流儀があると思います．

#!/bin/sh
# 地震月報の震源ファイルをWIN final形式に変換する関数
jma2final () {
  awk '{
    printf "%04d %02d %02d %02d %02d %5.2f %7.4f %8.4f %6.2f %3.1f %s\n",\
    substr($0,2,4),substr($0,6,2),substr($0,8,2),\
    substr($0,10,2),substr($0,12,2),substr($0,14,2)+substr($0,16,2)/100,\
    substr($0,22,3)+(substr($0,25,2)+substr($0,27,2)/100)/60,\
    substr($0,33,4)+(substr($0,37,2)+substr($0,39,2)/100)/60,
    substr($0,45,3)+substr($0,48,2)/100,
    substr($0,53,2)/10,
    substr($0,69,24);
  }'
}

# 図のサイズ
MAPWIDTH=12
TIMEWIDTH=14

# 描画範囲
YMIN=32.0
YMAX=33.5
XMIN=130
XMAX=132
TICK=a1f0.25
DRY=white
WET=grey80

# WIN final形式に変換，必要な範囲のデータだけ抽出，一時ファイルに保存
cat h2016 | \
jma2final | \
awk '$2==4{print $0}' | \
awk '$8>='$XMIN'&&$8<='$XMAX'&&$7>='$YMIN'&&$7<='$YMAX'{print $0}' > tmp.$$

gmt begin kumamoto ps,png

gmt set MAP_FRAME_TYPE plain MAP_DEGREE_SYMBOL degree \
FORMAT_DATE_MAP mm/dd PS_MEDIA a4 \
FONT_ANNOT_PRIMARY 12p,Helvetica,black \
FONT_LABEL 16p,GothicBBB-Medium-EUC-H,black

# 震央分布
gmt coast -JM$MAPWIDTH -R$XMIN/$XMAX/$YMIN/$YMAX -Dh -Wblack -B$TICK -G$DRY -S$WET -X3 -Y15

## 線分A-B（時系列を投影する）
gmt plot <<END -W
131.4 33.4
130.4 32.2
END

## M<5
cat tmp.$$ | \
awk '$10<5{print $0}' | \
awk '$7>='$YMIN'&&$7<='$YMAX'&&$8>='$XMIN'&&$8<='$XMAX'&&$9<=20{print $8,$7,$10/(10-$10)/5}' | \
gmt plot -Sc -Wthinnest,grey50

## M>=5
cat tmp.$$ | \
awk '$10>=5{print $0}' | \
awk '$7>='$YMIN'&&$7<='$YMAX'&&$8>='$XMIN'&&$8<='$XMAX'&&$9<=20{print $8,$7,$10/(10-$10)/5}' | \
gmt plot -Sc -Wthinner,black

## ラベルA-B
gmt text << END -F+f12p,0,black+jLT
131.4 33.4 A
130.4 32.2 B
END

## 凡例
gmt plot <<END -Sc -Wthinner,black -N
132.2 32.7 0.467
132.2 32.5 0.3
132.2 32.3 0.2
END

gmt text << END -F+f12p,0,black+jCM -N
132.2 32.62 M7.0
132.2 32.42 M6.0
132.2 32.22 M5.0
END

# 時系列
## 線分A-Bの長さ
DIST=`gmt mapproject -G131.4/33.4+uk <<END | awk '{print $3}'
130.4 32.2
END
`
gmt basemap -JX$TIMEWIDTH/-8 -R2016-04-14T/2016-05-01T/0/$DIST -Bxafg+l"2016年" -Byafg+l"距離, km" -BWSne \
-X-1 -Y-10

## M<5
cat tmp.$$ | \
awk '$10<5{print $0}' | \
awk '$7>='$YMIN'&&$7<='$YMAX'&&$8>='$XMIN'&&$8<='$XMAX'&&$9<=20\
{print $8,$7,$10/(10-$10)/5,$1"-"$2"-"$3"T"$4":"$5":"$6}' | \
gmt project -C131.4/33.4 -E130.4/32.2 -W-25/25 -Lw -Q | \
awk '{print $4,$5,$3}' | \
gmt plot -Sc -Wthinnest,grey50

## M>=5
cat tmp.$$ | \
awk '$10>=5{print $0}' | \
awk '$7>='$YMIN'&&$7<='$YMAX'&&$8>='$XMIN'&&$8<='$XMAX'&&$9<=20\
{print $8,$7,$10/(10-$10)/5,$1"-"$2"-"$3"T"$4":"$5":"$6}' | \
gmt project -C131.4/33.4 -E130.4/32.2 -W-25/25 -Lw -Q | \
awk '{print $4,$5,$3}' | \
gmt plot -Sc -Wthinner,black

## ラベルAとB
gmt text << END -F+f12p,0,black+jLT -N
2016-5-1T12 0 A
2016-5-1T12 $DIST B
END

gmt end

# 一時ファイルを削除
rm tmp.$$

例によって，出来上がりは本のほうで確認してください．

マグニチュード5以上の震央を赤い◯でしめすと下のようになります．

別の図についても記事を書いています．合わせて参考になれば幸いです．