mapproject
- 官方文档:
- 简介:
地图投影的正反变换、数据转换以及大地测量计算
mapproject 可完成如下计算功能：
地图投影，读取经纬度数据，并使用指定的地图投影和比例尺计算投影坐标 (x,y)
地图投影逆变换，读取 (x,y), 计算经纬度；此功能可以将已知投影类型的地图数字化，再将其变换为地理坐标
坐标反算/大地主题反算，计算沿轨点列之间的距离与方位角，或到固定点的距离和方位角
计算点到最近线段的距离
基准转换，用于将某基准/椭球下的坐标转换为另一基准/椭球
不同纬度之间的转换，包括大地纬度、地心纬度等
大地坐标与空间直角坐标的转换
确定地图范围/绘图范围
语法
gmt mapproject
[ table ]
[ -Ab|f|o[lon0/lat0][+v] ]
[ -C[dx/dy][+m] ]
[ -Dc|i|p ]
[ -E[datum] ]
[ -F[e|f|k|M|n|u|c|i|p] ]
[ -G[lon0/lat0][+a][+i][+uunit][+v] ]
[ -I ]
[ -Jparameters ]
[ -Ltable[+p][+uunit] ]
[ -N[a|c|g|m] ]
[ -Q[d|e] ]
[ -Rregion ]
[ -S ]
[ -T[h]from[/to] ]
[ -V[level] ]
[ -W[b|B|e|E|g|h|j|m|M|n|o|O|r|R|w|s|x][+n[nx[/ny]]] ]
[ -Z[speed][+a][+i][+f][+tepoch] ]
[ -bibinary ]
[ -bobinary ]
[ -dnodata[+ccol] ]
[ -eregexp ]
[ -fflags ]
[ -ggaps ]
[ -hheaders ]
[ -iflags ]
[ -jflags ]
[ -oflags ]
[ -pflags ]
[ -qflags ]
[ -sflags ]
[ -:[i|o] ]
[ --PAR=value ]
输入数据
- table
一个或多个ASCII或二进制表数据。若不提供表数据，则会从标准输入中读取。
可选选项
- -A
- -Ab|f|o[lon0/lat0][+v]
计算沿着某轨迹，或到某固定点 (lon0/lat0) 的方位角。
-Af 计算每个数据点的（前向）方位角
-Ab 计算数据点到固定点的反方位角
-Ao 可获得 (-90/90) 方向，而非方位角 (0/360)
使用
-je 选项可在椭球体而非球面上计算方位角。 如果没有给出固定点，则从上一个点开始计算方位角或反方位角。 附加 +v 选项通过输入文件的第 3-4 列获取“可变”的第二个点 (lon0/lat0)。有关
-A如何影响输出记录，请参见 输出顺序 。 如果指定了-R和-J选项，则先进行投影，然后再计算投影后笛卡尔坐标下的方位角。
- -C
- -C[dx/dy][+m]
将投影坐标的中心设置为地图投影中心 [默认为左下角]。 也可以通过指定 dx/dy 来设置投影中心在东方向和北方向的偏移量。 投影结果将加上该偏移量，详见 示例 中的高斯投影和 UTM 投影。 如果和
-I选项 同时使用，则先减去该偏移量，然后再进行逆投影变换。 偏移量的单位是有效的绘图距离单位（请参阅 PROJ_LENGTH_UNIT ）， 但如果使用了-F选项，则偏移量单位为米。 对于 Mercator 投影，附加 +m 选项可将投影 y 坐标的原点设置为与标准纬线重合 [默认为赤道]。
- -D
- -Dc|i|p
临时覆盖 PROJ_LENGTH_UNIT，并改用 c (cm)，i (inch) 或 p (points) 为单位。 该选项不能与
-F一起使用。
- -E
- -E[datum]
将大地坐标 (lon,lat,height) 转换为以地球为中心的地球固定坐标 (ECEF) (x,y,z)，使用
-I选项 可实现逆向转换。追加基准的 ID datum （见-Qd 选项）或者通过 ellipsoid:dx,dy,dz 可指定基准，其中 ellipsoid 为椭球名称（见-Qe）。ellipsoid 除可用-Qe 列出 的椭球外，还可使用 a[,inv_f] 的形式通过给定椭球长半轴和扁率的倒数来指定椭球，其中 扁率的倒数默认为 0。不指定基准时，默认使用 WGS-84。
- -F
- -F[e|f|k|M|n|u|c|i|p]
强制使用 1:1 比例尺，即输出数据（或输入，参见
-I）为真实投影下的米数。 如果要指定其他单位，请附加所需的单位（请参阅 单位 ）。 在不使用-F且使用经典-J语法时， 输出（或输入，请参见-I）使用 PROJ_LENGTH_UNIT 指定的单位（或由-D选项指定）。 当在-J中使用 PROJ4 语法或 EPSG 时，无论是否使用-F，单位都会变为米。
- -G
- -G[lon0/lat0][+a][+i][+uunit][+v]
计算沿某轨迹的距离，或到某固定点 lon0/lat0 的距离。 如果没有给出固定点，默认将计算沿输入点定义的轨迹的累积距离。 如果给出了固定点，则默认计算增量距离。可以附加以下选项改变默认计算方式：
- -I
- -I
进行逆变换，例如：从 (x,y) 数据获取经纬度。
- -J
- -Jprojection
设置地图投影方式。 (参数详细介绍)
- -L
- -Ltable[+p][+uunit|c|C]
获取输入点到 ASCII 多段文件 table 中的线的最短距离。该距离和最近点的坐标将追加到输出的最后三列。
+u ：指定单位及距离计算方式，参见 单位 ，默认使用等面积半径的大圆路径。 其中 unit 可以设置为 c 表示使用输入坐标的笛卡尔距离， 设置为 C 表示使用投影坐标的笛卡尔距离。 C 需要设置
-R和-J，且所有输出坐标都将以投影坐标形式报告。+p ：报告线段 ID seg 和小数形式的最近点编号 pnr，不报告最近点的坐标。
-L选项对输出记录的信息顺序的影响，请参见 输出顺序 。 注 ：对于地理坐标来说，距离的计算使用球面模型，因此本选项不能与-je 选项同时使用。
- -N
- -N[a|c|g|m]
从大地纬度（使用当前的椭球，参阅 PROJ_ELLIPSOID）转换为其他辅助纬度，经度不受影响。
a - 转换为等面积纬度 (authalic latitudes)。
c - 转换为正形纬度 (conformal latitudes)。
g - 转换为地心纬度 (geocentric latitudes) [默认]。
m - 转换为子午线纬度 (meridional latitudes)。
与
-I选项同时使用，可实现反向变换。
- -Q
- -Q[d|e]
列出所有投影参数。 -Qd，仅列出基准， -Qe 仅列出椭球。
- -R
- -Rxmin/xmax/ymin/ymax[+r][+uunit]
指定数据范围。 (参数详细介绍)
- -S
- -S
不转换区域之外的点
- -T
- -T[h]from[/to]
使用标准的 Molodensky 转换将基准从 from 转换到 to。 如果使用 -Th ，则表明除坐标转换外，还进行大地高转换，若输入数据包含 3 列， 则最后一列为高度，若不含第三列，则默认高度为 0，即点位于椭球面上。 基准可通过基准 ID（参见 -Qd）指定或通过 ellipsoid:dx,dy,dz 指定；其中 ellipsoid 为椭球，dx，dy，dz为椭球中心和地心在三个方向的差异。 其中椭球可以是椭球ID（参见 -Qe），也可以指定为 a[,inv_f]，其中 a 是椭球长半轴，inv_f 是椭球扁率的倒数，如果省略，则为 0）。若不指定转换后的 基准，则假定为 WGS84。
-T和-R以及-J同时使用时，会先做基准转换，然后 投影。用户在使用时需注意 PROJ_ELLIPSOID 设置的椭球是否正确。
- -V
- -V[level]
设置 verbose 等级 [w]。 (参数详细介绍)
- -W
- -W[b|B|e|E|g|h|j|m|M|n|o|O|r|R|s|w|x][+n[nx[/ny]]]
报告以投影单位或地理单位表示的各种绘图尺寸或地图区域，不读取任何输入文件。 在没有参数的情况下，将报告地图的宽度和高度。可以通过
-D更改报告绘图尺寸所采用的单位。 可以使用以下附加选项：b - 获取以经纬度表示的边界框。
B - 同上，但以 -Rw/e/s/n 字符串格式作为尾随文本返回结果。
E - 同上，但以 -Rw/e/s/n 字符串格式作为尾随文本返回。
g - 输出附加地图点 lon/lat 的绘图坐标。
h - 仅输出地图的高度。
j - 输出参考点 code 的地图坐标，参考点 code 的定义请见 锚点 。
n - 同上，但参考点 rx/ry 以 0-1 范围内的归一化位置给出。
o - 如果通过
-Rxmin/xmax/ymin/ymax+uunit 设置了斜向定义域， 则返回以度为单位的对角顶点坐标（顺序为 llx urx lly ury）。O - 同上，但以等效的 -Rw/e/s/n 字符串格式作为尾随文本返回。
m - 改为获取投影绘图坐标下的矩形区域。
M - 同上，但以 -Rw/e/s/n 字符串格式作为尾随文本返回。
R - 同上，但以 -Rw/e/s/n 字符串格式作为尾随文本返回结果。
s - 以 1:xxxxxx 的形式输出地图比例尺。
w - 以当前绘图单位输出地图的宽度。
x - 输出特定绘图参考点 px/py 的地图坐标。
其中， e 或 r 选项后追加 +n 可设置斜向区域闭合多边形每条边上的点数 nx/ny [默认为 100]。 输出所有点构成多边形，而不是只输出四个角点。
倾斜投影区域（红框范围）和常规区域（黑框范围）。左图中红色区域为倾斜的，可以使用
-Wr|R 获取对应的整个范围的矩形经纬度范围。右图中黑色框是立体投影的范围，红色框为其对应的矩形范围，使用-We|E 可获取红色框范围
- -Z
- -Z[speed][+a][+i][+f][+tepoch]
指定行进速度 speed ，计算沿
-G指定轨迹的行进时间，如果不指定速度，则会从 第 3 列读取。速度为单位时间（单位由 TIME_UNIT [m/s] 指定）内通过的距离。+i ：输出每两个相邻点之间的行进时间。
+a ：输出累积的行进时间。同时使用则获取两种时间信息。
+f ：将累积的行进时间转换为 ISO 8601 规定格式。 时间中秒的小数位数，默认使用 FORMAT_CLOCK_OUT 设置的值。
+t：设置一个时间历元 epoch ，输出相邻点的绝对时间。 由于上述操作需要点之间的距离，因此需要在
-G选项中使用 +i 。
本选项对输出记录的信息顺序的影响，请参见 输出顺序 。
- -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]zgap[+n|p]
确定数据或线段的间断。 (参数详细介绍)
- -h
- -h[i|o][n][+c][+d][+msegheader][+rremark][+ttitle]
在读/写数据时跳过文件开头的若干个记录。 (参数详细介绍)
- -i
- -icols[+l][+sscale][+ooffset][,…][,t[word]]
对输入的数据进行列选择以及简单的代数运算。 (参数详细介绍)
- -j
- -je|f|g
设置球面距离的计算方式。 (参数详细介绍)
- -o
- -ocols[,…][,t[word]]
对输出的数据进行列选择以及简单的代数运算。 (参数详细介绍)
- -p
- -p[x|y|z]azim[/elev[/zlevel]][+wlon0/lat0[/z0]][+vx0/y0]
设置3D透视视角。 (参数详细介绍)
- -q
- -q[i|o][~]rows[+ccol][+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 指定更多的有效数字。
存在的限制
使用 -R 指定的矩形区域在经过投影后，范围变成了非矩形。除非使用 -C, 否则投影
后的网格的最左边的坐标将为 0，最下面的坐标也为 0。因此，在对地图数字化前，应该
使用 mapproject 对地图的几个角点的坐标极值进行一定的计算，以确定适合的比例尺
和地图左下角的坐标。在数字化时，使用这些值以保证正反算结果正确；或者，可以使用
类 Unix 系统中的 awk 命令将比例尺和偏移添加到坐标中。
对于某些投影，即使设置了参考椭球为某个椭球，仍然会使用球解法。一般只有 -R 设置
范围超出椭球展开的有效范围时会出现这种情况。包括下面几种投影：
Lambert Conformal Comic (-JN) 和 Albers Equal-Area (-JB) 将会在地图比例尺 超过 1.0E7 时使用球解法。
Transverse Mercator (-JT) 和 UTM (-JU) 当
-R指定的范围的东西边界 超过中央子午线 10 度以上时，使用球解法；Cassini (-JC) 与上面相同，但超过 4 度就使用球解法；
所有的圆锥投影，在经度范围小于 90 度时，使用球解法。
椭球和球
上述计算中，如果用户选择的椭球不为球，且基于椭球的公式已经实现，则 GMT 使用相应的 椭球公式，椭球选择见 PROJ_ELLIPSOID。用户需要注意的一些隐患如下：
对于某些投影，如横轴墨卡托投影，Albers，以及 Lambert 投影，当计算区域较小时， 使用椭球公式，当计算区域较大时，则使用球面公式，同时将纬度替换为对应的辅助 纬度；具体情况在 存在的限制 章节中已经介绍；
当用户想检查一些历史数据时，可能会发现 GMT 给出的计算结果和历史数据有轻微的不同， 其中一个原因是，早期的 GMT 计算使用的有效数字位数比较少，例如，Clarke 1866 椭球 的扁率为 1/294.98，GMT 对扁率做平方后使用，值为 0.00676862818，但是早期的版本为 0.00676866。这会导致十几个厘米的差异。如果用户需要复现以前的计算结果，则需要自己 指定椭球扁率。另一个原因可能是，早期的数据使用了和现在不同的基准面，这可能会导致 几十甚至几百米的差异。
PROJ_SCALE_FACTOR 对于某些投影来说有特定的默认值，用户若想实现某些设置需 手动设置该参数，见 示例 中的高斯投影。
输出顺序
选项 -A、-G、-L 和 -Z 生成的数据一律按这几个选项的字母顺序排序。
因此，这些选项在命令中出现的顺序无关紧要。当然，输出顺序也可以用 -o 选项调整。
示例
Mercator 投影
将经纬度坐标使用 Mercator 投影转换为平面坐标，1 度对应画布上 0.5 cm:
gmt mapproject @waypoints.txt -R-180/180/-72/72 -Jm0.5c -:
UTM 投影与高斯投影及其逆变换
将经纬度坐标使用 UTM 投影转换为平面坐标，下面经纬度坐标位于 UTM 的 51 投影带， 投影中心的坐标为 (0,0)，投影后坐标单位为米:
$ echo 121 32 | gmt mapproject -Ju51/1:1 -C -F
311072.361931 3542183.49115
GMT 中并无高斯投影（Gauss-Kruger），但其与 UTM 投影都是横轴墨卡托投影的变种，因此 可用 UTM 投影经过一定的设置近似实现，其中 UTM 和高斯投影的主要不同在于带号和中央经线投影的比例因子，高斯投影为切椭圆柱投影， 因此其比例因子为 1， UTM 为割椭圆柱投影，其比例因子为 0.9996，因此上述经纬度的 6 度带高斯投影为:
$ echo 121 32 | gmt mapproject -Ju51/1:1 -C -F --PROJ_SCALE_FACTOR=1
310996.760635 3543600.93152
此外，可以使用横轴墨卡托投用实现 UTM 和高斯投影:
$ echo 121 32 | gmt mapproject -Jt123/0/1:1 -C500000/0 -R120/126/30/40 -F --PROJ_SCALE_FACTOR=0.9996
311072.361931 3542183.49115
$ echo 121 32 | gmt mapproject -Jt123/0/1:1 -C500000/0 -R120/126/30/40 -F --PROJ_SCALE_FACTOR=1
310996.760635 3543600.93152
需要注意的是，上述的代码均基于默认椭球（WGS-84），若需其他椭球，需手动设置，也可 通过 PROJ 的形式指定 EPSG 代码实现，详见 -J 选项。
将平面坐标使用 UTM 投影逆变换转换为经纬度坐标，平面坐标单位为 m，对应的带号为 51
$ echo 311072.4 3542183.5 | gmt mapproject -Ju51/1:1 -C -F -I
121.000000401 32.0000000858
# 将高斯坐标转换为经纬度
$ echo 310996.760635 3543600.93152 | gmt mapproject -Ju51/1:1 -C -F -I --PROJ_SCALE_FACTOR=1
121 31.9999999997
不同纬度的转换
# 将大地纬度转换为地心纬度
$ echo 0 45 | gmt mapproject -Ng
0 44.807576784
# 将地心纬度转换为大地纬度
$ echo 0 44.807576784 | gmt mapproject -Ng -I
0 45
笛卡尔坐标反算与大地主题反算
假定一平面坐标（0,0），计算其与坐标 (1,1) 之间的距离和方位角:
$ echo 0 0 | gmt mapproject -Ab1/1 -G1/1+uc
0 0 45 1.41421356237
输出信息分别为横纵坐标，方位角以及距离；设置输出单位为 km，距离计算方式 (这里使用椭球面距离)，可计算两个经纬度坐标之间的方位角和距离:
$ echo 0 0 | gmt mapproject -Ab1/1 -G1/1+uk -jg
0 0 44.9956364553 156.898594962
基准转换
将点 (0,0) 从香港使用的 1924 年国际椭球转换为 WGS-84 椭球:
$ echo 0 0 | gmt mapproject -T
-0.00243433862099 -0.00170923972985
对于带有大地高的坐标 (0,0,0)，将其从 CGCS2000 椭球转换为 WGS-84 椭球，由于 GMT 中没有 CGCS2000 的 ID，因此需通过参数指定:
$ echo 0 0 0 | gmt mapproject -Th6378137,298.257222101:0,0,0/219
0 0 0
由于 CSCS2000 椭球与 WGS-84 椭球极为接近，因此其转换后的经纬度和大地高均未有 明显变化
大地坐标和经纬度的转换
$ echo 10 10 10 | gmt mapproject -E219
6186446.76449 1090837.4793 1100250.28422
$ echo 6186446.76449 1090837.4793 1100250.28422 | gmt mapproject -E219 -I
10 10 9.99999782909
点到线段的距离
计算 quakes.txt 中的点到 coastline.txt 中每段线的最短距离，
距离单位为 km
gmt mapproject quakes.txt -Lcoastline.txt+uk > quake_dist.txt
计算沿轨累计距离和时间
假定船速度为 12 Knots，计算沿轨迹累计距离和航行时间:
gmt mapproject track.txt -G+un+a+i -Z12+a --TIME_UNIT=h > elapsed_time.txt
:option:`-W` 选项的使用
确定横轴墨卡托投影地图的中点的地图坐标:
gmt mapproject -R-80/-70/20/40 -Jt-75/1:500000 -WjCM > mid_point.txt
确定倾斜墨卡托投影覆盖的矩形区域:
gmt mapproject -R270/20/305/25+r -JOc280/25.5/22/69/2c -WR
确定矩形绘图范围对应的实际地图坐标范围:
gmt mapproject -R-2800/2400/-570/630+uk -Joc190/25/266/68/1:1 -WO
获取矩形绘图范围对应的实际地图范围的多边形，坐标为地理坐标:
gmt mapproject -R-2800/2400/-570/630+uk -Joc190/25/266/68/1:1 -Wr+n > polygon.txt
确定立体投影确定的非矩形区域对应的矩形范围:
gmt mapproject -JS36/90/30c -R-15/60/68/90 -WE
参考文献
Bomford, G., 1952, Geodesy, Oxford U. Press.
Snyder, J. P., 1987, Map Projections - A Working Manual, U.S. Geological Survey Prof. Paper 1395.
Vanicek, P. and Krakiwsky, E, 1982, Geodesy - The Concepts, North-Holland Publ., ISBN: 0 444 86149 1.