3216 words
16 minutes
用GMT绘制极坐标剖面
2020-01-19
Mapping
GMT

本文展示了在极坐标系下绘制地理剖面的示例。对于大尺度的地理剖面非常实用。本实例中GMT的难点主要在：

  • 极坐标投影与剖面长度的计算。
  • 在极坐标系下绘制410 km和600 km弧线。
  • 彩色图像绘制。
  • 速度矢量的求和与绘制。

1. 读取剖面经纬度#

剖面列表如下：

la1 lo1 la2 lo2 a
-- -- -- -- b
-- -- -- --  c

我们可以输入测线编号让脚本获取测线位置。

Terminal window
la1=`grep $line_name $lines|awk '{print $1}'`
lo1=`grep $line_name $lines|awk '{print $2}'`
la2=`grep $line_name $lines|awk '{print $3}'`
lo2=`grep $line_name $lines|awk '{print $4}'`

其中

  • line_name为测线编号
  • lines为剖面列表的路径名（相对或绝对路径）

2. 计算剖面长度#

由于每条剖面的长度不一样，我们通常指定剖面两个端点来确定剖面位置。在这个脚本的开始我们先通过project命令生成测线上的一系列网格点。

Terminal window
gmt project -C$lo1/$la1 -E$lo2/$la2 -G10 -Q > great_circle_points.xyp

其中

  • -C, -E分别表示测线起始点和终点坐标。
  • -G后的参数表示相邻两个点的距离（这里设置10 km）
  • -Q表示使用地图单位

生成的数据重定向到great_circle_points.xyp，后面还会用到:

-50     10      0
-49.9243897556  10.0504392132   10
-49.8487559173  10.1008612754   20
......
-10.1812240988  29.9436088093   4680
-10.08354655    29.974041736    4690
-10     30      4698.54844803

输出内容有三列，分别是每个点的经纬度和与第一个点的距离。最后一行的最后一列即为测线的总长（km）。

由于在极坐标投影时单位必须是度，我们要将上述测线单位转换为度，为了简化我们直接将公里数除以111。以下三行命令分别计算了剖面长度（设置绘图范围时使用）、剖面半长度（设置旋转角度是使用）和图幅大小（为了让不同长度剖面等比例）。由于在Shell语言中没有浮点型计算，所以我们用到了awk工具辅助

Terminal window
length=`awk '{print $NF/111}' great_circle_points.xyp|tail -1`
half_length=`awk 'BEGIN {print "'$length'"/2}'`
size=`awk 'BEGIN {print "'$length'"*2/3}'`

3. 根据计算的范围绘制极坐标底图#

绘制底图时最重要的两个参数是-R和-J，他们分别指定了绘制范围和投影方式。

绘制范围#

Terminal window
R=0/$length/5671/6371

其中

  • length表示上面我们计算的测线总长，单位是度。
  • 5671/6371表示上下两个边界的半径，5671意为下边界为700 km深

极坐标投影#

极坐标投影在GMT中用-JP表示。其标准格式为

Terminal window
-JP[a]<width>[/<theta0>][r][z]

其中

  • width表示图像的宽度
  • theta0表示旋转角度
  • 如果在P后面加a表示逆时针旋转，否则顺时针旋转
  • z 表示将 r 轴标记为深度而不是半径
  • r 表示将径向方向反转，此时r轴范围应在0到90之间

这里我们将投影方式设置为

Terminal window
J=Pa${size}c/${half_length}z
  • size为上面计算的图幅大小
  • half_length为旋转角度。

绘制底图#

有了-R和-J控制绘图范围和投影方式就可以绘制底图

Terminal window
gmt psbasemap -R$R -J$J -Bxa0f5 -Bya200f100 -BWSne -K -O -Yf7.7c >> $ps

需要说明的是

  • -Bx表示X轴的刻度，a0表示不画主刻度，f5表示没5度画一格副刻度。
  • -By表示Y轴刻度，每200 km画一格主刻度，每100 km画一格副刻度。
  • -BWSne表示只有左侧和下侧画刻度，上侧和右侧只画边框。
  • -O和-K在中文手册-K和-O选项的说明。
  • -Y表示相对与纸张的偏移量，也可以参见中文手册-X和-Y选项的说明.

4. 绘制彩色的网格文件#

下面需要在底图上绘制我们的数据内容（温度，速度等）

数据网格化#

虽然我们的数据已经是等间隔的网格数据，但是并不符合GMT绘图的格式标准。GMT要求网格文件是netCDF格式。所以这里使用xyz2grd命令只是将文本格式的数据转换成netCDF格式。当然如果想改变网格间隔或者输入数据是非等间隔的也可以surface、greenspline、nearneighbor 或 triangulate进行插值。

Terminal window
awk '{print $4,6371-$1,$8}' $sec_data | gmt xyz2grd -R -I0.1/10 -Gout.grd
  • 由于xyz2grd要求输入数据为三列（x,y,z 值）所以我们用awk命令从数据表里提取出绘图的三列
    • 第一列是横坐标，距离起始点的长度，单位是度。
    • 第二列是半径，所以用地球半径减去深度。
    • 第三列是数值。
  • -R选项表示继承上个命令设定的范围。
  • -I表示网格之间的间隔。0.1表示横向间隔为0.1度，10表示纵向间隔10 km。单位与数据表中的单位相同。
  • -G后面跟一个输出网格文件的文件名

生成一个色标文件#

这是所有脚本绘图工具的弱点，我们没有办法用GMT像ParaView一样拖动色块来自定义色标。所以有两种方式去生成色标。

  • 用GMT中的makecpt或grd2cpt命令生成等间隔色标文件。
  • 色标文件本质是一个文本文件，所以也可以手动修改文件中的数值来生成色标。

如果不需要刻意突出某个范围的数值，而且数值变化偏线性，那么使用makecpt命令即可

Terminal window
gmt makecpt -Crainbow -T200/1700/50 -Z > tmp.cpt
  • -C选项后跟输入的色标文件（CPT文件）。GMT提供了一些内置的色标（参见GMT中文手册内置CPT），在-C后加上这些内置色标名即可。色标也可以是一个自定义的色标文件，这时在-C后面加上文件的路径名即可。自定义色标可以在cpt-city中获取。
  • -T表示色标相对于数值的取值范围。200/1700/50意为从200至1700每50在输入色标中取一个值。
  • -Z表示使用连续色标。
  • 最后输出重定向至一个cpt文件

绘制网格文件#

有了网格文件和色标文件，就可以绘制图像了。

Terminal window
gmt grdimage out.grd -R -J -O -K -Ctmp.cpt >> $ps
  • out.grd是之前生成的网格文件名。
  • -C后加之前生成的色标文件名。

这时一副图的主体就画完了！

5. 在极坐标下画弧线#

这里我们想标注410 km和660 km的深度，在极坐标投影下是弧线，所以不能直接指定两个端点坐标，这样画出来是一条直线（在笛卡尔投影下这样做是没问题的）。我们的做法是利用之前生成的测线文件great_circle_points.xyp的点，把这些点连起来，来生成弧线。

Terminal window
awk '{print $3/111, 5961}' great_circle_points.xyp|gmt psxy -R -J -O -K -W0.7p,-+s  >> $ps
awk '{print $3/111, 5711}' great_circle_points.xyp|gmt psxy -R -J -O -K -W0.7p,-+s  >> $ps
  • awk用来生成每个点的坐标。$3/111是将每个点距起始点距离转化为度，5961和5711分别表示410 km和660 km对应的半径。
  • -W选项指定了线型。0.7p表示线的粗细，-表示虚线，+s表示两点之间以弧线连接（如果点数够密则影响不大，可以去掉）。

6. 绘制色标#

现在用刚才生成的色标文件绘制色标，用到psscale命令。

Terminal window
gmt psscale -DjML+w4c/0.4c+o`awk 'BEGIN {print "'$size'"+1}'`c/-0.5c -Ctmp.cpt -Bxafg+l"Temperature (K)" -R -J -O -K >> $ps
  • -D选项控制了色标的长宽和位置。
    • jML指定色标的锚点在中间左侧，参考绘制修饰物。
    • +w4c/0.4c表示色标的长为4厘米，宽为0.4厘米。
    • +o指定了色标相对于底图的位置。
      • `awk 'BEGIN {print "'$size'"+1}'`c表示向右偏移的距离为$size+1厘米。
      • -0.5c表示向下偏移1厘米
  • -C后面跟要绘制的色标文件
  • -Bx指定色标的刻度
    • afg表示自动绘制主刻度、副刻度和网格线。
    • +l"Temperature (K)"表示添加标注为Temperature (K)

7. 绘制另一个子图#

现在要在同一个文件中绘制另一个子图，所以要先画一个一样的底图

Terminal window
gmt psbasemap -R$R -J$J -Bxa0f5 -Bya200f100 -BWSne -K -O -Yf7.7c >> $ps

与上一个子图不同的是通过-Y对底图进行偏移，以免与之前的子图重叠，具体的偏移量只能手动慢慢调试。

8. 绘制速度矢量的模#

在数据表里速度有垂直和水平两个分量，我们要先求两个分量矢量和的模，然后绘制成彩色图像。

分别将两个分量网格化#

Terminal window
awk '{print $4,6371-$1,$5}' $sec_data | gmt xyz2grd -R -I0.1/10 -Gvr.grd
awk '{print $4,6371-$1,-1*$7}' $sec_data | gmt xyz2grd -R -I0.1/10 -Gvz.grd

用两个分量的网格文件求矢量和的模#

这里用GMT内置的grdmath命令，对两个分量的网格文件求矢量和的模。

∣vel∣=vx2+vz2|vel|=\sqrt{vx^2+vz^2}
Terminal window
gmt grdmath vr.grd vz.grd HYPOT = vel.grd
  • HYPOT表示运算符为上述公式
  • vel.grd为输出的网格文件名

下面用grdimage命令绘制网格文件即可。

Terminal window
gmt grdimage vel.grd -R -J -O -K -Cvel.cpt >> $ps

这里经过尝试makecpt和grd2cpt生成的色标都没办法满足我的要求。所以vel.cpt只能通过漫长的调试手动生成，其实没什么难度只是过程漫长。

9. 绘制矢量#

现在我们要将矢量方向以箭头的方式绘制在上一个图层上，同时箭头的长短表示大小。这里我们用到之前生成的两个分量的网格文件，以及grdvector命令。

Terminal window
gmt grdvector -R -J -O -K vr.grd vz.grd -W0.8p -Si9p -Ix10 -Q0.3c+b+jb+h1+p- -G0 >> $ps
  • -W指定矢量轮廓的粗细为0.8p
  • -S指定每个矢量数据单位在地图上为9p
  • -Q指定箭头的属性，参考绘制矢量/箭头
  • -Ix10表示每个每10个网格画一个箭头，以免箭头太密影响美观。
  • -G指定箭头填充颜色，0表示黑色

将箭头填充为渐变色#

为了突出数值大小，箭头填充的颜色可以按照其渐变，我们这里使用灰度渐变的色标。我们不知道具体数值的大小范围，所以先生成一个标准色标，在用grd2cpt命令自动生成。

  • 生成灰度值从0-200的连续色标（黑色至灰度值为200的灰色）

    Terminal window
    cat>tmp.cpt<<eof
    0   200  200  200  1  0   0   0
    eof

  • 按照数值大小自动生成色标文件

    Terminal window
    gmt grd2cpt vel.grd -Ctmp.cpt -Z > tmp.cpt1

  • 将上述grdvector中-G选项换成-C选项就可以把填充色改为渐变色。

    Terminal window
    gmt grdvector -R -J -O -K vr.grd vz.grd -W0.8p -Si9p -Ix10 -Q0.3c+b+jb+h1+p- -Ctmp.cpt1 >> $ps

只绘制部分箭头#

有时我们只想绘制部分数值较大的箭头，而省略数值较小的箭头。这是我们需要对数据表做一些筛选判断，并且重新网格化两个分量。

Terminal window
awk 'sqrt($5*$5+$7*$7)>0.9{print $4,6371-$1,$5}' $sec_data | gmt xyz2grd -R -I0.1/10 -Gvr.grd
awk 'sqrt($5*$5+$7*$7)>0.9{print $4,6371-$1,-1*$7}' $sec_data | gmt xyz2grd -R -I0.1/10 -Gvz.grd

与之前网格化命令不同的是我们在awk中多加了一行判断sqrt($5*$5+$7*$7)>0.9，来筛选出矢量和大小大于0.9的网格点。然后再用grdvector绘制矢量时则只绘制矢量和大小大于0.9的网格。

10. 绘制剖面的地形起伏#

我们还要在上幅图上方绘制地形起伏。由于地形的变化相对地球半径的尺度非常小，所以需要将地形起伏进行适当放大。

重新定义绘图范围#

Terminal window
R=0/$length/6371/6561

长度不变，高度进行放大。

Terminal window
gmt psbasemap -R$R -JPa${size}c/${half_length} -K -O -Bws -Yf12.3c >> $ps
  • 只标注-Bws可以去除四周的边框和刻度。

截取地形起伏数据#

这里使用grdtrack命令从地形起伏数据中截取该剖面的地形起伏。这条命令的输出有四列分别为经、纬度、距起始点的距离和高程。

Terminal window
gmt grdtrack great_circle_points.xyp -G$etopo | awk '{print $3/111,6371+($4/1000)*40}' > elev

在绘制地形时只需要横纵两个坐标，即距离和高程。我们对高程放大了40倍。

绘制地形起伏并填充颜色#

这时已经可以利用elev文件绘制高程了，但是要想填充颜色需要让图形闭合

Terminal window
tail -r great_circle_points.xyp | awk '{print $3/111,6371}' >> elev

这行命令中我们将测线上的点倒序追加到elev后面来形成闭合曲线。我们用下面的示意图来表述这一过程

最后我们用psxy命令绘制地形起伏

Terminal window
gmt psxy elev -R -J -W -Ggray -N -K -O>> $ps
  • -G表示填充颜色
  • -N表示超出$R范围的数据仍然绘制。

总结和问题#

  1. 同样是绘制剖面在极坐标下绘制比在笛卡尔坐标下绘制难得多。所以如果剖面不是很长就别用极坐标。
  2. 极坐标下没办法绘制Y轴的标签，通常情况下在-By选项后加+l"Depth (km)"即可，但是极坐标不行。需要的话只能用其他方式画了。
用GMT绘制极坐标剖面
https://fuwari.vercel.app/posts/gmt-section/
Author
Mijian Xu
Published at
2020-01-19
License
CC BY-NC-SA 4.0
用GMT绘制中国及周边构造分区与GPS速度场
判断地理坐标点是否在多边形内的Python实现
© 2025 Mijian Xu. All Rights Reserved. / RSS / Sitemap
Powered by Astro & Fuwari