绘制研究区域的地震射线空间分布图要求:

  1. 绘制所有地震到达所有观测台站的射线
  2. 绘制地震分布和台站分布

第一点在 gmt4 下不太好实现,因此采用 gmt5 进行绘制,脚本如下:

gmt set MAP_FRAME_TYPE plain
gmt set MAP_TICK_LENGTH -0.1c
gmt set FONT_ANNOT_PRIMARY=8p 
gmt set FORMAT_GEO_MAP ddd
gmt set MAP_FRAME_PEN 0.03c
set PS=seismicray.ps
set file=earthquake.txt
rem 提取地震经纬度
gawk "{if (NR>1) print $9,$8}" resultsn.dat > %file%
rem 设置图宽为10cm
gmt psbasemap -R115.5/120/38/41.5 -JM10c -Ba1 -BWS -K > %PS%
gmt psbasemap -R -J -B0 -Bne -K -O >> %PS%
rem 所有地震都与每个台站连线 gmt psxy -Fr
rem 依次提取文件中每个台站的经、纬度
for /f "tokens=1,2" %%i in (tstz.txt) do (
gmt psxy %file% -R -J -W0.3p -Fr%%i/%%j -O -K >> %PS%
)
gmt psxy %file% -R -J -Sc0.1c -Gblue -Wblack -K -O >> %PS% 
gmt psxy tstz.txt -R -J -St0.25c -Gred -Wblack -K -O >> %PS% 
echo H 10 39 图例 > legend.txt 
echo S 0.2c t 0.25c red black 0.5c @%%39%%台站@%%%% >> legend.txt 
echo S 0.2c c 0.15c blue black 0.5c @%%39%%地震@%%%% >> legend.txt
gmt pslegend legend.txt -R -J -F+p+gwhite -DjBL+w1.5c -L1.5 -O >> %PS%
gmt psconvert %PS% -Tg -A -P -C-sFONTPATH=C:\Windows\fonts -Z
del .gmt* gmt.* %file% legend.txt

效果见研究区域地震射线空间分布图。

研究区域地震射线空间分布图