请选择 进入手机版 | 继续访问电脑版
爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 44846|回复: 99

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

  [复制链接]

新浪微博达人勋

发表于 2015-6-30 10:38:37 | 显示全部楼层 |阅读模式

登录后查看更多精彩内容~

您需要 登录 才可以下载或查看,没有帐号?立即注册 新浪微博登陆

x
本帖最后由 MeteoInfo 于 2022-11-24 11:18 编辑

站点数据绘制等值线需要首先将站点数据插值为格点数据,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. #griddata function - interpolate
  8. x = arange(75, 135, 0.5)
  9. y = arange(18, 55, 0.5)
  10. prg, gx, gy = griddata((lon, lat), pr, xi=(x, y), method='cressman', radius=[7,4,2,1])

  11. #Plot
  12. figure(figsize=[700, 400], newfig=False)
  13. axesm()
  14. geoshow('cn_province', edgecolor='lightgray')
  15. bou1_layer = geoshow('cn_border', facecolor=(0,0,255))
  16. city_layer = geoshow('cn_cities', facecolor='r', size=4)
  17. city_layer.addlabels('NAME', fontname=u'楷体', fontsize=16, yoffset=15, avoidcoll=False)
  18. china_layer = geoshow('china', visible=False)
  19. city_layer.movelabel(u'西宁', -15)
  20. levs = [0.1, 1, 2, 5, 10, 20, 25, 50, 100]
  21. cols = [(255,255,255),(170,240,255),(120,230,240),(200,220,50),(240,220,20),(255,120,10),(255,90,10), \
  22.     (240,40,0),(180,10,0),(120,10,0)]
  23. layer = contourf(x, y, prg, levs, colors=cols)
  24. masklayer(china_layer, [layer])
  25. colorbar(layer)
  26. xlim(72, 136)
  27. ylim(16, 55)
  28. text(95, 52, u'全国降水量实况图', fontname=u'黑体', fontsize=16)
  29. text(95, 50, u'(2010-10-14 08:00 至 2010-10-14 14:00)', fontname=u'黑体', fontsize=14)

  30. #Add south China Sea
  31. sc_layer = bou1_layer.clone()
  32. axesm(position=[0.14,0.18,0.15,0.2], axison=False, frameon=True)
  33. geoshow(sc_layer, facecolor=(0,0,255))
  34. xlim(106, 123)
  35. ylim(2, 23)
  36. #savefig('D:/Temp/test/test_pre.png')


脚本程序(Lambert投影):
  1. #Set data folders

  2. basedir = 'D:/MyProgram/Distribution/java/MeteoInfo/MeteoInfo'
  3. datadir = os.path.join(basedir, 'sample/MICAPS')

  4. #Read station data
  5. f = addfile_micaps(os.path.join(datadir, '10101414.000'))
  6. lon = f['Longitude'][:]
  7. lat = f['Latitude'][:]
  8. pr = f['Precipitation6h'][:]

  9. #griddata function - interpolate
  10. x = arange(75, 135, 0.5)
  11. y = arange(18, 55, 0.5)
  12. #prg = griddata((lon, lat), pr, xi=(x, y), method='idw', radius=3)[0]
  13. prg = griddata((lon, lat), pr, xi=(x, y), method='cressman', radius=[10,8,6,4])[0]

  14. #Plot
  15. proj = projinfo(proj='lcc', lon_0=105, lat_1=25, lat_2=47)
  16. axesm(projinfo=proj, position=[0, 0, 0.9, 1], axison=False, gridlabel=False, frameon=False)
  17. geoshow('cn_province', edgecolor='lightgray')
  18. bou1_layer = geoshow('cn_border', facecolor=(0,0,255))
  19. geoshow('cn_cities', facecolor='r', size=4, labelfield='NAME', fontname=u'楷体', fontsize=16, yoffset=15)
  20. china_layer = geoshow('china', visible=False)
  21. levs = [0.1, 1, 2, 5, 10, 20, 25, 50, 100]
  22. cols = [(255,255,255),(170,240,255),(120,230,240),(200,220,50),(240,220,20),(255,120,10),(255,90,10), \
  23.     (240,40,0),(180,10,0),(120,10,0)]
  24. layer = contourfm(x, y, prg, levs, colors=cols)
  25. masklayer(china_layer, [layer])
  26. colorbar(layer, shrink=0.5, aspect=15)
  27. axism([78, 130, 14, 53])
  28. text(95, 53, u'全国降水量实况图', fontname=u'黑体', fontsize=18)
  29. text(95, 51, u'(2010-10-14 08:00 至 2010-10-14 14:00)', fontname=u'黑体', fontsize=16)

  30. #Add south China Sea
  31. sc_layer = bou1_layer.clone()
  32. axesm(position=[0.1,0.05,0.15,0.2], axison=False, frameon=True)
  33. geoshow(sc_layer, facecolor=(0,0,255))
  34. xlim(106, 123)
  35. ylim(2, 23)
  36. #savefig('D:/Temp/test/rain_test.png', 800, 800)


运行结果:
Image00861.png

Image00865.png

密码修改失败请联系微信:mofangbao

新浪微博达人勋

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

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

# coding=utf-8

试试
密码修改失败请联系微信:mofangbao
回复 支持 1 反对 0

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2015-7-11 18:12:42 | 显示全部楼层
气象万千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主界面。
密码修改失败请联系微信:mofangbao
回复 支持 1 反对 0

使用道具 举报

新浪微博达人勋

发表于 2015-6-30 11:52:43 | 显示全部楼层
MeteoInfo功能真是越来越强大了,发展方向也多元化,希望有一天能成为国际主流的绘图软件
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2015-6-30 14:34:26 | 显示全部楼层
谢谢楼主,辛苦了,保存了
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2015-7-1 18:21:24 | 显示全部楼层
好东西,收藏,谢谢楼主

评分

参与人数 1金钱 +1 收起 理由
凉生初雨 + 1 步步生莲 www.7455.com/play/117862.html

查看全部评分

密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2015-7-10 16:01:59 | 显示全部楼层
王老师,我把您的这个例子运行时出现了这个错误,请问怎么解决
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
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2015-7-10 17:30:37 | 显示全部楼层
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


密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2015-7-10 17:45:00 | 显示全部楼层
气象万千2011 发表于 2015-7-10 17:30
# coding=utf-8
#Set data folders
basedir = 'E:/MeteoInfo'

你是在什么环境下运行的?这个需要在MeteoInfoLab里运行。
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2015-7-11 17:24:55 | 显示全部楼层
MeteoInfo 发表于 2015-7-10 17:45
你是在什么环境下运行的?这个需要在MeteoInfoLab里运行。

谢谢老师,可以了。我想问下脚本都要遭MeteoInfoLab里运行吗?因为我先前把有些脚本在工具--脚本编辑器里面运行也可以,直到出现这个脚本有很多问题
密码修改失败请联系微信:mofangbao
您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表,不代表气象家园立场,禁止在本站发表与国家法律相抵触言论

快速回复 返回顶部 返回列表