南海トラフの巨大地震モデル検討会の公開データを使ってみた

内閣府の南海トラフの巨大地震モデル検討会において検討されたデータを使ってお絵かきしてみます．researchmap のほうに書いていたのですが，noteも試してみます．

『歴史のなかの地震・噴火――過去がしめす未来』に最大クラスの想定の図を掲載したかったのですが，同検討会の報告書の図はカラーで，白黒だとうまく色が出ないと思われたため，グレースケールで作りなおすことにしました．

自分でデータを見てみたいという方のお役に立てば幸いです．さまざまなケースについて，選んだ地域を拡大して作図するなどできると思います．あくまである（いくつかの）シナリオに沿った想定ですが，さまざまにプロットしてみることで，想定への理解が深まると思います．想定に用いられた手法や仮定にも興味を持っていただくとさらによいと思います．

データの詳細に興味はないという方も，行政の判断に使われるようなデータが，このような形で公開されている（オープンデータ・オープンガバメント）ということを知っていただければと思います．

対象とするデータ

強震断層モデル(1)データセットA https://www.geospatial.jp/ckan/dataset/1201 ※ダウンロードにはG空間情報センターWebサイトでのユーザ登録が必要です
国土数値情報・行政区域データ https://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-N03-v2_4.html#prefecture00

使ったツール・環境

gmt6 https://www.generic-mapping-tools.org/
ogr2ogr (GDAL 3.0.4, released 2020/01/28) https://gdal.org/programs/ogr2ogr.html
awk
WSL2上のubuntu Ubuntu 20.04.2 LTS 「Windows 10 用 Windows Subsystem for Linux のインストール ガイド」

もちろん，QGISでなどのGISソフトウェアで描画する，というのもありです．

描き方

震度分布データ

強震断層モデル(1)データセットAを https://www.geospatial.jp/ckan/dataset/1201 からダウンロードします．

地震のケース1から5のそれぞれと，地震のケース1から5までの震度のうち最大となる震度（最大クラス）のデータがあります．まずは最大クラスのデータを使ってみましょう．

必要なデータ列（経度，緯度と震度）を抜きだし，GMTのグリッドファイルに変換します．このデータの緯度，経度はメッシュの南西端と北東端の経度（6-7列目と8-9列目），緯度（日本測地系JGD2000）で表現されています．足して2で割ったメッシュの中央を代表点とします．震度（小数点以下第一位まで表示）が4列目です．地域メッシュについては総務省統計局の「地域メッシュ統計について」を参照してください．

合わせてデータを間引きします．強震断層モデル(1)データセットAは250 mメッシュでの計算値で，日本全国を描画するときは，そんなに細かいものは必要ないと思います．とりあえず1 kmメッシュ（緯度、経度方向の間隔がそれぞれ30秒、45秒）にしてみます．もっと間引いてもよいかもしれません．-rオプションでpixel registration（グリッド＝メッシュのなかに代表点がある），-Auオプションでグリッド中の最大値を採用することを指定します．

tail +2 Eq_champ.csv | awk -F, '{printf "%.6f %.6f %3.1f\n",($6+$8)/2,($7+$9)/2,$4}' | gmt xyz2grd -GEq_champ_1km.grd -I45s/30s -R128:0/142:0/30:0/38:0 -fg -r -Au -D+tEq_champ+rCAO-Nankai

都道府県境データ

都道府県境のデータを https://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-N03-v2_4.html#prefecture00 からダウンロードします．SHAPE形式データが含まれています．

SHAPE形式データをGMTで扱えるデータ形式に変換します．GDALを含むようにインストールしたGMTだと，直接読めますが，何度も描きなおす場合は変換しておいたほうが早いと思います．GMTのマニュアルには，さらにバイナリデータに変換してデータファイルを小さくする方法も書かれています．

ogr2ogr -f "GMT" N03-19_190101.gmt N03-19_190101.shp

県境の経度，緯度が並んだテキストデータになります．「>」を区切り文字とするmulti segmentデータです．

描画

データの準備ができたら，GMTで描画します．今回は，震度分布をメッシュの色であらわして，その上に県境を重ねて描いています．凡例はまあお好みで．

日本語でラベルなどを書きたい場合，シェルスクリプトの文字コードはEUC-JPで保存しましょう．

#!/bin/sh
# 描画範囲
YMIN=30
YMAX=38
XMIN=128
XMAX=142
gmt begin nankai_soutei_modern
gmt set MAP_FRAME_TYPE fancy MAP_DEGREE_SYMBOL degree PS_MEDIA a4
# 震度のカラーパレット
cat > intensity.cpt << END
0.0 grey100 3.5 grey100
3.5 grey90 4.5 grey90
4.5 grey80 5.0 grey80
5.0 grey70 5.5 grey70
5.5 grey50 6.0 grey50
6.0 grey30 6.5 grey30
6.5 grey0 9.0 grey0
B 255
F 0
N grey100
END
# 震度分布
gmt grdimage Eq_champ_1km.grd -JM15 -R$XMIN/$XMAX/$YMIN/$YMAX -Cintensity.cpt -Y10 -Vv
# 県境
cat N03-19_190101.gmt | \
gmt plot -Wfaint,black -Vv
# 凡例
gmt text << END -D0/-0.25 -F+f16p,37,black+jML -V
136.5 33.5 震度
136.5 33   7
136.5 32.5 6強
136.5 32   6弱
136.5 31.5 5強
136.5 31   5弱
136.5 30.5 4
END
# 7
gmt plot << END -Ggray0 -L
137.5 33
138.0 33
138.0 32.5
137.5 32.5
END
# 6+
gmt plot << END -Ggray30 -L
137.5 32.5
138.0 32.5
138.0 32
137.5 32
END
# 6-
gmt plot << END -Ggray50 -L
137.5 32
138.0 32
138.0 31.5
137.5 31.5
END
# 5+
gmt plot << END -Ggray70 -L
137.5 31.5
138.0 31.5
138.0 31
137.5 31
END
# 5-
gmt plot << END -Ggray80 -L
137.5 31
138.0 31
138.0 30.5
137.5 30.5
END
# 4
gmt plot << END -Ggray90 -L
137.5 30.5
138.0 30.5
138.0 30
137.5 30
END
gmt end

できあがりは，こちらの書籍で確認してみてください．誤植というかコピペの消し忘れがあります．お恥ずかしい．

カラーの図

カラーパレットのところを以下のように変更すれば，カラーの図が作れます．配色は，気象庁の震度の色にしました．南海トラフの巨大地震モデル検討会の報告書の図とは色が違います．

# 震度のカラーパレット（気象庁カラーコード）
cat > intensity.cpt << END
0.0 0 65 255 3.5 0 65 255
3.5 250 230 150 4.5 250 230 150
4.5 255 230 0 5.0 255 230 0
5.0 255 253 0 5.5 255 153 0
5.5 254 40 0 6.0 254 40 0
6.0 165 0 33 6.5 165 0 33
6.5 180 0 104 9.0 180 0 104
B 255
F 0
N grey100
END

こんな感じになります．