-J 选项

-J 选项用于指定坐标变换方式或地图投影方式,即将数据投影到画布上所采用的函数。

-J 选项后接投影代码以及投影参数。GMT 目前支持两种不同的方式指定投影代码和投影方式:

  • GMT 投影代码

  • GMT + PROJ 投影代码

GMT 投影代码

-J 选项有两种写法:

-J\(\delta\)[parameters/]scale

-J\(\Delta\)[parameters/]width

其中,δΔ 用于指定投影代码,前者为小写字母,后者为大写字母。 parameters 是零个或多个由斜杠分隔的投影参数,参数数目由投影方式决定。

投影代码使用小写字母时,-J 的最后一个参数 scale 表示底图比例尺, 即图上距离与真实地球距离之间的换算关系。scale 可以有两种格式:

  • 单个数字加单位,例如 2c,表示真实地球距离的 1 度投影到画布上为 2 厘米

  • 1:xxxx 格式,例如 1:10000000 表示画布上的 1 厘米对应真实地球距离的 10000000 厘米

投影代码为大写字母时,-J 的最后一个参数 width 表示底图宽度。还可以加上 子选项还进一步控制 width 的含义:

  • +dw 表示 width 为底图的宽度 [默认值]

  • +dh 表示 width 为底图的高度

  • +du 表示 width 为底图维度中最大的那个维度的长度

  • +dl 表示 width 为底图维度中最小的那个维度的长度

备注

几乎所有地图投影方式,都只能指定底图宽度或高度中的任一个,而不能同时指定二者, 因为地图高度会由地图宽度和投影方式自动决定。

画图时通常建议使用大写投影代码以直接指定图片宽度,除非需要明确指定比例尺。

例如:

  • -Jm1c 表示使用墨卡托投影,地图上的 1 度距离投影到画布上为 1 厘米

  • -Jm1:10000000 表示使用墨卡托投影,画布上的 1 cm代表实际距离中的 10000000 cm, 即 100 km

  • -JM15c 也表示使用墨卡托投影,整个地图的宽度是 15 厘米,地图的高度由 -R-J 自动确定

  • -JM15ch 表示整个地图的高度是 15 厘米

  • -JX10c/5c 使用线性投影,地图的宽度是 10 厘米,高度为 5 厘米

下表列出了 GMT 所支持的全部投影方式,详细介绍见 地图投影

GMT 投影代码

-J 代码

说明

-JA|lon0|/|lat0|[/horizon]/width

Lambert azimuthal equal area

-JB|lon0|/|lat0|/|lat1|/|lat2|/width

Albers conic equal area

-JC|lon0|/|lat0|/width

Cassini cylindrical

-JCyl_stere[|lon0|/[|lat0|/]]width

Cylindrical stereographic

-JD|lon0|/|lat0|/|lat1|/|lat2|/width

Equidistant conic

-JE|lon0|/|lat0|[/horizon]/width

Azimuthal equidistant

-JF|lon0|/|lat0|[/horizon]/width

Azimuthal gnomonic

-JG|lon0|/|lat0|[/horizon]/width

Azimuthal orthographic

-JG|lon0|/|lat0|/alt/azim/tilt/twist/W/H/width

General perspective

-JH[|lon0|/]width

Hammer equal area

-JI[|lon0|/]width

Sinusoidal equal area

-JJ[|lon0|/]width

Miller cylindrical

-JKf[|lon0|/]width

Eckert IV equal area

-JKs[|lon0|/]width

Eckert VI equal area

-JL|lon0|/|lat0|/|lat1|/|lat2|/width

Lambert conic conformal

-JM[|lon0|/[|lat0|/]]width

Mercator cylindrical

-JN[|lon0|/]width

Robinson

-JOa|lon0|/|lat0|/azim/width[+v]

Oblique Mercator, 1: origin and azim

-JOb|lon0|/|lat0|/|lon1|/|lat1|/width[+v]

Oblique Mercator, 2: two points

-JOc|lon0|/|lat0|/|lonp|/|latp|/width[+v]

Oblique Mercator, 3: origin and pole

-JPwidth[+a][+f[e|p|radius]][+kkind][+roffset][+torigin][+z[p|radius]]

Polar/Cylindrical

-JPoly[|lon0|/[|lat0|/]]width

(American) polyconic

-JQ[|lon0|/[|lat0|/]]width

Equidistant cylindrical

-JR[|lon0|/]width

Winkel Tripel

-JS|lon0|/|lat0|[/horizon]/width

General stereographic

-JT[|lon0|/[|lat0|/]]width

Transverse Mercator

-JUzone/width

Universal Transverse Mercator (UTM)

-JV[|lon0|/]width

Van der Grinten

-JW[|lon0|/]width

Mollweide

-JXwidth[l|pexp|T|t][/height[l|pexp|T|t]][d]

Linear, logarithmic, power, and time

-JY|lon0|/|lat0|/width

Cylindrical equal area

GMT + PROJ

从 GMT6 开始,GMT 支持使用 PROJ 库来实现坐标和基准面的转换。这一特性是通过 GDAL 实现的。详细的 PROJ 语法见 https://proj.org/en/latest/apps/index.html

在 GMT 中使用 PROJ 和单独使用 PROJ 提供的 projcs2cs 命令非常相似。 在 PROJ 中,投影一般有很多参数,多个参数之间用空格分隔。在 GMT 中,可以将所有参数 用双引号括起来:

-J"+proj=merc +ellps=WGS84 +units=m"

或直接将所有参数连在一起:

-J+proj=merc+ellps=WGS84+units=m

也可以直接使用 EPSG codes。 例如 -JEPSG:4326 表示使用 WGS-84 系统。

对于 mapprojectgrdproject 模块,可以直接使用 +to 关键字直接指定要将 A 参考系统转换为 B 参考系统,而不需要中间步骤。例如:

-JEPSG:4326+to+proj=aeqd+ellps=WGS84+units=m

对于使用 mapprojectgrdproject 进行点和网格 文件的转换,GMT 可以使用 所有的 PROJ 投影。 但对于绘图而言,其用处却很有限。一方面,只有一部分 PROJ 的投影方式可以被 映射到 GMT 的投影语法中, 只能使用这些投影绘图;另一方面,由于 PROJ 不是一个绘图库, 其不支持设置地图比例尺或地图大小,因而,GMT 为 PROJ 语法引入了两个扩展:+width=size+scale=1:xxxx 使得其与经典的 GMT 中的工作方式相似。 此外,也可以在投影参数的最后加上字符串 /1:xxx 来指定比例尺。

下面将给出一些示例:

高斯投影

虽然高斯投影可以用 GMT 自带的 UTM 投影做一定的修改近似实现,但由于分带和轴平移等设置 的不同,使用多有不变,同时近似导致在某些情况下精度不足。下面介绍在 GMT 中使用 PROJ 做高斯投影。

假定存在一点经纬度坐标 (87,0),现将其转换为西安 80 坐标系,西安 80 坐标系 为投影坐标系,为减小变形,使用 3 度或者 6 度高斯投影分带分别用于大比例尺和中比例尺 应用。上述坐标在 3 度分带中的代号为 29,通过 EPSG codes 可以查询到其对应的 EPSG 代码为 2353, 因此,可用如下代码投影:

$ echo 87 0 | gmt mapproject -J"EPSG:2353"
29500000    0

结果中横轴坐标 2950000 中 29 为带号,根据投影规则,为防止投影后的横轴坐标出现负值,投影 后需将横坐标加 500 km,由于经度 87 度刚好位于该投影带的中央经线上,所以其横坐标为 500000; 纬度为 0 度,位于赤道上,其纵坐标也为 0。在上述 EPSG 代码查询中,同时可获得其 PROJ 投影 语法为 +proj=tmerc +lat_0=0 +lon_0=87 +k=1 +x_0=29500000 +y_0=0 +a=6378140 +b=6356755.288157528 +units=m +no_defs, 因此,将有下面的等价语句:

$ echo 87 0 | gmt mapproject -J"+proj=tmerc +lat_0=0 +lon_0=87 +k=1 +x_0=29500000 +y_0=0 +a=6378140 +b=6356755.288157528 +units=m +no_defs"
29500000    0

上述仅给出 3 度带投影的例子,6 度带投影类似:

$ echo 87 0 | gmt mapproject -J"EPSG:2329"
15500000    0

投影转换

下面使用 +to 来实现投影转换,将上述 3 度带投影转换为 6 度带投影:

$ echo 29500000    0 | gmt mapproject -J"EPSG:2353 +to EPSG:2329"
15500000    0

将其中的 EPSG 代码替换为对应的代码则可实现任意投影的转换。