做中尺度海洋过程可视化的时候我经常需要把海平面异常SLA和地转流速叠加到一起做成一段能直接看出演变过程的动图。这个需求听起来简单真正做起来却有不少门道数据源怎么选、单位怎么统一、地转流怎么算、m_map工具箱怎么配、最后动图怎么输出才能不闪烁。这篇文章就把我实际操作下来的整套流程整理出来从NetCDF文件读取到最终GIF/MP4所有步骤都能直接复现。特别适合初次接触高度计数器、想快速出图的人也适合那些已经能画静态图、但还没跨出“动图”这一步的朋友。1. 项目背景与整体思路拆解1.1 为什么非要做成动图卫星高度计提供的海面高度是海洋中尺度过程研究的核心数据。所谓SLA是瞬时海面动力地形相对多年平均海面的偏差它把大尺度背景场去掉之后留下了涡旋、流轴摆动、锋面生消这些中尺度信号。单独看某一时刻的SLA填色图只能看到“这一刻”哪里有正异常、哪里有负异常可是涡旋怎么移动、强度怎么变化、流场怎么响应静态图根本讲不清楚。动图的意义就在这里把时间维展开人眼对运动最敏感一张动画几十帧看下来涡旋的旋转方向、传播路径、能量耗散过程一目了然。叠加u分量地转流矢量之后还能同时看到高度场和流场之间的配置关系比如正SLA中心对应顺时针流场负SLA中心对应逆时针流场。这种“高度场流场叠加动画”在论文补充材料、学术汇报、课题评审里都是非常直观的展示形式。1.2 技术路线拆解从原始数据到可发布的动画我把整条流程拆成了五个环节数据获取、预处理、物理量计算、单帧绘图、动图合成。网上很多教程把这五步混在一起讲读者看完一头雾水一旦某一步报错也不知道是哪里出的问题。拆开之后每一步都独立排查方便得多。数据获取环节解决的是“从哪里拿SLA和流速数据、怎么拿到”。预处理环节解决的是“把NetCDF里的数组变成我们需要的经纬度、时间、数值都正确的矩阵”。物理量计算是整条链路里最容易被忽略的一步SLA换算成地转流速u需要做梯度运算单位、纬向距离、科氏参数这些细节错一个后面图全是废的。单帧绘图解决的是“静态图怎么画得既专业又好看”动图合成则是把几十帧静态图串起来输出成GIF或MP4。1.3 为什么是MATLAB而不是别的工具我平时也经常用Python处理海洋数据但单独做这个“SLAu动图”任务我依然首选MATLAB。原因有三一是m_map工具箱在海洋地球物理圈子里非常成熟投影转换、海岸线、等值线都封装好了省掉自己写投影算法的麻烦二是MATLAB的矩阵操作对netCDF数组特别顺手切片、索引、循环都直观三是getframe配合VideoWriter可以一行代码把图形窗口变成视频帧这个工作流我用得很熟。当然Python也能做但要么依赖Basemap这类已经停止维护的工具要么得花时间配Cartopy环境中文化字体对只想快速出图的人来说门槛略高。本文所有代码都以MATLAB为基准版本不限我用的是日常更新的某个近几代版本较早版本只要支持ncread和VideoWriter都能跑通。2. 数据获取与预处理把NetCDF整理成干净场2.1 高度计数据里到底有哪些变量AVISO这条数据线的产品目前统一通过数据分发平台获取注册后按指引下载即可。我常用的产品是融合多颗卫星数据的网格化延迟处理产品时间分辨率为1天空间分辨率约0.25度。NetCDF文件的变量不多但个个都有用变量名含义常用单位sla海平面异常米adt绝对动力地形米ugos地转流纬向分量米/秒vgos地转流经向分量米/秒longitude经度度latitude纬度度time时间天相对参考日期这里需要特别说明如果你只想要SLA和u完全可以只读sla、ugos、vgos三个变量。不过很多老版本产品并不直接提供ugos和vgos或者你手里的数据已经经过了某种自定义处理这时候就需要自己从SLA场推导地转流。我建议两条腿走路原产品有流速就直接用没有就自己算。后面的第3节会专门讲怎么算。2.2 下载规划与网格选择下载前先想清楚三件事时间范围、空间范围、网格版本。时间范围很好选取决于你想展示的物理过程持续多久。比如一个中尺度涡旋的生命周期通常是一个月到几个月那就下载对应长度的日数据。空间范围建议圈一个适中的矩形区域不要贪大因为地图投影比例、矢量密度、动图文件大小都跟区域范围直接相关。网格版本值得多说一句。AVISO系列有一个“延迟处理”版本和一个“近实时”版本前者精度更高、校准更好适合做研究后处理后者时效性强但质量稍逊做过程分析我从不拿近实时版本出图。另一个容易忽略的点是经纬度网格的定义有些产品lon范围是0到360有些是-180到180有些网格中心在0度起始有些在0.125度起始。这些元数据信息在NetCDF文件的global attributes里都有说明下载后第一件事永远是用ncdisp把文件头看一遍。2.3 MATLAB读取和裁剪MATLAB读取NetCDF数据非常简单。先用ncdisp查看变量名和维度再用ncread读取。% 查看文件信息 ncdisp(dt_global_merged_sla.nc); % 读取经纬度和时间 lon ncread(dt_global_merged_sla.nc, longitude); lat ncread(dt_global_merged_sla.nc, latitude); time ncread(dt_global_merged_sla.nc, time); % 读取SLA和地转流分量注意维度顺序 % 不同文件维度可能是(lon, lat, time)也可能是(time, lat, lon) % 一切以ncdisp显示为准 sla ncread(dt_global_merged_sla.nc, sla); ugos ncread(dt_global_merged_sla.nc, ugos); vgos ncread(dt_global_merged_sla.nc, vgos);一个最常见的坑就在这里不同来源的NetCDF文件维度顺序五花八门。有的文件第一维是lon有的第一维是time。我习惯在读取之后立刻统一成(lon, lat, time)顺序全部用permute处理后面所有代码都不再关心原始维度。裁剪则用逻辑索引实现% 设定目标海区范围 min_lon 100; max_lon 140; min_lat 10; max_lat 40; % 经纬度筛选索引 lon_idx find(lon min_lon lon max_lon); lat_idx find(lat min_lat lat max_lat); % 裁剪数据 sla_sub sla(lon_idx, lat_idx, :); ugos_sub ugos(lon_idx, lat_idx, :); vgos_sub vgos(lon_idx, lat_idx, :); lon_sub lon(lon_idx); lat_sub lat(lat_idx);如果原始lon是0到360范围需要先转成-180到180再裁剪lon mod(lon, 360); lon(lon 180) lon(lon 180) - 360; [lon, sort_idx] sort(lon); sla sla(sort_idx, :, :);这一步很容易被忽略坐标没转一致的话后面m_proj里会直接报错或者画出奇怪的位置偏移。2.4 单位、NaN与坐标的几处细节预处理阶段有几个细节我吃过亏每一个都足以让最终图完全不能用。单位问题绝大多数新版高度计产品SLA单位是米但早期版本或某些再分析产品可能是厘米。判断方法看ncdisp里的units属性如果是cm务必除以100转成米。地转流速单位一般是米/秒但有些产品导出时是厘米/秒也需要仔细核对。这种低级错误一旦混进去画出来的色标量级会谬之千里。NaN问题陆地格点的SLA值是NaN近岸、冰区的数据也可能有空洞。如果直接用带NaN的矩阵去算梯度gradient会把NaN向周围扩散产生大片虚假的异常值。保险做法是先做一个陆点掩膜算完梯度之后再掩掉边缘区域。时间问题NetCDF里的time通常是相对某个参考日期的天数比如“days since 1950-01-01”。转换成实际日期用以下代码time_units ncreadatt(dt_global_merged_sla.nc, time, units); ref_str time_units(12:end); % 提取参考日期字符串 ref_date datetime(ref_str, InputFormat, yyyy-MM-dd); time_abs ref_date days(time);转换出来的时间做两件事一是排序确保时间维度是递增的二是给后续每帧图打时间标签动画里没有时间戳会显得很不专业。3. 核心物理量SLA与地转流u的计算3.1 从海面高度异常到地转流的物理关系中尺度过程的时间尺度和空间尺度都远大于惯性周期运动近似满足地转平衡。所谓地转平衡就是水平压强梯度力与科氏力相互抵消。海表高度异常直接对应海表压强异常因此SLA的水平梯度就决定了地转流异常。公式如下u -(g/f) * ∂η/∂y v (g/f) * ∂η/∂x其中η是SLAg是重力加速度约9.81 m/s²f是科氏参数。f随纬度变化f 2Ω sin(φ)Ω是地球自转角速度7.2921e-5 rad/sφ是纬度。在中纬度f量级在1e-4左右在赤道附近f趋近于0地转关系不再适用。所以u不是直接从SLA数值就能看出来的必须做空间梯度。这也解释了为什么需要精确的距离换算。3.2 用gradient数值微分MATLAB的gradient函数用于数值差分。关键问题是如何把经纬度网格转换到实际距离。纬度方向每1度约等于111公里这个换算相对固定经度方向每1度对应的距离则随纬度收缩是111公里乘以cos(φ)。下面这段代码是我反复验证过的直接可用g 9.81; Omega 7.2921e-5; % 目标海区的纬度向量 % 注意gradient的第一维对应行方向第二维对应列方向 % 行方向沿纬度列方向沿经度 dy_m 111e3 * ones(size(lat_sub)); % 每纬度间距单位米 dx_m 111e3 * cosd(lat_sub); % 每经度间距随纬度变化单位米 % 构造二维网格间距矩阵 % 每个格点的纬度对应不同的经度方向距离 DY repmat(dy_m(:), 1, length(lon_sub)); DX repmat(dx_m(:), 1, length(lon_sub)); % 对单帧SLA计算梯度得到的是(dsla/dy, dsla/dx) [dsla_dy, dsla_dx] gradient(squeeze(sla_sub(:,:,t)), DY, DX); % 科氏参数与纬度有关 f 2 * Omega * sind(lat_sub); F repmat(f(:), 1, length(lon_sub)); % 地转流分量 u -g ./ F .* dsla_dy; v g ./ F .* dsla_dx;需要注意gradient的两个输出第一个是沿第一维y方向的偏导第二个是沿第二维x方向的偏导。很多人在这一步搞反画出来的矢量场整体旋转了90度物理上完全错误。另一个关键细节DX必须用与纬度相关的矢量。如果全程用固定距离高纬度地区的经向距离被高估u分量会系统性偏小。虽然对于0.25度网格的中纬度小区域影响不算致命但在北极圈附近或跨大纬度范围时误差非常大。3.3 直接读取原生产品ugos作为替代方案如果不想自己算梯度AVISO产品里直接有ugos和vgos变量这是官方用地转关系处理过的流速场我在实际对比中发现它和用SLA梯度算出来的结果高度一致差异主要来自边界处理和网格差分方式的细节。直接用ugos/ugos省事且少一个出错环节。但自己算梯度也有不可替代的价值一是可以完全控制平滑程度二是在数据产品里缺少流速变量时仍然能出图。实操中我更推荐一个混合策略用ugos/vgos作为主数据画矢量用SLA作为填色背景。这样高度场和流场都来自官方处理视觉上更干净也不用担心gradient引入的NaN扩散问题。平滑是一个可选项。直接数值差分会放大SLA场中的高频噪声导致矢量箭头方向凌乱动图看起来像一锅乱炖。我通常会对SLA场做一个简单的空间滤波比如用3x3窗口的中值滤波或者直接在gradient之前对SLA做一次轻度平滑。滤波强度要克制否则涡旋的真实信号也会被抹掉。比较好的经验是先用原始场看一眼如果矢量箭头方向连续、不存在明显棋盘状跳跃就不要加平滑只有箭头乱得没法看时才引入平滑。4. 动态绘图实现从单帧图到GIF/MP44.1 m_map工具箱准备m_map是UTM类地图投影工具箱在MATLAB界经久不衰。它其实不是MathWorks官方的而是一套免费开源的工具箱。准备过程三步下载m_map压缩包解压到某个目录在MATLAB里addpath并保存路径。验证是否安装成功运行m_proj(mercator, lon, [100 140], lat, [10 40])如果出现地图坐标轴说明已经可用。m_map提供的投影方式很多做SLA图我一般用mercator墨卡托投影适合中低纬度的小区域。如果研究区域纬度跨度大可以考虑mollweide或等距圆柱投影但墨卡托在日常海洋可视化里最稳妥经纬网正交涡旋形状不容易被夸张扭曲。4.2 单帧图布局SLA填色u矢量叠加先说单帧图怎么画它是动图的地基。我的布局是上层用m_pcolor画SLA填色下层叠加等值线和流速矢量箭头底图加海岸线和经纬网。核心代码如下% 选定投影和区域 m_proj(mercator, lon, [min(lon_sub) max(lon_sub)], ... lat, [min(lat_sub) max(lat_sub)]); % SLA填色 m_pcolor(lon_sub, lat_sub, sla_frame); shading interp; % 加SLA等值线 [C, h_contour] m_contour(lon_sub, lat_sub, sla_frame, -0.5:0.1:0.5); clabel(C, h_contour, FontSize, 8); % 加海岸线 m_coast(patch, [0.7 0.7 0.7], edgecolor, k); % 叠加流速矢量 m_quiver(lon_q, lat_q, u_q, v_q, color, k, AutoScale, on);几个要点m_pcolor和普通pcolor不同它会把经纬度坐标转换到地图投影坐标所以传入的lon/lat是原始经纬度即可。shading interp是必须的否则填色格点锯齿明显。等值线条数不要太多SLA范围通常只有±0.5米左右间隔0.1米是一个适中的密度。矢量稠密度需要控制。0.25度网格全画上去箭头会糊成一片。我的做法是每隔2个格点画一个箭头也就是4倍网格间距。矢量颜色建议用黑色或深灰色在蓝红色SLA背景上辨识度最高。如果需要箭头颜色表达速度大小也可以用速度幅值着色的彩色箭头但动图中信息过载反而影响观看我一般保持单色。时间戳用text或m_text标注在图右上角m_text(0.98, 1.02, datestr(time_abs(t), yyyy-mm-dd), ... Units, normalized, HorizontalAlignment, right, FontSize, 12);色标范围一定要固定。SLA数据的季节变化导致不同时刻的min/max不同如果每帧都用自动色标范围动图播放时色标不断变化视觉上就像曝光在闪完全无法看。我通常取整个时间序列的全局5%和95%分位数作为色标上下限sla_all sla_sub(~isnan(sla_sub)); caxis_lim [prctile(sla_all, 5), prctile(sla_all, 95)]; caxis(caxis_lim); colorbar;4.3 循环成动画的两种出口单帧图画好之后动图就是一个循环。两种输出方式我都在用根据场景选择MP4用于汇报放映GIF用于插入PPT、网页或微信聊天记录。MP4方式用的是VideoWriterv VideoWriter(sla_uv.mp4, MPEG-4); v.FrameRate 8; % 8帧/秒涡旋演变速度适中 v.Quality 90; open(v); figure(Position, [100 100 900 700]); for t 1:nt % 每帧重新绘图 draw_sla_frame(t); % 封装好的绘图函数 frame getframe(gcf); writeVideo(v, frame); end close(v);GIF方式用imwrite追加写figure(Position, [100 100 900 700]); delay 0.15; % 每帧间隔0.15秒 for t 1:nt draw_sla_frame(t); [imind, cm] rgb2ind(frame2im(getframe(gcf)), 256); if t 1 imwrite(imind, cm, sla_uv.gif, Loopcount, inf, DelayTime, delay); else imwrite(imind, cm, sla_uv.gif, WriteMode, append, DelayTime, delay); end endGIF的注意点帧数太多会导致文件巨大。我做过一个60帧、900像素宽的GIF文件能到几十MB加载时明显卡顿。如果需要长期存档或插入论文补充材料MP4是更好的选择如果必须用GIF可以降低figure尺寸到600像素宽或者减少颜色索引到128色。4.4 让动图不闪烁、不变形的几个细节动图看似只是循环实际上有很多隐蔽的坑。首当其冲的是图形句柄复用问题。如果每帧都重新调用m_pcolor创建新对象旧对象不删除内存疯涨图也会越来越卡。正确做法是在第一帧创建句柄后续帧只更新CData% 第一帧创建 h_p m_pcolor(lon_sub, lat_sub, sla_frame); shading interp; % 后续帧只更新 set(h_p, CData, sla_frame);但m_pcolor在投影坐标下更新CData有一个局限如果坐标网格不变只换CData完全没问题。如果区域固定就可以用这个方式。矢量箭头同理用set更新UData、VData、XData、YData或者干脆每帧delete再重建只要对象数量不大性能影响可接受。第二是坐标轴范围必须锁死。ax gca; set(ax, XLim, xlim_lim, YLim, ylim_lim); 这句话加到每帧绘图函数的末尾防止个别帧因为NaN范围变化或投影转换异常导致坐标轴抖动。第三是图形窗口的位置和大小保持一致。getframe捕获的是当前figure的可视区域如果窗口被移动或遮挡帧画面会产生跳动。循环期间不要手动拖窗口或者设置figure的Position后就不再变动。第四是渲染性能。SLA填色图如果格点数很大每帧m_pcolor加m_contour会偏慢。我的经验是把矢量抽取步长加大到4个格点等值线只在每3帧更新一次。对日数据、60帧动图来说渲染速度完全够用。5. 常见问题与排查实录5.1 问题速查表我把实际运行中遇到的典型问题整理成一张速查表代码报错时直接对照能省很多时间。问题现象可能原因处理方法sla读出来全是NaN陆地点掩膜未处理或变量名引用错误ncdisp检查变量名和维度读取后用nansum做全局测试经度范围不对图上海区位置偏移lon是0到360格式未转换成-180到180mod加sort重排转换后再裁剪矢量箭头方向整体旋转了90度gradient的两个输出含义用反了明确第一输出对应dy行第二输出对应dx列矢量箭头乱飞、没有规律SLA噪声被数值差分放大先对SLA做3x3中值滤波或改用产品自带ugos/vgos动图播放时色标闪烁caxis范围没有固定用全局分位数统一色标范围循环内禁止用autoGIF文件过大、播放卡顿帧数多、尺寸大、颜色多降低figure宽度到600颜色索引128色或改用MP4m_quiver绘制报错m_map版本太老或投影未初始化确保在m_quiver前调用过m_proj升级m_map到最新版动图画到一半MATLAB内存耗尽图形句柄未复用每帧新建对象第一帧创建句柄后续只set更新CData日期标注位置漂移text用了data坐标经度纬度范围变化改用normalized坐标定位日期标签SLA与u符号关系看起来不对科氏参数f在程序里用了负值或单位没统一确认f符号用正数u/v单位都是米/秒5.2 我踩过的坑希望你别再踩第一个印象深刻的教训来自经度转换。当时下载的数据lon范围写的是0.125到359.875我在裁剪时直接写了min_lon 100max_lon 140却忘了目标海区有一部分在180度以东、对应lon是负数。结果画出来的区域包含了完全错误的海域SLA填色图跟真实海区对不上。排查了整整一个下午才发现是经度格式问题。第二个教训是单位。某次换了一个再分析产品从文件头看units标注的是cm但我惯性思维以为是米直接画图出来的SLA量级是±5米色标被拉得无比宽涡旋信号全被压成一片蓝色。后来认真ncdisp才发现单位是厘米。从那以后我每拿到新数据第一件事就是打印units绝不跳过。第三个教训是动图闪烁。我第一次做60帧GIF时色标范围用的caxis auto结果播放时色标每帧都在跳动涡旋区域的颜色也在不断变亮变暗组里同学看了说晃眼睛。改成全局分位数固定色标后整段动画才真正清晰起来。这个坑不踩一次很难意识到因为静态单帧图自动色标完全看不出来问题。第四个教训是gradient的方向问题。我第一次算u分量时直接把gradient第一个输出当成了∂SLA/∂x结果把所有箭头转了90度看起来像流场整体穿越等值线物理上非常离谱。后来在图上叠了一条等值线检查箭头是否沿等势线流动才发现问题。强烈建议画完后做一个物理一致性检查正SLA中心周边流场应当是顺时针旋转也就是在中心北侧u向东、南侧u向西方向反了就是梯度维度搞错了。最后一个小经验绘图函数一定要封装。把“读数据、裁剪、算流速、画单帧、写动图”各封装成独立函数后续换数据源或换区域时只需要改参数不用翻代码。我现在有一个make_sla_movie函数输入文件路径、经纬度范围、时间范围直接输出动图十几秒就能处理出一个新区域的动画效率非常高。这个流程本身还有扩展空间。比如把变量从SLA换成绝对动力地形ADT就能从“异常场视角”切换成“全流场视角”把u/v按模长着色能直观展示流速能量分布如果要展示涡旋追踪结果还可以把涡旋中心轨迹散点叠加在动画上。只要把数据读取、物理量计算、绘图这三层解耦清楚这些扩展都只是改函数参数的事。