设置绘制区域及投影方式
投影方式选择立体等角投影，在GMT6中的命令是-Js
# 定义区域变量和投影变量，纬度从北纬30度到极点
region='-180/180/30/90'
projection='0/90/1:60000000'
gmt set PROJ_ELLIPSOID WGS-84
定义CPT及地形展示
现在定义一个CPT用于显示地形的颜色
# 定义一个CPT
gmt makecpt -Cearth -T-7000/4000 > land.cpt
# 绘制地形并添加海岸线
gmt grdimage -R$region -Js$projection @earth_relief_01m -I+d -Cland.cpt
gmt coast -R$region -Js$projection -Bafg -Di -W0.25p -A10000
完整代码：
#!/usr/bin/env -S bash -e
# GMT modern mode bash template
# Date: 2024-06-13T20:17:44
# User: br
# Purpose: Purpose of this script
set -e
export GMT_SESSION_NAME=$$ # Set a unique session name
gmt makecpt -Cearth -T-7000/4000 > land.cpt
gmt begin North_Hemispher png
# Place modern session commands here
region='-180/180/30/90'
projection='0/90/1:60000000'
gmt set PROJ_ELLIPSOID WGS-84
gmt grdimage -R$region -Js$projection @earth_relief_01m -I+d -Cland.cpt
gmt coast -R$region -Js$projection -Bafg -Di -W0.25p -A10000
gmt end show
结果展示：
只展示陆地地形
现在我有个需求只想展示陆地的地形，因为后面还要叠加海洋温度的数据
此时就需要使用GMT中提供的掩膜方法将陆地的地形数据提取出来，遮住海洋地形的数据
使用grdlandmask模块创建分辨率为1m的陆地掩膜，使用grdcut模块将陆地掩膜的范围裁剪到指定区域方便后续的数学运算，使用grdmath模块中的乘法操作将裁剪后的地形数据与陆地掩膜相乘。这样就得到了陆地的地形数据
gmt grdlandmask -R-180/180/30/90 -I1m -Dl -Gland_mask.nc
gmt grdcut @earth_relief_01m -R-180/180/30/90 -Gearth_relief_cut.nc
gmt grdmath earth_relief_cut.nc land_mask.nc MUL = masked_land.nc
然后将上面的完整代码稍作修改即可
#!/usr/bin/env -S bash -e
# GMT modern mode bash template
# Date: 2024-06-13T20:17:44
# User: br
# Purpose: Purpose of this script
set -e
export GMT_SESSION_NAME=$$ # Set a unique session name
gmt makecpt -Cearth -T-7000/4000 > land.cpt
gmt begin North_Hemispher png
# Place modern session commands here
region='-180/180/30/90'
projection='0/90/1:60000000'
gmt set PROJ_ELLIPSOID WGS-84
gmt grdimage -R$region -Js$projection masked_land.nc -I+d -Cland.cpt
gmt coast -R$region -Js$projection -Bafg -Di -W0.25p -A10000
gmt end show
结果展示：
海洋温度数据
我这里有一个海温的数据，是ERA5月温度平均，是一张721*1440的0.25度的TIF文件
海温的范围大致在270-300K。陆地上没有数据，用白色表示
QGIS展示：
海温数据叠加到北半球
我们可以使用gdalinfo工具查看海温数据的相关信息：
gdalinfo SST.tif
相关信息：
Driver: GTiff/GeoTIFF
Files: SST.tif
SST.tif.aux.xml
Size is 1440, 721
Origin = (-179.875000000000000,90.000000000000000)
Pixel Size = (0.250000000000000,-0.250000000000000)
Image Structure Metadata:
INTERLEAVE=BAND
Corner Coordinates:
Upper Left (-179.8750000, 90.0000000)
Lower Left (-179.8750000, -90.2500000)
Upper Right ( 180.1250000, 90.0000000)
Lower Right ( 180.1250000, -90.2500000)
Center ( 0.1250000, -0.1250000)
Band 1 Block=1440x1 Type=Float32, ColorInterp=Gray
Min=269.884 Max=304.888
Minimum=269.884, Maximum=304.888, Mean=286.909, StdDev=11.739
Metadata:
STATISTICS_MAXIMUM=304.88809204102
STATISTICS_MEAN=286.90908880145
STATISTICS_MINIMUM=269.88391113281
STATISTICS_STDDEV=11.739334479632
STATISTICS_VALID_PERCENT=66.11
显然直接将海温数据叠加到北半球上是不行的，首先需要对海温数据的范围进行裁剪，裁剪到指定区域
gmt grdcut SST.tif -R-180/180/30/90 -Gcut_SST.tif
接下来就可以进行叠加了。需要注意的是在GMT中后面的grdimage图层会覆盖掉前面的图层，所以需要使用-t选项设置一定的透明度。同时给海洋温度添加rainbow的CPT。完整代码：
#!/usr/bin/env -S bash -e
# GMT modern mode bash template
# Date: 2024-06-13T20:17:44
# User: br
# Purpose: Purpose of this script
set -e
export GMT_SESSION_NAME=$$ # Set a unique session name
gmt makecpt -Cearth -T-7000/4000 > land.cpt
gmt makecpt -Crainbow -T270/300/1 > ocean.cpt
gmt begin North_Hemispher png
# Place modern session commands here
region='-180/180/30/90'
projection='0/90/1:60000000'
gmt set PROJ_ELLIPSOID WGS-84
gmt grdimage -R$region -Js$projection masked_land.nc -I+d -Cland.cpt
gmt grdimage -R$region -Js$projection cut_SST.tif -Cocean.cpt -t30 # -t设置透明度为30%
gmt coast -R$region -Js$projection -Bafg -Di -W0.25p -A10000
gmt colorbar -DJRM+w6.5c/0.5c+o1c/0+mc -Bxa -By -G270/300
gmt end show
结果展示：
可以看到北极圈附近的海温相对较低，往南海温逐渐升高
