MeteoinfoLab脚本示例：中国地形数据提取及3维地形图

︶ㄣ安定■丶

基于之前网上分享的一个全球地形数据（忘哪里下的了），现将这个数据用Meteoinfola脚本ncwrite函数提取出纬度在[15,55]，经度在[70,137]范围内的数据，用于绘制中国地形图。

脚本如下：

1. fn = 'H:/Meteoinfo/dixing.nc'
   
2. f = addfile(fn)
   
3. tdata = f['z']['15:55','70:137']
   
4. ncwrite('H:/Meteoinfo/dixing-china.nc',tdata,'z')
   
5. print 'Finish...'

fn = 'H:/Meteoinfo/dixing.nc'<br /> f = addfile(fn)<br /> tdata = f['z']['15:55','70:137']<br /> ncwrite('H:/Meteoinfo/dixing-china.nc',tdata,'z')<br /> print 'Finish...'

提取出来的地形数据，请从附件下载 dixing-china.nc (4.62 MB, 下载次数: 1231)

2016-12-6 18:15 上传

点击文件名下载附件


另外这个数据也能通过如下脚本绘制中国地区地形图：

1. fn = 'H:/Meteoinfo/dixing-china.nc'
   
2. f = addfile(fn)
   
3. data = f['z'][:,:]
   
4. 
   
5. #Plot
   
6. axesm(tickfontsize=16)
   
7. world = shaperead('D:/Temp/map/country1.shp')
   
8. lchina = shaperead('D:/Temp/map/bou2_4p.shp')
   
9. geoshow(world)
   
10. geoshow(lchina,edgecolor='k',size=0.3)
    
11. levels = arange(0, 5200, 400)
    
12. layer = imshowm(data,levels,cmap='temp_diff_18lev')
    
13. masklayer(lchina, [layer])
    
14. xlim(70,136)
    
15. ylim(15,55)
    
16. xlabel('Longitude',fontsize=17)
    
17. ylabel('Latitude',fontsize=17)
    
18. colorbar(layer,shrink=1,label='Height(m)')

fn = 'H:/Meteoinfo/dixing-china.nc'<br /> f = addfile(fn)<br /> data = f['z'][:,:] <br /> <br /> #Plot<br /> axesm(tickfontsize=16) <br /> world = shaperead('D:/Temp/map/country1.shp') <br /> lchina = shaperead('D:/Temp/map/bou2_4p.shp')<br /> geoshow(world) <br /> geoshow(lchina,edgecolor='k',size=0.3)<br /> levels = arange(0, 5200, 400)<br /> layer = imshowm(data,levels,cmap='temp_diff_18lev')<br /> masklayer(lchina, [layer])<br /> xlim(70,136)<br /> ylim(15,55)<br /> xlabel('Longitude',fontsize=17)<br /> ylabel('Latitude',fontsize=17)<br /> colorbar(layer,shrink=1,label='Height(m)')

绘制地形图如下所示：

也可以直接通过meteoinfo直接打开绘制。

基于MIlab的强大的axes3d函数，使用提取出来的中国地形数据，简单的绘制了一下四川盆地的三维地形图：

脚本如下：

1. #read data
   
2. fn = r'D:\test/dixing-china.nc'
   
3. f = addfile(fn)
   
4. elev = f['z']['26:34','102:108']/1000.0
   
5. lchina = shaperead('D:/Temp/map/bou2_4p.shp')
   
6. 
   
7. #Plot
   
8. figure(figsize=[1400,850],newfig=False)
   
9. ax = axes3d(tickfontsize=20,position=[0.1, 0.13, 0.75, 0.85],bbox=False)
   
10. levs = arange(0,6.1,0.1)
    
11. cols = makecolors(len(levs) + 1, cmap='temp_diff_18lev',reverse=False)
    
12. ls = ax.plot_surface(elev, levs, colors=cols, edge=False)
    
13. colorbar(ls,ticks=arange(0,6.1,1),orientation='horizontal',label='Altitude(km)',labelloc='right',fontsize=24,shrink=0.8,aspect=40,yshift=35)
    
14. ylabel('Latitude (degrees)',fontsize=20)
    
15. xlabel('Longitude (degrees)',fontsize=20)
    
16. zlabel('Altitude (km)',fontsize=20)
    
17. zlim(0, 6)
    
18. xlim(102, 108)
    
19. ylim(26,34)

#read data<br /> fn = r'D:\test/dixing-china.nc'<br /> f = addfile(fn)<br /> elev = f['z']['26:34','102:108']/1000.0<br /> lchina = shaperead('D:/Temp/map/bou2_4p.shp')<br /> <br /> #Plot<br /> figure(figsize=[1400,850],newfig=False)<br /> ax = axes3d(tickfontsize=20,position=[0.1, 0.13, 0.75, 0.85],bbox=False)<br /> levs = arange(0,6.1,0.1)<br /> cols = makecolors(len(levs) + 1, cmap='temp_diff_18lev',reverse=False)<br /> ls = ax.plot_surface(elev, levs, colors=cols, edge=False)<br /> colorbar(ls,ticks=arange(0,6.1,1),orientation='horizontal',label='Altitude(km)',labelloc='right',fontsize=24,shrink=0.8,aspect=40,yshift=35)<br /> ylabel('Latitude (degrees)',fontsize=20)<br /> xlabel('Longitude (degrees)',fontsize=20)<br /> zlabel('Altitude (km)',fontsize=20)<br /> zlim(0, 6)<br /> xlim(102, 108)<br /> ylim(26,34)

MeteoInfo

lyh0205 发表于 2020-4-30 09:10
还没有弄清楚

MeteoInfoLab软件Figure区域有专门的Rotate工具按钮，选中后可以手动旋转三维图形。或者用三维图形坐标系的set_elevation和set_rotation方法设置：

ax.set_elevation(-30)
ax.set_rotation(70)

寻路人

{:eb502:}{:eb502:}

andrewsoong

这个精细程度，不像一度的分辨率的啊

︶ㄣ安定■丶

andrewsoong 发表于 2016-12-6 20:52
这个精细程度，不像一度的分辨率的啊

哈哈，我看大概是0.033度的样子。。不过画一般的区域性地形是没有问题的

andrewsoong

︶ㄣ安定■丶 发表于 2016-12-6 21:01
哈哈，我看大概是0.033度的样子。。不过画一般的区域性地形是没有问题的

是的，很好了

贫道敬孔

多谢楼主分享{:5_213:}

悠哉~鱿鱼~

很棒！学习了！希望楼主多分享学习经验~

空谷幽竹

多谢楼主分享，很有用

寻路人

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

fjtzf

{:eb502:}{:eb502:}