使用GMT实现SHP矢量裁剪栅格与山体阴影地形图制作-平芜编程栈

简介：本资源面向地理信息科学、地球物理及遥感方向的科研人员与高年级本科生，提供一套基于GMT（Generic Mapping Tools）实现“本地Shapefile裁剪栅格+叠加山体阴影”的完整实操方案。资源解决实际制图中常遇的区域限定与地形可视化难题，适用于地形图制作、论文插图生成及教学演示等场景。压缩包共11个文件（289KB），涵盖shp矢量数据（含.shp/.shx/.dbf/.prj等标准组件）、GMT脚本（.gmt）、批处理命令（.bat）、中间结果图（.png）及元数据文件（.xml/.cpg等），结构规范，可直接运行调试。已有82人学习下载，用户可即刻获取可复用的裁剪流程脚本、光照参数配置模板、shading透明度控制技巧及配套地理数据示例，显著降低GMT高级制图的学习门槛与试错成本。

1. 项目概述：当本地矢量遇上全球栅格

在地球科学、地理信息或者遥感数据处理领域，我们常常会遇到一个非常具体的需求：手里有一份描述特定区域的矢量边界文件（比如一个县的行政边界、一个流域的范围，或者一个研究区的轮廓），同时还有一份覆盖范围更大的栅格数据（比如全球地形数据、卫星影像或者气候模型输出）。我们的目标很明确，就是“按图索骥”，用这个矢量边界像 cookie cutter（饼干模具）一样，从大块的栅格“面团”里，精准地切出我们关心的那一部分。切出来之后，往往还不够直观，我们还想给它“上个色”，通过山体阴影（Hillshade）这种经典的地形渲染手法，让高程起伏跃然纸上，生成一幅既专业又美观的地形图。

这个需求听起来简单直白，但实操起来，从工具选择、数据准备、命令执行到结果美化，每一步都有不少门道。你可能会想到用 ArcGIS 或者 QGIS 这类图形化桌面软件，拖拖拽拽也能完成。但对于需要批量处理、流程自动化，或者追求更高可重复性和计算效率的场景，命令行工具的优势就凸显出来了。GMT（Generic Mapping Tools）正是这样一个在学术界和工业界被广泛信赖的“瑞士军刀”。它没有花哨的界面，全靠一行行命令驱动，正因如此，它极其强大、灵活，且能无缝嵌入到脚本中实现流程化作业。

本项目要解决的，就是如何纯粹使用 GMT，利用本地的 Shapefile（.shp）矢量文件去裁剪一个栅格文件（如 GeoTIFF 格式的 DEM 数据），并为裁剪后的结果生成高质量的山体阴影图。这个过程不依赖任何其他图形化 GIS 软件作为前置或后置处理，完全在 GMT 的生态内完成。无论是制作区域地质图、规划图，还是进行局地地形分析，这都是一个基础且核心的技能栈。

2. 核心工具与数据准备

2.1 为什么选择 GMT？

在开始动手之前，我们先聊聊为什么是 GMT，而不是其他工具。首先，跨平台与一致性：GMT 在 Windows、macOS 和 Linux 上表现一致，编写的脚本可以几乎无修改地在不同系统间迁移，这对于团队协作和服务器部署至关重要。其次，出版级绘图质量：GMT 诞生于地球物理学领域，其默认的绘图美学和精细的可控性（如刻度、字体、颜色表）使其产出直接可用于学术出版物，这是很多工具的短板。再者，强大的数据处理能力：GMT 并非单纯的绘图工具，它内置了海量的数据操作模块，如网格处理、滤波、投影变换等，裁剪（grdcut）和山体阴影计算（grdgradient）只是其冰山一角。最后，脚本化与自动化：所有操作通过命令记录，形成脚本后，可以轻松修改参数、批量处理数据，确保了研究的可重复性。

当然，它的学习曲线相对陡峭，需要记忆命令和参数。但一旦掌握，效率的提升是巨大的，尤其当你需要处理成百上千个区域和数据集时。

2.2 理解我们的“原料”：SHP 与栅格文件

1. Shapefile 文件 (.shp)：你提供的“本地 shp 文件”通常不仅仅是一个 .shp 文件。一个完整的 Shapefile 由至少三个文件构成：

.shp：存储几何图形（点、线、面）的主体文件。
.shx：几何图形的索引文件，用于快速定位。
.dbf：属性表文件，存储每个几何图形对应的属性信息（如名称、编码）。有时还会有.prj（投影信息文件）、.cpg（字符编码文件）等。GMT 在读取 SHP 文件时，主要需要.shp和.dbf，并且非常依赖.prj来正确理解数据的坐标系统。如果缺失.prj，你可能需要手动通过-f参数指定坐标格式。

2. 栅格文件：常见的栅格文件格式有 GeoTIFF (.tif)、NetCDF (.nc)、GMT 自定义的 .grd 格式等。它们本质上是一个数值矩阵，每个像素（栅格单元）都有一个值（如高程、温度、浓度），并带有地理参考信息（范围、像素大小、投影）。本项目假设我们有一个全球或大范围的 DEM（数字高程模型）文件，例如 SRTM 数据或 ASTER GDEM 数据，格式为 GeoTIFF。

3. 山体阴影（Hillshade）：这是一种通过模拟太阳光照在地形表面的效果来增强地形立体感的可视化技术。它需要两个核心参数：太阳方位角（Azimuth，光源方向，0°为正北，顺时针增加）和太阳高度角（Altitude，光源高度，0°为地平线，90°为正头顶）。计算原理是根据 DEM 的坡度（Slope）和坡向（Aspect）来模拟每个栅格单元的亮度。GMT 的grdgradient命令就是专门干这个的。

2.3 环境搭建与数据检查

在运行任何命令之前，确保你的 GMT 版本在 6.0 以上（推荐使用最新版）。你可以通过终端运行gmt --version来检查。

数据检查是关键的第一步：

检查 SHP 文件的投影：用ogrinfo（GDAL 工具包的一部分，通常随 GMT 或 QGIS 安装）查看。
```
ogrinfo -al -so your_boundary.shp
```
重点关注输出中的Layer SRS WKT:部分。确认其投影是否是地理坐标系（如 WGS84，EPSG:4326），或者是与你的栅格数据相匹配的投影坐标系。如果 SHP 和栅格的投影不同，我们需要在 GMT 中进行投影转换或统一。
检查栅格文件信息：使用 GMT 的gmt grdinfo命令。
```
gmt grdinfo your_dem.tif
```
这会输出栅格的范围（xmin, xmax, ymin, ymax）、像素大小、数据范围、投影等信息。同样，记下它的投影。

注意：如果 SHP 文件和栅格文件的投影不一致，直接裁剪会导致错误或结果错位。我们必须在同一个坐标参考系下工作。通常的作法是将所有数据统一到地理坐标系（如 WGS84）或某个适合区域的投影坐标系（如 UTM）。GMT 的gmt grdproject和gmt grdsample可以用于栅格的重投影和重采样，而 SHP 文件可以在 QGIS 中预先重投影，或者利用 GMT 的gmt convert配合-R和-J参数在读取时进行实时转换（这需要更高级的技巧）。为了简化，本项目假设两者投影一致，或 SHP 文件已预先转换为与栅格相同的投影。

3. 核心步骤拆解与命令详解

整个流程可以清晰地分为四个阶段：数据准备、区域裁剪、山体阴影计算、地图绘制与合成。

3.1 第一阶段：提取并格式化裁剪区域

GMT 的裁剪命令grdcut需要一个明确的区域范围，通常由-R参数指定（如-Rminlon/maxlon/minlat/maxlat）。但我们的 SHP 文件可能是不规则的多边形，直接用其外包矩形（Bounding Box）裁剪可能会带入无关区域。更精确的做法是使用 SHP 文件本身作为掩膜（Mask）。GMT 6 提供了更强大的栅格操作能力，我们可以结合gmt grdclip、gmt grdmask和gmt grdmath来实现。

步骤 3.1.1：将 SHP 文件转换为 GMT 兼容的矢量格式虽然 GMT 可以直接读取.shp，但在复杂流程中，先将其转换为 GMT 更“原生”的格式（如 OGR_GMT 格式或简单的 ASCII 表）会更稳定。这里我们使用gmt convert。

gmt convert your_boundary.shp -fg > boundary.xy

-fg参数表示输入数据是地理坐标（经度/纬度）。如果你的 SHP 是投影坐标，则需要使用对应的-f参数，例如-fx表示横坐标/纵坐标。
boundary.xy是一个文本文件，里面包含了多边形的顶点坐标。如果 SHP 有多个多边形，它们会被记录为多个段（Segment），由>符号分隔。

步骤 3.1.2：创建二值掩膜网格（Mask Grid）这是实现精确裁剪的核心。我们创建一个和目标裁剪区域分辨率一致的空白网格，然后在 SHP 多边形内的区域赋值为 1，外的区域赋值为 NaN（非数字）或 0。

# 首先，从原始DEM获取我们大致需要的网格参数（范围、间距） gmt grdinfo your_dem.tif -Ie > grid_info.txt # 假设 grid_info.txt 内容为：-R-180/180/-90/90 -I0.00833333333333/0.00833333333333 # 我们可以用这些参数来创建掩膜网格，但为了更高效，我们通常用原DEM的网格来初始化。 # 方法：使用 grdmask 命令 gmt grdmask boundary.xy -Gmask.grd -Ryour_dem.tif -I+ -N0/1/1 -V

-Gmask.grd: 输出掩膜网格文件。
-Ryour_dem.tif: 设置网格范围与原始 DEM 一致。你也可以用-Rminlon/maxlon/minlat/maxlat手动指定一个稍大于 SHP 边界的范围以节省计算量。
-I+: 这是一个关键参数！它告诉 GMT 使用与-R所关联的网格（这里是your_dem.tif）相同的像素间距（Increment）。这确保了掩膜网格和原 DEM 在空间上完全对齐，避免后续操作出现像素错位。
-N0/1/1: 设置掩膜值。这里意味着：多边形外的值为 0，多边形内的值为 1，多边形边界上的值也为 1。
-V: 显示详细过程信息，便于调试。

执行后，你会得到mask.grd。你可以用gmt grdimage mask.grd -JM10c -B -C -P > mask.png快速查看，它应该是一个在你 SHP 区域内为白色（1），区域外为黑色（0）的图片。

3.2 第二阶段：执行栅格裁剪

有了精确对齐的掩膜，裁剪就变成了一个网格间的数学运算。

gmt grdmath your_dem.tif mask.grd MUL = dem_cropped.grd

grdmath: GMT 的网格计算器。
your_dem.tif mask.grd MUL: 将原始 DEM 网格与掩膜网格逐像素相乘。在掩膜为 1 的区域（内部），DEM 值保持不变；在掩膜为 0 的区域（外部），结果变为 0。但地形数据外部我们通常希望是 NaN（透明），所以更好的掩膜值设置是内部为 1，外部为 NaN。
= dem_cropped.grd: 输出结果。

更优的掩膜与裁剪方案：为了让裁剪区域外的部分在绘图时完全透明（而非 0 值），我们调整grdmask参数，并组合使用grdclip。

# 创建掩膜，外部为NaN，内部为1 gmt grdmask boundary.xy -Gmask_nan.grd -Ryour_dem.tif -I+ -NNaN/1/1 -V # 使用 grdclip 进行裁剪 gmt grdclip your_dem.tif -Gdem_raw_cut.grd -Smask_nan.grd/1/NaN -V

-NNaN/1/1: 现在多边形外为 NaN，内为 1。
grdclip ... -Smask_nan.grd/1/NaN:-S参数指定一个“掩膜网格”和阈值。这里的意思是：对于your_dem.tif的每个像素，查看mask_nan.grd中对应位置的值。如果该值等于1，则保留原 DEM 值；否则，将输出值设置为 NaN。

这样得到的dem_raw_cut.grd就是一个在 SHP 区域内保留原始高程值，区域外全部为 NaN 的完美裁剪网格。NaN 值在 GMT 绘图时会被自动忽略，实现“透明”效果。

3.3 第三阶段：生成山体阴影

现在我们对裁剪后的 DEM (dem_raw_cut.grd) 计算山体阴影。

gmt grdgradient dem_raw_cut.grd -Gdem_hillshade.grd -A315 -Nt0.8 -V

-Gdem_hillshade.grd: 输出山体阴影强度网格。
-A315: 设置太阳方位角为 315°，即西北方向照射。这是最常用的角度之一，能产生良好的立体感。
-Nt0.8:-N指定归一化方式。t0.8是一个组合参数：t表示使用 “Lambertian 反射 + 强度调整” 的算法，0.8是增强因子（1.0 为标准）。-Nt1是默认值，-Nt0.8或-Nt0.6可以使阴影对比度更柔和，避免过强的“浮雕感”，视觉效果更自然。这是个人偏好和数据类型决定的，需要多尝试。
为什么不直接用-Ne0.8？e表示简单的强度缩放，而t算法结合了朗伯反射模型，考虑了坡度，效果更符合物理光照，通常视觉效果更好。

得到的dem_hillshade.grd是一个灰度网格，值在 -1 到 1 之间（经过归一化），亮部接近 1，暗部接近 -1。

3.4 第四阶段：绘制合成地图

这是展现成果的一步。我们将山体阴影作为底图（营造立体感），再将裁剪后的 DEM 以颜色映射（CPT）的方式半透明覆盖在上面（显示高程分层）。

# 1. 为DEM创建颜色梯度表（CPT） gmt makecpt -Cdem2 -T0/3000/100 -Z > elevation.cpt # 解释：-Cdem2 使用GMT内置的‘dem2’配色（适合地形），-T0/3000/100 表示颜色从0米到3000米，每100米一个间隔，-Z 表示创建连续的颜色表。 # 2. 开始绘制 gmt begin final_map png # 设置绘图区域和投影。这里使用UTM投影，假设区域中心经度是105E，比例尺为1:500000。 gmt basemap -Rdem_raw_cut.grd -JUtm105/1:500000 -Baf -BWSen+t"区域地形图" # 首先，绘制山体阴影底图。使用 -I 指定强度网格，-M 表示使用该网格调制亮度。 gmt grdimage dem_hillshade.grd -I -M # 然后，半透明叠加高程填色图。-t50 表示设置50%的透明度。 gmt grdimage dem_raw_cut.grd -Celevation.cpt -t50 # 可选：叠加SHP边界线，使其更突出 gmt plot boundary.xy -W1p,black -V # 添加比例尺和图例 gmt colorbar -Celevation.cpt -Baf+l"Elevation (m)" -DJMR+o0.8c/0+w10c/0.3c gmt end

gmt begin/end: GMT 6 的现代模式，用于管理一个绘图会话，比古典模式更简洁。
-Rdem_raw_cut.grd: 自动使用裁剪后网格的范围作为绘图范围，非常方便。
-JUtm105/1:500000: 使用 UTM 投影，中央经线 105°，比例尺 1:500000。你需要根据你的数据位置选择合适的投影和参数。对于小区域，UTM 是常用选择；对于大范围或特定需求，可能需要-JM（墨卡托）、-JX（线性）等。
grdimage ... -I -M:-I选项后面没有直接跟文件名，是因为在gmt begin/end块中，上一个grdimage生成的“强度网格”上下文会被记住。-M表示使用强度调制（即山体阴影效果）。
-t50: 透明度参数是关键。它让上层的颜色填色图与下层的山体阴影融合，既显示了高程色彩，又保留了地形纹理，效果远好于单独使用任何一种。

执行后，将生成final_map.png，一幅具有专业质感的地形图就诞生了。

4. 常见问题、技巧与深度优化

4.1 踩坑记录与问题排查

问题1：执行grdmask或grdclip时速度极慢，或内存溢出。

原因：原始 DEM 分辨率过高（如 30米 SRTM），覆盖范围很大，而-I+参数会创建一个与原始 DEM 全范围同分辨率的掩膜网格，导致网格像元数爆炸。
解决：
1. 限制范围 (-R)：不要用-Ryour_dem.tif，而是先用gmt grdinfo your_dem.tif查看范围，然后用gmt minmax boundary.xy查看 SHP 的大致范围。取两者的交集并适当外扩，手动定义一个更小的-R范围。例如-R99.5/101.2/24.8/26.5。
2. 降低分辨率：如果最终成图比例尺不大，可以考虑对原始 DEM 和 SHP 边界进行重采样到较低分辨率。使用gmt grdsample降低 DEM 分辨率，用gmt grdtrack或gmt sample1d对边界线进行抽稀（但需小心保持多边形闭合性）。
3. 分块处理：对于超大区域，可以考虑用gmt grdcut先粗剪一个大致区域，再在此基础上进行精细掩膜。

问题2：生成的山体阴影看起来非常“脏”，有大量斑驳的噪声。

原因：DEM 数据本身可能存在噪声（特别是 SRTM 数据在水体或平坦区域），或者太阳高度角 (-E) 设置过低（默认 45°），放大了微小起伏。
解决：
1. 对 DEM 进行平滑滤波：在计算山体阴影前，先使用gmt grdfilter进行轻度高斯平滑。
```
gmt grdfilter dem_raw_cut.grd -Gdem_smooth.grd -Fg1 -D2
```
  -Fg1表示使用高斯滤波器，标准差为 1 个像素。-D2表示在滤波时忽略 NaN 值。然后用dem_smooth.grd去计算山体阴影。
2. 调整光照参数：尝试提高太阳高度角-E（如-E60），光照更“正”，阴影更短，噪声不明显。也可以尝试不同的-N参数，如-Nt0.4获得更柔和的效果。

问题3：SHP 文件包含多个多边形（如群岛），裁剪后只保留了其中一个。

原因：gmt convert导出的boundary.xy文件虽然包含多个段（由>分隔），但grdmask默认可能只处理第一个多边形，或者需要特殊参数处理多部件（Multi-part）多边形。
解决：确保你的 GMT 版本支持多部件多边形。一种更可靠的方法是使用 OGR 直接读取 SHP，并让grdmask处理.shp文件本身，而不是转换后的文本。grdmask命令可以直接读取.shp文件。如果问题依旧，可以在 QGIS 中将多部件多边形“分解为单部件”，然后保存为一个新的 SHP 再使用。

问题4：最终图上，SHP 边界外的区域不是透明，而是白色或黑色。

原因：裁剪后的 DEM 外部值不是 NaN，而是 0 或其他数值。或者在绘图时，-R范围设置得比裁剪网格大，GMT 会自动用-N参数指定的背景色（默认白色）填充无数据区域。
解决：
1. 确保按照3.2节的“更优方案”使用-NNaN/1/1和grdclip，得到外部为 NaN 的网格。
2. 在gmt grdimage绘制 DEM 填色图时，使用-N参数。例如gmt grdimage dem_raw_cut.grd -Celevation.cpt -t50 -N，其中的-N会将 NaN 像素设置为完全透明，这是叠加绘图时的关键。

4.2 高级技巧与性能优化

1. 使用gmt grdcut进行快速粗剪：在创建高分辨率掩膜之前，先用grdcut以 SHP 的外包矩形（Bounding Box）快速切出一块稍大的 DEM 子集，能极大减少后续grdmask和grdmath处理的数据量。

# 获取SHP文件的范围 gmt minmax boundary.xy -I0.1 > region.txt # 假设 region.txt 内容为：-R99.8/100.5/25.1/25.8 # 使用该范围裁剪原始DEM gmt grdcut your_dem.tif -Gdem_rough_cut.grd -R$(cat region.txt)

后续所有操作都基于dem_rough_cut.grd进行，速度会快很多。

2. 制作多光源融合的山体阴影：单一光源的山体阴影可能会有明显的方向性阴影死角。可以采用“多光源加权平均”的方法来生成更均衡、细节更丰富的山体阴影，模拟“天空光”效果。

# 生成四个方向的光源 gmt grdgradient dem_raw_cut.grd -Ghill_north.grd -A0 -Nt1 gmt grdgradient dem_raw_cut.grd -Ghill_east.grd -A90 -Nt1 gmt grdgradient dem_raw_cut.grd -Ghill_south.grd -A180 -Nt1 gmt grdgradient dem_raw_cut.grd -Ghill_west.grd -A270 -Nt1 # 合并（平均） gmt grdmath hill_north.grd hill_east.grd ADD hill_south.grd ADD hill_west.grd ADD 4 DIV = hill_multi.grd

使用hill_multi.grd作为底图，地形细节表现力更强。

3. 色彩与光照的精细调节：

CPT 定制：不要局限于内置的dem2。使用gmt makecpt或gmt grd2cpt基于你的 DEM 数据范围生成 CPT。gmt grd2cpt dem_raw_cut.grd -Crelief -Z > my.cpt可以创建与地形匹配的连续色标。
光照强度混合：在grdimage绘制山体阴影时，可以尝试-I+a<factor>参数，其中<factor>是一个加权因子，用于混合朗伯反射和强度调制，有时能获得更自然的效果。例如-I+a0.5。

4. 自动化脚本编写：将上述所有步骤写进一个 Bash 脚本（.sh）或批处理文件（.bat），并参数化关键变量（如输入文件路径、输出文件名、投影参数、颜色表、光照角度等）。这样，下次处理另一个区域时，只需修改几个参数即可运行，极大提升效率，也保证了处理流程的一致性。这是从“会操作”到“专业化”的关键一步。

5. 从成果到发布：流程整合与质量检查

当你走完整个流程，得到满意的final_map.png后，工作并未完全结束。一个专业的流程还包括整合与检查。

流程整合脚本示例：创建一个名为clip_and_hillshade.sh的脚本。

#!/bin/bash # 描述：使用GMT和本地SHP裁剪栅格并生成山体阴影图 # 用法：./clip_and_hillshade.sh <input_dem> <input_shp> <output_prefix> set -e # 遇到错误即停止 INPUT_DEM=$1 INPUT_SHP=$2 OUTPUT_PREFIX=$3 echo “步骤1：转换SHP边界...” gmt convert $INPUT_SHP -fg > ${OUTPUT_PREFIX}_boundary.xy echo “步骤2：创建NaN掩膜...” # 先获取SHP的大致范围，外扩0.1度，用于创建掩膜，提升速度 REGION=$(gmt minmax ${OUTPUT_PREFIX}_boundary.xy -I0.1) gmt grdmask ${OUTPUT_PREFIX}_boundary.xy -G${OUTPUT_PREFIX}_mask.grd $REGION -I+ -NNaN/1/1 -V echo “步骤3：精确裁剪DEM...” gmt grdclip $INPUT_DEM -G${OUTPUT_PREFIX}_dem_cropped.grd -S${OUTPUT_PREFIX}_mask.grd/1/NaN -V echo “步骤4：生成山体阴影...” gmt grdgradient ${OUTPUT_PREFIX}_dem_cropped.grd -G${OUTPUT_PREFIX}_hillshade.grd -A315 -Nt0.7 -V echo “步骤5：创建颜色表...” gmt grd2cpt ${OUTPUT_PREFIX}_dem_cropped.grd -Crelief -Z > ${OUTPUT_PREFIX}_color.cpt echo “步骤6：绘制最终地图...” gmt begin ${OUTPUT_PREFIX}_final_map png gmt basemap -R${OUTPUT_PREFIX}_dem_cropped.grd -JM15c -Baf -BWSen+t“${OUTPUT_PREFIX} 地形图” gmt grdimage ${OUTPUT_PREFIX}_hillshade.grd -I -M gmt grdimage ${OUTPUT_PREFIX}_dem_cropped.grd -C${OUTPUT_PREFIX}_color.cpt -t60 gmt plot ${OUTPUT_PREFIX}_boundary.xy -W0.8p,black gmt colorbar -C${OUTPUT_PREFIX}_color.cpt -Baf+l“高程 (米)” -DJMR+o0.8c/0+w12c/0.4c gmt end echo “处理完成！输出文件：” echo “- 裁剪后DEM: ${OUTPUT_PREFIX}_dem_cropped.grd” echo “- 山体阴影: ${OUTPUT_PREFIX}_hillshade.grd” echo “- 最终地图: ${OUTPUT_PREFIX}_final_map.png”

最终质量检查清单：

边界对齐：放大检查 SHP 矢量边界与地形裁剪边缘是否完全重合，有无错位或锯齿。
数据完整性：裁剪区域内的地形特征是否连续、自然？有无因掩膜错误导致的异常条带或孔洞？
视觉效果：山体阴影的明暗对比是否适中？颜色表是否清晰反映了高程分布？叠加透明度（-t）是否使色彩和阴影融合得当？
地图元素：比例尺、指北针、图例、标题等是否齐全、清晰、美观？坐标轴标注是否合理？
文件管理：中间过程生成的网格文件（如*_mask.grd,*_hillshade.grd）通常很大，在确认最终成果后，可以考虑将其删除或归档，只保留源数据、最终成果图和运行脚本。

通过这样一套从原理到命令，从步骤到脚本，从操作到排错的完整流程，你不仅能够完成一次“裁剪+山体阴影”的任务，更建立起了一套可重复、可优化、可扩展的 GMT 空间数据处理工作流。这套方法同样适用于其他类型的栅格数据（如温度、降水、重力异常等）与矢量数据的裁剪与可视化，只需替换数据源和调整可视化参数即可。

本文还有配套的精品资源，点击获取