GMT6绘制北半球

news2025/2/23 10:33:25

设置绘制区域及投影方式

投影方式选择立体等角投影，在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