绘制任意两点连线的剖面

MeteoInfo

MeteoInfo 3.9.1版本在meteolib包中增加了cross_section函数，可以通过起始和结束点的经纬度获取三维或者二维格点数据两点连线处的剖面数据。

详见MeteoInfo公众号：https://mp.weixin.qq.com/s/_AkarrVz4lnP5cdj3l1QLg

1. #Read data
   
2. fn = os.path.join(migl.get_sample_folder(), 'GrADS', 'model.ctl')
   
3. f = addfile(fn)
   
4. data = f['U'][0]
   
5. lev1 = data.dimvalue(0)
   
6. lev2 = meteo.pressure_to_height_std(lev1)
   
7. lev2 = lev2 / 1000
   
8. data.setdimvalue(0, lev2)
   
9. 
   
10. f1 = addfile('D:/Temp/nc/elev.0.25-deg.nc')
    
11. elev = f1['data'][0,::-1]
    
12. 
    
13. lon1 = 0; lat1 = -30; lon2 = 180; lat2 = 60
    
14. n = 73
    
15. terrain = meteo.cross_section(elev, start=(lon1, lat1), end=(lon2, lat2),
    
16.     steps=n)[0]
    
17. terrain[terrain<0] = 0
    
18. terrain = terrain / 1000.0
    
19. csdata, x, y = meteo.cross_section(data, start=(lon1, lat1),
    
20.     end=(lon2, lat2), steps=n)
    
21. 
    
22. #Plot
    
23. subplot(2,1,1)
    
24. yaxis(tickin=False)
    
25. xaxis(tickin=False)
    
26. yaxis(location='right', tickvisible=False)
    
27. xaxis(location='top', tickvisible=False)
    
28. levs = arange(-10, 81, 10)
    
29. layer = contourf(x, lev2, csdata, levs)
    
30. fill_between(x, terrain, color='gray')
    
31. plot(x, terrain, color='k')
    
32. ylabel('Hight(km)')
    
33. xlabel('Longitude')
    
34. ylim(0, 16)
    
35. colorbar(layer, aspect=25)
    
36. 
    
37. subplot(2,1,2,axestype='map')
    
38. lat = [lat1, lat2]
    
39. lon = [lon1, lon2]
    
40. geoshow(lat, lon, size=2, color='b')
    
41. geoshow('country')
    
42. xlim(-10, 190)
    
43. ylim(-55, 65)

#Read data<br /> fn = os.path.join(migl.get_sample_folder(), 'GrADS', 'model.ctl')<br /> f = addfile(fn)<br /> data = f['U'][0]<br /> lev1 = data.dimvalue(0)<br /> lev2 = meteo.pressure_to_height_std(lev1)<br /> lev2 = lev2 / 1000<br /> data.setdimvalue(0, lev2)<br /> <br /> f1 = addfile('D:/Temp/nc/elev.0.25-deg.nc')<br /> elev = f1['data'][0,::-1]<br /> <br /> lon1 = 0; lat1 = -30; lon2 = 180; lat2 = 60<br /> n = 73<br /> terrain = meteo.cross_section(elev, start=(lon1, lat1), end=(lon2, lat2), <br />     steps=n)[0]<br /> terrain[terrain<0] = 0<br /> terrain = terrain / 1000.0<br /> csdata, x, y = meteo.cross_section(data, start=(lon1, lat1), <br />     end=(lon2, lat2), steps=n)<br /> <br /> #Plot<br /> subplot(2,1,1)<br /> yaxis(tickin=False)<br /> xaxis(tickin=False)<br /> yaxis(location='right', tickvisible=False)<br /> xaxis(location='top', tickvisible=False)<br /> levs = arange(-10, 81, 10)<br /> layer = contourf(x, lev2, csdata, levs)<br /> fill_between(x, terrain, color='gray')<br /> plot(x, terrain, color='k')<br /> ylabel('Hight(km)')<br /> xlabel('Longitude')<br /> ylim(0, 16)<br /> colorbar(layer, aspect=25)<br /> <br /> subplot(2,1,2,axestype='map')<br /> lat = [lat1, lat2]<br /> lon = [lon1, lon2]<br /> geoshow(lat, lon, size=2, color='b')<br /> geoshow('country')<br /> xlim(-10, 190)<br /> ylim(-55, 65)