数据下载: earthquakes.dat

本文展示了如何绘制选定测线的地震剖面图。 earthquakes.dat 为地震数据文件,包含三列数据,分别是经度、纬度、深度。

图(a)绘制了区域范围内的地震震中分布图,蓝色的线表示选定的测线;图(b)绘制了选定测线两侧0.2°范围内地震的震中分布图;图(c)绘制了选定地震的经度/深度剖面图;图(d)绘制了选定地震的纬度/深度剖面图。

本文主要使用了 project 模块:

  • -C100/22 :指定测线起点为(100°, 22°);
  • -E106/30 :指定测线终点为(106°, 30°);
  • -G0.1 :指定了测线上测点的间隔;
  • -W-0.2/0.2 :指定测线两侧各0.2°范围内的地震;
  • -F :该选项指定期望输出的数据列
    • -Fxy 表示输出数据的前两列,即经度和纬度;
    • -Fxz 表示输出数据的第一、三列,即经度和深度;
    • -Fyz 表示输出数据的后两列,即纬度和深度。

示例中 –W-G 选项使用的是和输入文件相同的单位(即度),如果使用不同单位(比如km),则需要添加-Q选项。

绘图效果如下:

地震剖面图

bash绘图脚本如下,也提供了等效的 bat脚本

#!/bin/bash
PS=seismic-profile.ps

gmt set MAP_FRAME_TYPE=plain
# === (a) 绘制区域震中分布图和选定的测线 ===
gmt psbasemap -R100/110/20/30 -JM3i  -Ba5f2.5 -BWeSn -P -Y7i -K > $PS
# 绘制所有地震位置
gmt psxy earthquakes.dat -R -J -Sc0.1i -Gred -W0.25p -O -K >> $PS
# 沿着测线每隔0.1度生成一个测点并绘制测线
gmt project -C100/22 -E106/30 -G0.1 | gmt psxy -R -J -W2p,blue -O -K >> $PS
echo '(a)' | gmt pstext -R -J -F+cBR+f12p,1 -Dj0.2c/0.2c -O -K >> $PS

# === (b) 绘制选定测线0.2度范围内地震的震中分布图 ===
gmt psbasemap -R100/110/20/30 -JM3i  -Ba5f2.5 -BWeSn -X4i -O -K >> $PS
# 筛选出测线周围0.2度范围内的地震并绘制
gmt project earthquakes.dat -C100/22 -E106/30 -Fxy -W-0.2/0.2 | gmt psxy -R -J -Sc0.1i -Gred -W0.25p -O -K >> $PS
# 沿着测线每隔0.1度生成一个测点并绘制测线
gmt project -C100/22 -E106/30 -G0.1 | gmt psxy -R -J -W2p,blue -O -K >> $PS
echo '(b)' | gmt pstext -R -J -F+cBR+f12p,1 -Dj0.2c/0.2c -O -K >> $PS

# === (c) 绘制选定测线0.2度范围内地震的经度/深度剖面图 ===
gmt psbasemap -R100/106/0/70 -JX7i/-1.0i  -Bxa1+lLongitude -Bya20f10+l"Depth(km)" -BWeSn --MAP_TICK_LENGTH_PRIMARY=-0.03i -X-4i -Y-2.0i -O -K >> $PS
# 筛选出测线周围0.2度范围内的地震并输出地震的经度和深度
gmt project earthquakes.dat -C100/22 -E106/30 -Fxz -W-0.2/0.2 | gmt psxy -R -J -Sc0.1i -Gred -W0.25p -O -K >> $PS
echo '(c)' | gmt pstext -R -J -F+cBR+f12p,1 -Dj0.2c/0.2c -O -K >> $PS

# === (d) 绘制选定测线0.2度范围内地震的纬度/深度剖面图 ===
gmt psbasemap -R22/30/0/70 -JX7i/-1.0i  -Bxa1+lLatitude -Bya20f10+l"Depth(km)" -BWeSn --MAP_TICK_LENGTH_PRIMARY=-0.03i -Y-2.0i -O -K >> $PS
# 筛选出测线周围0.2度范围内的地震并输出地震的纬度和深度
gmt project earthquakes.dat -C100/22 -E106/30 -Fyz -W-0.2/0.2 | gmt psxy -R -J -Sc0.1i -Gred -W0.25p -O -K >> $PS
echo '(d)' | gmt pstext -R -J -F+cBR+f12p,1 -Dj0.2c/0.2c -O -K >> $PS

gmt psxy -R -J -T -O >> $PS
gmt psconvert $PS -A -P -Tg
rm gmt.*