MATLAB画风场图

目录

1.风场图简介

2.数据来源

3.代码实现

3.1 读取数据

---

1.风场图简介

如下是我用MATLAB绘制的一张风场图，它包含以下几个要素：

1.需要有正确的地图投影

2.需要有箭头表示风向

3.需要用色阶来区分风力大小

接下来本文将着手介绍，如何使用MATLAB画出上图这样的风场图。

2.数据来源

风场数据来源于 NOAA 网站提供的混合风场数据（NOAA NCEI Blended Seawinds，NBS)

下载链接：NOAA NCEI Blended Seawinds (NBS v2) | NOAA CoastWatch

我是在这个模块下载的数据

 

 本文使用的数据如下：

NBSv02_wind_monthly_201508.nc

NBSv02_wind_monthly_201601.nc

3.代码实现

3.1 读取数据

filepath = 'C:\Users\98485\Desktop\Academic\研一课程\海洋观测\海洋观测作业\DATA';
filename = 'NBSv02_wind_monthly_201508.nc';
path = fullfile(filepath,filename);
ncdisp(path);

ncdisp（）命令的作用是展示nc文件的内容，通过此命令我们可以看到nc文件中包含哪些变量，以便于下一步的变量提取。下图是该命令的执行结果：

可以看到风场数据被分解成了：u方向分量、v方向分量、风速，接下来我们用ncread（）函数将这些变量读取出来，代码如下：

u = ncread(path,'u_wind');u = u';
v = ncread(path,'v_wind');v = v';
speed = ncread(path,'windspeed');speed = speed';

注意：一定要将u、v、speed转置，才能和经纬度网格对应上！

接下来读取经纬度信息，并生成网格：

lon = ncread(path,'lon'); 
lat = ncread(path,'lat');
[LON LAT] =meshgrid(lon,lat);

3.2 绘制风速的伪彩图

figure();
subplot(121);
m_proj('Miller','lon',[100  300],'lat',[-70 70]);
m_pcolor(lon,lat,speed);
shading flat;colormap(parula);

首先选择地图投影为“米勒投影”，并指定投影的经度范围是100°E~60°W，纬度范围是70°S~70°N

注意：MATLAB中的经度取值范围0~360°，纬度-90~90°，经度0~180°是东经，180~360°是西经，例如要确认MATLAB中的300°是西经多少度，只需要将300-360= - 60，即西经60°

m_pcolor（x,y,z）函数用来绘制伪彩图，其中x，y参数就像直角坐标系中的x，y轴一样，用于确定色彩点的位置，z参数的值决定这个色彩点的颜色深浅。x，y，z一般都是矩阵，因此里面有很多个数，也就是说能点出很多个色彩点，这就是伪彩图的绘制原理，由此可见，x，y，z三个矩阵的维度必须一一对应。

h=colorbar('h');
set(get(h,'xlabel'),'string','Wind speed(m/s)','fontsize',14);

最后再加上色阶，就得到了一张这样的伪彩图：

3.3 绘制风向箭头

上文中我们已经完成了风速伪彩图的绘制，接下来我们需要在这张图上加上箭头来表示风向。

hold on;
m_quiver(LON(1:100:end),LAT(1:100:end),u(1:100:end),v(1:100:end),'k','linewidth',0.5,'autoscalefactor',1.6);

首先，我们需要先用hold on语句，来表示接下来的操作（即绘制风向箭头）需要覆盖在刚刚画好的风速伪彩图上。如果不使用hold on语句，MATLAB就会将伪彩图擦除，然后在新的空白画布上绘制风向箭头图（见下图）。

m_quiver(x,y,u,v)函数用来绘制箭头，x，y参数用于指定箭头起点位置，u，v参数用于确定箭头的方向，'linewidth'，0.5表示箭头线宽0.5，'autoscalefactor',1.6用于调整箭头自动缩放的比例，即调整箭头长度。

值得注意的是，我将经度、纬度、u分量、v分量的步长都设置成了100，这是因为我不想箭头过密，以确保图的美观性。

最后加上地形投影、网格、标题，就得到了完整的作品，代码如下：

m_gshhs_c('patch',[0.7 0.7 0.7],'Edgecolor','none');
m_grid('box','off','tickdir','in'); 
title('10m Wind speed in Equatorial Pacific,Jan 2023','fontsize',18);

这样我们就得到了一幅风场图

3.3 地图标点 

最后分享一段在地图上标点的代码，当我们需要说明采样点位于地图何处时，我们可以通过以下代码来实现：

ax = worldmap('world');   %选择区域为‘世界’
setm(ax,'Origin',[0 180 0]);
land = shaperead('landareas', 'UseGeoCoords', true); %导入陆地框架
geoshow(ax, land, 'FaceColor', [0.5 0.7 0.5]); %展示地图
plotm(13,200,'Marker','.','markersize',20,'markeredgecolor','r')

plotm（x，y）函数：x，y代表采样点的经纬度，'marker'，‘.’ 参数代表用“点”表示，'markersize',20 表示“点”的大小为20，'markeredgecolor','r' 表示点的颜色为红色，这样就可以得到清晰的地图标点：

4.完整代码

filepath = 'C:\Users\98485\Desktop\Academic\研一课程\海洋观测\海洋观测作业\DATA';
filename = 'NBSv02_wind_monthly_201508.nc';
path = fullfile(filepath,filename);
ncdisp(path);
lon = ncread(path,'lon'); 
lat = ncread(path,'lat');
[LON LAT] =meshgrid(lon,lat);
u = ncread(path,'u_wind');u = u';
v = ncread(path,'v_wind');v = v';
speed = ncread(path,'windspeed');speed = speed';
figure();
m_proj('Miller','lon',[100  300],'lat',[-70 70]);
m_pcolor(lon,lat,speed);
shading flat;colormap(parula);
h=colorbar('h');
set(get(h,'xlabel'),'string','Wind speed(m/s)','fontsize',14);
hold on;
m_quiver(LON(1:100:end),LAT(1:100:end),u(1:100:end),v(1:100:end),'k','linewidth',0.5,'autoscalefactor',1.6);
m_gshhs_c('patch',[0.7 0.7 0.7],'Edgecolor','none');
m_grid('box','off','tickdir','in'); 
title('10m Wind speed in Equatorial Pacific,Jan 2023','fontsize',18);

本文感谢江源同志的细致指导！