爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
楼主: MeteoInfo

MeteoInfoLab脚本示例:闪电位置图

[复制链接]

新浪微博达人勋

发表于 2016-10-18 18:41:30 | 显示全部楼层
本帖最后由 wuwei2163 于 2016-10-18 21:04 编辑
MeteoInfo 发表于 2016-10-18 17:25
读取每个闪电的经纬度,判断是否在区域内,就这么简单

嗯,若是个不规则的区域,我想的是先定义个对应的shape图层文件,然后判断区域用inpollygon函数进行判断,然后在统计分析。但是脚本有问题,分析其原因:一是我这样的思路根本就不对;二是对inpollygon函数不熟悉,进行判断时出错了。本人这个菜鸟尝试了,编写的脚本如下,请王老师和各位老师不要见笑,多多指导:
#Read lighting data
fn = 'e:/2016_07_09.txt'
tf = open(fn)
times = []
lats = []
lons = []
nums = []
for aline in tf:   
    datalist = aline.split()
    lat = float(datalist[3].split('=')[1])
    lon = float(datalist[4].split('=')[1])
    lats.append(lat)
    lons.append(lon)
    nums.append(0)
lon = array(lons)
lat = array(lats)  
for h in range(0, 24):
    times.append(datetime.datetime(2016,7,19,h))
    h = int(datalist[2].split(':')[0])
#Read shape file
m_henan = shaperead('D:/slh.shp')
#inpolygon function
isin = inpolygon(lon,lat, m_henan)
print isin
nums[h] = nums[h] + 1
#Plot
plot(times, nums, '-ro')
xlim(times[0], times[-1])

nmmm.PNG 麻烦王老师您看看,帮着分析分析,判断是否在区域内的脚本该怎样写呢?
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2016-10-18 21:08:41 | 显示全部楼层
wuwei2163 发表于 2016-10-18 18:41
嗯,若是个不规则的区域,我想的是先定义个对应的shape图层文件,然后判断区域用inpollygon函数进行判断 ...

两种方式都可以
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2016-10-18 21:55:27 | 显示全部楼层
本帖最后由 wuwei2163 于 2016-10-18 22:36 编辑

我选择第一种方式:先定义个对应的shape图层文件,然后判断区域用inpollygon函数进行判断,然后在统计分析。但是脚本有问题,分析其原因:一是这种思路根本就不对;二是对inpollygon函数不熟悉,进行判断时出错了。我编写的脚本如下,请您看看:
#Read lighting data
fn = 'e:/2016_07_09.txt'
tf = open(fn)
times = []
lats = []
lons = []
nums = []
for aline in tf:   
     datalist = aline.split()
     lat = float(datalist[3].split('=')[1])
     lon = float(datalist[4].split('=')[1])
     lats.append(lat)
     lons.append(lon)
     nums.append(0)
lon = array(lons)
lat = array(lats)  
for h in range(0, 24):
     times.append(datetime.datetime(2016,7,19,h))
     h = int(datalist[2].split(':')[0])
#Read shape file
m_henan = shaperead('D:/slh.shp')
#inpolygon function
isin = inpolygon(lon,lat, m_henan)
print isin
nums[h] = nums[h] + 1
#Plot
plot(times, nums, '-ro')
xlim(times[0], times[-1])










麻烦王老师您看看,关于错误在哪儿呢?求正确的判断是否在区域内的代码?
nmmm.PNG
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2016-10-18 22:35:19 | 显示全部楼层
wuwei2163 发表于 2016-10-18 21:55
先定义个对应的shape图层文件,然后判断区域用inpollygon函数进行判断,然后在统计分析。但是脚本有问题 ...

参考此脚本:
  1. fn = 'D:/Temp/ascii/lighting/2009_06_06.txt'
  2. tf = open(fn)
  3. m_henan = shaperead('D:/Temp/Map/henan.shp')
  4. lats = []
  5. lons = []
  6. vs = []
  7. for aline in tf:   
  8.     datalist = aline.split()
  9.     lat = float(datalist[3].split('=')[1])
  10.     lon = float(datalist[4].split('=')[1])
  11.     isin = inpolygon(lon, lat, m_henan.shapes()[0])
  12.     if isin:
  13.         v = float(datalist[5].split('=')[1])
  14.         lats.append(lat)
  15.         lons.append(lon)
  16.         vs.append(v)
  17. lon = array(lons)
  18. lat = array(lats)  
  19. v = array(vs)
  20. #Plot
  21. axesm()
  22. mlayer = shaperead('D:/Temp/map/bou2_4p.shp')
  23. geoshow(mlayer)
  24. ss = makesymbolspec('point', {'value':(-10000,0), 'color':'b', 'marker':'m', 'size':6, 'caption':'Negative'}, \
  25.     {'value':(0,10000), 'color':'r', 'marker':'+', 'size':6, 'caption':'Positive'})
  26. layer = scatterm(lon, lat, v, symbolspec=ss)
  27. legend(legend=layer.legend(), loc='lower left')
  28. xlim(90, 130)
  29. ylim(20, 50)
  30. title('Lighting locations')


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

新浪微博达人勋

发表于 2016-10-18 22:43:16 | 显示全部楼层

谢谢王老师的指导,我一而再再而三的麻烦您,您还很耐心的帮助我这个新手,我再好好学习和摸索摸索~
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2016-10-18 22:45:49 | 显示全部楼层
wuwei2163 发表于 2016-10-18 22:43
谢谢王老师的指导,我一而再再而三的麻烦您,您还很耐心的帮助我这个新手,我再好好学习和摸索摸索~

建议仔细看看User's guide:http://www.meteothinker.com/docs/meteoinfolab/user_guide.html
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2016-10-18 22:58:04 | 显示全部楼层
MeteoInfo 发表于 2016-10-18 22:45
建议仔细看看User's guide:http://www.meteothinker.com/docs/meteoinfolab/user_guide.html

好的,王老师,这个教程真是宝典啊,等我把省台布置的活忙完后,一定会从零基础开始一点一点的学习这本教程宝典的~~~
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2016-10-19 09:32:12 | 显示全部楼层

王老师,根据您的指导,我脚本重新编写之后运行,同一个脚本出现两种情况:一是脚本没有错误,但不出图;二是程序出错,出来的图也不正常。 1020.PNG 1019.PNG
具体的脚本如下:
fn = 'F:/shandianshujuchuli/2016_07_09.txt'
tf = open(fn)
m_henan = shaperead('D:/MeteoInfo1.36/map/henan.shp')
lats = []
lons = []
times = []
nums = []
for aline in tf:   
    datalist = aline.split()
    lat = float(datalist[3].split('=')[1])
    lon = float(datalist[4].split('=')[1])  
    isin = inpolygon(lon, lat, m_henan.shapes()[0])
    if isin:
        lats.append(lat)
        lons.append(lon)
for h in range(0, 24):
    times.append(datetime.datetime(2016,7,9,h))
    nums.append(0)
for aline in tf:   
    datalist = aline.split()
    h = int(datalist[2].split(':')[0])
    nums[h] = nums[h] + 1
#Plot
plot(times, nums, '-ro')
xaxis(axistype='time')
ylabel('Lighting number')
xlim(times[0], times[-1])
title('Lighting number hourly variation in ' + times[0].strftime('%Y/%m/%d')),估计还是脚本有问题吧,但我一直没有查出来,请您指导,学生不胜感激~
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2016-10-19 10:47:25 | 显示全部楼层
wuwei2163 发表于 2016-10-19 09:32
王老师,根据您的指导,我脚本重新编写之后运行,同一个脚本出现两种情况:一是脚本没有错误,但不出图; ...

你是要做河南省某日的逐时闪电次数变化图吗?
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2016-10-19 10:53:04 | 显示全部楼层
MeteoInfo 发表于 2016-10-19 10:47
你是要做河南省某日的逐时闪电次数变化图吗?

嗯,是的,王老师。我上面的脚本就是先对河南区域内的闪电判断,然后再统计变化。但运行程序后不出图或者报错,不知道什么原因?
密码修改失败请联系微信:mofangbao
您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

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

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

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