刚接触海洋数据处理的时候很多人第一步就卡在怎么把地形水深数据画成一张像样的图上。直接用MATLAB自带的pcolor、contourf画经纬度网格出来的图要么海岸线缺失要么经纬度变形严重要么投影一塌糊涂。后来我换成了m_map工具箱彻底解决了这些问题。m_map是海洋科学和大气科学领域用得最广的MATLAB地图绘制工具包专门用来处理带经纬度坐标的空间数据从地形水深图、海表温度图到航迹叠加图一套代码逻辑可以通吃。这篇博文就围绕用m_map画地形水深图这条主线从数据准备、工具箱配置到完整绘图流程把每个关键步骤和坑点都拆开讲清楚。1. 为什么画水深图要选m_map而不是MATLAB自带函数1.1 MATLAB自带绘图的三个老大难问题很多人一开始都是这么干的拿到lon、lat、depth三个变量直接pcolor(lon,lat,depth)或者contourf(lon,lat,depth)然后发现图出来不对味。问题主要集中在三个方面。第一投影变形。地球是个球面而绘图区域是平面如果直接把经纬度当作平面坐标去画高纬度地区的东西向距离会被严重拉长。比如画北纬60度以上的海域直接用pcolor画出的图看起来一个网格是正方形实际在地球上东西方向距离只是南北方向的一半。这种变形对于定性看图可能影响不大但一旦要叠加航迹、测站、等深线问题就很明显。第二海岸线缺失。MATLAB自带的worldmap函数也能画海岸线但它的海岸线数据分辨率较低放大到局部海域后岸线轮廓非常粗糙很多岛屿细节直接消失。第三坐标轴与图框控制太麻烦。用pcolor画完以后想让经纬度刻度变成度分秒格式或者让图框变成自己想要的形状MATLAB原生功能能实现但步骤繁琐而且不同版本兼容性一般。1.2 m_map到底解决了什么m_map是一套基于MATLAB的地图投影工具箱核心思路是先设定投影和区域范围再把所有空间数据通过投影变换绘制在图上。你不需要自己关心投影公式怎么算只需要告诉它你用哪种投影、想画哪个经纬度范围它会把海岸线、等深线、数据网格全部统一投影到目标坐标系里。和MATLAB自带的Mapping Toolbox相比m_map最大的优势是免费、轻量、社区成熟。Mapping Toolbox虽然功能更强但需要额外购买授权很多课题组并没有安装。m_map体积小解压即用只需要把文件夹加入MATLAB路径。在海洋科学领域大量公开发表的论文附图都是用m_map画的这意味着遇到问题很容易搜索到解决方案这一点在实际工作中非常重要。2. 数据准备地形水深数据从哪来怎么整理2.1 常用的公开水深数据源画水深图的前提是有数据。公开的全球地形水深数据主要有这么几个我按推荐程度排个序。GEBCO 2023 Grid是我最推荐的数据集分辨率15弧秒大约是450米左右全球覆盖陆地高程和海底水深都有。数据从GEBCO官网免费下载输出为netCDF格式包含elevation变量正值代表陆地海拔负值代表海水深度。这个产品更新及时很多海洋研究都用它。ETOPO1也是经典选择有ice surface和bedrock两个版本分辨率1弧分全全球范围文件相对小适合快速出图。缺点是分辨率低一些近岸细节不足。**SRTM15**是海底地形数据里分辨率较高的约15弧秒融合了卫星测高和船测数据比较适合深海区域。不过它在极浅水区的表现不一定比GEBCO好。如果你是做局部海域研究比如某个海湾、某段陆架还可以下载局地高精度数据比如各个国家的海道测量部门发布的DEM产品。不过通用流程还是推荐GEBCO起步后续按需替换。2.2 从netCDF读取并裁剪数据GEBCO下载下来是一个netCDF文件用MATLAB读取非常简单。如果你没有netCDF相关工具箱MATLAB基础环境就自带ncread和ncdisp函数。% 查看netCDF文件的变量结构 ncdisp(GEBCO_2023.nc); % 读取经纬度与水深变量 lon ncread(GEBCO_2023.nc, lon); lat ncread(GEBCO_2023.nc, lat); elev ncread(GEBCO_2023.nc, elevation);注意GEBCO的lon是从-180到180lat是从-90到90elevation是一个二维矩阵第一维对应lon第二维对应lat。读取后往往需要裁剪到你关心的海域范围否则一个全球数据几GBMATLAB处理起来会卡到怀疑人生。裁剪的核心就是做索引筛选% 设定目标区域以一片海域为例 lon_min 2; lon_max 12; lat_min 56; lat_max 60; % 找到对应索引 lon_idx find(lon lon_min lon lon_max); lat_idx find(lat lat_min lat lat_max); % 裁剪 lon_sub lon(lon_idx); lat_sub lat(lat_idx); elev_sub elev(lon_idx, lat_idx); % 转置让矩阵维度变成 (lat, lon)方便后续绘图 elev_sub elev_sub;这里有一个非常容易踩的坑GEBCO的elevation矩阵维度是(lon, lat)很多人在裁剪之后忘了调整维度顺序结果后面画图的时候发现矩阵尺寸对不上。我在代码里直接做了转置统一变成(lat, lon)这样后面无论是pcolor还是contourf参数的传入逻辑都清晰。2.3 经纬度网格与坐标类型的统一水深数据里的经纬度可能是double类型数组也可能是single这本身不影响绘图。但要注意网格的排列方向有的数据lat是从北往南排的也就是降序。这种情况下如果直接丢给m_map画图出来的图上下颠倒海岸线也对不上。判断方法很简单% 判断lat是升序还是降序 if lat_sub(1) lat_sub(end) lat_sub flipud(lat_sub); elev_sub flipud(elev_sub); end另外如果你的水深数据里有NaN值陆地或者数据缺口不处理也能画但会出现颜色空洞。一般处理方式是用shadem或者插值填掉或者直接把NaN设置为一个固定值。比如让陆地全部显示为同一颜色可以在绘图时用contourf的Level设置来控制。3. 用m_map画水深图的完整实操3.1 m_map的安装与路径配置m_map没有安装程序解压即可用。从官网下载最新版本解压后把整个文件夹放到MATLAB的toolbox目录下或者任意你方便管理的目录然后在MATLAB里执行% 将m_map路径添加到MATLAB搜索路径 addpath(genpath(D:\Tools\m_map)); % 保存路径避免下次启动又要重新添加 savepath;如果你不想永久保存路径也可以用pathtool图形界面手动添加。注意genpath会递归添加所有子文件夹m_map里有些示例代码子目录加进去没有副作用懒人可以直接genpath。安装完成后检查是否成功% 检查m_proj函数是否可以调用 which m_proj如果能返回完整路径说明安装成功。如果提示找不到函数多半是路径没加对重新检查一下目录层级。3.2 设置投影与绘图区域m_map把投影方式和区域范围绑定在m_proj函数上格式非常固定m_proj(mercator, long, [lon_min lon_max], lat, [lat_min lat_max]);表示使用墨卡托投影经度范围从2度到12度纬度范围从56度到60度。实际使用时要看工作海域的特征选择投影方式。墨卡托投影Mercator用得最多适合低纬和中纬度海域形状保持好且经纬度网格正交看起来舒服。高纬度地区例如极地附近应该选极方位立体投影Stereographic避免变形严重。等距圆锥投影Lambert Conformal Conic适合中纬度东西跨度大的区域。通用原则是中低纬度海域用mercator极地用stereographic区域不规则时先用mercator试试不满意再换。m_proj设置的区域范围会成为整个图的坐标系基准在这之后所有m_map相关函数m_pcolor,m_contour,m_line都会自动进行投影变换。需要说明的是这个投影设置是针对当前图窗的一幅图只能有一个m_proj设置重复设置会覆盖前一次。3.3 绘制底图与水深填色设置好投影后先画海岸线再画水深数据。海岸线用m_gshhs或者m_coast区别是海岸线数据分辨率。m_gshhs需要额外下载GSHHS全球海岸线数据集但效果好m_coast用的是内置数据胜在方便。实际使用中我默认用m_gshhs的high分辨率近岸细节完整但第一次绘制可能要联网下载数据文件。绘制水深的经典代码figure(Color, w); m_proj(mercator, long, [lon_min lon_max], lat, [lat_min lat_max]); % 先画海岸线 m_gshhs(high, patch, [0.7 0.7 0.7]); hold on; % 画水深填色 m_pcolor(lon_sub, lat_sub, elev_sub); shading flat; % 叠加等深线 [c, h] m_contour(lon_sub, lat_sub, elev_sub, [-500 -1000 -2000 -3000], k); clabel(c, h, fontsize, 8); % 添加colorbar caxis([-3500 0]); % 根据数据范围调整 colorbar; colormap(flipud(jet(256))); % 或者用自定义色带这里m_pcolor和普通的pcolor逻辑一样只是内部做了投影变换。shading flat最好紧跟其后否则网格线会影响观感。如果你希望填色过渡平滑可以用shading interp但注意使用interp后离散网格之间的颜色会插值遇到NaN边界可能出现淡淡的颜色条纹。3.4 等深线参数的选择逻辑等深线不是随便选的它要服务于你的展示目的。比如你关心的是陆架坡折那等深线的间隔就应该在200米、500米、1000米附近加密如果你要看深海平原则要重点画3000米、4000米等深线。我这里给出一个实用的选择策略先用contourc计算数据范围内的等深线值看直方图再决定画哪些层级。简单来说浅水区等深线可以画得密一些深水区等深线可以画得疏一些让图面信息量均衡。% 查看水深数据的分布 histogram(elev_sub(-elev_sub0), 50);根据分布情况再设定合适的等深线数组。比如某个区域水深主要分布在-50到-500米之间那等深线选[-50, -100, -200, -300, -500]就比较合理而不是硬套固定的[-500, -1000]。等深线的标注也要注意别让文字互相重叠如果标注太多可以设置clabel(c, h, manual)手动选择标注位置。3.5 自定义色带让深水浅水一眼区分MATLAB内置的jet色带虽然鲜艳但有一个致命问题浅水区和深水区之间的色差不够直观而且色带本身有彩虹纹。海洋领域比较推荐的是蓝色系渐变色带从浅蓝到深蓝深度越大颜色越深。也可以用cmocean工具箱Name, 2023里的deep、deep_r、matter这类专为海洋数据设计的色带。如果你不想额外装cmocean自己手动构造一个蓝色渐变色带也很简单% 自定义一个从浅蓝到深蓝的色带 mycol [0.9 0.95 1; 0.6 0.8 1; 0.3 0.6 0.9; 0 0.3 0.7; 0 0 0.4; 0 0 0.2]; colormap(mycol);注意colormap要在colorbar之前调用否则色带的颜色映射可能不对。另外设置色带范围用caxis([cmin cmax])在新版MATLAB中部分场景已被clim取代建议写clim([cmin cmax])兼容性更好。如果想让地形和水深在一张图里同时显示可以做一个正负值分离的双色色带比如陆地用棕色系海底用蓝色系。这样图的信息量更大但代码要多几行后面在进阶部分展开。3.6 经纬度网格、标注与图框美化m_map可视化非常容易出效果的另一个原因是m_grid一行代码就能生成符合学术规范的经纬度网格与边框标注。m_grid(box, fancy, tickdir, out, ... fontsize, 10, xtick, 2:2:12, ytick, 56:1:60, ... xlabeldir, end, ylabeldir, end);几个参数说明一下。box设置边框样式fancy会带一个自然的刻度副框on代表标准方框。tickdir控制刻度线方向out是向外学术图常用。xtick和ytick手动指定刻度位置避免默认密度太密或者太疏。xlabeldir设置为end可以让经度标注只出现在图的右侧末端更清爽。需要提醒的是m_grid必须在所有绘图元素绘制完成后调用否则网格线会覆盖在数据之上或者被数据覆盖。正确顺序是先画数据、等深线、航迹最后调用m_grid加框加刻度。3.7 导出高清图的几个细节论文投稿时图的分辨率有硬性要求。MATLAB导出图片我最常用的是exportgraphics函数它比print更稳定对中文字体和fig大小控制友好。exportgraphics(gcf, bathymetry_map.png, Resolution, 600);如果要导出矢量图投稿经常需要推荐导出为PDF或EPSexportgraphics(gcf, bathymetry_map.pdf, ContentType, vector);导出前把图形窗口尺寸设置好避免输出尺寸不对。一般设置为横向构图比如set(gcf, Position, [100 100 1000 700]);注意exportgraphics在MATLAB R2020a及以后版本支持如果还在用老版本退而求其次用print(gcf, -dpng, -r600, bathymetry_map.png)。4. 常见报错与排查技巧实录4.1 m_proj投影设置报错或图形空白最常见的问题是调用m_map系列函数时提示Projection not initialized。原因非常直接你还没有调用m_proj就想画图。m_map所有绘制函数都依赖当前的投影信息必须先设置投影再调用其他m开头函数。另一种情况是已经调用了m_proj但图形窗口里什么也没有。这时候先检查两点第一m_proj指定的经纬度范围里有没有实际数据第二数据范围是否与投影区域重叠。比如你下载的是全球数据但目标海域在大西洋东侧而lon变量范围是-180到180你设置的经度是2到12这些经度在数据中存在那就是绘图函数的问题。反过来如果你的数据lon范围是0到360而m_proj设置的是-180到180那所有数据点都落在投影区域之外图自然空白。解决方案是统一经纬度体系。如果数据是0到360可以用wrapTo180函数MATLAB Mapping Toolbox提供或者手动转换lon(lon 180) lon(lon 180) - 360;4.2 数据维度不匹配的提示m_pcolor(lon_sub, lat_sub, elev_sub)报维度错误几乎都是因为elev_sub的尺寸和lon、lat的网格不匹配。pcolor要求elev_sub是length(lat_sub) x length(lon_sub)的二维矩阵少一维都不行。我习惯在绘图前加一行断言检查方便快速定位问题assert(size(elev_sub, 1) length(lat_sub), 行数必须等于lat长度); assert(size(elev_sub, 2) length(lon_sub), 列数必须等于lon长度);还有个容易忽略的点如果lon_sub、lat_sub是从netCDF直接读出来的原始数组它们是列向量还是行向量都可能但m_pcolor不会帮你自动转置。统一在数据准备阶段就把维度理清楚形成固定套路比每次画图前临时去试效率高得多。4.3 海岸线绘制太慢或数据缺失m_gshhs(high)第一次使用时会尝试下载GSHHS数据库如果网络不稳定可能卡住或者报错。解决方法是手动从GSHHS官网下载GSHHS数据文件放到m_map的gshhs目录下或者改用m_coast绘图。m_coast的优点是快内置了基础海岸线数据不需要额外下载。缺点是分辨率不够高放大到局部海域后岸线会显得棱角分明。如果你以近岸细节为主要展示内容还是建议用gshhs的高分辨率版本。一个实用技巧先m_coast快速预览区域范围和数据叠加效果确认无误后再切换到m_gshhs(high)出终稿节省等待时间。4.4 colorbar显示异常和坐标刻度冲突有时你画完图colorbar显示出来了但颜色范围和数据的实际范围对不上或者colorbar标签堆叠混乱。这通常是因为caxis/clim设置和colormap调用顺序问题。记住一个原则先设定colormap再设colorbar。同时如果你用m_grid设置了边框colorbar默认出现在右侧如果改用了colorbar(southoutside)之类的位置注意别和m_grid的纬度标注重叠。还有一个隐蔽问题m_grid画出来的经纬度刻度文字如果用了默认字体在Linux下可能显示为方框。解决方法是统一设置字体set(gca, FontName, Helvetica);或者中文环境里设置SimHei。4.5 负值水深和陆地掩膜问题如果你的数据里陆地高程是正数你在画水深图时直接把所有值都pcolor出来陆地区域会被画成浅蓝看起来像浅滩不符合习惯。处理办法是给陆地一个单独的掩膜让它显示为灰色或者不显示。% 把陆地高程设为NaN这样陆地区域显示为透明 elev_plot elev_sub; elev_plot(elev_sub 0) NaN; % 绘图 m_pcolor(lon_sub, lat_sub, elev_plot);不过设成NaN后陆地区域会露出底下的白色和海岸线填充的灰色可能不一致。更好的办法是先把海岸线填充成陆地色再画水深数据并且水深数据中陆地部分设为NaN。这样水层颜色和陆地颜色完全分离。m_gshhs(high, patch, [0.7 0.7 0.7]); hold on; elev_plot(elev_sub 0) NaN; m_pcolor(lon_sub, lat_sub, elev_plot);由于m_pcolor的数据会在陆地填充之后绘制陆地从视觉上被覆盖在灰色面之上。如果你希望陆地色不透明则需要把透明度和绘制顺序配合好。5. 进阶让水深图不只是一张图5.1 在地形图上叠加测站与航迹用m_map画地形底图最重要的实际场景是作为其他海洋数据的底图。比如你在研究海域布设了CTD测站或者船载走航观测的航迹都可以叠加在这张水深图上。% 假设有测站经纬度 sta_lon [3.2 4.1 5.5 7.8]; sta_lat [57.1 57.8 58.2 59.0]; % 叠加测站符号 m_plot(sta_lon, sta_lat, r^, MarkerSize, 8, LineWidth, 1.2); % 叠加航迹线 track_lon 2.5:0.1:10; track_lat 56.5 0.2 * sin(0.3 * (1:length(track_lon))); m_plot(track_lon, track_lat, b-, LineWidth, 1.5); % 标注站点名称 for i 1:length(sta_lon) m_text(sta_lon(i)0.1, sta_lat(i)0.1, sprintf(C%d, i), FontSize, 8); end这里m_plot和普通plot的用法几乎一致区别是它会做投影变换。测站符号和航迹线会随着投影正确变形这是m_map相比手动plot经纬度数据的绝对优势。5.2 做一版陆地水陆分离的双色图如果你想一张图同时看清陆地高程和海底地形双色色带是很好的选择。核心逻辑是把高程和水深分成两个数据集分别用两个colormap绘制。% 复制数据 elev_land elev_sub; elev_land(elev_sub 0) NaN; elev_ocean elev_sub; elev_ocean(elev_sub 0) NaN; % 先画陆地 m_pcolor(lon_sub, lat_sub, elev_land); colormap(gca, [0.8 0.7 0.6; 0.6 0.5 0.3]); % 简化示意 hold on; % 再画海洋 m_pcolor(lon_sub, lat_sub, elev_ocean); colormap(gca, [0.9 0.95 1; 0 0.2 0.5]); % 简化示意用两个colormap同图时需要技巧因为MATLAB的figure只能有一个colormap一般做法是用colormap结合数据范围的偏移把陆地和水深映射到同一颜色轴的不同区段来实现。这个过程相对复杂如果只是做初步看图的图件可以用透明度叠加近似实现。不过我认为对这个需求更高效的方案是用两个axes叠加或者在同一个axes里把色带合并。这里不展开太多但建议优先用透明度叠加m_pcolor(lon_sub, lat_sub, elev_ocean); shading flat; alpha(0.9); hold on; m_pcolor(lon_sub, lat_sub, elev_land); shading flat;因为透明度的关系陆地和水深会自然融合视觉上够用代码简单适合快速中间图。5.3 把图画成动图或者多子图水深图经常要配合时间序列数据用。比如你要画连续多天的海流模拟结果用同一张水深图当底图循环输出多帧图片最后合成GIF。% 预先设置投影与底图 figure; m_proj(mercator, long, [lon_min lon_max], lat, [lat_min lat_max]); m_gshhs(high, patch, [0.7 0.7 0.7]); hold on; % 循环绘制 for t 1:10 cla; % 清空当前axes % 重新绘制水深作为背景 m_pcolor(lon_sub, lat_sub, elev_sub); shading flat; hold on; % 在这里绘制动态数据例如流速场 m_quiver(lon_sub(1:5:end), lat_sub(1:5:end), u_t(:,:,t), v_t(:,:,t)); m_grid(...); exportgraphics(gcf, sprintf(frame_%02d.png, t), Resolution, 150); endcla清空的是当前axes的内容但m_proj的投影设置是否保留需要验证。稳妥的做法是在循环外只设置一次投影循环内不要再次调用m_proj否则图像位置可能会跳动。5.4 图例、比例尺与指北针论文图件里比例尺和指北针是加分项。m_map提供了比例尺函数m_scale可以在图内添加一个地理比例尺。m_scale(location, southwest, width, 20, fontsize, 8);width表示比例尺代表的实际距离单位是kmlocation控制位置。指北针可以用m_northarrow部分版本有或者直接m_text手动写一个N加箭头。如果m_map版本没有现成指北针也可以用MATLAB的annotation函数叠加一个简单的箭头。annotations, arrow, text;比例尺在局部海域图上特别重要能让读者直观感知尺度。别嫌麻烦加上之后图件的专业度立刻上一个档次。6. 实操心得与踩坑总结6.1 画水深图的黄金三步工作流经过大量实践我整理了一个适合绝大多数海域的工作流照着走基本不会出错。第一步数据准备。用ncdisp查看netCDF结构确认维度顺序和变量名按目标海域裁剪处理NaN和陆地正值统一lon/lat方向。第二步投影与底图。根据纬度选投影方式中低纬度mercator高纬度stereographic调用m_proj设置区域绘制海岸线并填充陆地色绘制水深填色设置色带范围和colormap。第三步完善与导出。调用m_grid美化边框刻度添加等深线、测站、航迹等要素添加colorbar、比例尺设置图形尺寸后导出高清图。这三步看起来简单但每一步都有细节。最怕的是流程颠倒先画数据再设置投影、先colorbar再colormap这种顺序问题改起来费时。6.2 控制图形变量别在内存里裸奔GEBCO全球数据量不小如果你一次读入全球分辨率15弧秒的数据内存占用可能达到几百MB甚至更多。实际绘图前先裁剪到目标区域再传进绘图函数这是最省心的内存管理方式。同时绘图时尽量用elev_sub或者elev_plot这样的局部变量不要反复把整个global数据传来传去。MATLAB本身有写时复制机制但大数据在循环里处理时提前预分配和裁剪能减少内存碎片减少卡顿。6.3 配色好图就成功了一半很多同学把重心全放在代码和图层上忽视了配色。实际上水深图好不好看80%取决于色带选得对不对。深海蓝、浅海青、陆地棕灰是比较经典的三段式配色。用cmocean工具箱的话直接调用cmocean(deep)就是专门为水深设计的色带。如果一定要用jet建议至少把jet翻转一下让深水区是深色而不是红色否则视觉上深水区像高温区非常容易误导读者。6.4 做任何图之前先画一张快速预览图实际工作中我几乎不会一开始就用高分辨率海岸线、精致色带和所有标注画终稿。那样如果数据有问题光等gshhs下载和渲染就浪费不少时间。我的习惯是先用m_coast配合pcolor快速出一张预览图检查数据范围、投影方向、水深分布是否符合预期确认无误后再换高分辨率出正式图。这个习惯帮我省了很多时间。水深处处理最折磨人的就是看着好像哪里不对但说不上来是什么问题一张快速预览图往往能立刻暴露问题。6.5 最后分享一个填充海岸线的小细节m_gshhs(patch, color)填充海岸线时patch里传入的颜色会被用于所有陆地多边形但如果你的数据范围里包含大的湖泊这些湖泊也会被当作陆地填充。这在某些局部海域图里可能是堆积的大量湖泊色块很碍眼。解决方法是使用GSHHS数据的l湖泊属性区分或者在m_gshhs调用时传入一个更细心构造的patch结构对land和lake分别指定颜色。m_gshhs(high, patch, struct(land, [0.7 0.7 0.7], lake, [0.8 0.9 1]));这样陆地和湖泊的填充色分开图面更自然。类似地你还可以把river的属性标出来画河流线条。不过这一步按需使用不要为了炫技反而把图弄得太复杂。