MeteoinfoLab脚本示例：3维水平、经向剖面图

︶ㄣ安定■丶

分享一个使用MiLab绘制3维水平、经向剖面图的脚本
(1) 3维水平剖面图

脚本如下：

1. # add data
   
2. f = addfile(r'D:\test\MERRA2_100.instM_3d_asm_Np.198001.SUB.nc')
   
3. RH = f['RH'][0,:,:,:]*100.0
   
4. lev1 = f['lev'][:]
   
5. 
   
6. lev2 = meteo.pressure_to_height_std(lev1)
   
7. lev2 = lev2/1000.0
   
8. RH.setdimvalue(0, lev2)
   
9. 
   
10. tdata_0km = RH['0',:,:]
    
11. tdata_2km = RH['2',:,:]
    
12. tdata_4km = RH['4',:,:]
    
13. tdata_6km = RH['6',:,:]
    
14. 
    
15. #plot
    
16. levs = arange(0,100,1)
    
17. cols=makecolors(len(levs)+1,cmap='NCV_jet')
    
18. ax = axes3d(position=[0.11, 0.11, 0.775, 0.88],tickfontsize=24,bbox=True)
    
19. china = shaperead('D:/Temp/map/china.shp')
    
20. world = shaperead('D:/Temp/map/country1.shp')
    
21. layer=ax.contourf(tdata_0km,levs,colors=cols,offset=0,alpha=1,zdir='z',smooth=True)
    
22. ax.plot_layer(china,edgecolor='gray',offset=0)
    
23. layer=ax.contourf(tdata_2km,levs,colors=cols,offset=2,alpha=1,zdir='z',smooth=True)
    
24. ax.plot_layer(china,edgecolor='gray',offset=2)
    
25. layer=ax.contourf(tdata_4km,levs,colors=cols,offset=4,alpha=1,zdir='z',smooth=True)
    
26. ax.plot_layer(china,edgecolor='gray',offset=4)
    
27. layer=ax.contourf(tdata_6km,levs,colors=cols,offset=6,alpha=1,zdir='z',smooth=True)
    
28. ax.plot_layer(china,edgecolor='gray',offset=6)
    
29. xlim(70,140)
    
30. ylim(15,55)
    
31. zlim(0,6)
    
32. ylabel('Latitude (degrees)',fontsize=24)
    
33. xlabel('Longitude (degrees)',fontsize=24)
    
34. zlabel('Altitude (km)',fontsize=24)
    
35. colorbar(layer,labelloc='right',ticks=arange(0,101,20),orientation='horizontal',yshift=25,label='RH (%)',fontsize=25,shrink=0.8,aspect=45)
    
36. antialias(True)
    

# add data<br /> f = addfile(r'D:\test\MERRA2_100.instM_3d_asm_Np.198001.SUB.nc')<br /> RH = f['RH'][0,:,:,:]*100.0<br /> lev1 = f['lev'][:]<br /> <br /> lev2 = meteo.pressure_to_height_std(lev1)<br /> lev2 = lev2/1000.0<br /> RH.setdimvalue(0, lev2)<br /> <br /> tdata_0km = RH['0',:,:]<br /> tdata_2km = RH['2',:,:]<br /> tdata_4km = RH['4',:,:]<br /> tdata_6km = RH['6',:,:]<br /> <br /> #plot<br /> levs = arange(0,100,1)<br /> cols=makecolors(len(levs)+1,cmap='NCV_jet')<br /> ax = axes3d(position=[0.11, 0.11, 0.775, 0.88],tickfontsize=24,bbox=True)<br /> china = shaperead('D:/Temp/map/china.shp')<br /> world = shaperead('D:/Temp/map/country1.shp')<br /> layer=ax.contourf(tdata_0km,levs,colors=cols,offset=0,alpha=1,zdir='z',smooth=True)<br /> ax.plot_layer(china,edgecolor='gray',offset=0)<br /> layer=ax.contourf(tdata_2km,levs,colors=cols,offset=2,alpha=1,zdir='z',smooth=True)<br /> ax.plot_layer(china,edgecolor='gray',offset=2)<br /> layer=ax.contourf(tdata_4km,levs,colors=cols,offset=4,alpha=1,zdir='z',smooth=True)<br /> ax.plot_layer(china,edgecolor='gray',offset=4)<br /> layer=ax.contourf(tdata_6km,levs,colors=cols,offset=6,alpha=1,zdir='z',smooth=True)<br /> ax.plot_layer(china,edgecolor='gray',offset=6)<br /> xlim(70,140)<br /> ylim(15,55)<br /> zlim(0,6)<br /> ylabel('Latitude (degrees)',fontsize=24)<br /> xlabel('Longitude (degrees)',fontsize=24)<br /> zlabel('Altitude (km)',fontsize=24)<br /> colorbar(layer,labelloc='right',ticks=arange(0,101,20),orientation='horizontal',yshift=25,label='RH (%)',fontsize=25,shrink=0.8,aspect=45)<br /> antialias(True)<br />

(2) 3维经向剖面图


脚本如下：

1. # add data
   
2. f = addfile(r'D:\test\MERRA2_100.instM_3d_asm_Np.198001.SUB.nc')
   
3. RH = f['RH'][0,:,:,:]*100.0
   
4. 
   
5. lev1 = f['lev'][:]
   
6. lev2 = meteo.pressure_to_height_std(lev1)
   
7. lev2 = lev2/1000.0
   
8. RH.setdimvalue(0, lev2)
   
9. 
   
10. #plot
    
11. levs = arange(0,100,1)
    
12. cols=makecolors(len(levs)+1,cmap='NCV_jet')
    
13. ax = axes3d(position=[0.11, 0.11, 0.775, 0.88],tickfontsize=24,bbox=True)
    
14. china = shaperead('D:/Temp/map/china.shp')
    
15. world = shaperead('D:/Temp/map/country1.shp')
    
16. ax.plot_layer(world,edgecolor='gray')
    
17. ax.plot_layer(china,edgecolor='k')
    
18. layer = ax.contourf(RH[:,:,'70'],levs,colors=cols,offset=70,alpha=0.8,zdir='x',smooth=True)
    
19. layer = ax.contourf(RH[:,:,'90'],levs,colors=cols,offset=90,alpha=0.8,zdir='x',smooth=True)
    
20. layer = ax.contourf(RH[:,:,'110'],levs,colors=cols,offset=110,alpha=0.8,zdir='x',smooth=True)
    
21. layer = ax.contourf(RH[:,:,'130'],levs,colors=cols,offset=130,alpha=0.8,zdir='x',smooth=True)
    
22. xlim(70,140)
    
23. ylim(15,55)
    
24. zlim(0,15)
    
25. ylabel('Latitude (degrees)',fontsize=24)
    
26. xlabel('Longitude (degrees)',fontsize=24)
    
27. zlabel('Altitude (km)',fontsize=24)
    
28. colorbar(layer,labelloc='right',ticks=arange(0,101,20),orientation='horizontal',yshift=25,label='RH (%)',fontsize=25,shrink=0.8,aspect=45)
    
29. antialias(True)
    

# add data<br /> f = addfile(r'D:\test\MERRA2_100.instM_3d_asm_Np.198001.SUB.nc')<br /> RH = f['RH'][0,:,:,:]*100.0<br /> <br /> lev1 = f['lev'][:]<br /> lev2 = meteo.pressure_to_height_std(lev1)<br /> lev2 = lev2/1000.0<br /> RH.setdimvalue(0, lev2)<br /> <br /> #plot<br /> levs = arange(0,100,1)<br /> cols=makecolors(len(levs)+1,cmap='NCV_jet')<br /> ax = axes3d(position=[0.11, 0.11, 0.775, 0.88],tickfontsize=24,bbox=True)<br /> china = shaperead('D:/Temp/map/china.shp')<br /> world = shaperead('D:/Temp/map/country1.shp')<br /> ax.plot_layer(world,edgecolor='gray')<br /> ax.plot_layer(china,edgecolor='k')<br /> layer = ax.contourf(RH[:,:,'70'],levs,colors=cols,offset=70,alpha=0.8,zdir='x',smooth=True)<br /> layer = ax.contourf(RH[:,:,'90'],levs,colors=cols,offset=90,alpha=0.8,zdir='x',smooth=True)<br /> layer = ax.contourf(RH[:,:,'110'],levs,colors=cols,offset=110,alpha=0.8,zdir='x',smooth=True)<br /> layer = ax.contourf(RH[:,:,'130'],levs,colors=cols,offset=130,alpha=0.8,zdir='x',smooth=True)<br /> xlim(70,140)<br /> ylim(15,55)<br /> zlim(0,15)<br /> ylabel('Latitude (degrees)',fontsize=24)<br /> xlabel('Longitude (degrees)',fontsize=24)<br /> zlabel('Altitude (km)',fontsize=24)<br /> colorbar(layer,labelloc='right',ticks=arange(0,101,20),orientation='horizontal',yshift=25,label='RH (%)',fontsize=25,shrink=0.8,aspect=45)<br /> antialias(True)<br />

总结： 纬向剖面和径向类似，这里就不再展示，大家可以根据示例文件尝试，欢迎跟帖。
附件为示例文件：Merra-2 的RH数据

独孤酒见

可以占沙发了，顶鸭，酷毙了！

lu19320

顶，感谢分享！！！

咖喱嫑洋葱

支持一下。。

咖喱嫑洋葱

支持一下。。

兔子爱吃萝卜

{:eb502:}{:eb502:}{:eb502:}

遥感卫星图像

感谢老师分享

wwggiser

支持，向老师学习

MeteoInfo

已加入脚本汇总帖。

看云淡风清

楼主赞一个，谢谢分享