MATLAB画断面图：数据来源HYCOM，数据格式NetCDF

最新推荐文章于 2024-01-04 11:50:35 发布

玛特莱布DUMMY

最新推荐文章于 2024-01-04 11:50:35 发布

阅读量1k

点赞数 1

文章标签： matlab 开发语言

本文链接：https://blog.csdn.net/weixin_47011923/article/details/134269645

版权

1.简介

在阅读文献时，我们常常能看见这样的断面图：

要想画出这样的断面图，我们需要完成以下几个步骤：

（1）重构数据；

（2）生成网格；

（3）画断面图；

2.重构数据

2.1 HYCOM数据结构

我们读取HYCOM数据后，其数据结构是这样的：

可以看到，HYCOM数据里的水温变量，其结构为114×172×38×90，即经度×纬度×水层×时间。解释一下为什么水层数的含义。我在下载HYCOM数据时，下载了0 - 3000 m的数据，HYCOM将3000 m分成了很多个水层，如0 m，2 m，4 m，6 m，8m，...，3000 m（见下图）。在MATLAB中，四维数据是无法画图的，因此我们需要进行数据降维。

2.2 数据降维

我们先提取第一个水层，第一天的水温数据，需要通过以下代码实现：

clc;clear all;close all;
filename = 'C:\Users\98485\Desktop\Academic\工作文件\HYCOM模式\HYCOM数据\HYCOM_data\temp_0_3000_20201201_20210228.nc'
ncdisp(filename)
temp_all = ncread(filename,'water_temp');
temp = temp_all(:,:,1,1);

如果不会读取NetCDF格式数据，可以参照我以前的文章：

MATLAB读取nc文件，并绘制南海水深图_m_contourf-CSDN博客https://blog.csdn.net/weixin_47011923/article/details/129830867?spm=1001.2014.3001.5502 第一个水层，第一天的水温数据是一个114×172的double（下图的temp）。可以看到，水温数据是按照经纬度排列的，即温度数组行和列的数量刚好和经纬度对应。其中包含了整个南海的水温数据。我们需要提取出120.72 °E经度断面上，纬度范围为18.5 - 22 °N的水温数据。

2.3 数据提取

我们找到120.72 °E在lon数组中的序号，通过下图可以看到，序号是52。再在lat数组中找到纬度的起止序号，分别是31和118，不再赘述。

通过以下代码，将120.72 °E经度断面上，纬度范围为18.5 - 22 °N的水温数据提取出来：

temp = temp(52,:);
temp = temp(:,31:118);

这样我们就得到了第一个水层，第一天的，在120.72 °E经度断面上，纬度范围为18.5 - 22 °N的水温数据。见下图的temp，1×88 double。仅有一层的水温数据是画不出断面图的，我们下一步需要写一个循环，将38个水层的指定水温数据都读取出来，再按深度排列，组装在同一个数组中。

2.4 数据组装

这一步的循环非常简单，具体解释见代码内的注释。

for j = 1:38                  % j = levels（水层数目）
    temp = temp_all(:,:,j,1); %提取第j个水层，第1天的水温数据
    temp = temp(52,:);        %提取170.72°E经度断面的所有水温数据
    temp = temp(:,31:118);    %提取纬度范围为18.5 - 22 °N的水温数据
    temp_whole(j,:) = temp;   %把提取完成的水温数据，赋值给temp_whole数组的第j行
end

执行上述循环后，得到的数组——temp_whole，为38 × 88, 38代表38个水层，见下图：

3. 生成网格

生成网格通过以下代码实现：

depth = [0,2,4,6,8,10,12,15,20,25,30,35,40,45,50,60,70,80,90,100,125,150,200,250,300,350,400,500,600,700,800,900,1000,1250,1500,2000,2500,3000];
lat_range = lat(31:118);
[lat_range,depth] = meshgrid(lat_range,depth);

首先，创建深度数组（depth），在里面储存水层的深度数值，这就是断面图的纵坐标；

其次，创建纬度数组（lat_range），在里面存储18.5 - 22 °N范围内的纬度值，这就是断面图的横坐标；

最后，使用meshgrid（）函数，生成网格。

meshgrid（）函数的作用：depth是1 × 38的数组，lat_range是88 × 1的数组，使用meshgrid（）函数后，depth和lat_range都变成了38 × 88的数组。做这一步的目的是后续使用的contourf（a,b,c）函数需要a、b、c三个数组的维度一致。

4. 画断面图

画断面图需要通过contourf（）函数来实现，见下列代码：

figure(1)
set(0,'defaultfigurecolor','w')   %%% 图片底色为白色
contourf(lat_range,depth,temp_whole,18,'LineStyle','none');
ylabel('Depth(m)');
set(gca, 'YDir', 'reverse', 'FontSize', 14, 'FontName', 'Times New Roman','FontWeight', 'bold');
colorbar;
caxis([0,30])
title('Mean Temperature in winter at 120.72 E','fontname','Times New Roman','fontsize',18);

contourf(lat_range,depth,temp_whole,18,'LineStyle','none')：以lat_range为横坐标，以depth为纵坐标，绘制关于数组temp_whole的等值线图，总共绘制18根等值线。画出来的图见下图：

需要注意的是，物理海洋学画的图y轴是从上向下递增的，所以需要将y轴翻转180°，实现y轴翻转180°是通过这条代码的红字部分实现的：set(gca, 'YDir', 'reverse', 'FontSize', 14, 'FontName', 'Times New Roman','FontWeight', 'bold'); y轴翻转180°后效果见下图：

5. 附录

5.1 完整代码

clc;clear all;close all;
filename = 'C:\Users\98485\Desktop\Academic\工作文件\HYCOM模式\HYCOM数据\HYCOM_data\temp_0_3000_20201201_20210228.nc'
ncdisp(filename)
temp_all = ncread(filename,'water_temp');
lon = ncread(filename,'lon');
lat = ncread(filename,'lat');
for j = 1:38 % j = levels 先找一个水层，提取所有的日期的数据，求平均值
    temp = temp_all(:,:,j,1);
    temp = temp(52,:);
    temp = temp(:,31:118);
    temp_whole(j,:) = temp;
end
depth = [0,2,4,6,8,10,12,15,20,25,30,35,40,45,50,60,70,80,90,100,125,150,200,250,300,350,400,500,600,700,800,900,1000,1250,1500,2000,2500,3000];
lat_range = lat(31:118);
[lat_range,depth] = meshgrid(lat_range,depth);
% %%
figure(1)
set(0,'defaultfigurecolor','w')   %%% 图片底色为白色
contourf(lat_range,depth,temp_whole,18,'LineStyle','none');
ylabel('Depth(m)');
set(gca, 'YDir', 'reverse', 'FontSize', 14, 'FontName', 'Times New Roman','FontWeight', 'bold');
colorbar;
caxis([0,30])
title('Mean Temperature in winter at 120.72 E','fontname','Times New Roman','fontsize',18);