MeteoInfoLab脚本示例：斜剖面（2D、3D显示）

︶ㄣ安定■丶

结果2D图绘制如下（叠加了地形）：


具体脚本如下：

1. figure(figsize=[1500,400],newfig=False)
   
2. fig,((ax1,ax2)) = subplots(1,2,position=[0,0.08,1,0.89])
   
3. f = addfile(r'H:\test\data\test.nc')
   
4. fn1 = addfile(r'H:\alldata\dixing/elev.0.5-deg.nc')
   
5. lev1 = f['level'][:]
   
6. lev2 = meteo.pressure_to_height_std(lev1)
   
7. lev2 = lev2[::-1]/1000
   
8. 
   
9. lon1 = 80
   
10. lon2 = 120
    
11. lat1 = 32
    
12. lat2 = 40
    
13. lon = lon1
    
14. lat_c = []
    
15. tdata = []
    
16. tdata_terrain = []
    
17. while lon<=lon2:
    
18.     lat = lat1+(lat2-lat1)*(lon-lon1)/(lon2-lon1)
    
19.     lat_c.append(lat)
    
20.     print lat
    
21.     tdata.append(f['w'][0,:,'%f'%lat,'%f'%lon])
    
22.     tdata_terrain.append(fn1['data']['%f'%lat,'%f'%lon])
    
23.     lon = lon+1
    
24.    
    
25. tdata_terrain=array(tdata_terrain)   
    
26. tdata_terrain=tdata_terrain/1000.0
    
27. alldata=concatenate(tdata,axis=0)
    
28. alldata=reshape(alldata,(len(tdata),len(lev1)))   
    
29. alldata=alldata.T
    
30. alldata = alldata[::-1,:]
    
31. 
    
32. axes(newaxes=False)   
    
33. x=arange(lon1,lon2+1,1)
    
34. levs = arange(-0.2,0.22,0.02)
    
35. yaxis(tickin=False,tickfontsize=17)
    
36. xaxis(tickin=False,tickfontsize=17)
    
37. yaxis(location='right',tickin=True,tickfontsize=18)
    
38. layer = contourf(x,lev2,alldata,levs,cmap='BlRe')
    
39. fill_between(x,tdata_terrain,color='gray')
    
40. plot(x,tdata_terrain,color='k')
    
41. ylabel('Hight(km)',fontsize=18,bold=False)
    
42. xlabel('Longitude',fontsize=18)
    
43. ylim(lev2.min(),12.001)
    
44. colorbar(layer,aspect=25)
    
45. 
    
46. caxes(ax2)
    
47. axesm(newaxes=False,tickfontsize=17)
    
48. lat = [32, 40]
    
49. lon = [100, 120]
    
50. geoshow(lat, lon, size=2, color='b')
    
51. lchina = shaperead('D:/Temp/map/bou2_4p.shp')
    
52. geoshow(lchina,edgecolor='k',size=1)
    
53. xlim(70,140)
    
54. ylim(15,55)
    
55. antialias(True)
    

figure(figsize=[1500,400],newfig=False)<br /> fig,((ax1,ax2)) = subplots(1,2,position=[0,0.08,1,0.89])<br /> f = addfile(r'H:\test\data\test.nc')<br /> fn1 = addfile(r'H:\alldata\dixing/elev.0.5-deg.nc')<br /> lev1 = f['level'][:]<br /> lev2 = meteo.pressure_to_height_std(lev1)<br /> lev2 = lev2[::-1]/1000<br /> <br /> lon1 = 80<br /> lon2 = 120<br /> lat1 = 32<br /> lat2 = 40<br /> lon = lon1<br /> lat_c = []<br /> tdata = []<br /> tdata_terrain = []<br /> while lon<=lon2:<br />     lat = lat1+(lat2-lat1)*(lon-lon1)/(lon2-lon1)<br />     lat_c.append(lat)<br />     print lat<br />     tdata.append(f['w'][0,:,'%f'%lat,'%f'%lon])<br />     tdata_terrain.append(fn1['data']['%f'%lat,'%f'%lon])<br />     lon = lon+1<br />     <br /> tdata_terrain=array(tdata_terrain)    <br /> tdata_terrain=tdata_terrain/1000.0<br /> alldata=concatenate(tdata,axis=0)<br /> alldata=reshape(alldata,(len(tdata),len(lev1)))    <br /> alldata=alldata.T<br /> alldata = alldata[::-1,:]<br /> <br /> axes(newaxes=False)    <br /> x=arange(lon1,lon2+1,1) <br /> levs = arange(-0.2,0.22,0.02)<br /> yaxis(tickin=False,tickfontsize=17)<br /> xaxis(tickin=False,tickfontsize=17)<br /> yaxis(location='right',tickin=True,tickfontsize=18)<br /> layer = contourf(x,lev2,alldata,levs,cmap='BlRe')<br /> fill_between(x,tdata_terrain,color='gray')<br /> plot(x,tdata_terrain,color='k')<br /> ylabel('Hight(km)',fontsize=18,bold=False)<br /> xlabel('Longitude',fontsize=18)<br /> ylim(lev2.min(),12.001)<br /> colorbar(layer,aspect=25)<br /> <br /> caxes(ax2)<br /> axesm(newaxes=False,tickfontsize=17)<br /> lat = [32, 40]<br /> lon = [100, 120]<br /> geoshow(lat, lon, size=2, color='b')<br /> lchina = shaperead('D:/Temp/map/bou2_4p.shp')<br /> geoshow(lchina,edgecolor='k',size=1)<br /> xlim(70,140)<br /> ylim(15,55)<br /> antialias(True)<br />

看到有3D显示斜剖面的需求，感谢王老师及时进行了更新，现在可以将斜剖面显示在3D地图上了。
结果3D图绘制如下（叠加了地形）：


3D显示斜剖面具体脚本如下：

1. f = addfile(r'H:\test\data\test.nc')
   
2. fn1 = addfile(r'H:\alldata\dixing/elev.0.5-deg.nc')
   
3. lev1 = f['level'][:]
   
4. lev2 = meteo.pressure_to_height_std(lev1)
   
5. lev2 = lev2[::-1]/1000
   
6. 
   
7. lon1 = 80
   
8. lon2 = 120
   
9. lat1 = 32
   
10. lat2 = 40
    
11. lon = lon1
    
12. x = []
    
13. y = []
    
14. tdata = []
    
15. tdata_terrain = []
    
16. while lon<=lon2:
    
17.     lat = lat1+(lat2-lat1)*(lon-lon1)/(lon2-lon1)
    
18.     print lon
    
19.     tdata.append(f['w'][0,:,'%f'%lat,'%f'%lon])  
    
20.     tdata_terrain.append(fn1['data']['%f'%lat,'%f'%lon])
    
21.     x.append(lon)
    
22.     y.append(lat)
    
23.     lon = lon+1
    
24.    
    
25. tdata_terrain=array(tdata_terrain)   
    
26. tdata_terrain[tdata_terrain<0] = 0
    
27. tdata_terrain=tdata_terrain/1000.0
    
28. 
    
29. alldata=concatenate(tdata,axis=0)
    
30. alldata=reshape(alldata,(len(tdata),len(lev1)))   
    
31. alldata=alldata.T
    
32. alldata = alldata[::-1,:]
    
33. x=array(x)
    
34. y=array(y)
    
35. 
    
36. ax = axes3d(position=[0.11, 0.11, 0.775, 0.815],tickfontsize=17,bbox=False)
    
37. levs = arange(-0.2,0.22,0.02)
    
38. china = shaperead('D:/Temp/map/china.shp')
    
39. world = shaperead('D:/Temp/map/country1.shp')
    
40. ax.plot_layer(china,edgecolor='gray')
    
41. ax.plot_layer(world,edgecolor='gray')
    
42. layer=ax.contourf(x,lev2,alldata,levs,zdir='xy',alpha=0.8,sepoint=[lon1,lat1,lon2,lat2],cmap='BlRe')
    
43. k =(lat2-lat1)/float((lon2-lon1))
    
44. y1=k*(x-lon1)+lat1
    
45. ax.plot(x,y1,tdata_terrain,'-k',linewidth=1)
    
46. ax.fill_between(x,tdata_terrain, zdir='xy',color='gray',y=y1)
    
47. 
    
48. xlim(70,140)
    
49. ylim(15,55)
    
50. zlim(0,12)
    
51. colorbar(layer,aspect=30)
    
52. ylabel('Latitude (degrees)',fontsize=17)
    
53. xlabel('Longitude (degrees)',fontsize=17)
    
54. zlabel(label='Altitude (km)',fontsize=17)
    
55. antialias(True)
    

f = addfile(r'H:\test\data\test.nc')<br /> fn1 = addfile(r'H:\alldata\dixing/elev.0.5-deg.nc')<br /> lev1 = f['level'][:]<br /> lev2 = meteo.pressure_to_height_std(lev1)<br /> lev2 = lev2[::-1]/1000<br /> <br /> lon1 = 80<br /> lon2 = 120<br /> lat1 = 32<br /> lat2 = 40<br /> lon = lon1<br /> x = []<br /> y = []<br /> tdata = []<br /> tdata_terrain = []<br /> while lon<=lon2:<br />     lat = lat1+(lat2-lat1)*(lon-lon1)/(lon2-lon1)<br />     print lon<br />     tdata.append(f['w'][0,:,'%f'%lat,'%f'%lon])  <br />     tdata_terrain.append(fn1['data']['%f'%lat,'%f'%lon])<br />     x.append(lon)<br />     y.append(lat)<br />     lon = lon+1<br />     <br /> tdata_terrain=array(tdata_terrain)    <br /> tdata_terrain[tdata_terrain<0] = 0<br /> tdata_terrain=tdata_terrain/1000.0<br /> <br /> alldata=concatenate(tdata,axis=0)<br /> alldata=reshape(alldata,(len(tdata),len(lev1)))   <br /> alldata=alldata.T<br /> alldata = alldata[::-1,:]<br /> x=array(x)<br /> y=array(y)<br /> <br /> ax = axes3d(position=[0.11, 0.11, 0.775, 0.815],tickfontsize=17,bbox=False) <br /> levs = arange(-0.2,0.22,0.02)<br /> china = shaperead('D:/Temp/map/china.shp')<br /> world = shaperead('D:/Temp/map/country1.shp')<br /> ax.plot_layer(china,edgecolor='gray')<br /> ax.plot_layer(world,edgecolor='gray')<br /> layer=ax.contourf(x,lev2,alldata,levs,zdir='xy',alpha=0.8,sepoint=[lon1,lat1,lon2,lat2],cmap='BlRe')<br /> k =(lat2-lat1)/float((lon2-lon1))<br /> y1=k*(x-lon1)+lat1<br /> ax.plot(x,y1,tdata_terrain,'-k',linewidth=1)<br /> ax.fill_between(x,tdata_terrain, zdir='xy',color='gray',y=y1)<br /> <br /> xlim(70,140)<br /> ylim(15,55)<br /> zlim(0,12)<br /> colorbar(layer,aspect=30)<br /> ylabel('Latitude (degrees)',fontsize=17)<br /> xlabel('Longitude (degrees)',fontsize=17)<br /> zlabel(label='Altitude (km)',fontsize=17)<br /> antialias(True)<br />

woiaml

这个比较厉害！

xyan88

6666666 能不能把剖面直接加到地图上面啊

︶ㄣ安定■丶

xyan88 发表于 2017-11-19 11:44
6666666 能不能把剖面直接加到地图上面啊

这个比较简单，你参考这个帖子应该就能加上http://bbs.06climate.com/forum.p ... mp;page=1#pid843774

︶ㄣ安定■丶

︶ㄣ安定■丶 发表于 2017-11-19 11:56
这个比较简单，你参考这个帖子应该就能加上http://bbs.06climate.com/forum.php?mod=viewthread&tid=5690 ...

已经在楼上脚本中更新了，叠加地形

︶ㄣ安定■丶

xyan88 发表于 2017-11-19 11:44
6666666 能不能把剖面直接加到地图上面啊

看错了，你指的3D的斜剖面，暂时还不行

独孤酒见

安定老师和王老师真是高人，且是热心人。好好学习一下。

river

感谢分享，meteoinfo用法和脚本示例越来越丰富了！

︶ㄣ安定■丶

xyan88 发表于 2017-11-19 11:44
6666666 能不能把剖面直接加到地图上面啊

已经在楼顶加了3D显示斜剖面的功能，请参考

小1龙2虾3

安老师，腻害哦！