mapproject

贡献者:

@cugliming周茂

最近更新日期:

2022-12-05


官方文档:

mapproject

简介:

地图投影的正反变换、数据转换以及大地测量计算

mapproject 可完成如下计算功能:

  1. 地图投影,读取经纬度数据,并使用指定的地图投影和比例尺计算投影坐标 (x,y)

  2. 地图投影逆变换,读取 (x,y), 计算经纬度;此功能可以将已知投影类型的地图数字化,再将其变换为地理坐标

  3. 坐标反算/大地主题反算,计算沿轨点列之间的距离与方位角,或到固定点的距离和方位角

  4. 计算点到最近线段的距离

  5. 基准转换,用于将某基准/椭球下的坐标转换为另一基准/椭球

  6. 不同纬度之间的转换,包括大地纬度、地心纬度等

  7. 大地坐标与空间直角坐标的转换

  8. 确定地图范围/绘图范围

语法

gmt mapproject [ table ] -Jparameters -Rregion [ -Ab|B|f|F|o|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 ] [ -Ltable[+p][+uunit] ] [ -N[a|c|g|m] ] [ -Q[d|e] ] [ -S ] [ -T[h]from[/to] ] [ -V[level] ] [ -W[e|E|g|h|j|n|o|O|r|R|w|x][+n[nx[/ny]]] ] [ -Z[speed][+a][+i][+f][+tepoch] ] [ -bbinary ] [ -dnodata[+ccol] ] [ -eregexp ] [ -fflags ] [ -ggaps ] [ -hheaders ] [ -iflags ] [ -jflags ] [ -oflags ] [ -pflags ] [ -qflags ] [ -sflags ] [ -:[i|o] ] [ --PAR=value ]

必须选项

table

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

-Jprojection (more …)

设置地图投影方式

-Rxmin/xmax/ymin/ymax[+r][+uunit] (more …)

指定数据范围

UTM 投影比较特殊:如果使用 -C 选项且未指定 -R 范围,则将区域设置为给定的 UTM 区域,以保留完整的椭球解

可选选项

-Ab|B|f|F|o|O[lon0/lat0][+v]

计算沿着某轨迹,或到某固定点 (lon0/lat0) 的方位角。

  • -Af 计算每个数据点的(前向)方位角

  • -Ab 计算数据点到固定点的反方位角

  • -Ao 可获得 (-90/90) 方向,而非方位角 (0/360)

大写的 FBO,表示从大地纬度转换为地心纬度,然后再计算方位角(假设当前的椭球体不是球体)。如果没有给出固定点,则从上一个点开始计算方位角(或反方位角)。或者,通过输入文件中的 3-4 列,利用 +v 选项以获得可变的第二个点 (lon0/lat0)。有关 -A 如何影响输出记录,请参见 输出顺序 。如果指定了 -R-J 选项,则先进行投影,然后再计算投影后笛卡尔坐标下的方位角。

-C[dx/dy][+m]

将投影坐标的中心设置为地图投影中心 [默认为左下角]。也可以通过指定 dx/dy 来设置投影中心在东方向和北方向的偏移量,投影结果将加上该偏移量,详见 示例 中的高斯投影和 UTM 投影;如果和 -I 选项同时使用,则先减去该偏移量,然后再进行逆投影变换。偏移量的单位是有效的绘图距离单位(请参阅 PROJ_LENGTH_UNIT ),但如果指定了 -F 选项,则偏移量单位为米。对于 Mercator 投影,追加 +m 选项可设置投影坐标的原点与标准纬线重合 [默认为赤道]。

-Dc|i|p

临时覆盖 PROJ_LENGTH_UNIT,并改用 c (cm),i (inch) 或 p (points) 为单位。该选项不能与 -F 一起使用。

-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[e|f|k|M|n|u|c|i|p]

强制使用 1:1 比例尺,即输出数据(或输入,参见 -I )为真实投影下的米数。如果要指定其他单位,请附加所需的单位(请参阅 单位 )。如果没有 -F ,则输出(或输入,请参见 -I )使用 PROJ_LENGTH_UNIT 指定的单位(参见 -D )。

-G[lon0/lat0][+a][+i][+uunit][+v]

计算沿某轨迹的距离,或到某固定点 lon0/lat0 的距离。在距离后可通过 +u 指定单位及距离计算方式,参见 单位 ,默认使用大圆距离。例如,其中 unit 也可使用如下字符,c 表示使用输入坐标的笛卡尔距离,C 表示使用投影坐标的笛卡尔距离,C 需要设置 -R-J 。如果没有给出固定点,将计算沿输入点定义的轨迹的累积距离,指定 +a 也表示计算累积距离; 附加 +i 可获取连续点之间的增量距离。或者同时使用 +a+i,以同时获得累积距离及增量距离。若指定 +v,则通过输入文件中的 3-4 列,以获得可变的第二个点 lon0/lat0 。关于 -G 如何影响输出记录的信息,参见 输出顺序

-I

进行逆变换,例如:从 (x,y) 数据获取经纬度。

-Ltable[+p][+uunit|c|C]

获取输入点到 table 中的线的最短距离。距离和最近点的坐标将追加到输出的最后三列。追加 +u 可指定距离的单位和距离计算方式,见 单位 ,默认使用大圆距离。例如,其中 unit 也可使用如下字符,c 表示使用输入坐标的笛卡尔距离,C 表示使用投影坐标的笛卡尔距离,C 需要设置 -R-J 。追加 +p 选项可以报告线段 ID 和最近点编号,而不是最近点的坐标。 -L 选项对输出记录的信息顺序的影响,请参见 输出顺序 :对于地理坐标来说,距离的计算使用球模型,因此不能与 -je 选项同时使用。

-N[a|c|g|m]

从大地纬度(使用当前的椭球,参阅 PROJ_ELLIPSOID)转换为其他辅助纬度,经度不受影响。

  • a authalic 纬度

  • c conformal 纬度

  • g 地心纬度, 地心坐标系下的纬度, 默认

  • m meridional

-I 选项同时使用,可实现反向变换。

-Q[d|e]

列出所有投影参数。-Qd,仅列出基准; -Qe 仅列出椭球

-S

不转换区域之外的点

-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 设置的椭球是否正确。

-W[e|E|g|h|j|n|o|O|r|R|w|x][+n[nx[/ny]]]

将地图的宽度和高度等信息输出到标准输出。该选项不需要输入文件。输出的绘图的维度的单位可通过 -D 修改。

  • w 输出地图高度

  • h 输出地图宽度

  • glon/lat 输出 lon/lat 对应的绘图坐标

  • jcode 输出参考点(绘图坐标)对应的地图坐标,参考点请见 锚点

  • r 在使用 -R-J 定义倾斜区域时,输出其覆盖的矩形范围

  • Rr 相同,但输出 -Rw/e/s/n 形式

  • o 在使用 -Rxmin/xmax/ymin/ymax+uunit 将地理坐标投影为笛卡尔坐标时,返回对角点地图坐标,单位为度,输出顺序为 llx urx lly ury

  • Oo 相同,但以 -R 的形式输出

  • e 返回由 -R-J 定义的非矩形区域覆盖的矩形区域

  • Ee 相同,但以 -R 的形式输出

其中,er 选项后追加 +n 时,可以设置输出区域的每条边的点数 nx/ny ,然后输出所有点构成多边形,而不是只输出四个角点

Source Code

../../_images/261dea605b4d8c499c8368fa83fe3ffd.png

倾斜投影区域(红框范围)和常规区域(黑框范围)。

左图中红色区域为倾斜的,可以使用 -Wr|R 获取对应的整个范围的矩形经纬度范围。右图中黑色框是立体投影的范围,红色框为其对应的矩形范围,使用 -We|E 可获取红色框范围

-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

-V[level] (more …)

设置 verbose 等级 [w]

-bi[ncols][type][w][+l|b] (more …)

设置二进制输入数据的格式

-bo[ncols][type][w][+l|b] (more …)

设置二进制输出的数据格式

-d[i|o]nodata (more …)

将输入数据中等于 nodata 的记录替换为 NaN,或将输出数据中值为 NaN 的记录替换为 nodata

-e[~]“pattern” | -e[~]/regexp/[i] (more …)

筛选或剔除匹配指定模式的数据记录

-f[i|o]colinfo (more …)

指定输入或输出列的数据类型

-g[a]x|y|d|X|Y|D|[col]zgap[+n|p] (more …)

确定数据或线段的间断

-h[i|o][n][+c][+d][+msegheader][+rremark][+ttitle] (more …)

跳过或生成指定数目的头段记录

-icols[+l][+sscale][+ooffset][,][,t[word]] (more …)

设置输入数据列及简单变换(0表示第一列,t 表示文本列)

-je|f|g (more …)

设置球面距离的计算方式

-ocols[,…][,t[word]] (more …)

设置输出数据列(0表示第一列,t 表示文本列)

-p[x|y|z]azim[/elev[/zlevel]][+wlon0/lat0[/z0]][+vx0/y0] (more …)

设置3D透视视角

-s[cols][+a|+r] (more …)

设置 NaN 记录的处理方式

-:[i|o] (more …)

交换输入或输出中的第一和第二列

-^-

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

-++

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

-? 或无参数

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

--PAR=value

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

存在的限制

使用 -R 指定的矩形区域在经过投影后,范围变成了非矩形。除非使用 -C, 否则投影后的网格的最左边的坐标将为 0,最下面的坐标也为 0。因此,在对地图数字化前,应该使用 mapproject 对地图的几个角点的坐标极值进行一定的计算,以确定适合的比例尺和地图左下角的坐标。在数字化时,使用这些值以保证正反算结果正确;或者,可以使用类 Unix 系统中的 awk 命令将比例尺和偏移添加到坐标中。

对于某些投影,即使设置了参考椭球为某个椭球,仍然会使用球解法。一般只有 -R 设置范围超出椭球展开的有效范围时会出现这种情况。包括下面几种投影:

  1. Lambert Conformal Comic (-JN) 和 Albers Equal-Area (-JB) 将会在地图比例尺超过 1.0E7 时使用球解法。

  2. Transverse Mercator (-JT) 和 UTM (-JU) 当 -R 指定的范围的东西边界超过中央子午线 10 度以上时,使用球解法;

  3. Cassini (-JC) 与上面相同,但超过 4 度就使用球解法;

  4. 所有的圆锥投影,在经度范围小于 90 度时,使用球解法。

椭球和球

上述计算中,如果用户选择的椭球不为球,且基于椭球的公式已经实现,则 GMT 使用相应的椭球公式,椭球选择见 PROJ_ELLIPSOID。用户需要注意的一些隐患如下:

  1. 对于某些投影,如横轴墨卡托投影,Albers,以及 Lambert 投影,当计算区域较小时,使用椭球公式,当计算区域较大时,则使用球面公式,同时将纬度替换为对应的辅助纬度;具体情况在 存在的限制 章节中已经介绍;

  2. 当用户想检查一些历史数据时,可能会发现 GMT 给出的计算结果和历史数据有轻微的不同,其中一个原因是,早期的 GMT 计算使用的有效数字位数比较少,例如,Clarke 1866 椭球的扁率为 1/294.98,GMT 对扁率做平方后使用,值为 0.00676862818,但是早期的版本为 0.00676866。这会导致十几个厘米的差异。如果用户需要复现以前的计算结果,则需要自己指定椭球扁率。另一个原因可能是,早期的数据使用了和现在不同的基准面,这可能会导致几十甚至几百米的差异。

  3. 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

|-W| 选项的使用

  1. 确定横轴墨卡托投影地图的中点的地图坐标:

    gmt mapproject -R-80/-70/20/40 -Jt-75/1:500000 -WjCM > mid_point.txt
    
  2. 确定倾斜墨卡托投影覆盖的矩形区域:

    gmt mapproject -R270/20/305/25+r -JOc280/25.5/22/69/2c -WR
    
  3. 确定矩形绘图范围对应的实际地图坐标范围:

    gmt mapproject -R-2800/2400/-570/630+uk -Joc190/25/266/68/1:1 -WO
    
  4. 获取矩形绘图范围对应的实际地图范围的多边形,坐标为地理坐标:

    gmt mapproject -R-2800/2400/-570/630+uk -Joc190/25/266/68/1:1 -Wr+n > polygon.txt
    
  5. 确定立体投影确定的非矩形区域对应的矩形范围:

    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.

相关模块

project grdproject