想问一下王老师，wrfout数据里的500hpa风场要怎么画

James_z

之前用python画wrfout，有个wrf-python库可以给u,v插值到500hpa，现在想用meteoinfo画画wrfout图，但是不知道怎么插值画500hpa的风场，有人画过可以教教我吗

人品守恒

MeteoInfo 发表于 2021-6-20 23:56
wrfout里有些变量输出在交错格点上（stag），为了计算中保持变量维一致，需要去stag。destagger函数第二 ...

王老师您好，我最近使用WRF时也刚好遇到这种问题，如果要绘制标量，例如湿度、位势高度该怎么修改呢

MeteoInfo

参考下面的代码：

1. fn = 'D:/Temp/nc/wrfout_d01_2005-08-28_00-00-00'
   
2. f = addfile(fn)
   
3. 
   
4. #Calculate pressure array
   
5. pp = f['P'][:]
   
6. psb = f['PB'][:]
   
7. ptop = f['P_TOP'][0]
   
8. eta = f['ZNU'][0]
   
9. nt,nz,ny,nx = pp.shape
   
10. eta = eta.reshape(1,nz,1,1)
    
11. pres = (psb * eta + ptop + pp) * 0.01
    
12. 
    
13. #Read and de-stagger u and v
    
14. u = f['U'][:]
    
15. v = f['V'][:]
    
16. destag_u = meteo.wrf.destagger(u, -1)
    
17. destag_v = meteo.wrf.destagger(v, -2)
    
18. y = destag_u.dimvalue(-2)
    
19. x = destag_u.dimvalue(-1)
    
20. 
    
21. #Interpolate u and v to pressure level
    
22. levels = [700, 500]
    
23. u1 = meteo.log_interpolate_1d(levels, pres, destag_u, axis=1)
    
24. v1 = meteo.log_interpolate_1d(levels, pres, destag_v, axis=1)
    
25. 
    
26. #Plot
    
27. axesm(projinfo=f.proj)
    
28. geoshow('country')
    
29. quiver(x[::4], y[::4], u1[0,1,::4,::4], v1[0,1,::4,::4], color='b', proj=f.proj)
    
30. title('Pressure 500 hPa')

fn = 'D:/Temp/nc/wrfout_d01_2005-08-28_00-00-00'<br /> f = addfile(fn)<br /> <br /> #Calculate pressure array<br /> pp = f['P'][:]<br /> psb = f['PB'][:]<br /> ptop = f['P_TOP'][0]<br /> eta = f['ZNU'][0]<br /> nt,nz,ny,nx = pp.shape<br /> eta = eta.reshape(1,nz,1,1)<br /> pres = (psb * eta + ptop + pp) * 0.01<br /> <br /> #Read and de-stagger u and v<br /> u = f['U'][:]<br /> v = f['V'][:]<br /> destag_u = meteo.wrf.destagger(u, -1)<br /> destag_v = meteo.wrf.destagger(v, -2)<br /> y = destag_u.dimvalue(-2)<br /> x = destag_u.dimvalue(-1)<br /> <br /> #Interpolate u and v to pressure level<br /> levels = [700, 500]<br /> u1 = meteo.log_interpolate_1d(levels, pres, destag_u, axis=1)<br /> v1 = meteo.log_interpolate_1d(levels, pres, destag_v, axis=1)<br /> <br /> #Plot<br /> axesm(projinfo=f.proj)<br /> geoshow('country')<br /> quiver(x[::4], y[::4], u1[0,1,::4,::4], v1[0,1,::4,::4], color='b', proj=f.proj)<br /> title('Pressure 500 hPa')

James_z

MeteoInfo 发表于 2021-6-20 21:49
参考下面的代码：

王老师您好好，我这里还有些问题不懂，我参考了你的代码，uv那个插值完后，quiver绘制流场图那里u1[0,1,::4,::4]，v1[0,1,::4,::4],里面的1是什么意思，1是指500hpa吗，uv插值完后还有0这层代表什么

MeteoInfo

James_z 发表于 2021-6-20 22:14
王老师您好好，我这里还有些问题不懂，我参考了你的代码，uv那个插值完后，quiver绘制流场图那里u1[0,1,: ...

levels = [700, 500]
这个脚本插值到了两个等压面。0和1分布是700和500

James_z

MeteoInfo 发表于 2021-6-20 22:27
levels = [700, 500]
这个脚本插值到了两个等压面。0和1分布是700和500

好的弄懂了！谢谢王老师

James_z

James_z 发表于 2021-6-20 22:49
好的弄懂了！谢谢王老师

destag_u = meteo.wrf.destagger(u, -1)
destag_v = meteo.wrf.destagger(v, -2)
老师我还有一点点问题，我在你那本书上找不到关于destagger的说明，那个(u, -1)和(v, -2)后面的-1和-2是啥意思，是跟矢量有关吗，如果我想绘制500hpa温度这种标量变量该怎么用这个函数

MeteoInfo

James_z 发表于 2021-6-20 23:34
destag_u = meteo.wrf.destagger(u, -1)
destag_v = meteo.wrf.destagger(v, -2)
老师我还有一点点问题 ...

wrfout里有些变量输出在交错格点上（stag），为了计算中保持变量维一致，需要去stag。destagger函数第二个参数是stag所在的维，-1是倒数第一个维（也就是最后一个维）。

哈库拉玛塔塔000

MeteoInfo 发表于 2021-6-20 21:49
参考下面的代码：

老师，我使用了这个代码，但是得出的风场维度不对

MeteoInfo

哈库拉玛塔塔000 发表于 2022-1-26 11:21
老师，我使用了这个代码，但是报了这样的错误。

示例程序里 destag_u 变量是4维（time, level, y, x），log_interpolate_1d 函数中 axis=1 就制定了要对 level 维插值。如果你的 destag_u 变量是3维（没有时间维），那么 axis=0 。原则是 axis 参数要指定高度维的序号。