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

thinkfunny

& }# }5 K! _! c- B( c7 v% I8 U一、matlab读取NCEP再分析数据并绘制风场
  p) @4 S7 k1 |* L4 C: ]
%该程序用于求水汽通量散度
. W5 U) o  @) |- o%注意！我们这里读到的u是四维矩阵，分别是lon*lat*level*time，
# ^9 H2 \7 V3 W6 Fclc;clear;close all
! B+ I6 n' \& E* V; Wf_hgt = 'ps_level_20170121_0130.nc';
% ncdisp(f_hgt);
time=ncread(f_hgt,'time');
: j2 k& @) `0 i1 U6 N2 N" Slevel=ncread(f_hgt,'level');
lon=ncread(f_hgt,'longitude');
lat=ncread(f_hgt,'latitude');
%%%%%%时间转换
time  = double(time);
( W5 I# D6 I* G1 g- X8 \& ~# d: hformat = 'mm dd, yyyy HH:MM:SS.FFF AM';%转换格式
dstr = datestr((datenum('1900-01-01') + time./24),format);%转换后时间字符串存储
TM = datevec(dstr);%将时间字符数组转化为数值数组
2 O# @7 a2 e; B; `$ jtidx=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];%所指定变量的每一维的开始读取的位置
2 M! N: d+ R/ ]7 }8 ^count=[41,31,1,1];%从start指定的开始位置算起，一共读取的每一维要素的数目
strip=[1,1,1,1];%从start开始，每一维读取的数目为count时，每一维的读取的步长
hgt=ncread(f_hgt,'z',start,count,strip);%读取温度值，单位K
% N5 [4 C6 `! m# tu=ncread(f_hgt,'u',start,count,strip);%读取温度值，单位K
v=ncread(f_hgt,'v',start,count,strip);%读取温度值，单位K
6 Y% v. A2 V6 L. C, j  h[X,Y]=meshgrid(lon,lat);
/ H2 F% C7 |2 I  z3 p( K2 s) j' Zfigure(1)
9 a" g, d. ^3 d- h0 Xm_proj('Mercator','lat',[25,35],'lon',[100,115]);
% m_grid('linestyle','none','tickdir','out','fontsize',12,'fontname','Times New Roman');
6 O3 L4 c, A! t$ X4 Q. e9 om_grid('linestyle','none','box','fancy','fontsize',11,'tickdir','in','xtick',[100:3:115],'ytick',[25:2:35]);
5 d6 v9 e) Z7 k5 D$ Phold on
m_windbarb(X',Y',u,v,'color','k')
6 l7 G! Y/ J/ X: B/ p- X5 s%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]
+ @6 k2 }6 n5 U$ ^3 u0 M%set(b,'ShowText','on','TextStep',get(b,'LevelStep'),'fontsize',12,'fontname','Times New Roman');   %在等高线上叠加数值（文后详情）
: K7 j: H- u# b  xclabel(C,h,'FontSize',13,'fontname','Times New Roman');
$ f# N- o$ S7 M1 W, ?# E6 ]7 ima=shaperead('F:/RMeteoInfo/data/map/bou2_4l.shp');
- v7 e; q* b8 Z& o& B% 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')
5 d1 E3 @% L5 o2 W3 {
3 D" m8 w( z; b0 z" w3 G' M
# u) a3 O1 [5 r9 v

二、matlab读取ERA5并绘制全球风场图
3 S7 w  |$ ^( e! I
clc;clear;close all
u0=ncread('202008muwind.nc','uas');    %读取其中一项
v0=ncread('202008mvwind.nc','vas');
& `  k2 {. q0 H3 K: N( `0 v, atime0=ncread('202008muwind.nc','time');
lat0=ncread('202008muwind.nc','lat');
5 V( q  w3 A: s- G4 q. Flon0=ncread('202008muwind.nc','lon');      
! L, C  C/ `, k( d) ]# {2 D, ?, ][lat,lon]=meshgrid(lat0,lon0);
' X, o3 D' z6 `+ Y: W/ \. w( Lu=u0(:,:,2);
# C2 ^+ M9 d+ N( o6 [. X- dv=v0(:,:,2);
b=sqrt(u.^2+v.^2);
figure(1);
  ^/ d$ ?4 D" {1 i: x3 s4 Um_proj('Equidistant Cylindrical','long',[-180 180],'lat',[-90 90]);%矩形投影；取区域观察
[lon1,lat1]=meshgrid([-180:2.5:179.75],[-90:2.5:90]);
! D) c% c  v) ]" V; Wu2=u(1:10:end,1:10:end)';
% u3=[u2,u2(:,end)];
( m7 ^" O+ j' V0 ?' z; G. }v2=v(1:gap:end,1:gap:end)';
' m6 Y! x3 s' g2 {( W: i& E, H! o% v3=[v2,v2(:,end)];
m_quiver(lon1,lat1,u2,v2,0);
' `; E6 @7 }8 s* W: s$ @hold off
: m6 [9 B+ e- y9 G$ ]& e7 r  nm_coast('patch',[1 .85 .7]);
' v% V$ G5 C' Q1 X  f0 Dm_grid('box','on','tickdir','out');
" Q1 h5 S* Q8 @% H


  ]( [& n5 @5 z0 {- a6 f2 \
5 c9 l2 r4 c0 Q: x' V
m_contourf(lon,lat,b);%在mmap基础上的画
2 z* Z9 ^: D0 }, F9 Q: Bshading interp;%使数据插值
hold on;
2 ]9 N1 I# v7 l: f. rm_quiver(lon1,lat1,u2,v2,0);
hold off;


: _, }" B; i& T. ~( @+ x0 M3 G( j
9 Z# V7 N! H& a& j
" V9 C. \$ J( A+ `- {$ g4 u
m_quiver(lon1,lat1,u2,v2,0);
( O! o1 \" q% z7 a% k+ T+ T& ]# b$ ]# ^
, z3 {" R1 O( c" |


9 T: O% m5 a* A- ?- L
m_contourf(lon,lat,b,500,'linestyle','none')   %在mmap基础上的画
shading interp;  %使数据插值
' A5 Y4 A' v1 x; {
+ l, {, H- B) I' D% K4 Q
4 k5 @, W/ s$ u& U& C

SsaaM7

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