找回密码
 注册
关于网站域名变更的通知
查看: 2771|回复: 1
打印 上一主题 下一主题

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

[复制链接]

该用户从未签到

跳转到指定楼层
1#
发表于 2021-3-3 09:50 | 只看该作者 |只看大图 回帖奖励 |倒序浏览 |阅读模式

EDA365欢迎您登录!

您需要 登录 才可以下载或查看,没有帐号?注册

x
: ]" z! ]0 J" d. v. u: _
一、matlab读取NCEP再分析数据并绘制风场0 O) j( p+ y; t% Q" u1 o/ N

# 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');1 N' v( {  |5 G
lat=ncread(f_hgt,'latitude');
3 K0 w4 h6 H* v) S. O& L%%%%%%时间转换0 j5 K9 k3 x% h, q
time  = double(time);; z1 x5 P. d$ b2 f+ z
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);%将时间字符数组转化为数值数组! g, a$ m5 e1 L2 z, ?* ~
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指定的开始位置算起,一共读取的每一维要素的数目" p( Q# h1 U2 [5 v
strip=[1,1,1,1];%从start开始,每一维读取的数目为count时,每一维的读取的步长5 H$ G1 }1 {7 I$ A8 J- V4 b
hgt=ncread(f_hgt,'z',start,count,strip);%读取温度值,单位K
4 [: s3 f9 D) S. U' lu=ncread(f_hgt,'u',start,count,strip);%读取温度值,单位K: P5 x5 z. c# h. _
v=ncread(f_hgt,'v',start,count,strip);%读取温度值,单位K/ }* @+ R2 M) F  t0 S
[X,Y]=meshgrid(lon,lat);. r$ J9 ^9 O- E; R/ D
figure(1)- B" K, z  h) G
m_proj('Mercator','lat',[25,35],'lon',[100,115]);; {% y' x/ A) d1 r+ ~+ V/ N# ?
% m_grid('linestyle','none','tickdir','out','fontsize',12,'fontname','Times New Roman');7 Q% l% U* V( \8 o0 ?
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')  ?1 n  Q# Q7 E! Z) c+ P$ D
%m_coast('patch',[.9 .9 .9],'edgecolor','none');" h) P  e/ z* e/ y5 w
[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]);%绘制范围内的地图2 \! @; r5 `* u9 G( G1 r
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')( d/ g  R: z; A
' b8 ?+ C8 d* m4 D: f5 Q+ U

9 u1 X! Q( y6 I8 A" s+ N- a0 @
! F$ H0 A8 g! |" H2 d
二、matlab读取ERA5并绘制全球风场图3 {# I, [% q! Y5 C$ ~3 I

% 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');! D& j1 K! v; X) j& W
time0=ncread('202008muwind.nc','time');' \5 f$ J! l' i# V% ]2 `. v
lat0=ncread('202008muwind.nc','lat');& E8 |+ q5 y2 D7 J5 r. p
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);$ i' k: v6 }# \6 c, `0 v7 Z) f
b=sqrt(u.^2+v.^2);/ ]. \5 s2 `2 ~* |- I" o) J9 c
figure(1);: u: x7 V# m* Q. y6 s
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)];, f) P0 w( V; @0 l1 q" @
v2=v(1:gap:end,1:gap:end)';. E7 M' E- S+ j4 h( T
% v3=[v2,v2(:,end)];
' n* g9 f* |, z: C3 V6 l( Wm_quiver(lon1,lat1,u2,v2,0);; U6 S0 x2 O7 b/ B4 Q% y- `7 A
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 q5 c5 g$ _2 V& z5 u& _) s
$ L8 i" w. r+ A

. 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基础上的画& T4 [' r: |/ G7 C8 ^7 W
shading interp;%使数据插值
8 Q+ l6 u. u6 Q( w' l: o* thold on;& Y0 m* Q6 w- ^/ E. I& ^! w9 r# K
m_quiver(lon1,lat1,u2,v2,0);- Z6 @5 k: `4 Y# \/ j) @/ J/ K. C
hold off;% b1 M6 h( V. v2 D

3 O( f( J, h( g0 i% X 7 u: z5 B% H* _3 U  i" S

. @5 }5 V; w, W' p9 r) O5 X2 T9 u" W) D5 o

: k; @6 b9 c9 f) r9 mm_quiver(lon1,lat1,u2,v2,0);- x" w- ]2 g; I5 `0 J2 ~

6 m8 V1 x' }8 j: T6 c
, M4 L- y4 M$ Z  P+ Y. i6 h9 c+ a" l7 V; k5 y/ p+ }+ |

# b5 P3 j. `5 Z- ~* X/ w4 z' H. O" B( q8 r$ T8 T% [& y! [5 i. y
m_contourf(lon,lat,b,500,'linestyle','none')   %在mmap基础上的画
# T$ Y8 B- |: C/ q9 g, Eshading interp;  %使数据插值
. M& u( }8 ~+ A2 f3 n( ^0 [+ [, c* ~9 j  f& _
1 k) J: d3 N- X2 a

2 v8 S# b% A: ?$ m2 A( o5 n0 u2 y/ n' u# x

该用户从未签到

2#
发表于 2021-3-3 10:52 | 只看该作者
Matlab全球风场图—以ERA5、NCEP和ECMWF再分析数据为例
您需要登录后才可以回帖 登录 | 注册

本版积分规则

关闭

推荐内容上一条 /1 下一条

EDA365公众号

关于我们|手机版|EDA365电子论坛网 ( 粤ICP备18020198号-1 )

GMT+8, 2025-10-30 15:29 , Processed in 0.171875 second(s), 26 queries , Gzip On.

深圳市墨知创新科技有限公司

地址:深圳市南山区科技生态园2栋A座805 电话:19926409050

快速回复 返回顶部 返回列表