[转载]GMT绘制不同投影方式的地形图

以ETOPO5为底图数据，使用GMT绘制多种投影的全球地形图。

ETOPO有多种类型的数据，之前有ETOPO5和ETOPO2版本，最新的是ETOPO1版本，其分辨率是1′，虽然不如SRTM的分辨率，但是ETOPO数据集成了海底地形和极地冰盖高程数据，对于绘制一些海洋和极地区域的地形图很有价值。

(下载：ETOPO1 Global Relief | NCEI )

（1）Robinson投影

首先以Robinson投影为例子，绘制全球陆地地形和海底地形图。

下面是bat代码：

::xyz2grd ETOPO5.DAT -Getopo5.grd -I5m -Rg -ZTLxyhw

::makecpt -Cglobe -T-10000/10000/1000 -Z > colors.cpt


grdimage ETOPO1_Ice_g_gmt4.grd -Ba60g30 -Rg -JN180/25c  -Ccolors.cpt -K > world.ps

psscale -Ba2500f500::/:"m": -Ccolors.cpt -D12.5c/0.3c/15c/.35ch -Y-0.5i -O -U"Author:YangLei">> world.ps

psconvert world.ps -A0.5c -E600 -Tj

pause

GMT图片格式转换--psconvert

1. 转换为普通图片格式

gmt psconvert input.ps -A0.5c -E600 -Tp -V
解释：
input.ps是待转换的PS图件
-A是切除画板（A4页面）上的留白。不带任何参数值是切除所有的留白，0.5c是在四个边框都留出0.5c（这样好看一点），也可以用0.5c/0.4c/0.3c/0.2c这种方式分别指定四边的多余留白
-E是设置分辨率，默认300dpi，这里设置为600
-T是设置输出图片的格式：
- g表示png
- j表示jpg
- f表示pdf
- t表示tiff
详情参见官方文档，或在命令行输入：gmt psconvert --help
-V输出所有运行的信息。

2. 转换为KML

要点有两个：
1. 投影方式必须为笛卡尔，即-JX
猜测原因，笛卡尔投影相当于线性投影，基本没有转换，而kml格式是直接覆盖在三维地球上，不需要投影，相当于要把GMT做的投影反变换回来，所以用笛卡尔投影才能反变换？？
2. 设置图片边框 在开头加入：

gmt gmtset MAP_FRAME_TYPE=inside

这样转为kml后图件的边框才和google earth的坐标对的上
如果是plain的话，则是图件的边缘和真实坐标对的上。注意下图海南岛偏了，就是这个缘故。

1. gmt psconvert input.ps -W+k
   -W+k即表示输出kml

3. 转换为geotiff

需要安装：gdal_translate,

安装方式：

- Windows下：64位下载，默认会安装在：C:\OSGeo4W64\bin这个地址，所以要把这个地址加入系统环境变量中，如何加参见：windows下添加PATH环境变量

- Linux/ubuntu下： sudo apt install gdal-bin

在GMT中输出geotiff要点也是两个：
1. 投影方式不能为笛卡尔，可以为墨卡托或其他，如：-JM 2. gmt psconvert input.ps -W+g
注意，输出的带地理坐标的图件后缀是.tiff，
用envi打开tiff发现出现了偏色情况，因为导出来的tiff图件不是单波段的，在导入envi是做了彩色变换。
技巧:如果tif图件是单波段的，可以用gdal_translate转换为grd文件，再作为grdimage的输入数据。

GMT图片格式转换--psconvert - 知乎

（2） Eckert IV projection[equal-area] 

该投影也专门用于绘制全球地图，属于伪圆柱等积投影。其特征是中央子午线为直线，最外侧子午线为半圆形，Scale is true along the parallel at 40:30 North and South. 

需要修改grdimage 程序中的投影为-JKf180/5c，GMT绘制的全球地图如下：

（3）Sinusoidal projection

修改投影为 -JI0/25c，结果如下：

（4）Interrupted Sinusoidal projection

其绘制方法有点区别，需要分次绘制

grdimage etopo5.grd -Ba60g30 -R-160/-20/-90/90 -Ji-90/0.02i -Ccolors.cpt -K > world.ps

grdimage etopo5.grd -Ba60g30 -R-20/60/-90/90 -Ji20/0.02i -X2.8i -Ccolors.cpt -K -O>> world.ps

grdimage etopo5.grd -Ba60g30 -R60/200/-90/90 -Ji130/0.02i -X1.6i -Ccolors.cpt -K -O>> world.ps

（5）墨卡托投影

投影方式设置为-JM25c，还要注意-R的设置，墨卡托投影不能设置为-Rg。将-R设置为不同区域时，因为投影中心变化大，所以需要自己注意设置。 投影纬度区域如果设置为-90/90无法成图，如果设置为-85/85之间图件出幅。采用墨卡托投影距离极点跃进，变形误差越大，无法采用。

grdgradient的-Ne后的数值调节幅度。数值越大越明显。

使用-X -Y选项 ，设置图幅上下移动位置。

makecpt -Cglobe -T-10000/10000/1000 -Z > colors.cpt

grdgradient ETOPO1_Ice_g_gmt4.grd -A0 -Ne1 -fg -Gibcso_v1_hd.nc
grdimage ETOPO1_Ice_g_gmt4.grd  -Iibcso_v1_hd.nc -R-180/180/-75/75   -Ba30g30/a20g20 -JM0/0/25c -Ccolors.cpt -K > world.ps

pscoast -JM0/0/25c -R-180/180/-75/75 -N1  -O -K >>world.ps

psscale -Ba2500f500::/:"m": -Ccolors.cpt -D12.5c/0.3c/15c/.35ch -Y-0.5i -O -U"Author:YangLei">> world.ps

psconvert world.ps -A0.5c -E600 -Tj

pause

分享：