MeteoInfoLab脚本示例：站点数据散点图

MeteoInfo

这里演示从micaps第一类数据（地面全要素观测）中读取一个变量（用DimDataFile类的stationdata方法），然后maskout掉中国区域之外的数据，利用scatterm函数绘制散点图。


脚本程序：

1. fn = os.path.join(migl.get_sample_folder(), 'MICAPS', '10101414.000')
   
2. f = addfile(fn)
   
3. data = f['Visibility'][:]
   
4. lon = f['Longitude'][:]
   
5. lat = f['Latitude'][:]
   
6. lchina = shaperead('china')
   
7. data, lon, lat = rmaskout(data, lon, lat, lchina.shapes())
   
8. t = f.gettime(0)
   
9. 
   
10. #Plot
    
11. geoshow('cn_province', edgecolor='gray')
    
12. geoshow('country')
    
13. levs = arange(3, 31, 3)
    
14. layer = scatter(lon, lat, data, levs, cmap='matlab_jet_r')
    
15. colorbar(layer)
    
16. yticks([20,30,40,50])
    
17. title('Visibility (' + t.strftime('%Y-%m-%d %H:00') + ')')
    
18. xlim(72, 136)
    
19. ylim(16, 55)
    
20. 
    
21. #Add south China Sea
    
22. axesm(position=[0.105,0.12,0.18,0.24], bgcolor='w', axison=False)
    
23. geoshow('cn_border', facecolor=(0,0,255))
    
24. xlim(106, 123)
    
25. ylim(2, 23)
    
26. 
    
27. #savefig('D:/Temp/figure/micaps1_points.png', 600, 350, dpi=300)

fn = os.path.join(migl.get_sample_folder(), 'MICAPS', '10101414.000')<br /> f = addfile(fn)<br /> data = f['Visibility'][:]<br /> lon = f['Longitude'][:]<br /> lat = f['Latitude'][:]<br /> lchina = shaperead('china')<br /> data, lon, lat = rmaskout(data, lon, lat, lchina.shapes())<br /> t = f.gettime(0)<br /> <br /> #Plot<br /> geoshow('cn_province', edgecolor='gray')<br /> geoshow('country')<br /> levs = arange(3, 31, 3)<br /> layer = scatter(lon, lat, data, levs, cmap='matlab_jet_r')<br /> colorbar(layer)<br /> yticks([20,30,40,50])<br /> title('Visibility (' + t.strftime('%Y-%m-%d %H:00') + ')')<br /> xlim(72, 136)<br /> ylim(16, 55)<br /> <br /> #Add south China Sea<br /> axesm(position=[0.105,0.12,0.18,0.24], bgcolor='w', axison=False)<br /> geoshow('cn_border', facecolor=(0,0,255))<br /> xlim(106, 123)<br /> ylim(2, 23)<br /> <br /> #savefig('D:/Temp/figure/micaps1_points.png', 600, 350, dpi=300)

平流层的萝卜

真给力，赞一个

wsjqxjy

谢谢楼主共享.

rceclx

和王老师一起进步！

半颗麦芽糖

王老师，如果我想把这个图画成填色的该怎么画呢？

MeteoInfo

半颗麦芽糖 发表于 2016-12-14 16:30
王老师，如果我想把这个图画成填色的该怎么画呢？

参考此贴：MeteoInfoLab脚本示例：站点数据绘制等值线
http://bbs.06climate.com/forum.p ... 422&fromuid=106
(出处: 气象家园)

半颗麦芽糖

王老师，我参照您的帖子，改了脚本成功出图，但是当我改色标的颜色和大小的时候就报错了，不知道为什么，附脚本和图：

1. 
   
2. <P>#用这个脚本可以出图</P>
   
3. <P>f = addfile_micaps('E:/test/14102708.000')
   
4. pr = f.stationdata('Var1')
   
5. #griddata function - interpolate
   
6. x = arange(70, 120, 0.5)
   
7. y = arange(20, 50, 0.5)
   
8. prg = pr.griddata((x, y), method='idw', radius=3)
   
9. #Plot
   
10. axesm()
    
11. lchina = shaperead('E:/MeteoInfo/map/bou2_4p.shp')
    
12. ltibet = shaperead('E:/MeteoInfo/map/DBATP/tibet_ASCII.wmp')
    
13. geoshow(lchina, edgecolor='k')
    
14. geoshow(ltibet, edgecolor='r',size=2)
    
15. levs = [0.1, 1, 2, 5, 10, 20, 25, 50, 100]
    
16. cols = [(255,255,255),(170,240,255),(120,230,240),(200,220,50),(240,220,20),(255,120,10),(255,90,10), \
    
17.     (240,40,0),(180,10,0),(120,10,0)]
    
18. layer = contourfm(prg, levs, colors=cols)
    
19. masklayer(lchina, [layer])
    
20. colorbar(layer)
    
21. xlim(70, 120)
    
22. ylim(20, 50)</P>

复制代码

当时当我把
levs = [0.1, 1, 2, 5, 10, 20, 25, 50, 100]
cols = [(255,255,255),(170,240,255),(120,230,240),(200,220,50),(240,220,20),(255,120,10),(255,90,10), \
    (240,40,0),(180,10,0),(120,10,0)]
改成这两句
levs = [0,10,20,30,40,50,60,70]
cols = makecolors(len(levs)+1)
或者任意修改levs的时候，报错，请王老师指导，谢谢，最后贴上数据
14102708.000 (2.89 KB, 下载次数: 4)

2016-12-14 19:30 上传

点击文件名下载附件

MeteoInfo

半颗麦芽糖 发表于 2016-12-14 19:28
王老师，我参照您的帖子，改了脚本成功出图，但是当我改色标的颜色和大小的时候就报错了，不知道为什么，附 ...

是程序里的bug，会尽快修正。

半颗麦芽糖

MeteoInfo 发表于 2016-12-14 22:16
是程序里的bug，会尽快修正。

恩恩，等待王老师修复。

MeteoInfo

半颗麦芽糖 发表于 2016-12-14 22:36
恩恩，等待王老师修复。

发布了新版本，你下载试试吧