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

thinkfunny

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

# x+ e4 V- k! }: p6 b$ l%该程序用于求水汽通量散度
0 O1 |8 K/ ]9 r1 e" J- Z5 v7 B%注意！我们这里读到的u是四维矩阵，分别是lon*lat*level*time，
; R$ \. I4 D" |8 e8 ?! Jclc;clear;close all
" _8 L# b- c4 U4 [f_hgt = 'ps_level_20170121_0130.nc';
% K% r8 d0 v$ X: l$ n' e  K% ncdisp(f_hgt);
% O- [" U# l7 `' T1 ytime=ncread(f_hgt,'time');
. u% ^& N$ b. [5 U, ?6 ?' Q+ y  plevel=ncread(f_hgt,'level');
8 C; S& \8 T& U% R* ]. Blon=ncread(f_hgt,'longitude');
lat=ncread(f_hgt,'latitude');
3 K0 w4 h6 H* v) S. O& L%%%%%%时间转换
time  = double(time);
format = 'mm dd, yyyy HH:MM:SS.FFF AM';%转换格式
* q! G4 b! L' ?3 b9 b2 l4 e  bdstr = datestr((datenum('1900-01-01') + time./24),format);%转换后时间字符串存储
/ j; o, w# ?3 X! a+ STM = datevec(dstr);%将时间字符数组转化为数值数组
tidx=find(TM(:,2)==1 & TM(:,3)==28 & TM(:,4)==00);%筛选7月25日08时（世界时加8）
: J& G( t8 L" Q! M7 d; R+ `! \ps_lev=find(level ==850);%%删选出850hPa高度
& t5 G; A2 o: A4 s! W8 hstart=[1,1,ps_lev,tidx];%所指定变量的每一维的开始读取的位置
: @. \7 E) V" }4 V! Ccount=[41,31,1,1];%从start指定的开始位置算起，一共读取的每一维要素的数目
strip=[1,1,1,1];%从start开始，每一维读取的数目为count时，每一维的读取的步长
hgt=ncread(f_hgt,'z',start,count,strip);%读取温度值，单位K
4 [: s3 f9 D) S. U' lu=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]);
- z, [; w! w2 ohold on
$ C& j- s/ N5 E/ Q8 H. wm_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]
  j* z/ s6 Y/ ?% c% {4 N0 N2 \- N$ P%set(b,'ShowText','on','TextStep',get(b,'LevelStep'),'fontsize',12,'fontname','Times New Roman');   %在等高线上叠加数值（文后详情）
; D# n# o. D6 q  h# P# f3 r, S7 p/ Jclabel(C,h,'FontSize',13,'fontname','Times New Roman');
3 D3 G+ d0 K$ \" B: s2 n' X( Z& C) lma=shaperead('F:/RMeteoInfo/data/map/bou2_4l.shp');
& |0 Q; l6 E4 u* r% m_line( [ma(：).X], [ma(：).Y],'color',[0.5,0.5,0.5]);%绘制范围内的地图
m_line([ma(：).X],[ma(：).Y],'color','r');%绘制范围内的地图
* I7 |5 X( ~+ q5 Y" B% Em_plot(105.5,29.43,'marker','^','MarkerSize',7,'color','k','MarkeRFaceColor','k')


9 u1 X! Q( y6 I

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

% T* n4 t/ E9 Z: r2 Eclc;clear;close all
0 f' h8 T& {5 y# {3 m  [/ Mu0=ncread('202008muwind.nc','uas');    %读取其中一项
( g, \7 f0 V5 @8 j" u1 j; n' Hv0=ncread('202008mvwind.nc','vas');
time0=ncread('202008muwind.nc','time');
lat0=ncread('202008muwind.nc','lat');
lon0=ncread('202008muwind.nc','lon');      
8 c" n( B. j6 A[lat,lon]=meshgrid(lat0,lon0);
- H  m- S; H# D+ O- Y- @% Fu=u0(:,:,2);
4 H0 T. G* p# u2 E9 n/ mv=v0(:,:,2);
b=sqrt(u.^2+v.^2);
figure(1);
m_proj('Equidistant Cylindrical','long',[-180 180],'lat',[-90 90]);%矩形投影；取区域观察
" W5 \2 T# x5 {" `[lon1,lat1]=meshgrid([-180:2.5:179.75],[-90:2.5:90]);
0 J* i8 p& W$ H9 u6 f7 Ju2=u(1:10:end,1:10:end)';
% h, U: R4 p9 P9 U% u3=[u2,u2(:,end)];
v2=v(1:gap:end,1:gap:end)';
% v3=[v2,v2(:,end)];
' n* g9 f* |, z: C3 V6 l( Wm_quiver(lon1,lat1,u2,v2,0);
hold off
! y# g* `2 S- S8 a6 S! ?' v* Mm_coast('patch',[1 .85 .7]);
2 f% ]1 F" ?! H* c) P8 }  B1 s# Am_grid('box','on','tickdir','out');
9 y/ A1 Z" d) `7 ?3 |5 u1 G4 q


. W% ^! p+ y$ \, W% w9 D) z
6 x( \3 L; k% O  L2 q
" p9 {9 t6 Q# N3 Lm_contourf(lon,lat,b);%在mmap基础上的画
shading interp;%使数据插值
8 Q+ l6 u. u6 Q( w' l: o* thold on;
m_quiver(lon1,lat1,u2,v2,0);
hold off;

3 O( f( J, h( g0 i% X

. @5 }5 V; w, W' p

: k; @6 b9 c9 f) r9 mm_quiver(lon1,lat1,u2,v2,0);

6 m8 V1 x' }8 j: T6 c
, M4 L- y4 M$ Z  P+ Y. i

# b5 P3 j. `5 Z- ~* X/ w4 z
m_contourf(lon,lat,b,500,'linestyle','none')   %在mmap基础上的画
# T$ Y8 B- |: C/ q9 g, Eshading interp;  %使数据插值
. M& u( }8 ~+ A2 f


2 v8 S# b% A: ?$ m2 A

SsaaM7

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