首先解决几个基本的问题。
提问：我有一个tif格式的文件，可以在GMT中直接使用吗？
回答：可以直接使用。因为GMT已经引用了GDAL。
提问：我的tif坐标是UTM，能转成经纬度吗？
回答：可以，有两种方法可以转化，并互相检验。一种是grdproject，一种是先通过grd2xyz转成文本，然后通过mapproject转换为经纬度，在通过xyz2grd转为nc。经检验，两种方法一致，因此直接采用第一种即可。需要注意的是选择UTM投影代号。
其中grdproject用法：
1
| gmt grdproject -Ju51S/1:1 -C -F -I $tif2 -Ggeo2.grd
|
mapproject用法：
1
2
| $ echo 274980 4039220 |gmt mapproject -Ju51S/1:1 -C -I -F
120.488507928 36.4718659261
|
mapproject反转换的用法：
1
2
| $ echo 274980 4039220 |gmt mapproject -Ju51S/1:1 -C -I -F | gmt mapproject -Ju51S/1:1 -C -F
274979.999879 4039220
|
转高斯克吕格，使用tm投影，加带号。或者是UTM加PROJ_SCALE_FACTOR。
1
2
| echo 121 32 | gmt mapproject -Jt123/1:1 -R120/126/0/90 -C -F
echo 274980 4039220 |gmt mapproject -Ju51S/1:1 -C -I -F --PROJ_SCALE_FACTOR=1
|
提问：我有一个nc文件，包含陆地和海洋区域，现在我想把陆地填充为0或者无效值，保留海洋上的数据，怎么办？
回答：可以先通过grdlandmask做一个陆地掩模，然后通过grdmath进行运算。
示例：
1
| gmt grdlandmask -Rgeo2.grd -Df -N1/NaN -Gland_mask.nc -V
|
设置陆地区域的数据值为NaN，海洋区域的数据为1，并生成等网格间隔的nc文件。
1
| gmt grdmath geo2.grd land_mask.nc MUL = geo2.grd
|
进行网格文件的数学运算，两个网格文件相乘，NaN相乘得到NaN，海洋数据和1相乘还是海洋数据。
绘图和数据处理的完整代码，本数据为天宫2数据，点击下载，提取码cdl3，也可以替换成任意tif数据，但需要根据坐标修改UTM的代号。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
|
gmt gmtset FORMAT_GEO_MAP = ddd:mm:ssF
gmt gmtset MAP_FRAME_WIDTH=2p PS_CHAR_ENCODING = ISOLatin1+
gmt gmtset FONT_ANNOT_PRIMARY 7p,Helvetica,black FONT_LABEL 7p,Helvetica,black
tif2=T2_IALT_SRRH_SCI_20170806150925_20170806152709_L4_1_58_20180713110926_V200.tif
ps=tg2_test.ps
gmt grdproject -Ju51S/1:1 -C -F -I $tif2 -Ggeo2.grd
gmt grdlandmask -Rgeo2.grd -Df -N1/NaN -Gland_mask.nc -V
gmt grdmath geo2.grd land_mask.nc MUL = geo3.grd
gmt pscoast -Ju51S/1:1500000 -Df+ -W0.5p,blue -Ba12mg12m -BNESW -K --MAP_ANNOT_OBLIQUE=40 -R >$ps
gmt grd2cpt geo2.grd -Crainbow > mydata.cpt
gmt grdimage geo2.grd -R -J -Cmydata.cpt -K -Q -O >> $ps
gmt pscoast -Ju51S/1:1500000 -Df+ -W0.5p,blue -Ba12mg12m -BNESW -K --MAP_ANNOT_OBLIQUE=40 -R -O -X2.5i >>$ps
gmt grdimage geo3.grd -R -J -Cmydata.cpt -Q -O >> $ps
gmt psconvert $ps -A -P -Tg
rm *.nc *.cpt *.grd
|