matlab可视化综合实验

发表于 | 分类于 | 阅读次数

实验目的

1.了解科学数据netcdf的存储模式，学会读取和创建netcdf等科学数据。
2.学会应用Matlab对数据进行处理分析。
3.应用Matlab和Surfer等软件对已有数据进行可视化显示。
4.学会应用contour绘制等值线图和pcolor绘制伪彩色图。
5.掌握Matlab句柄操作和对图像显示的控制。

主要内容

1.NetCDF格式的读取
2.NetCDF格式的写入
3.Matlab绘图与格式化
4.Surfer绘图

实验主要仪器设备、软件

1.硬件准备：PC机；
2.操作系统：Linux系统或Windows 2k以上系统；
3.软件开发工具：Matlab 2013a 、Surfer 80
4.使用数据：ECMWF 1965年全年海表压强(Mean sea level pressure)和皮温(skin temperature)的netcdf文件。

NetCDF文件介绍

NetCDF简介

  NetCDF(network Common Data Form)网络通用数据格式是由美国大学大气研究协会（University Corporation for Atmospheric Research，UCAR)的Unidata项目科学家针对科学数据的特点开发的，是一种面向数组型并适于网络共享的数据的描述和编码标准。对程序员来说，它和zip、jpeg、bmp文件格式类似，都是一种文件格式的标准。netcdf文件开始的目的是用于存储气象科学中的数据，现在已经成为许多数据采集软件的生成文件的格式。利用NetCDF可以对网格数据进行高效地存储、管理、获取和分发等操作。由于其灵活性，能够传输海量的面向阵列（array-oriented）数据，目前广泛应用于大气科学、水文、海洋学、环境模拟、地球物理等诸多领域，例如，NCEP(美国国家环境预报中心)发布的再分析资料，NOAA的CDC(气候数据中心)发布的海洋与大气综合数据集(COADS)均采用NetCDF作为标准。

NetCDF的数据结构

  NetCDF数据集(文件名后缀为.nc) 的格式不是固定的，它是使用者根据需求 自己定义的。一个NetCDF数据集包含维(dimensions)、变量(variables)和属性(attributes)三种描述类型，每种类型都会被分配一个名字和一个ID，这些类型共同描述了一个数据集，NetCDF库可以同时访问多个数据集，用ID来识别不同数据集。变量存储实际数据，维给出了变量维度信息，属性则给出了变量或数据集本身的辅助信息属性，又可以分为适用于整个文件的全局属性和适用于特定变量的局部属性，全局属性则描述了数据集的基本属性以及数据集的来源。一个NetCDF文件的结构包括以下对象：

1
2
3
4
5
6
NetCDF name｛
Dimensions：…  //定义维数
Variables：… //定义变量
Attributes：… //属性
Data：…//数据
｝

NetCDF数据的主要特点

1) 自描述性：
  它是一种自描述的二进制数据格式，包含自身的描述信息；
2) 易用性：
  它是网络透明的，可以使用多种方式管理和操作这些数据；
3) 高可用性：
  可以高效访问该数据，在读取大数据集中的子数据集时不用按顺序读取，可以直接读取需要访问的数据；
4) 可追加性：
  对于新数据，可沿某一维进行追加，不用复制数据集和重新定义数据结构；
5) 平台无关性：
  NetCDF数据集支持在异构的网络平台间进行数据传输和数据共享。可以由多种软件读取并使用多种语言编写，其中包括C语言，C++，Fortran，IDL，Python，Perl和Java语言等。

实验内容

使用数据

1.ECMWF简介
  欧洲中期天气预报中心（European Centre for Medium-Range Weather Forecasts 简称 ECMWF）是一个包括24个欧盟成员国的国际性组织，是当今全球独树一帜的国际性天气预报研究和业务机构。其前身为欧洲的一个科学与技术合作项目。1973年有关国家召开大会宣布ECMWF正式成立，总部设在英国的Bracknell。
该中心于1979年6月首次做出了实时的中期天气预报。自1979年8月１日中心开始发布业务性中期天气预报，为其成员国提供实时的天气预报服务。
  此外，ECMWF还与克罗地亚、冰岛、匈牙利和斯洛文尼亚达成了合作协议，与世界气象组织（WMO）、欧洲气象卫星开发组织（EUMETSAT）和非洲气象应用发展中心（ACMAD）也有业务工作上的具体协议。同时ECMWF与世界各国气象预报机构在天气预报领域有广泛的联系。ECMWF主要提供10天的中期数值预报产品，各成员国通过专用的区域气象数据通信网络得到这些产品后做出各自的中期预报，同时ECMWF也通过由世界气象组织（WMO）维护的全球通信网络向世界所有国家发送部分有用的中期数值预报产品。其使用的模式充分利用四维同化资料，可提供全球在65公里高度内60层的40公里网格密度共20，911，680个点的风、温、湿预报。
2.数据下载
  本次实习使用数据是从ECMWF官网(http://data-portal.ecmwf.int/data/d/era40_daily/ )下载的NetCDF文件，主要内容为1965年全球海表压强(Mean sea level pressure)和皮温(skin temperature)。

并读取数据

1.获取文件信息
  matlab中使用ncdisp命令查看NetCDF文件信息。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
>> ncdisp('1965.nc')
Source:
           C:\Users\Jiangyayu\Documents\MATLAB\1965.nc
Format:
           classic
Global Attributes:
           Conventions = 'CF-1.0'
           history     = '2013-09-04 04:50:02 GMT by mars2netcdf-0.92'
Dimensions:
           longitude = 144
           latitude  = 73
           time      = 1460  (UNLIMITED)
Variables:
    longitude
           Size:       144x1
           Dimensions: longitude
           Datatype:   single
           Attributes:
                       units     = 'degrees_east'
                       long_name = 'longitude'
    latitude 
           Size:       73x1
           Dimensions: latitude
           Datatype:   single
           Attributes:
                       units     = 'degrees_north'
                       long_name = 'latitude'
    time     
           Size:       1460x1
           Dimensions: time
           Datatype:   int32
           Attributes:
                       units     = 'hours since 1900-01-01 00:00:0.0'
                       long_name = 'time'
    msl      
           Size:       144x73x1460
           Dimensions: longitude,latitude,time
           Datatype:   int16
           Attributes:
                       scale_factor  = 0.222
                       add_offset    = 9.95e+04
                       _FillValue    = -3.28e+04
                       missing_value = -3.28e+04
                       units         = 'Pa'
                       long_name     = 'Mean sea level pressure'
    skt      
           Size:       144x73x1460
           Dimensions: longitude,latitude,time
           Datatype:   int16
           Attributes:
                       scale_factor  = 0.00238
                       add_offset    = 267
                       _FillValue    = -3.28e+04
                       missing_value = -3.28e+04
                       units         = 'K'
                       long_name     = 'Skin temperature'

2.读取数据
（1）读取变量

1
2
3
4
5
6
7
%% Read data from the netcdf file
%Read Variables 
longitude_data = ncread('1965.nc','longitude');
latitude_data = ncread('1965.nc','latitude');
time_data = ncread('1965.nc','time');
msl = ncread('1965.nc','msl');
skt = ncread('1965.nc','skt');

（2）读取属性

1
2
3
4
5
6
7
8
9
10
% Read Attributes
msl_FillValue = ncreadatt('1965.nc','msl','_FillValue');
msl_missing_value = ncreadatt('1965.nc','msl','missing_value');
msl_scale_factor = ncreadatt('1965.nc','msl','scale_factor');
msl_add_offset = ncreadatt('1965.nc','msl','add_offset');
 
skt_FillValue = ncreadatt('1965.nc','skt','_FillValue');
skt_missing_value = ncreadatt('1965.nc','skt','missing_value');
skt_scale_factor = ncreadatt('1965.nc','skt','scale_factor');
skt_add_offset = ncreadatt('1965.nc','skt','add_offset');

（3）读取结果

（4）计算年平均msl和skt

绘制等值线图和伪彩色图

(1)contour绘制等值线图

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
%% figure
AX = figure(1);
% Figure Mean sea level pressure
AX1 = subplot(211);
v1 = min(min(ave_msl))/100:2.5:max(max(ave_msl))/100;
[c1,h1] = contour(longitude_data,latitude_data,ave_msl'/100,v1);
axis image
xlabel('longitude');
ylabel('latitude');
title('1965 Mean sea level pressure (+1000hpa)');
clabel(c1,h1,v1([1,2,4,6,9,11,13,14,15,16]),'LabelSpacing',1000);
ht = clabel(c1,h1);
for ii = 1:length(ht)
    num = get(ht(ii),'UserData');
    set(ht(ii),'String',num2str(num-1000,'%.1f'))
end

(2)pcolor绘制伪彩色图

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
%Figure Skin temperature
AX2 = subplot(212); 
h2 = pcolor(longitude_data,latitude_data,ave_skt');
axis image
xlabel('longitude');
ylabel('latitude');
title('1965 Average Skin temperature (K)');
shading interp;
3)	 绘制colorbar
% colorbar
h = colorbar;
set(h,'Position',[0.85,0.12,0.05,0.8]);
v2 = min(min(ave_skt)):10:max(max(ave_skt));
v3 = cell(1,length(v2));
for ii = 1:length(v2)
    v3(ii) = {num2str(v2(ii),'%.1f')};
end
set(h,'YTick',v2,'YTickLabel',v3)
set(get(h,'Title'),'string','K')

(4) 保存图片

1
2
% save the figure 
saveas(gcf,'msl&skt','jpg')

(5)结果图

将数据写入新文件

(1)创建NetCDF文件和变量

1
2
3
4
5
6
%% Rewrite the data to a new netcdf file
%Create File and Variables
nccreate('average_1965.nc','longitude','Dimensions',{'longitude',144},'Datatype','single','Format','classic');
nccreate('average_1965.nc','latitude','Dimensions',{'latitude',73},'Datatype','single','Format','classic');
nccreate('average_1965.nc','msl','Dimensions',{'longitude',144,'latitude',73},'Datatype','int16','Format','classic');
nccreate('average_1965.nc','skt','Dimensions',{'longitude',144,'latitude',73},'Datatype','int16','Format','classic');

(2)创建属性值

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
%Attributes
fileattrib('average_1965.nc','+w');
%Global Attributes
ncwriteatt('average_1965.nc','/','discription','Mean sea level pressure and skin temperature of 1965');
ncwriteatt('average_1965.nc','/','author','Yayu Jiang');
ncwriteatt('average_1965.nc','/','data','20130910');
 
% Variables Attributes
ncwriteatt('average_1965.nc','longitude', 'units','degrees_east');
ncwriteatt('average_1965.nc','longitude', 'long_name' ,'longitude');
 
ncwriteatt('average_1965.nc','latitude', 'units','degrees_north');
ncwriteatt('average_1965.nc','latitude', 'long_name' ,'latitude');
 
ncwriteatt('average_1965.nc','msl', 'scale_factor',msl_scale_factor);
ncwriteatt('average_1965.nc','msl', 'add_offset',msl_add_offset);
ncwriteatt('average_1965.nc','msl', '_FillValue',msl_FillValue);
ncwriteatt('average_1965.nc','msl', 'missing_value',msl_missing_value);
ncwriteatt('average_1965.nc','msl', 'units','K');
ncwriteatt('average_1965.nc','msl', 'long_name','Skin temperature');
 
ncwriteatt('average_1965.nc','skt', 'scale_factor',skt_scale_factor);
ncwriteatt('average_1965.nc','skt', 'add_offset',skt_add_offset);
ncwriteatt('average_1965.nc','skt', '_FillValue',msl_FillValue);
ncwriteatt('average_1965.nc','skt', 'missing_value',msl_missing_value);
ncwriteatt('average_1965.nc','skt', 'units','Pa');
ncwriteatt('average_1965.nc','skt', 'long_name','Mean sea level pressure');

(3) 写入数据

1
2
3
4
5
% input the new data
ncwrite('average_1965.nc','longitude',longitude_data);
ncwrite('average_1965.nc','latitude',latitude_data);
ncwrite('average_1965.nc','msl',ave_msl);
ncwrite('average_1965.nc','skt',ave_skt);

Surfer绘图

（1)创建dat文件

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
%% creat data for surfer
m = zeros(144*73,3);
s = zeros(144*73,3);
for ii = 1:144
    for jj = 1:73
        m((ii-1)*73+jj,1) = longitude_data(ii);
        m((ii-1)*73+jj,2) = latitude_data(jj);
        m((ii-1)*73+jj,3) = ave_msl(ii,jj);
        s((ii-1)*73+jj,1) = longitude_data(ii);
        s((ii-1)*73+jj,2) = latitude_data(jj);
        s((ii-1)*73+jj,3) = ave_skt(ii,jj);
    end
end
csvwrite('msl.dat',m);
csvwrite('skt.dat',s);

(2)网格化数据
  打开Surfer 8软件，通过菜单下网格->数据，将matlab生成的dat文件通过插值方式变成网格数据，由于在matlab中已考虑行列的关系，因此这里的数据列就不用考虑了，数据点数也比较多，网格化方式也可默认。

(3) 绘制等值线
  在菜单栏下选择地图->等值线图->新建等值线图，依次导入刚刚生成的msl.grd和skt.grd网格数据。
(4) 设置参数
  对skt等值线图作如下设置：
  勾选填充等值线选项；将等值线填充成由篮到红的渐变；设置线条为不可见；将右侧和上侧的刻度属性设置为无。


  最终效果图：

  对msl等值线进行颜色设置标注显示等最终效果：

  添加全球大陆的轮廓图，设置填充：

  选中这三个图层，在菜单下选择：地图->覆盖地图，将三幅图覆盖到一起，调整合适的图层顺序使显示正常。
  在菜单下选择：绘制->符号，找到指北针符号，为等值线图添加指北针。
  点击工具栏中的T，为地图添加标题、坐标轴名称等基本信息。

(5) 添加colorbar
  根据画skt等值线图的基本思想，要制作一个colorbar可以通过设置一个从左往右均匀渐变的网格，然后填充颜色的做法来绘制。这样绘制的colorbar其实就是一种特殊的等值线。
  同样用matlab生成这样的一个dat文件：
1
2
3
4
5
6
7
8
9
10
a = 0:2.55:255;
m = zeros(1010,3);
for ii = 1:10
    for jj = 210:310
        m((ii-1)*101+jj-209,1) = jj;
        m((ii-1)*101+jj-209,2) = ii;
        m((ii-1)*101+jj-209,3) = a(jj-209);
    end
end
csvwrite('colorbar.dat',m);

  然后像skt图一样设置填充、坐标轴、等值线、刻度线等，使之美观。

(6)保存文件

小结

matlab等值线标注问题

  用Matlab中contour命令画等值线的等值线虽然设置了每隔2.5hpa画一条，但min（min（ave_msl））和max（max（ave_msl））这中间还有16个数据，画出来的等压线十分密集。因此我们在使用clabel进行标注时可以有选择性地对等值线密集的地方少一些标注，稀疏的地方多一些标注，例如：

1
clabel(c1,h1,v1([1,2,4,6,9,11,13,14,15,16]),'LabelSpacing',1000);

  此外我们还可以设置标注数据的有效位数，并将其值共同减去一个较大的值来实现用较少位数的数进行标注，例如：

1
2
3
4
5
ht = clabel(c1,h1);
for ii = 1:length(ht)
    num = get(ht(ii),'UserData');
    set(ht(ii),'String',num2str(num-1000,'%.1f'))
end

matlab colorbar的标注问题

  matlab中的colorbar命令可以直接绘制出比较优美的colorbar，但是有时候为了是设置colorbar的标注，我们需要用到命令：

1
set(h,'YTick',v2,'YTickLabel',{’string1’,’strng2’,’string3’});

  当colorbar上的标注比较少的时候这种做法可行，但是当colorbar上的标注达到一定的数量时，这种做法很明显会比较耗时耗力，我们可以尝试用：

1
2
3
4
5
6
v2 = min(min(ave_skt)):10:max(max(ave_skt));
v3 = cell(1,length(v2));
for ii = 1:length(v2)
    v3(ii) = {num2str(v2(ii),'%.1f')};
end
set(h,'YTick',v2,'YTickLabel',v3);

  此时就能按照所要求的标注间隔和所要保留的位数进行自动标注了。

surfer 8绘制colorbar的方法

  相比较matlab而言，surfer中没有提供直接画colorbar的方便操作，但是我们还是可以通过已有的功能来实现colorbar的制作。

(1)首先我们要做的是生成用来制作colorbar的数据，我们可以在matlab中生成这样的一组数据：
  x：将来作为colorbar的刻度值，应设置为skt的整个范围
  y：colorbar的宽度，可以任意设置，将来仍然可以调节
  (x，y):0:255要保证维度和x，y个数，且不同值的数目与x相同

(2) 用matlab实现生成数据
  这里为了后面绘图方便，在生成数据的时候就生成了一个等价的二维数组colorbar.dat

(3) 转换成网格数据并导

(4) 设置参数
  分别设置填充颜色、线条颜色、坐标轴标注、坐标轴刻度等属性，一个美观的colorbar就绘制出来了。

若本文对您有帮助，请打赏鼓励本人！
---------------- End ----------------
扫二维码
扫一扫，使用手机查看

扫一扫，使用手机查看

发送邮件
QQ