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

thinkfunny

一、matlab读取NCEP再分析数据并绘制风场
# _( c  z1 J/ n5 j, T0 y
%该程序用于求水汽通量散度
%注意！我们这里读到的u是四维矩阵，分别是lon*lat*level*time，
clc;clear;close all
0 L4 {! v  q5 S# rf_hgt = 'ps_level_20170121_0130.nc';
% ncdisp(f_hgt);
" u' p! |- ~  l1 l5 \time=ncread(f_hgt,'time');
9 u! T2 D. W' F' B# m3 Dlevel=ncread(f_hgt,'level');
lon=ncread(f_hgt,'longitude');
( y, h( I! |- [; r8 ]lat=ncread(f_hgt,'latitude');
%%%%%%时间转换
time  = double(time);
% {; g  q; }8 I) q" j# [  D& f. a  \format = 'mm dd, yyyy HH:MM:SS.FFF AM';%转换格式
dstr = datestr((datenum('1900-01-01') + time./24),format);%转换后时间字符串存储
; @2 S, Q7 E5 Z5 |6 P, U0 A$ ?0 fTM = datevec(dstr);%将时间字符数组转化为数值数组
5 |/ E& i/ y4 a9 W) ]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];%所指定变量的每一维的开始读取的位置
/ j& h8 F2 u7 K& L) bcount=[41,31,1,1];%从start指定的开始位置算起，一共读取的每一维要素的数目
strip=[1,1,1,1];%从start开始，每一维读取的数目为count时，每一维的读取的步长
0 H, p6 z3 b( L+ e. g7 O0 ihgt=ncread(f_hgt,'z',start,count,strip);%读取温度值，单位K
( K7 e% r. F: P& |u=ncread(f_hgt,'u',start,count,strip);%读取温度值，单位K
v=ncread(f_hgt,'v',start,count,strip);%读取温度值，单位K
[X,Y]=meshgrid(lon,lat);
4 {8 U2 a! s. |/ `3 }figure(1)
m_proj('Mercator','lat',[25,35],'lon',[100,115]);
% m_grid('linestyle','none','tickdir','out','fontsize',12,'fontname','Times New Roman');
2 Z1 n# W1 s! I7 B+ D0 C2 T( z# H1 e* sm_grid('linestyle','none','box','fancy','fontsize',11,'tickdir','in','xtick',[100:3:115],'ytick',[25:2:35]);
, V' I# U2 U' y! e/ ehold on
0 P7 @& n! i  m. w( b7 K# }9 Wm_windbarb(X',Y',u,v,'color','k')
  f  C" z# L1 [%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');   %在等高线上叠加数值（文后详情）
6 E- @5 `0 e( P' Wclabel(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')
8 N8 F4 N# k/ j5 i, R


( b. V, j% y4 P* z3 {
+ ~) Z; H2 U& }" M二、matlab读取ERA5并绘制全球风场图
* v( P) b/ C# ^3 \# _7 l/ N( n
3 i7 M$ a: S' L' x2 aclc;clear;close all
u0=ncread('202008muwind.nc','uas');    %读取其中一项
/ O0 |& i0 b; c* E7 _- E" Fv0=ncread('202008mvwind.nc','vas');
+ X& O1 L" ~0 n$ ]  W- m- G% Etime0=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);
0 w3 J: s. _( km_proj('Equidistant Cylindrical','long',[-180 180],'lat',[-90 90]);%矩形投影；取区域观察
0 c; O/ _' e4 g9 P9 L' c6 t[lon1,lat1]=meshgrid([-180:2.5:179.75],[-90:2.5:90]);
u2=u(1:10:end,1:10:end)';
4 T" h7 c7 p" d) m; E: S& S% u3=[u2,u2(:,end)];
7 P9 u: _$ _8 [$ g0 e, Rv2=v(1:gap:end,1:gap:end)';
% v3=[v2,v2(:,end)];
8 w0 I) Y# ]* v5 g9 ?m_quiver(lon1,lat1,u2,v2,0);
hold off
; ^7 ?6 h% R5 q. {# @3 @m_coast('patch',[1 .85 .7]);
6 O5 m1 w% y. a+ K7 J! e$ qm_grid('box','on','tickdir','out');
/ T( _( [3 s& G- n
$ m! _0 g) ]- ?! w5 g+ x


- B( N! @4 o9 s) b5 ^0 m. f
m_contourf(lon,lat,b);%在mmap基础上的画
0 \; T2 P  h5 T3 Oshading interp;%使数据插值
hold on;
m_quiver(lon1,lat1,u2,v2,0);
hold off;
0 l* K4 D$ O& v: u
1 P7 l" c8 x- o

% g6 b2 v! h5 G1 p2 p! @
- M3 F, Q/ p  i8 U. J+ R
; ~* d' I5 z3 `0 N) Sm_quiver(lon1,lat1,u2,v2,0);
( n5 J$ S4 ]' T) g4 D



0 V. M# e" `, ~$ v
m_contourf(lon,lat,b,500,'linestyle','none')   %在mmap基础上的画
shading interp;  %使数据插值

! a" _. E9 |% A

! `% k6 l" h  Y
7 _% F  a# @+ I3 n5 H* q; I

SsaaM7

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