[Matlab] Matlab教程:nc文件的打开和使用m_map绘制海温图

[复制链接]

在大气科学中,matlab可以用于小规模的科学计算,也可以绘制各类气象图,做各种统计运算,功能强大。由于疫情关系,顺便写写使用m_map和matlab读取和绘制海温图的教程。笔者花费了一天的时间学习了m_map(入门,熟悉下语法),遇到一些难点,希望通过我对这个例子的分解,能化简入门难度,从而能帮助到有需要的人。

如果没有安装过m_map,则需要额外的安装,可以参考这个帖子:

游客,如果您要查看本帖隐藏内容请回复

总体思路是:输入路径,文件名,读取范围,时间层,变量名等,调用m_map绘制

1)输入路径,文件名,读取范围

input('string','s')输入函数,string的内容会打印在屏幕,有 's' 的情况下是输入字符串,没有的情况下输入数字

游客,如果您要查看本帖隐藏内容请回复

2)打开NC文件,显示NC文件信息

strcat(a, b, ......) 合并字符串a,b......,如果直接打上去的字符串记得用引号括起来;例如:strcat('app', 'ple')输出就是apple

ncdisp( 'source' ) 显示NC文件信息,其中source是文件位置

%文件打开与数据载入source1 = strcat(InPath,file_name);  %文件源,strcat函数:合并字符串ncdisp(source1);                     %查阅和显示NC文件信息
' p/ H3 Y- s. B0 v" t. t. h

3)读取NC文件的精度,时间层数等,方便后面的计算

同时需要注意的是,不同的NC文件对经纬度和时间层数的表述并不一致,一些NC文件会将longitude简写为lon;latitude简写为lat;时间层数为其他的表述,这个需要我们注意并适当修改;以下代码是ECMWF的描述:

length( 矩阵名 )   查询矩阵长度lon = ncread(source1,'longitude');   %查阅经度信息loncount = length(lon);              %查阅经度的精度(有多少格点)lat = ncread(source1,'latitude');    %查阅纬度信息latcount = length(lat);              %查阅纬度精度(有多少格点)time = ncread(source1,'time');       %查阅时间层数信息ticount = length(time);              %查阅时间层数
9 N0 L0 T! G# |& ^. g

4)显示时间层数,避免超出NC文件的范围,并且根据需要选择变量

disp(  )     把括号内的内容打印在屏幕上disp('时间层数为:')disp(ticount);                       %显示时间层数t = input('输入绘制的时间层:');varname = input('输入变量:','s');    %根据ncdisp显示的变量输入绘图6 |' E+ C% g6 C) ~$ N/ t

5)我们在前面输入了我们要绘制的经纬度,这个语句是对我们规定经纬度进行查询,把绘图区域规定在一定的范围内,得出对应的矩阵区域(matlab中数据可以看成一个n维矩阵)

find( 条件 ) 查询函数,查询对应的矩阵位置

%查找绘制范围对应的所在矩阵的位置,相当于截取一小段矩阵lon_scope = find(lon >= boundary(1) & lon<=boundary(2));lat_scope = find(lat >= boundary(3) & lat<=boundary(4));lon_number = length(lon_scope);     lat_number = length(lat_scope);
. c, D0 y8 Z0 @# T# S6 C: A# U

6)得出对应的矩阵,使用ncread语句读取NC文件,并储存在变量中

ncread(source,varname ,start,count,stride)

其中:source:文件路径

varname:要绘制的变量名

start:初始读取位置

count:读取范围

stride:读取步长

imrotate( 矩阵或图像名, 旋转角度n ) 将图像所对应的矩阵旋转n度,要注意的是正方向是逆时针方向;

start = [lon_scope(1),lat_scope(1),1];    %初始位置count = [lon_number,lat_number,ticount];  %读取范围 stride1 = [1,1,1];                        %读取步长sst1 = ncread(source1,varname ,start,count,stride1);sst_plot = imrotate(sst1(:,:,t), 90);    %旋转矩阵,因为matlab是列优先
7 B" p& L( H; \9 |4 R+ o

7)设置投影方式,绘制范围

m_proj('投影方式', 'lat', [纬度范围], 'lon', [经度范围])

其他的投影方式参见官方文档

m_proj('Mercator','lat',[boundary(3) boundary(4)],'lon',[boundary(1) boundary(2)]);7 T' ^* k$ F+ t* x: q

8)生成网格

linspace(x1, x2, N) x1 开始数值;x2 结束数值;N生成个数,生成从x1开始到x2结束的N个均匀分布的行向量

meshgrid( x, y ) 将向量x和向量y规定的区域(类似于经纬度)生成格点,比如x是经度向量,y是纬度向量的话,将会生成一个以x为经度,y为纬度的坐标网格。

lat_1=linspace(boundary(3),boundary(4),lat_number);lon_1=linspace(boundary(1),boundary(2),lon_number); [plon,plat]=meshgrid(lon_1,lat_1);hold on" P; S8 H4 l7 w  n! Y- Q3 s

9)绘制内容

m_pcolor(lon, lat, varname) 这个函数的功能是把要绘制的格点数据,绘制在地图上,样式是伪色彩图。其中:lon上个步骤中meshgrid函数生成的经度格点;lat上个步骤meshgrid函数生成的格点(应该说是句柄比较合适),也就是meshgrid等号前的2个句柄。

m_coast( ) 这个函数是绘制海岸线,并规定海岸线颜色,陆地填充颜色等。(这个代码我忘了添加陆地颜色了)

m_grid( ) 绘制边框,其具体设置见官方文档,代码中,('box','fancy')指的是绘制外边框,选择fancy样式,这时候默认是黑白交替的样式

m_pcolor(plon,plat,sst_lot)             %添加我们要画的内容m_coast('color',[0 0 0],'linewidth',2);  %绘制海岸线,填充陆地m_grid('box','fancy')                    %添加边框hold on
! M: v  y! t* `6 {, o) T

10)添加标题和色标

title('标题文字', 'fontsize', 15) fontsize后面的数字是设置字体大小

colorbar('h') 添加色标

下面的set 设置色标属性

title('SST','fontsize',15) %标题%添加色标h = colorbar('h');set(get(h,'title'),'string','摄氏度℃');hold on
% U/ c/ y7 R# i0 Z" o8 ~2 v

我们对其进行编译,运行:


0 u. B# w+ |1 W. l9 K7 N* ^                               
登录/注册后可看大图

输入路径时结尾记得加上“ \ ”符号,输入文件名记得加上后缀名,输入范围时经纬度用空格区分,输入时记得加上中括号[ ],我们绘制西太平洋的SST:

其中输入的变量名可以在上面输出的NC文件信息中找到

4 G% a' }0 H7 I( D" i. n6 P& i
                               
登录/注册后可看大图

输出结果


0 i( p% e# K1 d! y) @7 J! h                               
登录/注册后可看大图

由于本文较为仓促,仅学了一天m_map,难免会有疏漏;笔者代码水平还属于入门中的入门水平,有很多代码写的时候并不规范,希望引以为戒!水平较低,希望前辈们指出不足的地方

附完整代码:

游客,如果您要查看本帖隐藏内容请回复

. t; P5 u, k9 Q$ b( c- U) f
回复

举报 使用道具

相关帖子

全部回帖
感谢分享
发表于 2022-4-18 22:01:02

举报 回复 使用道具

懒得打字?点击右侧快捷回复 【吾爱海洋论坛发文有奖】
您需要登录后才可以回帖 登录 | 立即注册
杨小燕
活跃在2022-11-2
快速回复 返回顶部 返回列表