MeteoInfoLab脚本示例：站点数据绘制等值线

MeteoInfo

站点数据绘制等值线需要首先将站点数据插值为格点数据，MeteoInfo中提供了反距离权法（IDW）和cressman两个方法，其中IDW方法可以有插值半径的选项。这里示例读取一个MICAPS第一类数据（地面全要素观测），获取6小时累积降水数据（Precipitation6h），然后用站点数据的griddata函数将站点数据插值为格点数据，再利用contourfm函数创建等值线填色图层（等值线间隔和颜色可以自定义）。

脚本程序（经纬度投影）：

1. #Read station data
   
2. fn = os.path.join(migl.get_sample_folder(), 'MICAPS', '10101414.000')
   
3. f = addfile_micaps(fn)
   
4. pr = f['Precipitation6h'][:]
   
5. lon = f['Longitude'][:]
   
6. lat = f['Latitude'][:]
   
7. 
   
8. #griddata function - interpolate
   
9. x = arange(75, 135, 0.5)
   
10. y = arange(18, 55, 0.5)
    
11. prg, gx, gy = griddata((lon, lat), pr, xi=(x, y), method='cressman', radius=[7,4,2,1])
    
12. 
    
13. #Plot
    
14. figure(figsize=[700, 400], newfig=False)
    
15. axesm()
    
16. geoshow('cn_province', edgecolor='lightgray')
    
17. bou1_layer = geoshow('cn_border', facecolor=(0,0,255))
    
18. city_layer = geoshow('cn_cities', facecolor='r', size=4)
    
19. city_layer.addlabels('NAME', fontname=u'楷体', fontsize=16, yoffset=15, avoidcoll=False)
    
20. china_layer = geoshow('china', visible=False)
    
21. city_layer.movelabel(u'西宁', -15)
    
22. levs = [0.1, 1, 2, 5, 10, 20, 25, 50, 100]
    
23. cols = [(255,255,255),(170,240,255),(120,230,240),(200,220,50),(240,220,20),(255,120,10),(255,90,10), \
    
24.     (240,40,0),(180,10,0),(120,10,0)]
    
25. layer = contourf(x, y, prg, levs, colors=cols)
    
26. masklayer(china_layer, [layer])
    
27. colorbar(layer)
    
28. xlim(72, 136)
    
29. ylim(16, 55)
    
30. text(95, 52, u'全国降水量实况图', fontname=u'黑体', fontsize=16)
    
31. text(95, 50, u'(2010-10-14 08:00 至 2010-10-14 14:00)', fontname=u'黑体', fontsize=14)
    
32. 
    
33. #Add south China Sea
    
34. sc_layer = bou1_layer.clone()
    
35. axesm(position=[0.14,0.18,0.15,0.2], axison=False, frameon=True)
    
36. geoshow(sc_layer, facecolor=(0,0,255))
    
37. xlim(106, 123)
    
38. ylim(2, 23)
    
39. #savefig('D:/Temp/test/test_pre.png')

#Read station data<br /> fn = os.path.join(migl.get_sample_folder(), 'MICAPS', '10101414.000')<br /> f = addfile_micaps(fn)<br /> pr = f['Precipitation6h'][:]<br /> lon = f['Longitude'][:]<br /> lat = f['Latitude'][:]<br /> <br /> #griddata function - interpolate<br /> x = arange(75, 135, 0.5)<br /> y = arange(18, 55, 0.5)<br /> prg, gx, gy = griddata((lon, lat), pr, xi=(x, y), method='cressman', radius=[7,4,2,1])<br /> <br /> #Plot<br /> figure(figsize=[700, 400], newfig=False)<br /> axesm()<br /> geoshow('cn_province', edgecolor='lightgray')<br /> bou1_layer = geoshow('cn_border', facecolor=(0,0,255))<br /> city_layer = geoshow('cn_cities', facecolor='r', size=4)<br /> city_layer.addlabels('NAME', fontname=u'楷体', fontsize=16, yoffset=15, avoidcoll=False)<br /> china_layer = geoshow('china', visible=False)<br /> city_layer.movelabel(u'西宁', -15)<br /> levs = [0.1, 1, 2, 5, 10, 20, 25, 50, 100]<br /> cols = [(255,255,255),(170,240,255),(120,230,240),(200,220,50),(240,220,20),(255,120,10),(255,90,10), \<br />     (240,40,0),(180,10,0),(120,10,0)]<br /> layer = contourf(x, y, prg, levs, colors=cols)<br /> masklayer(china_layer, [layer])<br /> colorbar(layer)<br /> xlim(72, 136)<br /> ylim(16, 55)<br /> text(95, 52, u'全国降水量实况图', fontname=u'黑体', fontsize=16)<br /> text(95, 50, u'(2010-10-14 08:00 至 2010-10-14 14:00)', fontname=u'黑体', fontsize=14)<br /> <br /> #Add south China Sea<br /> sc_layer = bou1_layer.clone()<br /> axesm(position=[0.14,0.18,0.15,0.2], axison=False, frameon=True)<br /> geoshow(sc_layer, facecolor=(0,0,255))<br /> xlim(106, 123)<br /> ylim(2, 23)<br /> #savefig('D:/Temp/test/test_pre.png')

脚本程序（Lambert投影）：

1. #Set data folders
   
2. 
   
3. basedir = 'D:/MyProgram/Distribution/java/MeteoInfo/MeteoInfo'
   
4. datadir = os.path.join(basedir, 'sample/MICAPS')
   
5. 
   
6. #Read station data
   
7. f = addfile_micaps(os.path.join(datadir, '10101414.000'))
   
8. lon = f['Longitude'][:]
   
9. lat = f['Latitude'][:]
   
10. pr = f['Precipitation6h'][:]
    
11. 
    
12. #griddata function - interpolate
    
13. x = arange(75, 135, 0.5)
    
14. y = arange(18, 55, 0.5)
    
15. #prg = griddata((lon, lat), pr, xi=(x, y), method='idw', radius=3)[0]
    
16. prg = griddata((lon, lat), pr, xi=(x, y), method='cressman', radius=[10,8,6,4])[0]
    
17. 
    
18. #Plot
    
19. proj = projinfo(proj='lcc', lon_0=105, lat_1=25, lat_2=47)
    
20. axesm(projinfo=proj, position=[0, 0, 0.9, 1], axison=False, gridlabel=False, frameon=False)
    
21. geoshow('cn_province', edgecolor='lightgray')
    
22. bou1_layer = geoshow('cn_border', facecolor=(0,0,255))
    
23. geoshow('cn_cities', facecolor='r', size=4, labelfield='NAME', fontname=u'楷体', fontsize=16, yoffset=15)
    
24. china_layer = geoshow('china', visible=False)
    
25. levs = [0.1, 1, 2, 5, 10, 20, 25, 50, 100]
    
26. cols = [(255,255,255),(170,240,255),(120,230,240),(200,220,50),(240,220,20),(255,120,10),(255,90,10), \
    
27.     (240,40,0),(180,10,0),(120,10,0)]
    
28. layer = contourfm(x, y, prg, levs, colors=cols)
    
29. masklayer(china_layer, [layer])
    
30. colorbar(layer, shrink=0.5, aspect=15)
    
31. axism([78, 130, 14, 53])
    
32. text(95, 53, u'全国降水量实况图', fontname=u'黑体', fontsize=18)
    
33. text(95, 51, u'(2010-10-14 08:00 至 2010-10-14 14:00)', fontname=u'黑体', fontsize=16)
    
34. 
    
35. #Add south China Sea
    
36. sc_layer = bou1_layer.clone()
    
37. axesm(position=[0.1,0.05,0.15,0.2], axison=False, frameon=True)
    
38. geoshow(sc_layer, facecolor=(0,0,255))
    
39. xlim(106, 123)
    
40. ylim(2, 23)
    
41. #savefig('D:/Temp/test/rain_test.png', 800, 800)

#Set data folders<br /> <br /> basedir = 'D:/MyProgram/Distribution/java/MeteoInfo/MeteoInfo'<br /> datadir = os.path.join(basedir, 'sample/MICAPS')<br /> <br /> #Read station data<br /> f = addfile_micaps(os.path.join(datadir, '10101414.000'))<br /> lon = f['Longitude'][:]<br /> lat = f['Latitude'][:]<br /> pr = f['Precipitation6h'][:]<br /> <br /> #griddata function - interpolate<br /> x = arange(75, 135, 0.5)<br /> y = arange(18, 55, 0.5)<br /> #prg = griddata((lon, lat), pr, xi=(x, y), method='idw', radius=3)[0]<br /> prg = griddata((lon, lat), pr, xi=(x, y), method='cressman', radius=[10,8,6,4])[0]<br /> <br /> #Plot<br /> proj = projinfo(proj='lcc', lon_0=105, lat_1=25, lat_2=47)<br /> axesm(projinfo=proj, position=[0, 0, 0.9, 1], axison=False, gridlabel=False, frameon=False)<br /> geoshow('cn_province', edgecolor='lightgray')<br /> bou1_layer = geoshow('cn_border', facecolor=(0,0,255))<br /> geoshow('cn_cities', facecolor='r', size=4, labelfield='NAME', fontname=u'楷体', fontsize=16, yoffset=15)<br /> china_layer = geoshow('china', visible=False)<br /> levs = [0.1, 1, 2, 5, 10, 20, 25, 50, 100]<br /> cols = [(255,255,255),(170,240,255),(120,230,240),(200,220,50),(240,220,20),(255,120,10),(255,90,10), \<br />     (240,40,0),(180,10,0),(120,10,0)]<br /> layer = contourfm(x, y, prg, levs, colors=cols)<br /> masklayer(china_layer, [layer])<br /> colorbar(layer, shrink=0.5, aspect=15)<br /> axism([78, 130, 14, 53])<br /> text(95, 53, u'全国降水量实况图', fontname=u'黑体', fontsize=18)<br /> text(95, 51, u'(2010-10-14 08:00 至 2010-10-14 14:00)', fontname=u'黑体', fontsize=16)<br /> <br /> #Add south China Sea<br /> sc_layer = bou1_layer.clone()<br /> axesm(position=[0.1,0.05,0.15,0.2], axison=False, frameon=True)<br /> geoshow(sc_layer, facecolor=(0,0,255))<br /> xlim(106, 123)<br /> ylim(2, 23)<br /> #savefig('D:/Temp/test/rain_test.png', 800, 800)

运行结果：

MeteoInfo

气象万千2011 发表于 2015-7-10 16:01
王老师，我把您的这个例子运行时出现了这个错误，请问怎么解决
File "", line None
SyntaxError: Non-ASC ...

你在脚本程序最前面加上一句：

# coding=utf-8

试试

MeteoInfo

气象万千2011 发表于 2015-7-11 17:24
谢谢老师，可以了。我想问下脚本都要遭MeteoInfoLab里运行吗？因为我先前把有些脚本在工具--脚本编辑器里 ...

MeteoInfoLab是用Jython对MeteoInfo库的功能封装了很多类似NCL和MatLab函数，编写脚本更加简洁、方便。可以通过下面几点简单判断是否需要在MeteoInfoLab中运行。

1、在MeteoInfo中写脚本需要引入大量类（脚本开始的地方有很多 import语句），而MeteoInoLab的脚本通常不需要引入MeteoInfo类的语句（实际上在封装时已经提前引用过了）。

2、MeteoInfoLab中打开数据文件的命令通常为addfile。

其实绝大部分之前在MeteoInfo里运行的脚本都可以在MeteoInfoLab中运行，除非脚本中有miapp变量，该变量就是正在运行的MeteoInfo主界面。

传说中的谁

MeteoInfo功能真是越来越强大了，发展方向也多元化，希望有一天能成为国际主流的绘图软件

小船

谢谢楼主，辛苦了，保存了

小船

好东西，收藏，谢谢楼主

气象万千2011

王老师，我把您的这个例子运行时出现了这个错误，请问怎么解决
File "<string>", line None
SyntaxError: Non-ASCII character in file '<iostream>', but no encoding declared; see http://www.python.org/peps/pep-0263.html for details

气象万千2011

MeteoInfo 发表于 2015-7-10 17:03
你在脚本程序最前面加上一句：

# coding=utf-8

# coding=utf-8
#Set data folders
basedir = 'E:/MeteoInfo'
datadir = os.path.join(basedir, 'sample/MICAPS')
mapdir = os.path.join(basedir, 'map')
#Read shape files
bou2_layer = shaperead(os.path.join(mapdir, 'bou2_4p.shp'))
bou1_layer = shaperead(os.path.join(mapdir, 'bou1_4l.shp'))
china_layer = shaperead(os.path.join(mapdir, 'china.shp'))
city_layer = shaperead(os.path.join(mapdir, 'res1_4m.shp'))
#Read station data
f = addfile_micaps(os.path.join(datadir, '10101414.000'))
pr = f.stationdata('Precipitation6h')
#griddata function - interpolate
x = arange(75, 135, 0.5)
y = arange(18, 55, 0.5)
prg = pr.griddata((x, y), method='idw', radius=3)
#Plot
axesm()
geoshow(bou2_layer, edgecolor='lightgray')
geoshow(bou1_layer, facecolor=(0,0,255))
#geoshow(city_layer, facecolor='r', size=4, labelfield='NAME', fontname=u'楷体', fontsize=16, yoffset=15)
geoshow(china_layer, visible=False)
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)]
layer = contourfm(prg, levs, colors=cols)
masklayer(china_layer, [layer])
colorbar(layer)
xlim(72, 136)
ylim(16, 55)
text(95, 52, u'全国降水量实况图', fontname=u'黑体', fontsize=16)
text(95, 50, u'(2010-10-14 08:00 至 2010-10-14 14:00)', fontname=u'黑体', fontsize=14)
#Add south China Sea
sc_layer = bou1_layer.clone()
axesm(position=[0.14,0.18,0.15,0.2], axison=False)
geoshow(sc_layer, facecolor=(0,0,255))
xlim(106, 123)
ylim(2, 23)


王老师，以上是我的程序，还是会报错，如下：
Traceback (most recent call last):
  File "<iostream>", line 4, in <module>
NameError: name 'os' is not defined

MeteoInfo

气象万千2011 发表于 2015-7-10 17:30
# coding=utf-8
#Set data folders
basedir = 'E:/MeteoInfo'

你是在什么环境下运行的？这个需要在MeteoInfoLab里运行。

气象万千2011

MeteoInfo 发表于 2015-7-10 17:45
你是在什么环境下运行的？这个需要在MeteoInfoLab里运行。

谢谢老师，可以了。我想问下脚本都要遭MeteoInfoLab里运行吗？因为我先前把有些脚本在工具--脚本编辑器里面运行也可以，直到出现这个脚本有很多问题