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

thinkfunny

一、matlab读取NCEP再分析数据并绘制风场
; \" L2 g. D2 u2 V) y5 O
, T4 g% j$ R  o%该程序用于求水汽通量散度
%注意！我们这里读到的u是四维矩阵，分别是lon*lat*level*time，
; m3 ^. ^) A2 A6 M/ z7 Yclc;clear;close all
+ ]. @9 J7 |# b# t# jf_hgt = 'ps_level_20170121_0130.nc';
% ncdisp(f_hgt);
time=ncread(f_hgt,'time');
level=ncread(f_hgt,'level');
& W0 b" Y2 J' elon=ncread(f_hgt,'longitude');
lat=ncread(f_hgt,'latitude');
%%%%%%时间转换
time  = double(time);
, @+ v; f% a1 q6 G! W/ I5 Vformat = 'mm dd, yyyy HH:MM:SS.FFF AM';%转换格式
' G7 W$ l5 _/ c7 G& j* Bdstr = datestr((datenum('1900-01-01') + time./24),format);%转换后时间字符串存储
, {% s% f# t4 q6 f# e/ c$ e/ ITM = datevec(dstr);%将时间字符数组转化为数值数组
! ?: k# D0 }6 Wtidx=find(TM(:,2)==1 & TM(:,3)==28 & TM(:,4)==00);%筛选7月25日08时（世界时加8）
  b: i1 i8 _& j! vps_lev=find(level ==850);%%删选出850hPa高度
start=[1,1,ps_lev,tidx];%所指定变量的每一维的开始读取的位置
% Q& S0 b# q* K1 F6 W' x7 I& rcount=[41,31,1,1];%从start指定的开始位置算起，一共读取的每一维要素的数目
$ U3 J3 V) D2 |- xstrip=[1,1,1,1];%从start开始，每一维读取的数目为count时，每一维的读取的步长
5 G0 T# T' Z' P' z$ I4 G. x- {) \* Y, phgt=ncread(f_hgt,'z',start,count,strip);%读取温度值，单位K
0 S. S; C& y0 P1 A! v5 o, Eu=ncread(f_hgt,'u',start,count,strip);%读取温度值，单位K
v=ncread(f_hgt,'v',start,count,strip);%读取温度值，单位K
[X,Y]=meshgrid(lon,lat);
* v: c( A8 C8 ?( p& V4 y: afigure(1)
6 k( ^8 j; \) y7 I! F7 k% Q; im_proj('Mercator','lat',[25,35],'lon',[100,115]);
% m_grid('linestyle','none','tickdir','out','fontsize',12,'fontname','Times New Roman');
4 ~1 d2 J. v6 O- Bm_grid('linestyle','none','box','fancy','fontsize',11,'tickdir','in','xtick',[100:3:115],'ytick',[25:2:35]);
hold on
0 ^# ^9 D2 A- J! T5 _; N: D- Qm_windbarb(X',Y',u,v,'color','k')
$ B2 ]0 Q! }6 w- s/ y: v* K%m_coast('patch',[.9 .9 .9],'edgecolor','none');
0 ?4 @- D* z. l[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');   %在等高线上叠加数值（文后详情）
; ?' i$ F. U9 v' I9 Gclabel(C,h,'FontSize',13,'fontname','Times New Roman');
ma=shaperead('F:/RMeteoInfo/data/map/bou2_4l.shp');
# ]# t1 N( u, v' o( X; W% 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')
, a+ M! e  p( _* U

  G3 s- b  J' D! c' d

& g9 l4 r7 a. g二、matlab读取ERA5并绘制全球风场图
: p/ N! n. h4 S
clc;clear;close all
& u3 c0 j/ a$ H" yu0=ncread('202008muwind.nc','uas');    %读取其中一项
v0=ncread('202008mvwind.nc','vas');
3 j% m( `0 R* W3 D0 F* e5 ?1 Etime0=ncread('202008muwind.nc','time');
lat0=ncread('202008muwind.nc','lat');
lon0=ncread('202008muwind.nc','lon');      
# x3 L: y) d$ M; [[lat,lon]=meshgrid(lat0,lon0);
u=u0(:,:,2);
v=v0(:,:,2);
4 [# M' m$ L6 F) }" ~( Mb=sqrt(u.^2+v.^2);
figure(1);
0 |: ^2 ^  T/ |8 }$ M# Z' Gm_proj('Equidistant Cylindrical','long',[-180 180],'lat',[-90 90]);%矩形投影；取区域观察
1 z& C7 J# ~0 T9 ^0 y* q[lon1,lat1]=meshgrid([-180:2.5:179.75],[-90:2.5:90]);
; Q8 f- S* u1 Z8 J4 z5 C/ d& Zu2=u(1:10:end,1:10:end)';
% u3=[u2,u2(:,end)];
5 x5 m) K# ]) g% T5 Ev2=v(1:gap:end,1:gap:end)';
% v3=[v2,v2(:,end)];
m_quiver(lon1,lat1,u2,v2,0);
. Z+ d( F( i" n) s8 V# phold off
m_coast('patch',[1 .85 .7]);
  ^9 |; x6 P. hm_grid('box','on','tickdir','out');

3 `3 z$ A4 x  v& L6 Z
. Q6 t7 u* K$ U  F$ Q  M+ p7 R% A


5 F% _6 |8 D6 `3 k$ Q( ~9 ~m_contourf(lon,lat,b);%在mmap基础上的画
& N8 T( D  ^; |; X  Tshading interp;%使数据插值
hold on;
m_quiver(lon1,lat1,u2,v2,0);
. Z$ j$ x' j  E, l7 ohold off;

3 G4 U- I1 K9 P- t% u1 t

7 e; v7 ~- X8 `- z

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

2 C+ y" z" c* \  ^0 {
. q1 C; Q0 Z3 v* v" @- E5 b, ^
$ x, _) w( i/ H0 `. X! X

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

, Z' s$ Z4 ~3 y3 h. D

- g% u- K0 F8 U. o8 J0 s$ g* @" ^
' Z$ _2 L1 e8 x$ M) ]7 O( ~& P

SsaaM7

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