- 积分
- 55955
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2011-6-21
- 最后登录
- 1970-1-1
|
楼主 |
发表于 2021-9-3 16:57:34
|
显示全部楼层
可以,比较麻烦一些,有空了可以考虑做些简化。
- fn = os.path.join(migl.get_sample_folder(), 'GrADS', 'model.ctl')
- f = addfile(fn)
- u = f['U'][0,:,'20','0']
- lev1 = meteo.pressure_to_height_std(u.dimvalue(0)) / 1000
- lon1 = 0.
- lon2 = 180.
- lat1 = -50.
- lat2 = 50.
- lon = lon1
- lat_c = []
- tdata = []
- #tdata_terrain = []
- while lon <= lon2:
- lat = lat1 + (lat2 - lat1) * (lon - lon1) / (lon2 - lon1)
- lat_c.append(lat)
- print lon
- tdata.append(f['U'][0,:,'%f'%lat,'%f'%lon])
- lon = lon + 2.5
-
- alldata = concatenate(tdata, axis=0)
- alldata = reshape(alldata, (len(tdata), len(lev1)))
- alldata = alldata.T
- x = arange(lon1, lon2 + 1, 2.5)
- #Plot
- ax = axes3d(tickfontsize=14)
- geoshow('continent', color='c', edgecolor='b')
- ls = ax.contourf(x, lev1, alldata, 10, edgecolor='k', zdir='xy', alpha=0.8, \
- sepoint=[lon1,lat1,lon2,lat2])
- colorbar(ls)
- xlabel('Longitude (degrees)', fontsize=16)
- ylabel('Latitude (degrees)', fontsize=16)
- zlabel('Altitude (km)', fontsize=16)
- xlim(0, 180)
- title('vertical slice plot example')
|
|