ここでは最近GMTで作成した図を紹介します。
2024/4/28 オクラホマ
図
スクリプト
-
震源
#!/bin/bash gmt begin oklahoma_hypo png gmt basemap -JM12 -R-103/-93/33/41 -Bafg -BWSNE+t"Oklahoma_2001-2023_USGS" gmt coast -Df -W0.25 -Glightgray -Slightcyan -LJBR+jTR+o0/1.5+w100+f -N2 gmt makecpt -Cseis -T0/40/1 -Z awk -F"," '{print $3, $2, $4, $5*0.1}' query.csv | gmt plot -Sc -C -W0.25 gmt colorbar -DJBR+jBL+o1.5/0+w-5/0.2 -Ba10f10+l"Depth (km)" gmt inset begin -DjBR+w2.5/1.4+o0.1c -F+gwhite+p0.25 -M0.01 gmt set MAP_FRAME_TYPE = "plain" gmt basemap -JM2.5 -R-125/-65/24.6/50 -Bafg -BWSNE gmt coast -Df -W0.25 -Glightgray -Slightcyan -LJBR+jTR+o0/1.5+w100+f -N2 gmt plot -W0.25,blue -L <<END -103 33 -93 33 -93 41 -103 41 END gmt inset end gmt basemap -JX12/-12 -R-103/-93/0/40 -Bxafg+l"longitude" -Byafg+l"depth" -BWSne -Y-15 awk -F"," '{print $3, $4, $4, $5*0.1}' query.csv | gmt plot -Sc -W0.01 -C gmt colorbar -DJBR+jBL+o1.5/0+w-5/0.2 -Ba10f10+l"Depth (km)" gmt end
-
N-T図
#!/bin/bash input="query2.csv" # ${input}のうち2001から2023が含まれる行の数を数える gmt begin count png awk -F"," '{print $1}' ${input} | gmt histogram -JX12 -R2001/2023/0/4000 -Bxa1fg+a45+l"Year" -By1000a+l"Count" -BWSne+t"Oklahoma_2001-2023_USGS_M2.5+" -Glightgray -W -D -T1 gmt inset begin -DjTL+w5/2.8+o0.1c -F+gwhite+p0.25 -M0.01 gmt set MAP_FRAME_TYPE = "plain" gmt basemap -JM5 -R-125/-65/24.6/50 -Bafg -BWSNE gmt coast -Df -W0.25 -Glightgray -Slightcyan -LJBR+jTR+o0/1.5+w100+f -N2 gmt plot -W3.5,blue -L <<END -103 33 -93 33 -93 41 -103 41 END gmt inset end gmt end
謝辞
データ出典:USGS
記して感謝申し上げます。USGSはデータ取得しやすい!便利!!
2024/4/17 豊後水道
図
スクリプト
- 震源
#!/bin/bash
input="./h2023.csv"
gmt begin bungo_suido png
gmt basemap -JM12 -R132.2/132.6/33.0/33.4 -Bafg -BWSNE+t"2023/01/01-2023/12/31_JMA"
gmt coast -Df -W0.25 -Glightgray -Slightcyan
gmt makecpt -Cseis -T0/50/1 -Z
awk -F"," '{print $11,$9,$13,$15*0.1}' ${input} | gmt plot -Sc -W0.25 -C
gmt plot -Sa -W0.25 -C << END
132.4 33.2 41 0.62
END
gmt colorbar -DJBR+jBL+o1.5/0+w-5/0.2 -Baf+l"depth"
gmt basemap -JX12/-12 -R132.3/132.5/0/50 -Bxafg+l"longitude" -Byafg+l"depth" -BWSne -Y-15
awk -F"," '{print $11,$13,$13,$15*0.1}' ${input} | gmt plot -Sc -W0.01 -C
gmt plot -Sa -W0.25 -C -l"202404172314" << END
132.4 41 41 0.65
END
gmt colorbar -DJBR+jBL+o2.5/0+w-5/0.2 -Baf+l"Depth"
gmt end
- 発震機構解
#!/bin/bash
# ローカルのファイル先を指定したりURLを入力してください。
input="./mecha/mecha2011-2023.txt"
rm -f mecha.txt mecha_aki.txt
# データ収集(URLから取得する場合のみ)
# wget $url -O $input
## 安芸・リチャーズ
## 1,2行目は削除してカンマ区切りにする
sed -e 1,2d $input | sed -e 's/ */,/g' | awk -F"," '{print $4,$3,$5,$8,$9,$10,$6}' > mecha.txt
## 負の値の対応
sed -e 's/;[0-9]*//g' mecha.txt | sed -e 's/-([0-9]*)-([0-9]*)/-1/g' > mecha_aki.txt
# GMT
gmt begin bungo_suido_mecha png
gmt basemap -JM12 -R132.2/132.6/33.0/33.4 -Bafg -BWSNE+t"2011-2023_F-net"
gmt coast -Df -W0.25 -Glightgray -Slightcyan -LJBR+jTR+o0/1+c40+w10+f
gmt makecpt -Cseis -T0/50/1 -Z
gmt meca "mecha_aki.txt" -Sa0.5 -C
gmt meca -Sa0.5 -C << END
132.4 33.2 41 166.1 69.6 -114 6.2 0 0 202404172314
END
gmt colorbar -DJBR+jBL+o2/0+w-5/0.2 -Baf+l"Depth(km)"
gmt end
注意事項
- 震源の位置はAQUAシステム メカニズム解カタログに基づくため今後変わる可能性があります。
- 2023年の気象庁一元化震源データは今後見直される可能性があります。
謝辞
以下のデータを使用させていただきました。ここに記して感謝申し上げます。
- 震源データとして、気象庁の「地震カタログ」を一部改変したもの。
気象庁|震源データ
気象庁が提供するページです
- 発震機構解のデータとして、国立研究開発法人防災科学技術研究所のF-netデータ
防災科研機関リポジトリ
- HinetPy
Dongdong Tian, 2020 , HinetPy: A Python package to request and process seismic waveform data from Hi-net. , doi: 10.5281/zenodo.3885779
2024/8/8 日向灘
図
スクリプト
#!/bin/bash
# 131.8593,31.704,32,0.35
rm danmen.txt danmen.csv
magnitude="2"
input="./h1998-h2024_m${magnitude}over.csv"
lon1="130"
lon2="134"
lat1="31"
lat2="34"
output="1998-2024_m${magnitude}over"
title="1998-2024 Dep<100km M${magnitude}>"
depth_max="100"
gmt set GMT_VERBOSE="long"
gmt begin $output png
gmt basemap -R$lon1/$lon2/$lat1/$lat2 -JM12 -Bafg -BWSNE+t"$title"
gmt coast -W0.25 -Df -LJBR+jTR+w100+f+o0/1 -Slightcyan -Glightgray
gmt makecpt -Cseis -T0/$depth_max/1 -Z
# 9カラム目が緯度、11カラム目が経度、13カラム目が深度、15カラム目がマグニチュード
# awkに変数を渡して13カラム目が$depth_max以下の行を抽出
# 15カラム目が小さい順に並びかえる
awk -F, -v depth_max="${depth_max}" '$13 <= depth_max {print $11, $9, $13, $15*0.05}' "$input"| sort -t ' ' -n -k4 | gmt plot -Scc -W0.25 -C
# 20240808
echo "131.8593,31.704,32,0.35" | gmt plot -Ss -W0.25 -C
#
gmt colorbar -DJBR+jBL+o1.5/0+w-5/0.2 -Baf+l"Depth (km)"
# lon1 lon2 lat1 lat2の中の震源の個数を図の左下にN=で記載する
awk -F, -v lon1="${lon1}" -v lon2="${lon2}" -v lat1="${lat1}" -v lat2="${lat2}" '$11 >= lon1 && $11 <= lon2 && $9 >= lat1 && $9 <= lat2 {n++} END {print lon1, lat1, "N="n}' "$input" | gmt text -F+f12p,Helvetica-Bold+jBL -D0/0.1
gmt legend -DJTL+jTL+o0.2c -F+p2p,black+gwhite << END
N 1
S 0.3 c 0.10 red 0.5p,black 0.8 M2
G0.25
S 0.3 c 0.15 red 0.5p,black 0.8 M3
G 0.25
S 0.3 c 0.20 red 0.5p,black 0.8 M4
G 0.25
S 0.3 c 0.25 red 0.5p,black 0.8 M5
G 0.25
S 0.3 c 0.30 red 0.5p,black 0.8 M6
G 0.25
S 0.3 c 0.35 red 0.5p,black 0.8 M7
END
# 断面図の作成
start_lon="131"
start_lat="32.5"
end_lon="133"
end_lat="31.5"
width="15"
echo -e "$start_lon $start_lat n$end_lon $end_lat" | gmt plot -W1,blue
# startからendまでの距離をmapprojectで計算
distance=(`echo "$start_lon $start_lat" | gmt mapproject -G$end_lon/$end_lat+uk`)
awk -F, -v depth_max="${depth_max}" '$13 <= depth_max {print $11, $9, $13, $15*0.05}' "$input"| sort -t ' ' -n -k4 | gmt project -C$start_lon/$start_lat -E$end_lon/$end_lat -W-${width}/${width} -Lw -Q > danmen.txt
# danmen.txtをタブ区切りからカンマ区切りに変換
sed 's/t/,/g' danmen.txt > danmen.csv
gmt basemap -JX12/-12 -R${lon1}/${lon2}/0/${depth_max} -Bxafg+l"longitude" -Byafg+l"depth[km]" -BWSne -Y-14
awk -F, '{print $7,$3,$3,$4}' "danmen.csv" | gmt plot -Sc -W0.25 -C -N
# 20240808
echo "131.8593,32,32,0.35" | gmt plot -Ss -W0.25 -C
#
# danmen.csvの行数を数えて図の左下に"N="で記載する。また、widthを記載する
awk -F, 'END {print "N="NR, "Width="width}' width="${width}km" "danmen.csv" | gmt text -F+f12p,Helvetica-Bold+cBL -D0.1/0.1
gmt end
注意事項
- 震源の位置はAQUAシステム メカニズム解カタログに基づくため今後変わる可能性があります。
- 2024年の気象庁一元化震源データは今後見直される可能性があります。
謝辞
以下のデータを使用させていただきました。ここに記して感謝申し上げます。
- 震源データとして、気象庁の「地震カタログ」を一部改変したもの。
気象庁|震源データ
気象庁が提供するページです
- 発震機構解のデータとして、国立研究開発法人防災科学技術研究所のF-netデータ
防災科研機関リポジトリ
- HinetPy
Dongdong Tian, 2020 , HinetPy: A Python package to request and process seismic waveform data from Hi-net. , doi: 10.5281/zenodo.3885779
2024/9/1 富士山
高いぞ高いぞ…
図
スクリプト
- 震源分布図
#!/bin/bash
input="./h2014-h2024zan.csv"
gmt begin mt_fuji pdf,png
gmt set GMT_VERBOSE=long
gmt basemap -R138.5/139.5/35/36 -Baf -BWSNE+t"Mt. Fuji(2014-2024)"
gmt coast -Df -W0.25
gmt makecpt -Cgeo -T-8000/8000/200 -Z
gmt grdcut @earth_relief_15s -R138.5/139.5/34/36 -Gdem.nc
gmt grdgradient dem.nc -Ggrad.grd -A45 -Ne0.8
gmt grdimage dem.nc -Igrad.grd -C
gmt makecpt -Cseis -T0/100/1 -Z
# 11カラム目は138-140, 9カラム目は34-36, 15カラム目の大きい順に並び変える
awk -F"," '{if($11>=138.5 && $11<=139.5 && $9>=35 && $9<=36 && $13<=100) print $11,$9,$13,$15*0.1}' $input | sort -k4 -n | gmt plot -Scc -W0.25 -C
# 138.9647, 35.4896 , 24.0, 4.3
gmt plot -Sa -W0.25 -C << END
138.9647 35.4896 24 0.43
END
gmt coast -Df -W0.25 -LJRB+jTR+o0/1.5+w50+f
gmt colorbar -DJBR+jBL+w-5/0.2+o1.5/0 -Baf+l"Depth(km)"
gmt end
- N-T図
#!/bin/bash
input="./h2014-h2024zan.csv"
north=35.6
south=35.3
west=139.3
east=138.9
# 9カラム目が緯度、11カラム目が経度なので、絞り込む
awk -F"," '{if($11>=138.5 && $11<=139.5 && $9>=35 && $9<=36 && $13<=100) print $0}' $input > range_hypo.csv
# 2カラム目が年、3カラム目が月、4カラム目が日なのでyyyy-mm-ddに変換
awk -F"," '{print $2"-"$3"-"$4","$0}' range_hypo.csv > range_hypo_date.csv
# 横軸に日付、縦軸にその日に発生した地震の回数をプロットする
awk -F"," '{print $1}' range_hypo_date.csv | sort | uniq -c | awk '{print $2","$1}' > date_count.csv
# 2カラム目を大きい順に並び替えて上位10行のみを取得して、range_hypo_date.csvからその日付のデータを抽出する
sort -t"," -k2 -n -r date_count.csv | head -n 10 | awk -F"," 'NR==FNR{a[$1];next}($1 in a){print $0}' - range_hypo_date.csv > top10.csv
gmt begin n-t pdf,png
gmt set GMT_VERBOSE=long
gmt basemap -JX24/10 -R2014-01-01/2024-08-31/0/500 -Bxa6Ofg+a45+l"Date" -Byafg+l"Count" -BWSne+t"2014-2024"
awk -F"," '{print $1}' range_hypo_date.csv | gmt histogram -T0.5 -W1
gmt histogram top10.csv -T0.5 -W1 -D -N
gmt inset begin -DjTR+w4/4+o0.1/0.1 -F+gwhite+p0.25 -M0.01
gmt set MAP_FRAME_TYPE = "plain"
gmt basemap -JM? -R138.5/139.5/35/36 -Bxy
gmt coast -Df -W0.25 -Glightbrown
gmt makecpt -Cseis -T0/100/1 -Z
awk -F"," '{if($11>=138.5 && $11<=139.5 && $9>=35 && $9<=36 && $13<=100) print $11,$9,$13,$15*0.01}' $input | sort -k4 -n | gmt plot -Scc -W0.01 -C
gmt plot -W0.5,blue -L << END
${west} ${north}
${east} ${north}
${east} ${south}
${west} ${south}
END
gmt inset end
gmt end
注意事項
- 震源の位置はAQUAシステム メカニズム解カタログに基づくため今後変わる可能性があります。
- 2024年の気象庁一元化震源データは今後見直される可能性があります。
謝辞
以下のデータを使用させていただきました。ここに記して感謝申し上げます。
- 震源データとして、気象庁の「地震カタログ」を一部改変したもの。
気象庁|震源データ
気象庁が提供するページです
- 発震機構解のデータとして、国立研究開発法人防災科学技術研究所のF-netデータ
防災科研機関リポジトリ
- HinetPy
Dongdong Tian, 2020 , HinetPy: A Python package to request and process seismic waveform data from Hi-net. , doi: 10.5281/zenodo.3885779
コメント