MeteoInfoLab脚本示例：闪电位置图

wuwei2163

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])

麻烦王老师您看看，帮着分析分析，判断是否在区域内的脚本该怎样写呢？

MeteoInfo

wuwei2163 发表于 2016-10-18 18:41
嗯，若是个不规则的区域，我想的是先定义个对应的shape图层文件，然后判断区域用inpollygon函数进行判断 ...

两种方式都可以

wuwei2163

MeteoInfo 发表于 2016-10-18 21:08
两种方式都可以

我选择第一种方式：先定义个对应的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])










麻烦王老师您看看，关于错误在哪儿呢？求正确的判断是否在区域内的代码？

MeteoInfo

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')

fn = 'D:/Temp/ascii/lighting/2009_06_06.txt'<br /> tf = open(fn)<br /> m_henan = shaperead('D:/Temp/Map/henan.shp')<br /> lats = []<br /> lons = []<br /> vs = []<br /> for aline in tf:    <br />     datalist = aline.split()<br />     lat = float(datalist[3].split('=')[1])<br />     lon = float(datalist[4].split('=')[1])<br />     isin = inpolygon(lon, lat, m_henan.shapes()[0])<br />     if isin:<br />         v = float(datalist[5].split('=')[1])<br />         lats.append(lat)<br />         lons.append(lon)<br />         vs.append(v)<br /> lon = array(lons)<br /> lat = array(lats)  <br /> v = array(vs)<br /> #Plot<br /> axesm()<br /> mlayer = shaperead('D:/Temp/map/bou2_4p.shp')<br /> geoshow(mlayer)<br /> ss = makesymbolspec('point', {'value':(-10000,0), 'color':'b', 'marker':'m', 'size':6, 'caption':'Negative'}, \<br />     {'value':(0,10000), 'color':'r', 'marker':'+', 'size':6, 'caption':'Positive'})<br /> layer = scatterm(lon, lat, v, symbolspec=ss)<br /> legend(legend=layer.legend(), loc='lower left')<br /> xlim(90, 130)<br /> ylim(20, 50) <br /> title('Lighting locations')

wuwei2163

MeteoInfo 发表于 2016-10-18 22:35
参考此脚本：

谢谢王老师的指导，我一而再再而三的麻烦您，您还很耐心的帮助我这个新手，我再好好学习和摸索摸索~

MeteoInfo

wuwei2163 发表于 2016-10-18 22:43
谢谢王老师的指导，我一而再再而三的麻烦您，您还很耐心的帮助我这个新手，我再好好学习和摸索摸索~

建议仔细看看User's guide：http://www.meteothinker.com/docs/meteoinfolab/user_guide.html

wuwei2163

MeteoInfo 发表于 2016-10-18 22:45
建议仔细看看User's guide：http://www.meteothinker.com/docs/meteoinfolab/user_guide.html

好的，王老师，这个教程真是宝典啊，等我把省台布置的活忙完后，一定会从零基础开始一点一点的学习这本教程宝典的~~~

wuwei2163

MeteoInfo 发表于 2016-10-18 22:35
参考此脚本：

王老师，根据您的指导，我脚本重新编写之后运行，同一个脚本出现两种情况：一是脚本没有错误，但不出图；二是程序出错，出来的图也不正常。
具体的脚本如下：
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'))，估计还是脚本有问题吧，但我一直没有查出来，请您指导，学生不胜感激~

MeteoInfo

wuwei2163 发表于 2016-10-19 09:32
王老师，根据您的指导，我脚本重新编写之后运行，同一个脚本出现两种情况：一是脚本没有错误，但不出图； ...

你是要做河南省某日的逐时闪电次数变化图吗？

wuwei2163

MeteoInfo 发表于 2016-10-19 10:47
你是要做河南省某日的逐时闪电次数变化图吗？

嗯，是的，王老师。我上面的脚本就是先对河南区域内的闪电判断，然后再统计变化。但运行程序后不出图或者报错，不知道什么原因？