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

thinkfunny

$ N1 s- }5 k% @  P) _4 d一、matlab读取NCEP再分析数据并绘制风场

%该程序用于求水汽通量散度
%注意！我们这里读到的u是四维矩阵，分别是lon*lat*level*time，
+ V5 e7 r9 T& r/ s8 b5 }6 p6 xclc;clear;close all
f_hgt = 'ps_level_20170121_0130.nc';
% ncdisp(f_hgt);
time=ncread(f_hgt,'time');
2 B4 v8 K3 W8 V9 b5 `$ t' glevel=ncread(f_hgt,'level');
lon=ncread(f_hgt,'longitude');
lat=ncread(f_hgt,'latitude');
+ j- ^5 O2 S' a+ W5 }%%%%%%时间转换
) i3 `( v8 I" c' i5 D7 Vtime  = double(time);
format = 'mm dd, yyyy HH:MM:SS.FFF AM';%转换格式
dstr = datestr((datenum('1900-01-01') + time./24),format);%转换后时间字符串存储
TM = datevec(dstr);%将时间字符数组转化为数值数组
" O# |# {" v$ \0 T7 B5 R- Mtidx=find(TM(:,2)==1 & TM(:,3)==28 & TM(:,4)==00);%筛选7月25日08时（世界时加8）
$ G# j$ n/ u' H8 K% c" Wps_lev=find(level ==850);%%删选出850hPa高度
: O5 s8 T' B) G+ o& @2 Sstart=[1,1,ps_lev,tidx];%所指定变量的每一维的开始读取的位置
; x6 Y5 O  ], Wcount=[41,31,1,1];%从start指定的开始位置算起，一共读取的每一维要素的数目
strip=[1,1,1,1];%从start开始，每一维读取的数目为count时，每一维的读取的步长
  A# U8 w4 W" F- `. c4 Hhgt=ncread(f_hgt,'z',start,count,strip);%读取温度值，单位K
4 v- K$ O/ x4 Ou=ncread(f_hgt,'u',start,count,strip);%读取温度值，单位K
4 |$ U6 }3 X5 Yv=ncread(f_hgt,'v',start,count,strip);%读取温度值，单位K
[X,Y]=meshgrid(lon,lat);
figure(1)
8 @- \3 ], z, B# R( Z0 E9 mm_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]);
6 u8 W% t3 q5 A* f% i( U8 o' ohold on
2 k+ `8 z" U) fm_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');
! u( o+ h1 C9 `4 R( \+ B% m_line( [ma(：).X], [ma(：).Y],'color',[0.5,0.5,0.5]);%绘制范围内的地图
m_line([ma(：).X],[ma(：).Y],'color','r');%绘制范围内的地图
4 a$ _1 d- k: K; B# P+ e1 Zm_plot(105.5,29.43,'marker','^','MarkerSize',7,'color','k','MarkeRFaceColor','k')



  W' T( P, V( p
二、matlab读取ERA5并绘制全球风场图
- w9 T2 j2 N' N  T* U
clc;clear;close all
* C: z- i/ c- E" V' `u0=ncread('202008muwind.nc','uas');    %读取其中一项
8 p4 c' p2 t0 }3 j9 zv0=ncread('202008mvwind.nc','vas');
time0=ncread('202008muwind.nc','time');
lat0=ncread('202008muwind.nc','lat');
8 ?1 S6 A3 [% J5 w& R9 O% tlon0=ncread('202008muwind.nc','lon');      
- v6 ]0 {$ ?  S0 a[lat,lon]=meshgrid(lat0,lon0);
u=u0(:,:,2);
9 g' ~9 S4 `. ]( c, {v=v0(:,:,2);
b=sqrt(u.^2+v.^2);
figure(1);
m_proj('Equidistant Cylindrical','long',[-180 180],'lat',[-90 90]);%矩形投影；取区域观察
; t, R) f) Y% s4 q[lon1,lat1]=meshgrid([-180:2.5:179.75],[-90:2.5:90]);
8 J9 t$ g$ X" S9 o: a) E/ O" s- hu2=u(1:10:end,1:10:end)';
* g+ e# z6 T! B+ ?- l) ~% u3=[u2,u2(:,end)];
$ O2 C/ r8 W6 M% _* n, v7 _. J4 m" Av2=v(1:gap:end,1:gap:end)';
4 k$ Z" ?6 ]$ k3 F' A+ q% v3=[v2,v2(:,end)];
m_quiver(lon1,lat1,u2,v2,0);
8 \/ f$ L6 l4 i  @# O6 Bhold off
+ p  G# K! N5 Q9 A/ ^* Km_coast('patch',[1 .85 .7]);
2 u' K4 c( L7 ]m_grid('box','on','tickdir','out');
& @) c* q% H% S

1 d7 D" \( X$ _3 K) w! X8 Y. O' Y  K$ a
# ?# \" r$ b1 `7 @2 f2 F

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

: {8 f4 `6 M  q* {& d  i
. F8 X- I+ z6 u5 Y" I. ?! |; W
% a# ?9 @& R! ]( L* i3 U( o" W+ U9 g

: f9 M3 K! ^  s0 H' vm_quiver(lon1,lat1,u2,v2,0);
# Y! n. h# {6 E: I1 ~4 r  P
7 y: B$ R, p+ n: ~5 w0 a


9 ^  l) y9 M5 r' S8 B
: o8 B! i$ d/ m9 q. S! U  om_contourf(lon,lat,b,500,'linestyle','none')   %在mmap基础上的画
' ~' i- ]1 S, v( M: Rshading interp;  %使数据插值

5 G( Z. M. `4 \4 _* @. ]& m6 ?/ O


2 l; ^9 @4 F. ?0 S# L' f6 _

SsaaM7

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