2024 能登震源地図作成備忘録

2024 能登震源地図作成備忘録．2024-01-28

１．震源データ作成

JMAデータは「2022_大阪北部地震震源地図作成備忘録.html」で作成した自動処理スクリプトが気象庁URL変更のため
動かなかったので，それを修正．動いたバージョンは

curl https://www.data.jma.go.jp/eqev/data/daily_map/202401[01-27].html|grep -e "2024 1" > test1.txt

ここで，|grep -e "2024 1"
を用いないと不要な行が残る．

これで作成したtest1.txtを元データとして用いる．
まず，plumaエディタの置換で度分記号を取り除く

2024 1 1 00:00 4.5 32°38.6'N 131°59.1'E   36     0.4 日向灘
2024 1 1 00:00 30.9 35°18.6'N 136°54.8'E    7     0.5 愛知県西部
2024 1 1 00:00 36.1 36° 3.0'N 137°32.9'E    6     0.5 岐阜県飛騨地方
2024 1 1 00:02 15.3 36° 0.7'N 137°31.4'E    8     0.0 岐阜県飛騨地方
2024 1 1 00:08 13.3 40° 3.1'N 140°38.9'E    8     0.3 秋田県内陸北部
2024 1 1 00:08 38.1 27°30.6'N 126°23.1'E 186     3.3 沖縄本島北西沖
2024 1 1 00:25 20.9 37°21.4'N 141°55.4'E   34     1.2 福島県沖

これをtest2.txtとする
2024 1 1 00:00 4.5 32 38.6 N 131 59.1 E   36     0.4 日向灘
2024 1 1 00:00 30.9 35 18.6 N 136 54.8 E    7     0.5 愛知県西部
2024 1 1 00:00 36.1 36 3.0 N 137 32.9 E    6     0.5 岐阜県飛騨地方
2024 1 1 00:02 15.3 36 0.7 N 137 31.4 E    8     0.0 岐阜県飛騨地方
2024 1 1 00:08 13.3 40 3.1 N 140 38.9 E    8     0.3 秋田県内陸北部
2024 1 1 00:08 38.1 27 30.6 N 126 23.1 E 186     3.3 沖縄本島北西沖
2024 1 1 00:25 20.9 37 21.4 N 141 55.4 E   34     1.2 福島県沖

これをさらに，10進数に直す．

awk '{ print $9+$10/60,$6+$7/60,$12,$13,$1,$2,$3,$4,$5 }' test2.txt > test3.txt

131.985 32.6433 36 0.4 2024 1 1 00:00 4.5
136.913 35.31 7 0.5 2024 1 1 00:00 30.9
137.548 36.05 6 0.5 2024 1 1 00:00 36.1
137.523 36.0117 8 0.0 2024 1 1 00:02 15.3
140.648 40.0517 8 0.3 2024 1 1 00:08 13.3
126.385 27.51 186 3.3 2024 1 1 00:08 38.1

次にMが"-"のデータを除く．
awk '{if($4 != "-" && $4<10) print $0}' test3.txt > test4.txt

最後に，深さでソートをかける．
sort -k 3 -n -r test4.txt > test_ds.dat

138.67 30.6467 460 3.5 2024 1 21 14:33 25.7
138.948 30.38 456 3.5 2024 1 7 03:42 40.9
139.912 29.1833 439 3.8 2024 1 21 10:13 51.4
137.99 32.6117 435 3.3 2024 1 25 23:01 49.3
146.228 46.325 434 3.6 2024 1 22 07:31 4.5
139.122 30.5483 428 4.4 2024 1 12 21:43 15.8
139.233 29.8783 422 3.1 2024 1 6 21:57 29.7
137.472 32.6483 406 3.1 2024 1 11 12:56 21.4
137.925 31.9883 399 3.6 2024 1 16 05:10 55.0

これでGMTに食わせる，震源ファイルは完成．

つぎにGMT（Ver3.4）での処理に移る．とりあえずの試作版．これで一応動く．

#!/bin/bash
# [Noto_AFS_JMA.sh]*****************************
#　GMTで能登地震余震地図を書くスクリプト
# data from JMA 日々の地震，震源リスト
# --data format--
# -> lon,    lat,    dep,mag (depth sorted)
#    132.4093 40.9786 681 3.4
#    130.5966 44.8036 644 6.3
#    35-color mode
# by Y.Okamoto 4/May 2008, 28th Jan.2024
# ************************************************
export PATH=$PATH:/home/seagull/GMT4/gmt-4.5.18/bin
#seis_file="JMA_20240116-20240126_ds.dat"
seis_file="test_ds.dat"
#out_file="JMA_20240116-20240126.eps"
out_file="test4.eps"        # temporary name

Pen_w=1     # pen width
cpt_file="34_cpt.csv"        # color code for Chroma Depth 3D Glasses
Area=134.0/140.0/35.5/39.0 # Draw area
scale=137.5/37.5/1:2300000   # center & scale
Mag_s=0.002                   # Eq plot scale factor

D_max=34
D_min=0
# int(34/34)
D_inc=1
i=0

#col 0   1   2   3   4   5   6   7   8   9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34
R=( 0   0   0   0   0   0   0   0   0   0   0   7 47 83 113 140 163 183 199 213 224 233 240 246 250 253 255 255 255 255 255 255 255 255 255)
G=( 0   4   9 16 24 34 45 57 71 87 104 123 142 164 187 211 237 248 227 207 188 169 151 134 118 102 87 73 59 47 35 23 13   3   0)
B=(255 255 255 254 250 243 234 222 209 193 175 155 133 109 82 53 22   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0)

echo "0 3 10 0 0 6 2024Noto Aftershock Map (data from JMA)" > title.txt
echo "0 0 10 0 0 6 script:Noto.sh by Yoshio Okamoto 2024" > credit.txt

# sea: water blue -S180/220/240
pscoast -R$Area -Jt$scale -Df -Gwhite -S180/220/240 -B1.0g1.0 -W$Pen_w/0/0/0 -K> $out_file

# Hypo plot for deeper earthquakes in deep blue
awk "{if (\$3 >= $D_max+100) { print \$1,\$2,\$4*\$4*\$4*$Mag_s }}" $seis_file | psxy -J -R -V -O -K -Sc -L -W$Pen_w/0/0/10 -N>> $out_file
awk "{if (\$3 >= $D_max+50) { print \$1,\$2,\$4*\$4*\$4*$Mag_s }}" $seis_file | psxy -J -R -V -O -K -Sc -L -W$Pen_w/0/0/40 -N>> $out_file
awk "{if (\$3 >= $D_max) { print \$1,\$2,\$4*\$4*\$4*$Mag_s }}" $seis_file | psxy -J -R -V -O -K -Sc -L -W$Pen_w/0/0/120 -N>> $out_file

let DD=D_max
let D=D_max-D_inc

# Hypo plot for each depth range
while [ $i -le 34 ]; do
echo The color counter is $i
echo The RGB values are ${R[$i]} ${G[$i]} ${B[$i]}
echo The Depth Range is $D km to $DD km
awk "{ if ((\$3 >= $D) && (\$3 < $DD)) { print \$1,\$2,\$4*\$4*\$4*$Mag_s }}" $seis_file | psxy -J -R -V -O -K -Sc -L -W$Pen_w -G${R[$i]}/${G[$i]}/${B[$i]} -N>> $out_file
let i=i+1
let D=D-D_inc
let DD=DD-D_inc
done

# 県境
psxy Japan_PB.dat -R -J -B -W1 -m -N -O -K >> $out_file

# Legend -R0/6/0/4:dummy
awk '{ print $1,$2,$3,$4,$5,$6,$7}' mag2.txt | pstext -Jx1.0c -R0/6/0/4 -D-0.15/1 -O -K -G0/0/0 -N>> $out_file
awk '{ print $1,$2+0.16,$2*$2*$2*0.002 }' mag2.txt | psxy -Jx1.0c -R -D0.8/1 -O -K -Sc -L -W1/0/0/0 -N>> $out_file

# scale                 x / y /length/width
psscale -C$cpt_file -D0.4/5.7/-0.6/0.10 -B10:km: -K -O >> $out_file
pstext title.txt -Jx1.0c -R -D9/20.8 -G0/220/0 -G0/220/0 -N -O -K>> $out_file
pstext credit.txt -Jx1.0c -R -D2.0/-0.65 -G0/220/0 -G0/220/0 -N -O >> $out_file

＜震源データ作成の古いバージョン＞
https://www.data.jma.go.jp/eqev/data/daily_map/20240126.html
などから日別にコピペ．

これをJMA_20240116-20240126.txtとする．下記のような感じ．

2024  1 26 23:47 38.6  33°58.0'N 134°31.9'E   37     1.0  徳島県北部                   
2024  1 26 23:47 50.7  24°26.6'N 124°14.7'E   14     0.5  石垣島近海                   
2024  1 26 23:50  8.2  37°11.1'N 136°42.9'E    4     1.4  石川県能登地方                 
2024  1 26 23:51 58.4  37°20.3'N 136°49.3'E   10     3.3  石川県能登地方                 
2024  1 26 23:53 42.7  37°41.3'N 141°33.0'E   75     1.1  福島県沖                    
2024  1 26 23:54  8.4  34°18.2'N 140° 1.4'E   19     1.7  房総半島南方沖                 
2024  1 26 23:54 43.8  32° 4.5'N 130°23.5'E   11     0.2  鹿児島県薩摩地方                
2024  1 26 23:54 58.7  36°59.7'N 136°26.1'E   11     2.4  石川県西方沖                  
2024  1 26 23:55 55.6  37°22.0'N 136°52.3'E    6     1.6  石川県能登地方                 
2024  1 26 23:56 32.6  37°12.2'N 136°48.0'E   11     1.3  石川県能登地方   

これをGMTで描ける形式に直す．

ここで，昔大阪北部地震備忘録で作ったスクリプトによる取得．はうまくできなかった．URLは合っているが
なぜか，震源リスト部分がうまく展開して取得できていないように思える．本当はこのリストの状態でデータを提供すべきだと思う．

How to compile the JMA Web site EQ_list for Pov-Ray  Y.Okamoto 2021.10.01　from Y.Okamoto 2016.07.27
The most recent data 2020 to 2021 are from

curl http://www.data.jma.go.jp/svd/eqev/data/daily_map/202001[01-31].html|grep -e 2020 > test1.txt
awk '{if ($1=="2024" && $2=="1") print $0}' test1.txt > 202001.dat
cat *.dat > test.dat
So, we can get the dara from 2020 to 2021 September.

仕方なく上記ファイルから，エディタで度，分を除く（これ本当にめんどくさい，頼むからJMAこの度分表記止めてくれ！あと地震の年号表記も！） 

2024  1 16 15:55 50.2  35 22.2 N 141  7.9 E   15     2.3  千葉県東方沖                  
2024  1 16 15:56 23.9  27 22.0 N 128 56.3 E   30     1.0  沖縄本島近海                  
2024  1 16 15:56 37.5  37 24.5 N 137  8.9 E    6     2.6  石川県能登地方                 
2024  1 16 15:58  2.6  33 12.8 N 132 24.2 E   38     2.2  豊後水道                    
2024  1 16 15:58 34.0  37 44.1 N 137 37.8 E   15     2.2  佐渡付近                    
2024  1 16 16:03 11.2  37 23.2 N 137 11.2 E    9     0.8  石川県能登地方                 
2024  1 16 16:05 34.0  37 23.0 N 141 49.3 E   37     2.1  福島県沖                    
2024  1 16 16:07  7.4  31 54.9 N 130  9.7 E    5     0.3  薩摩半島西方沖                 
2024  1 16 16:09 32.5  37 12.3 N 136 41.1 E    8     1.3  能登半島沖                   
2024  1 16 16:12 32.4  21 49.1 N 120 55.0 E   53     4.6  台湾付近   

これを，GMTに直すが，まず空行を取り除く．
awk '{if ($1=="2024" && $2=="1") print $0}' test.txt > 2024Noto_temp.dat
次に60度分秒->10進

awk '{ print $9+$10/60,$6+$7/60,$12,$13,$1,$2,$3,$4,$5 }' 2024Noto_temp.dat > JMA_20240116-20240126.dat

次にMが"-"のデータを除く．
awk '{if($4 != "-" && $4<10) print $0}' JMA_20240116-20240126.dat > JMA_20240116-20240126_v1.dat
さらに

でとりあえず震源データ作成が完了．次にGMTを使う．