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

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

[复制链接]

该用户从未签到

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

EDA365欢迎您登录!

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

x
- _4 @( i& o( w7 o/ N
一、matlab读取NCEP再分析数据并绘制风场
% A- N. o$ r, f- u. p8 [
; r8 p( R$ S3 C% _( U%该程序用于求水汽通量散度
  k8 D' ?/ n! T%注意!我们这里读到的u是四维矩阵,分别是lon*lat*level*time,, K/ ]3 V7 Q% }& F; e. H9 {
clc;clear;close all9 J8 F$ E+ c3 S
f_hgt = 'ps_level_20170121_0130.nc';
. y$ e# x1 z2 m! `2 O/ t% ncdisp(f_hgt);
6 o4 N9 f+ R2 u0 o5 b  ntime=ncread(f_hgt,'time');8 o" |" j. D) H% u3 o2 \
level=ncread(f_hgt,'level');4 u% D$ r: R# O! ?) \3 M# Q
lon=ncread(f_hgt,'longitude');
3 d1 `: i; O! Q: Blat=ncread(f_hgt,'latitude');5 i* z4 z( \! G/ K1 g2 }3 s6 V9 I
%%%%%%时间转换0 I2 y3 |& I* f; S
time  = double(time);' m8 \! T+ k) f4 q: ]( S( m
format = 'mm dd, yyyy HH:MM:SS.FFF AM';%转换格式6 \% ~7 |: z$ @. `. U  J2 ]+ O, A
dstr = datestr((datenum('1900-01-01') + time./24),format);%转换后时间字符串存储
. \+ V; g# O0 zTM = datevec(dstr);%将时间字符数组转化为数值数组
8 H, @7 Y" P, w7 ntidx=find(TM(:,2)==1 & TM(:,3)==28 & TM(:,4)==00);%筛选7月25日08时(世界时加8)
2 l2 g3 y: I# bps_lev=find(level ==850);%%删选出850hPa高度
4 n$ D3 D, m5 o# M/ j- [" Fstart=[1,1,ps_lev,tidx];%所指定变量的每一维的开始读取的位置7 @9 p% Y8 Q- U9 F: U% K
count=[41,31,1,1];%从start指定的开始位置算起,一共读取的每一维要素的数目
9 z! {0 }, ^6 C) ?' x0 U$ \& D5 y1 wstrip=[1,1,1,1];%从start开始,每一维读取的数目为count时,每一维的读取的步长$ @  c/ i/ o3 Z5 I/ V
hgt=ncread(f_hgt,'z',start,count,strip);%读取温度值,单位K. [  s, q, C- G; |1 z" W2 Q/ f
u=ncread(f_hgt,'u',start,count,strip);%读取温度值,单位K# ]) {  x# M1 n
v=ncread(f_hgt,'v',start,count,strip);%读取温度值,单位K( G! J% K+ S  B2 a
[X,Y]=meshgrid(lon,lat);
/ \  m9 Z9 Z1 ]  F/ L" }) X( Kfigure(1)
4 `: Y, U/ n* U: Q3 tm_proj('Mercator','lat',[25,35],'lon',[100,115]);
- R" _% M  D" t" e1 P1 q: |2 K0 N% m_grid('linestyle','none','tickdir','out','fontsize',12,'fontname','Times New Roman');
0 ^2 E) P5 P, u3 h" s6 }$ em_grid('linestyle','none','box','fancy','fontsize',11,'tickdir','in','xtick',[100:3:115],'ytick',[25:2:35]);
" w8 X3 {7 z2 D9 {; W1 bhold on
  w: N8 T4 H, H9 Dm_windbarb(X',Y',u,v,'color','k')
7 ^1 C/ p. \8 W0 ]%m_coast('patch',[.9 .9 .9],'edgecolor','none');
" }* r; i/ q0 I- r) y3 l5 d1 E0 S[C,h]=m_contour(X',Y',hgt/10 ,'color','k','LineWidth',1);%[5000:80:5900]%[1520:25:1680]! H5 k$ k# y' G- u5 w6 }
%set(b,'ShowText','on','TextStep',get(b,'LevelStep'),'fontsize',12,'fontname','Times New Roman');   %在等高线上叠加数值(文后详情)* s/ b8 x: B# O$ k9 ~
clabel(C,h,'FontSize',13,'fontname','Times New Roman');
8 h) j7 v5 R: k9 v2 r# l9 U) Sma=shaperead('F:/RMeteoInfo/data/map/bou2_4l.shp'); * P8 w& N& g& [  q5 W1 h; N# j
% m_line( [ma(:).X], [ma(:).Y],'color',[0.5,0.5,0.5]);%绘制范围内的地图/ c  w9 P  `  e; U
m_line([ma(:).X],[ma(:).Y],'color','r');%绘制范围内的地图
* q3 Q& Z7 _0 M8 \# X, zm_plot(105.5,29.43,'marker','^','MarkerSize',7,'color','k','MarkeRFaceColor','k')
( m" j/ X! U, h8 B3 w7 _* |- f3 R6 Y9 _$ M& r

( r  v5 v- O/ w, Q9 [# D
. C( F! x& Z( i5 X1 Z; D8 ?) Q% E& u
二、matlab读取ERA5并绘制全球风场图
# C5 e: ?5 x/ V1 ]" c% f& H% e' h3 W' l) L% w; t* x
clc;clear;close all
4 v$ J' R& s* l. m( }6 |! lu0=ncread('202008muwind.nc','uas');    %读取其中一项
1 G1 A) m( M- C3 \, K, Ov0=ncread('202008mvwind.nc','vas');
; O7 r7 c+ @) D0 F; }4 D* s, Z& |time0=ncread('202008muwind.nc','time');( p  d6 z9 x: @8 C/ }
lat0=ncread('202008muwind.nc','lat');4 Z- y1 m( P6 p1 `- H9 G' n! t) `2 Q+ P
lon0=ncread('202008muwind.nc','lon');      
9 @) m0 e* R  v# G8 p[lat,lon]=meshgrid(lat0,lon0);
: J0 L: W# H4 V" y4 ]u=u0(:,:,2);* }, G7 E! B( Q/ {$ j
v=v0(:,:,2);3 _& R' R$ _) Y. D) a/ k& Y
b=sqrt(u.^2+v.^2);
: }& J, G0 B0 @+ r/ _; d+ {& Bfigure(1);5 I- Q* K8 ]8 b$ R
m_proj('Equidistant Cylindrical','long',[-180 180],'lat',[-90 90]);%矩形投影;取区域观察. Q' e" C9 y9 m6 l# d
[lon1,lat1]=meshgrid([-180:2.5:179.75],[-90:2.5:90]);" i% r9 I( [' q* y3 t
u2=u(1:10:end,1:10:end)';; k' B+ g; ?, N, L( b) G. {8 \1 Q
% u3=[u2,u2(:,end)];. D5 q- x" D8 z1 g8 N0 I
v2=v(1:gap:end,1:gap:end)';
3 F! R5 w+ a0 b; ?! ?% v3=[v2,v2(:,end)];- ]( L, o) j) V) I1 {. K6 s* u
m_quiver(lon1,lat1,u2,v2,0);
0 z+ {% Q0 y+ D& R5 o% Xhold off
3 ]* u$ q% c/ ~" @- l, m) {m_coast('patch',[1 .85 .7]);
( M: N1 S3 ^" W) {% ^7 Jm_grid('box','on','tickdir','out');! e: Z& h' R: f+ U
. i& ^- N' I( F+ H- o0 W( z

3 M. x4 [2 [. P. x  U# R" U  q' _" r6 l4 [# P3 _1 P9 l. B, s) P
, X9 V. |9 s# f9 x& I$ N

: M1 E3 N% l! e4 k. w" {$ x$ Y9 Qm_contourf(lon,lat,b);%在mmap基础上的画
4 _. g( R1 D2 P3 m8 Zshading interp;%使数据插值" a, \5 @1 V7 o# `% ]
hold on;, F* ^9 a; J! U
m_quiver(lon1,lat1,u2,v2,0);% S6 R/ E" y# `) y/ T5 T2 A3 ]5 D
hold off;6 m# _7 r& w, Z, l3 J) k
  \  V& E9 q5 B( m) b! w# }3 ^# m: y
/ e9 y( ]) p# o/ d) h6 E" `

6 l5 l$ y- n; f. Q% L* L: R
+ ^) p2 s  _* s9 E9 J+ q
) m, c8 |# s/ B: S, Z' q4 Dm_quiver(lon1,lat1,u2,v2,0);0 t  d& j1 `* Y

: g3 z- @) I1 r8 k1 _' m . P" w8 J9 V5 J9 n3 Y3 m5 [

0 @0 J: F5 D  C( `' u1 R3 k. m" \1 c- P4 K7 u

2 J; v6 q$ h& A+ N( om_contourf(lon,lat,b,500,'linestyle','none')   %在mmap基础上的画- w0 ~5 {# I! [8 f4 q7 N# S
shading interp;  %使数据插值; J3 _5 M& m: F" j  o: q) C3 t

- M/ p6 A0 l: x5 p: M8 _& z' Z8 j! o; \& u
0 M) {2 [# [- i

  {# S0 R, K! k" T! N7 E

该用户从未签到

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

本版积分规则

关闭

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

EDA365公众号

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

GMT+8, 2025-10-6 05:30 , Processed in 0.125000 second(s), 26 queries , Gzip On.

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

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

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