做MATLAB可视化做到一定阶段“球面投影”是绕不开的一个坎。不管是地球物理数据展示、雷达覆盖范围模拟还是全景图像球面映射最终都要把三维球面上的信息“拍扁”到二维平面上。上一篇讲了球面投影的基本概念和坐标系变换的数学基础这篇直接把代码实现、可视化技巧和踩坑记录一起放出来。内容适合正在研究MATLAB球面投影、想在项目里直接落地的人。我会先对比几种常见投影的数学特点再给出可以直接运行的投影函数最后把视觉美化、动画和常见问题一并整理。不绕弯子直接开始。1. 球面投影的数学基础先把投影类型和公式吃透1.1 为什么要先算经纬度转笛卡尔坐标很多人在做球面投影时第一反应是直接套公式把经纬度代入某个投影方程结果经常在极点附近出错。原因很简单大多数投影公式本身就是在笛卡尔坐标系里推导出来的经纬度只是球面上的参数化表达直接代入会让数学关系绕路。正确的做法是先统一到笛卡尔坐标。球面上任意一点给定经度lon和纬度lat可以写成x R * cos(lat) .* cos(lon); y R * cos(lat) .* sin(lon); z R * sin(lat);这里的R是球半径。经纬度建议用弧度制MATLAB里直接调deg2rad转换不要手写pi/180减少出错机会。有了(x, y, z)投影的本质就清楚了找一个映射把三维坐标降到二维。正交投影直接丢一个轴球极投影做透视除法等距方位投影转极坐标——每个投影其实都是在处理这组笛卡尔坐标。先算笛卡尔坐标还有一个好处投影中心的旋转可以直接用矩阵乘法完成比在经纬度里反复做球面三角运算快得多也稳得多。1.2 四类常见投影的对比与选型方位投影Azimuthal Projection是球面投影里最常用的一类它的特点是所有点都投影到一个与球面相切或平行的平面上并且投影中心对称。做MATLAB球面投影时我接触最多的是下面四类投影类型基本原理核心性质投影范围典型场景正交投影 Orthographic从无穷远处平行投影视觉效果真实像卫星视角只能看到半个球地球外观模拟、气象展示球极投影 Stereographic从球面一点投影到对面切平面保角圆投影后仍是圆除投影点外全平面极区海图、全景图像拼接等距方位投影 Equidistant Azimuthal从球心投影到切平面到中心点的距离真实保持全局可用地震台网覆盖、飞行路径分析心射投影 Gnomonic从球心投影到切平面大圆投影为直线小于半球超出部分发散航海大圆航线规划选型时主要看需求。如果是做“从太空看地球”的效果正交投影最合适视觉上最直观如果要做极地附近的等角航线分析球极投影的保角性最好用如果要展示点与点之间的真实距离关系选等距方位如果只想看大圆航线心射投影直接把大圆画成直线比弧线直观很多。我自己的习惯是不确定用什么投影时先实现正交和等距方位两个做一个快速对比再按需求调整。因为正交投影的代码最简等距方位投影的物理意义最明显这两个摸透了另外两个只是换一下计算公式的事。1.3 旋转矩阵让投影中心不再是北极点很多教程讲投影公式都是以北极点为中心的。但实际项目里投影中心往往是任意经纬度比如东京、成都或者某个小岛。解决思路不是每次重新推导公式而是先做坐标旋转把想看的中心点转到“北极”位置再用标准公式投影。旋转的本质是把数据点从原始坐标系变换到以投影中心为原点的局部坐标系。具体分两步先绕z轴旋转-lon0让中心经度对准x轴再绕y轴旋转-(90-lat0)让中心点被抬到z轴顶端。MATLAB里可以写成function [xr, yr, zr] rotateToCenter(lon, lat, lon0, lat0, R) % 将经纬度旋转到以(lon0, lat0)为北极的新坐标系 lon_r deg2rad(lon); lat_r deg2rad(lat); lon0_r deg2rad(lon0); lat0_r deg2rad(lat0); x R * cos(lat_r) .* cos(lon_r); y R * cos(lat_r) .* sin(lon_r); z R * sin(lat_r); xr -sin(lon0_r) .* x cos(lon0_r) .* y; yr -sin(lat0_r) .* cos(lon0_r) .* x - sin(lat0_r) .* sin(lon0_r) .* y cos(lat0_r) .* z; zr cos(lat0_r) .* cos(lon0_r) .* x cos(lat0_r) .* sin(lon0_r) .* y sin(lat0_r) .* z; end旋转后(xr, yr, zr)中zr接近R的点就对应投影中心投影计算可以直接套标准公式。这里必须注意旋转矩阵是正交矩阵用它变换不会改变点与点之间的距离关系也不改变数据本身的形状所以投影结果不会引入额外畸变。2. 核心实现从经纬度网格到平面坐标2.1 准备投影数据网格生成与实际数据模拟手上没有现成数据时通常先生成一个全球经纬度网格来验证投影算法。MATLAB里用meshgrid生成网格很直接[lon, lat] meshgrid(-180:2:180, -90:2:90);这里lon对应当前网格点的经度lat对应纬度。注意维度顺序第一维是经度方向还是纬度方向会影响后面surf、pcolor的显示方向建议生成后马上用size检查一下。模拟数据我习惯用几个高斯函数叠加这样在投影图上能明显看出不同位置的变形。比如模拟两个异常体一个在北半球中纬度一个在赤道附近data 20 * exp(-((lat - 25).^2 (lon - 80).^2) / 2000) ... 15 * exp(-((lat 10).^2 (lon 120).^2) / 1500);这样做的原因是高斯函数的形态在投影前后变化直观投影中心附近基本不变形边缘区域会被拉伸或压缩一眼就能看出投影的性质是否符合预期。如果需要加载真实数据比如.nc文件或.mat文件里的全球气象场直接替换data变量就行。投影函数只关心经纬度网格和对应数值数据来源不影响后续处理。2.2 投影函数代码实现与参数说明把前面的旋转矩阵和各类投影公式封装成一个统一入口调用时就省心多了。我写了一个projectSphere函数支持ortho、stereo、eqdist、gnomonic四种方法function [X, Y] projectSphere(lon, lat, lon0, lat0, method, R) % 球面投影统一入口 % lon, lat: 经纬度网格(度) % lon0, lat0: 投影中心经纬度(度) % method: ortho正交, stereo球极, eqdist等距方位, gnomonic心射 % R: 球半径, 默认1 if nargin 6, R 1; end [xr, yr, zr] rotateToCenter(lon, lat, lon0, lat0, R); rho sqrt(xr.^2 yr.^2); switch lower(method) case ortho X xr; Y yr; invalid zr 0; case stereo denom R - zr; X R * xr ./ denom; Y R * yr ./ denom; invalid denom 0; case eqdist alpha acos(min(max(zr / R, -1), 1)); r R * alpha; safe rho 1e-12; X zeros(size(xr)); Y zeros(size(yr)); X(safe) r(safe) .* xr(safe) ./ rho(safe); Y(safe) r(safe) .* yr(safe) ./ rho(safe); invalid false(size(xr)); case gnomonic denom zr; X R * xr ./ denom; Y R * yr ./ denom; invalid zr 0; otherwise error(未知投影方法: %s, method); end X(invalid) NaN; Y(invalid) NaN; end几个值得注意的地方第一ortho投影的背面裁剪。zr 0意味着该点位于切平面后方物理上不可见如果不裁剪投影图上会出现内外重叠的混乱图形。这点很多新手会忽略。第二stereo投影的分母是R - zr。当zr接近R时也就是点接近投影中心旋转后的北极投影坐标会迅速发散到无穷大。这是球极投影的数学特性不是代码bug。如果渲染时发现中心附近有异常亮点多半是数据点刚好落在投影中心附近需要检查分母是否过小。第三eqdist投影在中心点处会出现0/0的情况。代码里用safe rho 1e-12先判断中心点附近直接置0这样既保证数学正确性也避免MATLAB输出NaN。第四gnomonic投影只有前半球的点有意义。zr 0的点需要裁剪否则投影结果会跑到“背后”去图形混乱且没有物理意义。2.3 背面裁剪与边界处理投影后的边界处理是最容易出细节问题的地方。拿正交投影来说理论上投影边界是一个圆圆外没有数据。但如果不裁剪直接用pcolor绘制会发现边界处颜色溢出或者出现细线。我的做法是裁剪后在X、Y中置入NaN。MATLAB的pcolor、patch等绘图函数会默认把NaN对应的面片跳过这样投影边界自然呈现出干净的圆形或扇形边界。还有一种情况需要单独处理当经纬度网格分辨率很粗时比如间隔10度投影边界会出现明显的锯齿。这时候不要在投影坐标里强行插值正确的做法是回到原始经纬度网格加密网格后再投影。比如把meshgrid的步长从10改成2边界就平滑很多。加密网格对投影算法的计算量影响不大因为整个投影过程都是向量化运算。如果投影数据本身是离散的站点数据不是规则网格那就要用scatteredInterpolant先插值到规则网格再走投影流程。这一步的插值方法选择后面第4章会专门讲。3. 可视化实战把投影结果画得专业又好看3.1 散点、曲面与patch三种绘制方式拿到投影坐标(X, Y)和对应的数据data后绘制方式有讲究。不同数据类型适合不同绘图函数。规则网格数据最推荐用pcolor配合shading interpfigure(Color, w); pcolor(X, Y, data); shading interp; colormap(parula); axis equal; axis off; colorbar;pcolor的优势在于它天然支持NaN裁剪。投影后不可见的区域会留白不需要额外画遮罩。如果是站点数据比如气象站、地震台网直接用scatter画散点图scatter(X(:), Y(:), 30, data(:), filled);这里X(:)把矩阵展开成向量scatter能正确识别每个点的坐标和数值。建议把散点大小调成20到50之间太小看不清颜色变化太大会互相覆盖。如果是球面上不均匀分布的多边形比如各国家的行政边界、板块边界用patch更合适。patch需要顶点坐标和面索引MATLAB自带海岸线数据可以直接用load coastlines; [coastX, coastY] projectSphere(coastlat, coastlon, lon0, lat0, ortho, R); plot(coastX, coastY, k-, LineWidth, 0.5);注意海岸线数据是经纬度折线不是网格投影后直接plot就行。3.2 光照、颜色映射与坐标轴控制投影到平面上以后大部分光照效果用不上但如果同时展示三维球面光照能给结果增色不少。做法很简单surf(x, y, z, data, EdgeColor, none); axis equal; camlight(headlight); lighting gouraud;camlight(headlight)表示光源放在观察者位置这样看到的地方基本都是亮的。lighting gouraud比flat更平滑适合连续场数据。如果要突出地形起伏可以加material dull降低高光反射避免数据被亮斑遮挡。颜色映射的选择也影响展示效果。MATLAB自带的parula是默认选择色盲友好连续场数据基本不会出错。jet色彩丰富但存在色带伪影建议在正式报告里少用。thermal、turbo这类新colormap也可以试试特别是数据跨度大、需要突出异常区域时。还有一个容易忽略的点设置颜色范围时不要用默认的全局极值因为投影裁剪后经常有异常大值或小值把色标拉伸变形。手动设置climclim([min(data(:)), max(data(:))]);如果想让某些区域突出可以在中位数附近做分位数截断这样色标分布更均匀。3.3 动画演示让投影中心绕着球面旋转球面投影的动画是很多人喜欢做的效果特别是展示全球数据时中心点沿经度旋转一圈整个球的数据仿佛滚动起来。实现思路很简单循环改变lon0重新调用projectSphere在循环里用drawnow刷新最后导出GIF。% 生成动画帧 F struct(cdata, [], colormap, []); k 1; for lon0 0:2:358 [X, Y] projectSphere(lon, lat, lon0, 15, ortho, R); pcolor(X, Y, data); shading interp; axis equal; axis off; title(sprintf(投影中心: %d°E, lon0)); drawnow; F(k) getframe(gcf); k k 1; end导出GIF时有一个容易踩的坑getframe抓到的RGB数据需要先做索引化再写入filename rotating_ortho.gif; for i 1:length(F) [A, map] rgb2ind(F(i).cdata, 256); if i 1 imwrite(A, map, filename, LoopCount, Inf, DelayTime, 0.05); else imwrite(A, map, filename, WriteMode, append, DelayTime, 0.05); end end如果帧数太多导致GIF文件体积过大可以每隔一帧保存或者把图像尺寸调小。另外MATLAB R2021a之后推荐用exportgraphics结合animatedline不过GIF格式还是imwrite这招最通用兼容性最好。4. 常见问题与实测排查4.1 极点附近为什么总是一片空白这是我在做球面投影时遇到最多的问题。正交投影和心射投影中极点附近出现大片空白是正常现象因为这两类投影本身只能显示半个球面位于背面或中心平面之后的点被裁剪掉了。但如果是等距方位投影或球极投影也出现极点空白那就要检查代码里有没有不小心把zr 0的点全裁剪掉。如果用的是我上面的projectSphere代码eqdist分支的invalid是false(size(xr))也就是不裁剪任何点正常情况下全球都能显示。还有一种情况数据本身的方案问题。比如全球网格里lat从-90到90如果中心设在lat030那么南极点对应的投影坐标在等距方位投影下会是一个很大的圆环边缘如果绘图窗口或者坐标轴范围被手动缩小就看不到边缘区域误以为数据空白。解决办法是用axis equal并手动放大坐标轴范围或者先绘制完整投影坐标确认范围后再调整视图。4.2 数据点太多绘图卡成PPT全球0.1度分辨率的网格有650万个点这种情况下pcolor直接绘制会非常卡甚至在低配电脑上直接内存溢出。我实测的优化方案有三个。第一个方案是降采样。投影前先对经纬度网格做粗化比如从0.1度抽稀到1度然后用griddata把原始数据插值到粗网格上。视觉差异在高倍放大时才明显整体展示效果几乎一样。第二个方案是改用texturemap纹理映射。把数据当成纹理贴到球面上MATLAB只需要处理一张二维纹理图和一个粗糙的球面网格效率比逐面片绘制高得多。但texturemap需要对图像做预处理数据本身必须是规则的适用于全球规则网格。第三个方案是分块绘制。如果计算量主要来自单个图形对象太大可以按经纬度分成几个子区域分别用surf绘制然后用hold on叠加。这样每个图形对象的数据量小一些交互响应更快。4.3 投影后的插值怎么处理才合理投影后坐标是非均匀分布的直接调用interp2做规则网格插值会得到错误结果因为在投影平面上“规则网格”已经不再是规则的。正确的做法是回到原始经纬度空间处理。如果手头是散点数据要投影到网格上用scatteredInterpolantFinterp scatteredInterpolant(lon_pts(:), lat_pts(:), data_pts(:), natural); data_grid Finterp(lon, lat);natural插值自然邻域插值在球面投影场景下表现最稳比linear更容易避开边界区域的震荡代价是计算稍慢。如果数据量特别大再换linear。如果手头是规则网格数据但投影后要重新采样的网格先插值到经纬度规则网格再投影不要先投影再插值。顺序反了会导致采样点在投影平面分布极不均匀边缘区域信息丢失。4.4 常见报错速查表报错信息可能原因快速解决Matrix dimensions must agree经纬度矩阵形状不一致或经纬度用了向量没转成网格用meshgrid生成网格检查size(lon) size(lat)NaN/Inf not supported投影分母为0比如stereo投影中心点检查denom用eps保护或裁剪参照projectSphere代码Index exceeds matrix dimensions动画帧索引越界或绘图数据维度不对检查F的索引和循环变量范围Error using surf: Z must be a matrixsurf的z参数不是二维矩阵投影坐标需要在网格下保持矩阵形状不要用X(:)展开Colormap must be a valid colormapcolormap名称拼写错误用doc colormap查看内置颜色表全名这里特意提一个隐藏问题rotateToCenter函数输入的是经纬度矩阵如果用户在调用时把经度向量和纬度向量传进去而没有经过meshgrid维度对不上就会报第一个错误。我自己的习惯是在函数入口加一行断言assert(isequal(size(lon), size(lat)), lon和lat必须同维度);5分钟的调试时间能省下不少。还有一个小坑在R2022a之后的版本里pcolor会默认把最后一行和最后一列丢弃。如果投影图的边缘缺了一圈检查一下是不是这个版本特性导致的。解决办法是在调用pcolor前给X、Y、data补一圈重复值或者改用surf(..., EdgeColor, none)代替。末尾补充这两个小习惯让我少踩很多坑文章写到这里再分享两个实际项目里的小经验。第一投影函数里所有角度统一用弧度制可以避免大量错误但输入输出接口保留角度制让调用方更顺手——用注释写清楚单位即可。第二写投影代码前先用一个分辨率很低的网格比如30度间隔快速画出投影轮廓确认投影方向和形状没有反再换成高分辨率。低分辨率下错误一眼可见高分辨率下往往被细节掩盖排查起来很费劲。希望这篇能帮到正在做球面投影的同行。如果遇到我这里没提到的奇怪现象建议从旋转矩阵和分母检查入手球面投影90%的异常都出在这两个地方。