✍️ 田冬冬 , 陈箫翰 , 王亮 , 姚家园 , 朱邓达 • 📅 2026-01-09

# project

官方文档:

  project

简介:

  生成测线、将数据点投影到测线上

该模块具有三个主要功能：

1. 生成测线

2. 指定测线和数据点，得到该点在测线坐标系下的坐标，即下文中提及的 p 和 q

3. 指定测线和数据点，得到该点在测线上的投影点的坐标，即下文中提及的 r 和 s

以上三个功能均要求用户首先定义测线，测线可以用如下三种方式中的任意一种来定义：

1. 用 `-C` 选项定义测线的起点，并用 `-A` 定义测线方位角。

2. 用 `-C` 选项定义测线的起点，并用 `-E` 选项定义测线的结束点。

3. 用 `-C` 选项定义测线的中心，并用 `-T` 选项定义旋转极点的位置。当使用 `-N` 设置笛卡尔坐标变换时，不允许使用此方式。

project 可以从输入数据中读取 (\(x, y [,\ z]\)) 坐标数据，并输出 (\(x, y, z, p, q, r, s\)) 的任意组合。其中：

- x 和 y 是数据在原坐标系下的坐标

- z 是输入数据中的其余所有列

- p 和 q 是数据点 (x,y) 在测线坐标系下的坐标

- r 和 s 是数据点 (x,y) 在测线上的投影点在原坐标系下的坐标

可以使用 `-F` 选项设置要输出哪些变量。

或者使用 `-G` 选项生成测线，输出测线上的点的坐标。 project 可以沿剖面以等间距 dist 输出 (\(r, s, p\)) 。在这种情况下，不需要读取任何输入数据。

为了将数据沿大圆路径进行球面投影，会创建一个斜坐标系，其赤道沿该路径，且零经度线通过 cx/cy。 此时，斜经度 (\(p\)) 对应于沿大圆路径距离 cx/cy 的距离，斜纬度 (\(q\)) 对应于垂直于大圆路径的距离。 当向 (\(p\)) 增加的方向（即由 `-A` azimuth 设定的方向）移动时，正 (\(q\)) 方向位于左侧。 如果通过 `-T` 指定了极点，则正 (\(q\)) 方向指向极点。

若要指定斜投影，请使用 `-T` 选项设置极点。此时投影的赤道已经确定，通过 `-C` 选项来定位 \(p = 0\) 的经线。 中心点 cx/cy 将被视为 \(p = 0\) 经线所通过的一个点。如果不需要选择特定的点，请使用南极点（cx = 0, cy = -90）。

可以使用 `-L` 和 `-W` 选项对数据进行过滤。 如果使用 `-W`，则仅输出 \(w_{min} < q < w_{max}\) 的点。 如果设置了 `-L`，则仅输出 \(l_{min} < p < l_{max}\) 的点。

下面详细解释一下这些变量的物理意义。

`Source Code`

project 示意图1

图中的红点就是给出的点(x,y)。绿色粗线为测线。 测线上有 3 个点，蓝色和绿色两个点为测线的起点和终点。 以测线起点为原点，以测线为P轴，在测线起点按右手螺旋法则做垂直于P轴的Q轴，构成测线坐标系。 点(x,y)在测线坐标系的坐标即为(p,q)。 紫色点为点(x,y)在测线上的投影点，坐标为(r,s)。

`Source Code`

project 示意图2

输入数据（红色圆圈）在原始的 \(x\text{-}y\) （或经纬度 \(lon\text{-}lat\) ）坐标系中给出， 并被投影到由起点 (C) 以及终点 (E) 或方位角 (\(\alpha\)) 定义的 \(p\text{-}q\) 测线坐标系中。 蓝色点的投影坐标为 \((p, 0)\)，在原始坐标系中为 \((r, s)\)。 选项 `-L` （限制 \(p\) 的范围）和 `-W` （限制 \(q\) 的范围）可用于排除指定范围之外的数据（浅灰色区域）。

对于特定的大圆和大地距离计算，或反方位角/方位角的计算，最好使用 mapproject，因为 project 仅严格遵循球面计算。

## 语法

gmt project [ table ] `-C` cx/cy [ `-A` azimuth ] [ `-E` bx/by ] [ `-F` flags ] [ `-G` dist[unit][/colat][+c][+h][+n] ] [ `-L`[w|lmin/lmax] ] [ `-N` ] [ `-Q` ] [ `-S` ] [ `-T` px/py ] [ `-V`[level] ] [ `-W` wmin/wmax ] [ `-Z` major[unit][/minor/azimuth][+e] ] [ `-bi` binary ] [ `-bo` binary ] [ `-d` nodata[+c col] ] [ `-e` regexp ] [ `-f` flags ] [ `-g` gaps ] [ `-h` headers ] [ `-i` flags ] [ `-o` flags ] [ `-q` flags ] [ `-s` flags ] [ `-:`[i|o] ] [ --PAR=value ]

## 输入数据

table

  一个或多个ASCII或二进制表数据。若不提供表数据，则会从标准输入中读取。

## 必须选项

-C

-C cx/cy

  当与 `-A` 或 `-E` 结合使用时，设置测线的起点 cx/cy。 当与 `-T` 结合使用时，设置斜零经线 (\(p = 0\)) 所通过的点 cx/cy 的坐标。cx/cy 并不要求必须位于距离极点 90 度的位置。

## 可选选项

-A

-A azimuth

  定义1中用于指定测线的方位角。参数 azimuth 定义为正北开始的顺时针旋转角。

-E

-E bx/by

  定义2中用于指定测线的终点

-F

-F flags

  指定输出格式 flags 。使用 xypqrsz 的任意组合按任意顺序指定所需的输出格式 [默认输出格式为 xypqrsz]。 如果输出格式为 ASCII，则 z 还包含任何尾随文本（无论 z 在 flags 中的顺序如何，文本均放置在记录的末尾）。必须使用小写字母且不要在字母之间添加空格。 注意：如果使用了 `-G`，则输出顺序固定为 rsp，并且此时 `-F` 选项无效。

-G

-G dist[unit][/colat][+c][+h][+n]

  该选项用于生成测线，此时不需要输入文件 table 。 输出数据有三列：经度、纬度以及当前点离测线起点的距离，即 (r, s, p) ，测点的距离间隔为 dist 。

  - 不使用 `-Q` 选项时， unit 的设置无效。 (r, s, p) 和 dist 的单位强制为弧度。

  - 而使用 `-Q` 选项时， r、s 的单位是弧度， p 和 dist 的单位默认为km。可以设置 unit 为其他长度单位，例如 e 表示为m。

  可以使用以下附加选项：

  - 附加 /colat 使用小圆路径 [默认余纬为 90，即按照大圆路径生成测线]。注意，当使用 `-C` 和 `-E` 指定测线的两个端点时，起点和终点之间的距离不能超过 \(2|colat|\)。

  - 在使用 `-T` 时附加 +c，以计算使小圆经过中心点 cx/cy 的余纬。

  - 在使用 `-T` 时附加 +h，以在段头中报告极点的位置 [默认无文件头]。

  - 附加 +n 则 dist 不再表示为距离间隔，而是测点的总数。需要 `-C` 和 `-E` 或 `-Z`，以便计算长度。

-L

-L[w|lmin/lmax]

  只输出坐标满足 \(l_{min} < p < l_{max}\) 的数据，单位规定见 `-Q` 选项。 若使用了 `-E` 选项，则可以简单使用 -Lw 表示只选取 p 在测线起点和终点之间的数据。

-N

-N

  指定展平地球，在平面内使用笛卡尔坐标变换。不使用本选项默认使用球面三角。 注意： 即使使用了本选项，方位角 azimuth 也是定义为从正北开始的顺时针旋转角， 不是 通常笛卡尔坐标系中的 theta（从 \(x\) 轴开始的逆时针旋转角）。 即：\(azimuth = 90 - theta\)。

-Q

-Q

  指定原坐标系下 x、y、r、s 的单位是弧度， 测线坐标系下 p、 q、 dist、 lmin、 lmax、 wmin、 wmax 的单位是km。 若不使用本选项，则以上所有量都默认为相同的单位（球面三角中默认全部为弧度）。

-S

-S

  将输出按照 p 增序排列。在原数据坐标随机分布时较为有用。

-T

-T px/py

  将投影旋转极点的位置设置为 \(px/py\)。

-V

-V[level]

  设置 verbose 等级 [w]。 (参数详细介绍)

-W

-W wmin/wmax

  只输出坐标满足 \(w_{min} < q < w_{max}\) 的数据，单位规定见 `-Q` 选项。

-Z

-Z major[unit][/minor/azimuth][+e]

  创建长轴为 \(major\)、短轴为 \(minor\) 的椭圆坐标。轴长单位默认为 km，长轴方位角 \(azimuth\) 的单位为度。 此功能需与 `-C` 和 `-G` 结合使用。

  注意：对于笛卡尔椭圆（使用 `-N` ），则应该给出相对于水平方向逆时针旋转的 \(direction\) （方向），而非 \(azimuth\) （方位角）。 地理椭圆的 \(major\) 长轴可以通过附加单位后缀来指定所需的单位 [默认为 km]。短轴 \(minor\) 和 `-G` 的 dist 也使用相同的单位。 对于退化为正圆的椭圆，可以仅提供单个 \(diameter\) （直径）。

  附加 +e 调整 `-G` 设置的增量 dist ，使椭圆具有相等的距离增量 [默认使用给定的增量并闭合椭圆]。

-bi

-bi[ncols][type][w][+l|b]

  控制二进制文件的输入格式。 (参数详细介绍)

-bo

-bo[ncols][type][w][+l|b]

  控制二进制文件的输出格式。 (参数详细介绍)

-d

-d[i|o]nodata

  将某些特定值当作 NaN。 (参数详细介绍)

-e

-e[~]"pattern" | -e[~]/regexp/[i]

  筛选或剔除匹配指定模式的数据记录。 (参数详细介绍)

-f

-f[i|o]colinfo

  显式指定当前输入或输出数据中每一列的数据类型。 (参数详细介绍)

-g

-g[a]x|y|d|X|Y|D|[col]z gap[+n|p]

  确定数据或线段的间断。 (参数详细介绍)

-h

-h[i|o][n][+c][+d][+m segheader][+r remark][+t title]

  在读/写数据时跳过文件开头的若干个记录。 (参数详细介绍)

-i

-i cols[+l][+s scale][+o offset][,...][,t[word]]

  对输入的数据进行列选择以及简单的代数运算。 (参数详细介绍)

-o

-o cols[,...][,t[word]]

  对输出的数据进行列选择以及简单的代数运算。 (参数详细介绍)

-q

-q[i|o][~]rows[+c col][+a|f|s]

  对输入或输出的行进行筛选，该选项在一定程度上可以代替 gawk 的某些功能。 (参数详细介绍)

-s

-s[cols][+a|+r]

  设置 NaN 记录的处理方式。 (参数详细介绍)

-:

-:[i|o]

  交换输入或输出数据的前两列。 (参数详细介绍)

-^ 或 -

  显示简短的帮助信息，包括模块简介和基本语法信息（Windows下只能使用 -）

-+ 或 +

  显示帮助信息，包括模块简介、基本语法以及模块特有选项的说明

-? 或无参数

  显示完整的帮助信息，包括模块简介、基本语法以及所有选项的说明

--PAR=value

  临时修改GMT参数的值，可重复多次使用。参数列表见 配置参数

## ASCII 格式精度

ASCII 格式输出数据通过 gmt.conf 配置文件控制。控制经纬度格式的参数为 FORMAT_GEO_OUT ；控制绝对时间的的参数包括 FORMAT_DATE_OUT 和 FORMAT_CLOCK_OUT ；普通浮点数通过参数 FORMAT_FLOAT_OUT 控制。上述格式控制可能会导致精度损失，这会在下游计算中导致一些问题。 如果用户需要保证数据精度，则应考虑将数据写为二进制文件，或者使用 FORMAT_FLOAT_OUT 指定更多的有效数字。

## 示例

把数据文件ship_03.txt（格式为lon,lat,depth）投影到由两点（330,-18）、（53,21）定义的大圆路径上，并进行排序，只输出距离 p 和深度。 p 的单位是千米:

    gmt project ship_03.txt -C330/-18 -T53/21 -S -Fpz -Q > ship_proj.txt

指定测线的起点和终点，在测线上每隔10千米生成一个点。距离测线起点的距离单位为千米:

    gmt project -C-50/10 -E-10/30 -G10 -Q > great_circle_points.xyp

指定测线的起点和终点，沿着colatitude=60的小圆上，每隔10千米生成一个点:

    gmt project -C-50/10 -E-10/30 -G10/60 -Q > small_circle_points.xyp

利用 `-F` 选项指定输出坐标 r, s 来得到某点在某测线上的投影点:

    echo 102 30 | gmt project -C103/31 -A225 -L0/500 -Frs -Q

已知某点，根据方位角和大圆距离计算另一点。已知一点(120, 25)，根据方位角 45 度，大圆距离 123 千米的点位置

    gmt project -C120/25 -A45 -L0/123 -G123 -Q | tail -1

根据地震目录数据，将地震事件投影到深度剖面并绘制：

    #!/usr/bin/env bash
    gmt begin ex
        a=122/22
        ap=148/46
        # 从GMT远程服务器下载示例地震目录文件
        gmt which -Gl @quakes_2018.txt

        gmt basemap -JM10c -R116/149/20/48 -Baf
        gmt grdimage @earth_relief_04m_p -Cgeo
        # 绘制 a-ap 测线。生成测线时以0.1度为间隔输出坐标点。
        gmt project -C${a} -E${ap} -G0.1 | gmt plot -W1p,cyan
        # 示例文件前三列分别为经度、纬度、深度
        # 根据深度绘制不同颜色的圆点
        gmt makecpt -Cmagma -T0/600/1 -I
        gmt colorbar -Bxaf+lDepth -C
        gmt plot quakes_2018.txt -Sc0.2c -C -W0.1

        # 绘制地表
        gmt basemap -JP10c+a+t18 -R0/36/5711/6371 -Byg6371 -BS -Y11c
        # 绘制 410 界面
        gmt basemap -Byg6371+5961 -BS
        # 绘制 660 界面
        gmt basemap -Byg6371+5711 -BS
        # 把数据投影到 a-ap 剖面，限制为剖面两侧30度范围的数据
        gmt project quakes_2018.txt -C${a} -E${ap} -Fpz -Lw -W-30/30 | gawk '{print $1,6371-$2,$2}' | gmt plot -Sc0.15c -C -W0.1p -N
    gmt end show

## 相关模块

fitcircle, vector, grdtrack, mapproject, grdproject