沙尘模式地球坐标系三维动图

MeteoInfo

今年沙尘暴目前发生比较频繁，很多沙尘暴模式能够比较好地模拟预报沙尘浓度的时空变化，这里利用中国气象局的沙尘暴模式 CUACE/Dust 预报结果绘制一个地球坐标系中的三维动图。在MeteoInfoLab中，axes3d(projection='earth') 创建一个三维地球坐标系，缺省绘制一个偏暗色的地形纹理图像，可以加上 image 参数指定地球表面的纹理图像（在MeteoInfo的map文件夹中有几个jpg文件可用）。模式预报结果是一个netCDF文件，读取沙尘浓度和风场U/V分量的三维数组。沙尘浓度可以用等值面（isosurface）来表示，这里用了两个浓度：100微克每立方米和500微克每立方米，高浓度用更深的颜色和更小的透明度可以展示出浓度的层次。风场用streamslice函数绘制850 hPa风场的流线图，geoshow函数绘制地图行政界线。

单个时次的示例脚本程序如下：

1. #Set date
   
2. sdate = datetime.datetime(2023, 4, 9, 0)
   
3. 
   
4. #Set directory
   
5. datadir = r'V:Asian-dustmodelCMA'
   
6. datadir = os.path.join(datadir, sdate.strftime('%Y%m%d'))
   
7. #datadir = r'D:Tempcuace‚21212'
   
8. 
   
9. #Read data
   
10. fn = os.path.join(datadir, 'CUACE-DUST_CMA_{}.nc'.format(sdate.strftime('%Y%m%d%H')))
    
11. f = addfile(fn)
    
12. st = f.gettime(0)
    
13. t = 3
    
14. dust = f['CONC_DUST'][t]
    
15. u = f['U'][t]
    
16. v = f['V'][t]
    
17. speed = sqrt(u*u + v*v)
    
18. levels = dust.dimvalue(0)
    
19. #dust[dust<5] = 0
    
20. height = meteolib.pressure_to_height_std(levels)
    
21. height = height / 10
    
22. lat = dust.dimvalue(1)
    
23. lon = dust.dimvalue(2)
    
24. 
    
25. #Plot
    
26. ax = axes3d(projection='earth', image='shadedrelief.jpg')
    
27. lighting(True)
    
28. ax.set_rotation(160)
    
29. ax.set_elevation(-50)
    
30. ax.set_head(0)
    
31. ax.set_pitch(-61.5)
    
32. 
    
33. geoshow('cn_province', edgecolor='gray', lighting=False)
    
34. geoshow('country', edgecolor='gray', lighting=False)
    
35. #Beijing location
    
36. plot3([116.39,116.39], [39.91,39.91], [0,1200], color='w', lighting=False)
    
37. #ax.add_zaxis(116.39, 39.91)
    
38. #streamslice
    
39. ss = slice(None, None, 4)
    
40. levs = arange(2, 20, 2)
    
41. hidx = 4    # 850hPa
    
42. #hidx = 11   # 500hPa
    
43. gss = streamslice(lon[ss], lat[ss], height, u[:,ss,ss], v[:,ss,ss], u[:,ss,ss],
    
44.     speed[:,ss,ss], levs=levs, zslice=[height[hidx]], interval=10)
    
45. colorbar(gss[0], aspect=30, tickcolor='w', xshift=150, label='m/s')
    
46. #isosurface
    
47. isosurface(lon, lat, height, dust, 500, facecolor='r', edgecolor=None,
    
48.     alpha=0.6, nthread=4)
    
49. isosurface(lon, lat, height, dust, 100, facecolor=[255,180,0],
    
50.     edgecolor=None, alpha=0.4, nthread=4)
    
51. 
    
52. v = 1200
    
53. axis([-v,v,-v,v,-v,v])
    
54. 
    
55. tt = st + datetime.timedelta(hours=t*3)
    
56. title('Dust concentration higher than 100/500 ug/m3 ({}UTC)'.format(tt.strftime('%Y-%m-%d %H:00')),
    
57.     color='w')
    
58. #savefig('D:/Temp/test/dust_3d_{}.png'.format(tt.strftime('%Y%m%d%H')), dpi=150)

#Set date<br /> sdate = datetime.datetime(2023, 4, 9, 0)<br /> <br /> #Set directory<br /> datadir = r'V:Asian-dustmodelCMA'<br /> datadir = os.path.join(datadir, sdate.strftime('%Y%m%d'))<br /> #datadir = r'D:Tempcuace&#8218;21212'<br /> <br /> #Read data<br /> fn = os.path.join(datadir, 'CUACE-DUST_CMA_{}.nc'.format(sdate.strftime('%Y%m%d%H')))<br /> f = addfile(fn)<br /> st = f.gettime(0)<br /> t = 3<br /> dust = f['CONC_DUST'][t]<br /> u = f['U'][t]<br /> v = f['V'][t]<br /> speed = sqrt(u*u + v*v)<br /> levels = dust.dimvalue(0)<br /> #dust[dust<5] = 0<br /> height = meteolib.pressure_to_height_std(levels)<br /> height = height / 10<br /> lat = dust.dimvalue(1)<br /> lon = dust.dimvalue(2)<br /> <br /> #Plot<br /> ax = axes3d(projection='earth', image='shadedrelief.jpg')<br /> lighting(True)<br /> ax.set_rotation(160)<br /> ax.set_elevation(-50)<br /> ax.set_head(0)<br /> ax.set_pitch(-61.5)<br /> <br /> geoshow('cn_province', edgecolor='gray', lighting=False)<br /> geoshow('country', edgecolor='gray', lighting=False)<br /> #Beijing location<br /> plot3([116.39,116.39], [39.91,39.91], [0,1200], color='w', lighting=False)<br /> #ax.add_zaxis(116.39, 39.91)<br /> #streamslice<br /> ss = slice(None, None, 4)<br /> levs = arange(2, 20, 2)<br /> hidx = 4    # 850hPa<br /> #hidx = 11   # 500hPa<br /> gss = streamslice(lon[ss], lat[ss], height, u[:,ss,ss], v[:,ss,ss], u[:,ss,ss],<br />     speed[:,ss,ss], levs=levs, zslice=[height[hidx]], interval=10)<br /> colorbar(gss[0], aspect=30, tickcolor='w', xshift=150, label='m/s')<br /> #isosurface<br /> isosurface(lon, lat, height, dust, 500, facecolor='r', edgecolor=None,<br />     alpha=0.6, nthread=4)<br /> isosurface(lon, lat, height, dust, 100, facecolor=[255,180,0],<br />     edgecolor=None, alpha=0.4, nthread=4)<br /> <br /> v = 1200<br /> axis([-v,v,-v,v,-v,v])<br /> <br /> tt = st + datetime.timedelta(hours=t*3)<br /> title('Dust concentration higher than 100/500 ug/m3 ({}UTC)'.format(tt.strftime('%Y-%m-%d %H:00')),<br />     color='w')<br /> #savefig('D:/Temp/test/dust_3d_{}.png'.format(tt.strftime('%Y%m%d%H')), dpi=150)

读取多个时次的数据可以绘制时间动画 gif 图，主要用到 gifanimation, gifaddframe 等函数。

1. #Set date
   
2. sdate = datetime.datetime(2023, 4, 9, 0)
   
3. 
   
4. #Set directory
   
5. datadir = r'V:Asian-dustmodelCMA'
   
6. datadir = os.path.join(datadir, sdate.strftime('%Y%m%d'))
   
7. outdir = 'V:/Asian-dust/figures'
   
8. outdir = os.path.join(outdir, 'CMA', sdate.strftime('%Y%m%d'))
   
9. if not os.path.exists(outdir):
   
10.     os.makedirs(outdir)
    
11. 
    
12. #Read data
    
13. fn = os.path.join(datadir, 'CUACE-DUST_CMA_{}.nc'.format(sdate.strftime('%Y%m%d%H')))
    
14. f = addfile(fn)
    
15. st = f.gettime(0)
    
16. t = 10
    
17. dust = f['CONC_DUST'][t]
    
18. levels = dust.dimvalue(0)
    
19. #dust[dust<5] = 0
    
20. height = meteolib.pressure_to_height_std(levels)
    
21. height = height / 10
    
22. lat = dust.dimvalue(1)
    
23. lon = dust.dimvalue(2)
    
24. 
    
25. #Plot
    
26. ax = axes3d(projection='earth', image='shadedrelief.jpg')
    
27. lighting(True)
    
28. ax.set_rotation(160)
    
29. ax.set_elevation(-50)
    
30. ax.set_head(0)
    
31. ax.set_pitch(-61.5)
    
32. 
    
33. geoshow('cn_province', edgecolor='gray', lighting=False)
    
34. geoshow('country', edgecolor='gray', lighting=False)
    
35. #Beijing location
    
36. plot3([116.39,116.39], [39.91,39.91], [0,1200], color='w', lighting=False)
    
37. #streamslice
    
38. ss = slice(None, None, 4)
    
39. levs = arange(2, 20, 2)
    
40. 
    
41. #Loop
    
42. giffn = os.path.join(outdir, 'Dust_3d_isosurface_relief_{}--loop_1.gif'.format(st.strftime('%Y%m%d%H')))
    
43. print giffn
    
44. animation = gifanimation(giffn)
    
45. tn = f.timenum()
    
46. tn = 24
    
47. for t in range(1, tn):
    
48.     print(t)
    
49.     if t > 1:
    
50.         cll()
    
51.         cll()
    
52.         cll()
    
53. 
    
54.     dust = f['CONC_DUST'][t]
    
55.     u = f['U'][t]
    
56.     v = f['V'][t]
    
57.     speed = sqrt(u*u + v*v)
    
58. 
    
59.     gss = streamslice(lon[ss], lat[ss], height, u[:,ss,ss], v[:,ss,ss], u[:,ss,ss],
    
60.         speed[:,ss,ss], levs=levs, zslice=[height[4]], interval=10)
    
61.     colorbar(gss[0], aspect=30, tickcolor='w', xshift=135, label='m/s')
    
62.     #isosurface
    
63.     isosurface(lon, lat, height, dust, 500, facecolor='r', edgecolor=None,
    
64.         alpha=0.6, nthread=4)
    
65.     isosurface(lon, lat, height, dust, 100, facecolor=[255,180,0],
    
66.         edgecolor=None, alpha=0.4, nthread=4)
    
67. 
    
68.     vv = 1200
    
69.     axis([-vv,vv,-vv,vv,-vv,vv])
    
70. 
    
71.     tt = f.gettime(t)
    
72.     title('Dust concentration higher than 100/500 ug/m3 ({}UTC)'.format(tt.strftime('%Y-%m-%d %H:00')),
    
73.         color='w')
    
74.     #Add frame to gif animation
    
75.     gifaddframe(animation, width=1095, height=647)
    
76. 
    
77. #Finish gif animation
    
78. animation.finish()
    
79. print 'Finished...'

#Set date<br /> sdate = datetime.datetime(2023, 4, 9, 0)<br /> <br /> #Set directory<br /> datadir = r'V:Asian-dustmodelCMA'<br /> datadir = os.path.join(datadir, sdate.strftime('%Y%m%d'))<br /> outdir = 'V:/Asian-dust/figures'<br /> outdir = os.path.join(outdir, 'CMA', sdate.strftime('%Y%m%d'))<br /> if not os.path.exists(outdir):<br />     os.makedirs(outdir)<br /> <br /> #Read data<br /> fn = os.path.join(datadir, 'CUACE-DUST_CMA_{}.nc'.format(sdate.strftime('%Y%m%d%H')))<br /> f = addfile(fn)<br /> st = f.gettime(0)<br /> t = 10<br /> dust = f['CONC_DUST'][t]<br /> levels = dust.dimvalue(0)<br /> #dust[dust<5] = 0<br /> height = meteolib.pressure_to_height_std(levels)<br /> height = height / 10<br /> lat = dust.dimvalue(1)<br /> lon = dust.dimvalue(2)<br /> <br /> #Plot<br /> ax = axes3d(projection='earth', image='shadedrelief.jpg')<br /> lighting(True)<br /> ax.set_rotation(160)<br /> ax.set_elevation(-50)<br /> ax.set_head(0)<br /> ax.set_pitch(-61.5)<br /> <br /> geoshow('cn_province', edgecolor='gray', lighting=False)<br /> geoshow('country', edgecolor='gray', lighting=False)<br /> #Beijing location<br /> plot3([116.39,116.39], [39.91,39.91], [0,1200], color='w', lighting=False)<br /> #streamslice<br /> ss = slice(None, None, 4)<br /> levs = arange(2, 20, 2)<br /> <br /> #Loop<br /> giffn = os.path.join(outdir, 'Dust_3d_isosurface_relief_{}--loop_1.gif'.format(st.strftime('%Y%m%d%H')))<br /> print giffn<br /> animation = gifanimation(giffn)<br /> tn = f.timenum()<br /> tn = 24<br /> for t in range(1, tn):<br />     print(t)<br />     if t > 1:<br />         cll()<br />         cll()<br />         cll()<br /> <br />     dust = f['CONC_DUST'][t]<br />     u = f['U'][t]<br />     v = f['V'][t]<br />     speed = sqrt(u*u + v*v)<br /> <br />     gss = streamslice(lon[ss], lat[ss], height, u[:,ss,ss], v[:,ss,ss], u[:,ss,ss],<br />         speed[:,ss,ss], levs=levs, zslice=[height[4]], interval=10)<br />     colorbar(gss[0], aspect=30, tickcolor='w', xshift=135, label='m/s')<br />     #isosurface<br />     isosurface(lon, lat, height, dust, 500, facecolor='r', edgecolor=None,<br />         alpha=0.6, nthread=4)<br />     isosurface(lon, lat, height, dust, 100, facecolor=[255,180,0],<br />         edgecolor=None, alpha=0.4, nthread=4)<br /> <br />     vv = 1200<br />     axis([-vv,vv,-vv,vv,-vv,vv])<br /> <br />     tt = f.gettime(t)<br />     title('Dust concentration higher than 100/500 ug/m3 ({}UTC)'.format(tt.strftime('%Y-%m-%d %H:00')),<br />         color='w')<br />     #Add frame to gif animation<br />     gifaddframe(animation, width=1095, height=647)<br /> <br /> #Finish gif animation<br /> animation.finish()<br /> print 'Finished...'

动画图效果如下：

https://i0.hdslb.com/bfs/article/daf0448192b0394a9c9e766143c68f8bf9b6f036.gif@942w_558h_progressive.webp

网哇哇哇哇

数据在哪里可以下载呢请问

MeteoInfo

网哇哇哇哇 发表于 2023-4-26 20:11
数据在哪里可以下载呢请问

可以用这个文件做测试：WMO_SDS-WAS_Asian_Center_Model_Forecasting_CUACE-DUST_CMA_2019041500.nc （http://www.meteothink.org/downloads/data.html）

xrpamhf

MeteoInfo 发表于 2023-4-27 16:54
可以用这个文件做测试：WMO_SDS-WAS_Asian_Center_Model_Forecasting_CUACE-DUST_CMA_2019041500.nc （ht ...

王老师您好，测试数据我看是19年4月15-18日的，要想获得其他时段的数据去哪里下载呢？比如今年的，谢谢

MeteoInfo

xrpamhf 发表于 2023-5-25 15:28
王老师您好，测试数据我看是19年4月15-18日的，要想获得其他时段的数据去哪里下载呢？比如今年的，谢谢

这是内部数据，没有公开下载渠道。

xrpamhf

MeteoInfo 发表于 2023-5-25 16:40
这是内部数据，没有公开下载渠道。

好的 谢谢王老师