企业官网建设流程全解析

matlab球面投影（二）来了。上一期做完了经纬度网格以后，留言区里问得最多的一类问题，基本都绕不开同一个小目标：怎么把一张平面的图像正确地贴到球面上？怎么把球面上的离散点投到平面里做密度统计？这两个问题背后的核心其实都是同一件事——matlab球面投影。只是有人做的是“球面→平面”的正向投影，有人做的是“平面→球面”的逆向贴图。这一期我就把这两条路线一次性讲透，从坐标换算、等距柱状投影、正射投影，到纹理贴图、常见坑，全部用MATLAB代码落地。适合做全景图、遥感可视化、物理场分布、雷达覆盖和地震点分布这类需求的朋友，新手也能跟着跑通。

1. 从球面到平面的坐标换算：先搞清楚投影模型

1.1 经纬度与三维坐标的对应关系

先别急着写投影公式，球面投影的第一步永远是“球面上的点”用什么样的坐标描述。绝大多数MATLAB可视化场景里，我们拿到的原始数据是经度lon和纬度lat，单位是度。要把它们画成三维球面，必须先转成直角坐标。

地理学里有一个约定俗成的转换关系，假设球的半径为R：

X = R * cosd(lat) * cosd(lon); Y = R * cosd(lat) * sind(lon); Z = R * sind(lat);

这里用cosd和sind是为了省去手动转弧度的麻烦，MATLAB里这两个函数接受度数，非常方便。注意纬度lat的范围是[-90, 90]，经度lon的范围是[-180, 180]。有些场景里经度会用0到360度表示，转换前建议统一成[-180, 180]，否则投影公式里会出现偏移，后面排查起来很头疼。

如果你想生成整个球面的网格点，经典做法是这样的：

R = 1; lon = linspace(-180, 180, 361); lat = linspace(-90, 90, 181); [LON, LAT] = meshgrid(lon, lat); X = R * cosd(LAT) .* cosd(LON); Y = R * cosd(LAT) .* sind(LON); Z = R * sind(LAT); figure; surf(X, Y, Z, 'EdgeColor', 'none', 'FaceAlpha', 0.8); axis equal;

我习惯把这段代码单独抽成一个函数，比如sph2cart_vec(lon, lat, R)，后面所有投影、贴图、绘制点分布都会用到它。球面投影项目里，坐标转换函数是地基，地基不稳，后面全是白搭。

1.2 常见球面投影方式的取舍

“球面投影”这四个字其实是一个很大的家族。不同场景选错投影方式，出来的图形会非常误导人。我这里列一个常用对照表，方便你选型。

这期我重点讲两个最常用的：等距柱状投影和正射投影。前者做全景图、纹理贴图特别方便，经纬度和像素坐标几乎一一对应；后者适合做“从太空看地球”的效果，肉眼看着最自然。其他投影方式如果以后有需要，可以在这一套框架基础上扩展。

1.3 把投影方式落实到代码的思路

我推荐你在一开始就建立两个对称的函数：一个是正向投影projectSphericalToPlane，输入经纬度和投影类型，输出平面像素坐标；另一个是逆向投影unprojectPlaneToSpherical，输入像素坐标，反算经纬度。这样无论后面是画经纬网、贴纹理还是做点云统计，都是调这两个函数，逻辑清晰，也不容易改崩。

一个最小可用的正向投影函数长这样：

function [u, v] = projectSphericalToPlane(lon, lat, projType, width, height) % 默认等距柱状投影 switch lower(projType) case 'equirectangular' u = (lon + 180) / 360 * width; v = (90 - lat) / 180 * height; case 'ortho' % 正射投影暂时只做沿Z轴方向的简单投影 R = 1; X = R * cosd(lat) .* cosd(lon); Y = R * cosd(lat) .* sind(lon); Z = R * sind(lat); mask = Z >= 0; u = nan(size(X)); v = nan(size(Y)); u(mask) = X(mask) * width / 2 + width / 2; v(mask) = -Y(mask) * height / 2 + height / 2; otherwise error('未知的投影类型'); end end

这里有个很多人容易忽略的细节：图像坐标系的v轴是向下为正的，屏幕上第一行是图像顶部。所以我用v = (90 - lat) / 180 * height，这样北纬90度对应v=0，南纬90度对应v=height。如果后面发现贴图上下颠倒，多半就是这里的方向选择反了。

2. 等距柱状投影的MATLAB实现：全景图映射的核心

2.1 正向投影：经纬度网格到平面

等距柱状投影的另一个名字叫“经纬度直投影”。它的核心思想很简单：把经度当作横坐标，纬度当作纵坐标，然后线性拉伸到目标图像尺寸。

假如我们要把全球的经纬度网格画到一张宽720、高360的图像上，那么：

width = 720; height = 360; [lonGrid, latGrid] = meshgrid(linspace(-180, 180, 37), linspace(-90, 90, 19)); [u, v] = projectSphericalToPlane(lonGrid, latGrid, 'equirectangular', width, height); figure; plot(u, v, 'k-', 'LineWidth', 0.8); hold on; plot(u', v', 'k-', 'LineWidth', 0.8); axis equal; axis([0 width 0 height]); set(gca, 'YDir', 'reverse'); xlabel('u / pixel'); ylabel('v / pixel');

运行之后你会看到，经纬线在平面上是互相垂直的直线网格。这看起来很规整，但它有个天生的毛病：高纬度地区的面积被剧烈放大。北极圈附近的一小格，实际面积可能只有赤道附近一小格的几分之一，但在等距柱状投影图上两者占据同样的像素数量。

所以如果你以后要做“全球密度统计图”，用等距柱状投影做底图没问题，但解读统计值时一定要小心，不能直接把每个像素的数值当作等权样本。最稳妥的办法是在统计时给每个像素乘一个cosd(lat)的权重，用来补偿高纬度的面积膨胀。

2.2 逆向投影：二维图像到球面纹理映射

等距柱状投影最大的价值，其实是它的逆过程：给一张全景图，怎么把它贴回球面。我以前第一次做的时候以为要自己算纹理坐标，研究半天，其实MATLAB里surf函数配合texturemap可以直接搞定。

假设你有一张全景图world.jpg，宽度是w，高度是h，图像第一行是北极，最后一行是南极，最左边是经度-180度，最右边是经度180度。那么逆向投影代码为：

img = imread('world.jpg'); [h, w, ~] = size(img); % 把图像坐标映射到经纬度 lon = linspace(-180, 180, w) * pi / 180; lat = linspace(90, -90, h) * pi / 180; % 注意是从90到-90 [LON, LAT] = meshgrid(lon, lat); R = 1; X = R * cos(LAT) .* cos(LON); Y = R * cos(LAT) .* sin(LON); Z = R * sin(LAT); figure; surf(X, Y, Z, 'CData', img, 'FaceColor', 'texturemap', 'EdgeColor', 'none'); axis equal; view([40, 20]);

这里有几个非常关键的细节。第一，lat是从90度到-90度，不能写反，否则图像上下颠倒。第二，meshgrid生成的LAT矩阵行方向对应纬度，列方向对应经度，所以矩阵尺寸是h x w，正好和图像h x w一致。第三，CData直接传图像矩阵，FaceColor必须设置为'texturemap'，否则MATLAB会把图像当作一张普通的平面色图，不做纹理映射。

我建议你先不要用真实的全景图测试，而是用一张自带明显方向感的图，比如左上角涂成红色、右下角涂成蓝色的测试图。贴上球面之后，红色应该出现在球面的左上方，蓝色在右下方。这样一测，有没有翻转、镜像，一眼就能看出来。

2.3 全景图旋转与视角调整

贴好球体之后，最常见的需求是把球体转一转，模拟视角变化。有些人会不自觉地调用rotate函数，但我想提醒一句：rotate修改的是图形对象的坐标数据，用不好会让纹理发生畸变。更优雅的做法是直接改变纹理图像的偏移，让全景图在球面上“滚动”。

比如你想把全景图水平旋转一个角度，只需要把图像的列做循环位移：

shiftPixels = 120; imgShift = circshift(img, shiftPixels, 2); % 重新用imgShift做上面同样的surf贴图

这样做的好处是纹理始终是同一个球面坐标下的自然映射，只是经度起点变了，不会产生任何几何变形。如果你做全景播放器或者交互式展馆漫游，这一招非常实用。

视角本身则用view函数控制，比如view([40, 20])表示从方位角40度、仰角20度观察球体。更进一步，可以设置campos和camtarget来模拟从外部相机看球体。记住一个原则：投影是数据本身的事，视角是相机的事，两者不要混在一起。

3. 正射投影与透视投影：模拟“从太空看地球”

3.1 正射投影公式与向量求交

正射投影的直观理解就是“从无穷远的地方垂直看向球体”。它不考虑近大远小，所有可见点都沿着平行线投影到平面上。因为透视被压平，所以它很适合用来做卫星云图、全球洋流、地震点分布这类需要“俯看半球”的图。

最简单的正射投影，视线方向沿Z轴，投影平面是XY平面，那么投影坐标就是：

u = X; v = Y;

同时只能看到面向观察者的半球，也就是Z大于等于0的那部分。如果用前面的投影函数，你会发现输出结果里出现一个圆形区域，圆外面全是NaN。

如果要把经纬网格画成正射投影的圆盘，记得不要直接用散点法去画，因为靠近边缘的经纬线会因Z<0被滤掉，导致线条断裂。更好的办法是把每一根经纬线当成独立的曲线，逐条投影。例如绘制经线：

figure; hold on; for lon0 = -180:30:180 latLine = linspace(-90, 90, 181); lonLine = lon0 * ones(size(latLine)); [u, v] = projectSphericalToPlane(lonLine, latLine, 'ortho', 800, 800); plot(u, v, 'b-'); end % 同理绘制纬线

这里投影函数内部已经把Z<0的位置置成了NaN，所以plot会自动断开不可见的部分，整根经线看起来会很干净，边缘正好贴合圆形边界。

3.2 透视投影下的球面边缘、光照与质感

正射投影虽然计算简单，但“从太空看地球”的真实观感还是透视投影更自然。透视投影和正射很大的区别在于：透视投影会近大远小，球面边缘部分会因为视线和表面相切而明显变形。

在MATLAB里，如果你只是在屏幕上看一个3D球体，最省事的做法不是自己写透视公式，而是直接设置相机位置：

hSurf = surf(X, Y, Z, 'CData', img, 'FaceColor', 'texturemap', 'EdgeColor', 'none'); axis equal; set(gca, 'CameraPosition', [0, 0, 3]); set(gca, 'CameraTarget', [0, 0, 0]); light('Position', [1, 0, 2]); material dull; set(hSurf, 'AmbientStrength', 0.7);

CameraPosition决定了相机离球体的距离，距离越近透视效果越强烈，距离越远画面越接近正射投影。如果你需要输出这种投影效果为一张二维图像，可以直接用exportgraphics(gca, 'spherical_view.png', 'Resolution', 300)把整个坐标区的内容导出来，效果比手动做几何投影更可控。

另外提醒一个光照相关的坑：贴图球体默认的FaceLighting会影响纹理亮度。如果发现贴图看起来发灰发暗，可以调高AmbientStrength，或者干脆设置lighting none。热词里经常有人搜“matlab亮度平衡”，在很多可视化里都是光照方向没调好造成的。

3.3 从离散点云到球面：一个实用案例

前面讲的都是连续曲面贴图，但球面投影还有一个高频应用：把离散的球面点云投影到二维平面做密度统计。比如全球地震事件分布、气象站点数据、卫星地面轨迹，这些数据本质上都是“经纬度 + 数值”。

假设你有一个数组lonData和latData，想统计全球分布密度，可以直接用等距柱状投影的二维直方图：

edgesLon = linspace(-180, 180, 181); edgesLat = linspace(-90, 90, 91); N = histcounts2(lonData, latData, edgesLon, edgesLat); figure; imagesc(edgesLon, edgesLat, N'); set(gca, 'YDir', 'normal'); axis xy; xlabel('经度'); ylabel('纬度'); colorbar; colormap(parula);

这里重点是histcounts2的输出矩阵维度顺序，矩阵行数等于edgesLon的区间数，列数等于edgesLat的区间数，所以显示时要用N'转置，并且YDir设为normal，让纬度从下到上递增。很多人第一次跑这里都会发现图像纬度方向反了，几乎成了必踩的坑。

如果你想把直方图结果再贴回球面，只需要把密度矩阵N作为CData传给surf，用texturemap方式渲染即可。这样一个“全球密度分布球”就出来了，既能转着看，又能导出平面图，效果相当能打。

4. 进阶：把任意图像“糊”到球面并导出高质量图片

4.1 常用纹理贴图注意事项

很多朋友第一次做纹理贴图，都会遇到“贴上去怎么是反的”或者“图怎么碎了”的问题。这里我把最常踩的雷一次性说清楚。

第一，图像上下颠倒。原因是网格的纬度方向与图像的Y方向不一致。解决办法很简单：在surf之前检查你的lat是不是从90到-90，如果不是，改成linspace(90, -90, h)通常就解决了。如果还颠倒，就flipud(img)，不要硬调坐标轴方向，容易引入其他问题。

第二，左右镜像。这个一般在图像处理软件里不容易发现，但贴到球面上，本来面朝左的文字会变成面朝右。检查起来比较麻烦。我的做法是准备一张带“L”和“R”标记的测试图，贴在球面上，根据实际显示结果决定是否fliplr。

第三，经度接缝。等距柱状投影图像最左边经度是-180度，最右边经度是180度。因为MATLAB把球面用网格表示，第一列和最后一列在球面上其实是同一个位置，但如果在纹理生成时不处理，接缝处常见一条明显的“裂缝”。处理方案是把经纬度网格多包一圈：

lonWrap = [lon, lon(:, 1) + 360]; latWrap = [lat, lat(:, 1)]; imgWrap = [img, img(:, 1, :)];

这样首尾自然相接，surf绘制时球面就是一个真正的闭环。

4.2 分辨率与计算精度的取舍

球面贴图的分辨率选择是一门学问。纹理图像分辨率太高，比如4000x2000的图片，如果你用同样尺寸的surf网格去承载，MATLAB会卡成幻灯片。其实surf的texturemap并不要求网格点数和图像像素数一致。我们可以用粗网格表示球面，用高分辨率图片作为纹理。

推荐的做法是这样：先用粗网格生成球面，比如纬向81点、经向161点，再把纹理图像用imresize缩到同样尺寸或者接近的尺寸，然后贴图。如果纹理图像和网格尺寸不匹配，可以这样：

targetR = 64; targetC = 128; imgResized = imresize(img, [targetR, targetC]);

imresize会自动做插值，不会让画面出现明显锯齿。反过来，如果网格比纹理更细，反而容易出现特效纹理的“摩尔纹”。所以网格点数和纹理分辨率不要差太多，一般长宽比保持2:1，分辨率适中即可。

导出的环节同样重要。我强烈建议用exportgraphics替代saveas和print，特别是导出PNG或EPS的时候。代码非常简单：

exportgraphics(gca, 'sphere_view.png', 'Resolution', 300); exportgraphics(gca, 'sphere_view.eps', 'ContentType', 'vector');

exportgraphics会严格按照当前坐标区的显示区域导出，没有多余白边，也支持300dpi以上的高分辨率输出。如果你需要写论文配图，这个函数比saveas省心太多了。

4.3 结合图像处理：色阶、叠加等值线

球面投影不只是几何映射，很多时候我们其实是想把“场”画到球面上。比如全球温度场、重力异常、电磁场覆盖范围。这些数据通常是经纬度网格上的数值矩阵，只需要把数值矩阵当成纹理，贴到球面就变成了一张立体物理场图。

我这里给一个可以直接抄的示例，模拟全球温度分布并叠加等值线：

lat = linspace(-90, 90, 181); lon = linspace(-180, 180, 361); [LON, LAT] = meshgrid(lon, lat); % 模拟温度：赤道高、两极低，再加上一些经向扰动 T = 30 * cosd(LAT) + 10 * sind(LON) .* cosd(3 * LAT); % 平面投影图 figure; imagesc(lon, lat, T); axis xy; hold on; contour(lon, lat, T, 'k', 'LineWidth', 0.5); colormap(jet); colorbar; xlabel('经度'); ylabel('纬度'); title('全球温度场等距柱状投影'); exportgraphics(gca, 'temperature_equirect.png', 'Resolution', 200); % 球面贴图 [LAT2, LON2] = meshgrid(lat * pi/180, lon * pi/180); X = cos(LAT2) .* cos(LON2); Y = cos(LAT2) .* sin(LON2); Z = sin(LAT2); figure; surf(X, Y, Z, 'CData', T, 'FaceColor', 'interp', 'EdgeColor', 'none'); axis equal; view([40, 20]); exportgraphics(gca, 'temperature_sphere.png', 'Resolution', 300);

注意这里第2段代码里meshgrid(lat*pi/180, lon*pi/180)生成的是lon x lat的矩阵，所以X的尺寸和T要一致。如果实际运行中尺寸对不上，检查一下是不是T的行列方向和X反了。这类问题往往很小，但特别耗时间，所以我习惯在surf之前加一行size(T), size(X), size(Y), size(Z)打印，先确认尺寸全对再渲染。

5. 常见问题与排查技巧实录

5.1 经纬线在平面投影图上为什么扭曲得那么夸张

等距柱状投影的“扭曲”是它本身的设计决定的，不是代码写错了。高纬度地区在真实球面上面积收缩，但等距柱状投影却把它们拉成和赤道一样宽，所以看起来所谓“变形”其实是为了保持经纬线横平竖直付出的代价。如果你做的是面积统计类工作，我建议底层数据不要直接用投影后的像素面积，要用cosd(lat)做面积权重。

5.2 surf贴纹理后图像上下颠倒或左右镜像

这是最典型的纹理映射问题。优先检查你的lat序列方向。正确的做法是从北极对应图像第一行，也就是lat = 90对应图像row=1。如果还是颠倒，就试flipud(img)。左右镜像则用fliplr(img)。不要盲目改AxisDirection，在3D视图下这个设置不如二维图直观，反而容易造成混乱。

5.3 导出后的PNG/EPS和屏幕显示不一样

这个我踩过太多次。屏幕上看好好的，导出之后图就拉扁了，或者边缘多了一圈空白。原因通常是坐标轴的DataAspectRatio没有设置。解决方法是：在exportgraphics之前，先执行axis equal，如果仍不满意，再执行daspect([1 1 1])。对于需要矢量导出的图，用ContentType','vector'导出的EPS会在放大后保持干净线条，但要注意有些复杂的纹理贴图导出矢量格式后体积巨大，不如直接导出高分辨率PNG。

5.4 画全球范围曲线时经度±180°处出现横穿整张图的斜线

这个一眼看上去像“bug”，其实是plot把不连续的点连起来了。比如一条从170度到-170度的曲线，中间隔着180度，直接连线就会横穿整个地图。解决办法是提前检测经度跳变，插入NaN断线：

lonJump = abs(diff(lon)) > 180; breakIdx = find(lonJump); lon(breakIdx + 1) = NaN; lat(breakIdx + 1) = NaN;

这样plot会在断点处停止画线，不会生成跨屏幕的斜线。如果你要处理的曲线不止一条，可以对每条线都做一遍这个处理。

5.5 正射投影边缘老是有缺口

正射投影的边缘就是球面的“地平圈”，也就是Z接近0的地方。如果只用Z >= 0过滤整根线，那么起点或者终点位于边缘附近的线会部分缺失。解决办法有两个：一是提高采样密度，让边缘附近有足够多的点；二是绘制每一根经纬线前，先把整条线按经度或纬度的可见范围切分成两段，分段绘制。对于大多数可视化需求，第一条就够了，别在缺口上死磕。

5.6 球面贴图有“接缝”怎么处理

接缝问题在前面已经提到，核心原因是首尾两列在球面上其实是同一个位置，但纹理坐标在边界处断开了。最实用的一招是给纹理图像和经纬度网格都“多包一圈”，让末尾和开头重合，这样MATLAB在绘制时就不会出现缝隙。具体代码在4.1节里已经给了，你直接用在lonWrap、latWrap、imgWrap之后，重新算出X,Y,Z再贴图即可。

6. 常用函数封装建议

走到这里，球面投影的主要流程你已经都接触过了。最后我想给你一个“项目健壮性”建议：把下面几个功能封装成独立脚本，之后做任何球面可视化都直接调用，不用每次重新写。

% 1. 经纬度转球面坐标 function [X, Y, Z] = ll2sphere(lon, lat, R) X = R * cosd(lat) .* cosd(lon); Y = R * cosd(lat) .* sind(lon); Z = R * sind(lat); end % 2. 球面坐标转经纬度 function [lon, lat] = sphere2ll(X, Y, Z) lon = atan2d(Y, X); lat = asind(Z ./ sqrt(X.^2 + Y.^2 + Z.^2)); end % 3. 等距柱状投影正反变换 function [u, v] = equirectProject(lon, lat, width, height) u = (lon + 180) / 360 * width; v = (90 - lat) / 180 * height; end function [lon, lat] = equirectUnproject(u, v, width, height) lon = u / width * 360 - 180; lat = 90 - v / height * 180; end

这些函数都很短，但胜在“一次定义，到处使用”。我自己的经验是，凡是在可视化脚本里反复出现的坐标转换代码，最好都抽成这种小函数。否则改一个方向符号，可能要把整个文件翻一遍，还容易漏。

说到个人体会，我做了这么多年MATLAB可视化，最深的感受是：球面投影的核心从来不是背公式，而是想清楚你正在做的是“正向投影”还是“逆向贴图”，以及屏幕坐标系的方向到底怎么定义。这两个问题一旦想明白，剩下的一切都是查API和调参数的事。希望这期内容能帮你把球面投影的套路彻底打通，后续再遇到全景图、地球纹理、空间点分布这类需求时，可以直接照葫芦画瓢。