【matlab】三维图画法-以卫星测高波形数据为例

1. 前言

matlab中画三维图的方法（函数）可以分为如下类别：

1.三维曲线

代表函数：plot3

2.三维曲面

代表函数：surf、mesh、pcolor

注：绘图函数有很多，上面未完全列出。

绘制三维曲面时一定要注意数据的准备：

确定自变量的取值范围和取值间隔；
构成自变量x,y的自变量“格点”矩阵（理解这个很重要！可以参考：【matlab】meshgrid生成网格原理1）；
计算在自变量采样“格点”上的函数值。

2. 任务

获取一些测高波形数据，作图显示。

3. 分析

作图其实与测高没有关系，作图方法都是通用的。由于测高波形数据是三维的，所以要用到matlab里面的三维绘图函数，可以选择的方法有很多，这里列出几个我试过的。

值得注意的是数据的正确使用，也就是上面提到的“格点”矩阵的生成，如果理解画图原理，就很容易理解。

3.1 一个不太正确的画法

先看两幅图：

图1

图2

这两幅图就可以看作是波形数据的3维图。一眼看过去没有什么大问题，图1是一幅3维图，图2也是一幅3维图，只不是过是从图1的顶上往下看，所以成一个二维的形状，下面是代码。如果一定要生成跟我一样的图，请去下载相应数据，数据下载方法可以参考：如何下载Jason2测高数据，还是提供一个数据：链接：http://pan.baidu.com/s/1jIdvc02 密码：862w

clc,clear,clf
filename  = 'JA2_GPS_2PdP117_240_20110914_142421_20110914_152034.nc';
lat       = ncread(filename,'lat');
lat_20hz  = ncread(filename,'lat_20hz');
lat_index = find(lat >= 22 & lat <= 25); % 选出纬度在该区域的数据
wv_all    = ncread(filename,'waveforms_20hz_ku');

cnt = 1;
for i = min(lat_index) : max(lat_index)
    new_wv(:,cnt) = wv_all(:,1,i);
    new_lat(cnt) = lat(i);
    cnt = cnt + 1;
end

bin = 1:104;                         % 采样数
[BIN, LAT] = meshgrid(new_lat, bin); % 生成网格
surf(BIN, LAT, new_wv);              % 画曲面

shading interp;         % 表面光滑
colormap(jet);          
% view(0,90) % 视角
colorbar
xlabel('lat');
ylabel('bin');
zlabel('wave poer');
print(gcf, '-djpeg', 'error2.jpg')

为什么说它不太正确？因为波形数据都是一条一条的，像下面这样：

图3

每条波形之间其实没有什么关系的，而图1、2画的是一整个曲面，因为插值的原因，每条波形相互影响，所以“不太正确”。

3.2 应该是正确的画法

“正确”的图看起来应该是这样的：

图4

图5

图4和图5才应该是一个正确的画法，可以看到每条波形都独立开了，下面是代码：

filename  = 'JA2_GPS_2PdP117_240_20110914_142421_20110914_152034.nc';
lat       = ncread(filename,'lat');
lat_20hz  = ncread(filename,'lat_20hz');
lat_index = find(lat >= 22 & lat <= 25); % 选出纬度在该区域的数据
wv_all    = ncread(filename,'waveforms_20hz_ku'); % 所有波形数据

figure;
% hold on; % 用pcolor时打开
for i = min(lat_index) : max(lat_index)
    x = [lat(i), lat(i + 1)];
    y = 1 : 104;
    [X, Y] = meshgrid(x, y);
    Z = [wv_all(:,1,i), wv_all(:,1,i)];
%     pcolor(X, Y, Z); % 平面图可以是surf从顶往下看，也可以用pcolor
    surf(X, Y, Z);
    hold on;
end
hold off;
view(0, 90)
shading interp; 
colormap(jet);
colorbar
xlabel('lat');
ylabel('bin');
axis([22, 25, 1, 104])
% axis([22, 25, 1, 104, 0, 250]) % 作平面图时需注释
print(gcf, '-djpeg', 'right3.jpg')

这里用到的技巧就是分开画，一条一条的画，看代码更直观。

ssh_alitheia

生而遗忘，文字挽留