MeteoInfoLab脚本示例：闪电位置图

wuwei2163 · 发表于 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])

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

MeteoInfo · 发表于 2016-10-18 21:08:41

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

两种方式都可以

wuwei2163 · 发表于 2016-10-18 21:55:27

本帖最后由 wuwei2163 于 2016-10-18 22:36 编辑

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 · 发表于 2016-10-18 22:35:19

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

参考此脚本：

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

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

wuwei2163 · 发表于 2016-10-18 22:43:16

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

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

MeteoInfo · 发表于 2016-10-18 22:45:49

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

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

wuwei2163 · 发表于 2016-10-18 22:58:04

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

好的，王老师，这个教程真是宝典啊

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

wuwei2163 · 发表于 2016-10-19 09:32:12

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 · 发表于 2016-10-19 10:47:25

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

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

wuwei2163 · 发表于 2016-10-19 10:53:04

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

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

		自动登录	找回密码
密码			立即注册