Matlab全球风场图—以ERA5、NCEP和ECMWF再分析数据为例

thinkfunny · 发表于 2021-3-3 09:50

EDA365欢迎您登录！

您需要登录才可以下载或查看，没有帐号？注册

x

一、matlab读取NCEP再分析数据并绘制风场

%该程序用于求水汽通量散度
%注意！我们这里读到的u是四维矩阵，分别是lon*lat*level*time，
clc;clear;close all
f_hgt = 'ps_level_20170121_0130.nc';
% ncdisp(f_hgt);
time=ncread(f_hgt,'time');
level=ncread(f_hgt,'level');
lon=ncread(f_hgt,'longitude');
lat=ncread(f_hgt,'latitude');
%%%%%%时间转换
time = double(time);
format = 'mm dd, yyyy HH:MM:SS.FFF AM';%转换格式
dstr = datestr((datenum('1900-01-01') + time./24),format);%转换后时间字符串存储
TM = datevec(dstr);%将时间字符数组转化为数值数组
tidx=find(TM(:,2)==1 & TM(:,3)==28 & TM(:,4)==00);%筛选7月25日08时（世界时加8）
ps_lev=find(level ==850);%%删选出850hPa高度
start=[1,1,ps_lev,tidx];%所指定变量的每一维的开始读取的位置
count=[41,31,1,1];%从start指定的开始位置算起，一共读取的每一维要素的数目
strip=[1,1,1,1];%从start开始，每一维读取的数目为count时，每一维的读取的步长
hgt=ncread(f_hgt,'z',start,count,strip);%读取温度值，单位K
u=ncread(f_hgt,'u',start,count,strip);%读取温度值，单位K
v=ncread(f_hgt,'v',start,count,strip);%读取温度值，单位K
[X,Y]=meshgrid(lon,lat);
figure(1)
m_proj('Mercator','lat',[25,35],'lon',[100,115]);
% m_grid('linestyle','none','tickdir','out','fontsize',12,'fontname','Times New Roman');
m_grid('linestyle','none','box','fancy','fontsize',11,'tickdir','in','xtick',[100:3:115],'ytick',[25:2:35]);
hold on
m_windbarb(X',Y',u,v,'color','k')
%m_coast('patch',[.9 .9 .9],'edgecolor','none');
[C,h]=m_contour(X',Y',hgt/10 ,'color','k','LineWidth',1);%[5000:80:5900]%[1520:25:1680]
%set(b,'ShowText','on','TextStep',get(b,'LevelStep'),'fontsize',12,'fontname','Times New Roman'); %在等高线上叠加数值（文后详情）
clabel(C,h,'FontSize',13,'fontname','Times New Roman');
ma=shaperead('F:/RMeteoInfo/data/map/bou2_4l.shp');
% m_line( [ma(：).X], [ma(：).Y],'color',[0.5,0.5,0.5]);%绘制范围内的地图
m_line([ma(：).X],[ma(：).Y],'color','r');%绘制范围内的地图
m_plot(105.5,29.43,'marker','^','MarkerSize',7,'color','k','MarkeRFaceColor','k')

二、matlab读取ERA5并绘制全球风场图

clc;clear;close all
u0=ncread('202008muwind.nc','uas'); %读取其中一项
v0=ncread('202008mvwind.nc','vas');
time0=ncread('202008muwind.nc','time');
lat0=ncread('202008muwind.nc','lat');
lon0=ncread('202008muwind.nc','lon');
[lat,lon]=meshgrid(lat0,lon0);
u=u0(:,:,2);
v=v0(:,:,2);
b=sqrt(u.^2+v.^2);
figure(1);
m_proj('Equidistant Cylindrical','long',[-180 180],'lat',[-90 90]);%矩形投影；取区域观察
[lon1,lat1]=meshgrid([-180:2.5:179.75],[-90:2.5:90]);
u2=u(1:10:end,1:10:end)';
% u3=[u2,u2(:,end)];
v2=v(1:gap:end,1:gap:end)';
% v3=[v2,v2(:,end)];
m_quiver(lon1,lat1,u2,v2,0);
hold off
m_coast('patch',[1 .85 .7]);
m_grid('box','on','tickdir','out');

m_contourf(lon,lat,b);%在mmap基础上的画
shading interp;%使数据插值
hold on;
m_quiver(lon1,lat1,u2,v2,0);
hold off;

m_quiver(lon1,lat1,u2,v2,0);

m_contourf(lon,lat,b,500,'linestyle','none') %在mmap基础上的画
shading interp; %使数据插值

SsaaM7 · 发表于 2021-3-3 10:52

Matlab全球风场图—以ERA5、NCEP和ECMWF再分析数据为例

帐号		自动登录	找回密码
密码			注册

Matlab全球风场图—以ERA5、NCEP和ECMWF再分析数据为例

EDA365欢迎您登录！

浏览过的版块

推荐内容 /1