MeteoInfoLab脚本示例：闪电位置图

MeteoInfo

这个脚本示例读取文本格式的闪电数据，读出每条闪电记录的经纬度和强度，在地图上绘制出每个闪电的位置，并用符号和颜色区分强度正负。

数据格式如下：
0 2009-06-06 00:01:16.6195722 纬度=35.0297  经度=114.8923  强度=23.4   陡度=9.1    误差=71.9   定位方式：四站算法
1 2009-06-06 01:34:58.5911760 纬度=34.4243  经度=115.6057  强度=-1.3   陡度=-0.4   误差=0.0    定位方式：二站振幅
2 2009-06-06 01:36:59.6324818 纬度=35.9314  经度=114.9388  强度=-35.2  陡度=-6.1   误差=0.0    定位方式：三站混合
3 2009-06-06 01:54:25.2727088 纬度=34.4163  经度=115.6121  强度=-1.9   陡度=-0.3   误差=0.0    定位方式：二站振幅
4 2009-06-06 01:58:35.1179896 纬度=35.3270  经度=106.5536  强度=-384.1 陡度=-80.6  误差=0.0    定位方式：二站振幅
5 2009-06-06 01:59:15.9764295 纬度=34.5339  经度=114.4466  强度=58.4   陡度=26.7   误差=0.0    定位方式：三站混合
7 2009-06-06 02:00:35.0869527 纬度=34.3934  经度=115.5976  强度=-3.3   陡度=-1.0   误差=0.0    定位方式：二站振幅
8 2009-06-06 02:00:43.0728269 纬度=34.4124  经度=115.7238  强度=-5.3   陡度=-3.1   误差=0.0    定位方式：二站混合
9 2009-06-06 02:00:43.2227403 纬度=34.4131  经度=115.7171  强度=-8.4   陡度=-4.1   误差=0.0    定位方式：二站混合
10 2009-06-06 02:02:21.5951061 纬度=34.3817  经度=115.6268  强度=-7.2   陡度=-1.6   误差=0.0    定位方式：二站振幅


脚本程序：

1. fn = 'D:/Temp/ascii/lighting/2009_06_06.txt'
   
2. tf = open(fn)
   
3. lats = []
   
4. lons = []
   
5. vs = []
   
6. for aline in tf:   
   
7.     datalist = aline.split()
   
8.     lat = float(datalist[3].split('=')[1])
   
9.     lon = float(datalist[4].split('=')[1])
   
10.     v = float(datalist[5].split('=')[1])
    
11.     lats.append(lat)
    
12.     lons.append(lon)
    
13.     vs.append(v)
    
14. lon = array(lons)
    
15. lat = array(lats)  
    
16. v = array(vs)
    
17. axesm()
    
18. mlayer = shaperead('D:/Temp/map/bou2_4p.shp')
    
19. geoshow(mlayer)
    
20. ss = makesymbolspec('point', {'value':(-10000,0), 'color':'b', 'marker':'m', 'size':6}, \
    
21.     {'value':(0,10000), 'color':'r', 'marker':'+', 'size':6})
    
22. layer = scatterm(lon, lat, v, symbolspec=ss)
    
23. xlim(90, 130)
    
24. ylim(20, 50)
    
25. title(u'闪电位置图', fontname=u'黑体', fontsize=18)

fn = 'D:/Temp/ascii/lighting/2009_06_06.txt'<br /> tf = open(fn)<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 />     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 /> 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}, \<br />     {'value':(0,10000), 'color':'r', 'marker':'+', 'size':6})<br /> layer = scatterm(lon, lat, v, symbolspec=ss)<br /> xlim(90, 130)<br /> ylim(20, 50) <br /> title(u'闪电位置图', fontname=u'黑体', fontsize=18)

计算1*1度格点闪电次数并绘图：

1. #Read lighting data
   
2. fn = 'D:/Temp/ascii/lighting/2009_06_06.txt'
   
3. tf = open(fn)
   
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.     v = float(datalist[5].split('=')[1])
    
12.     lats.append(lat)
    
13.     lons.append(lon)
    
14.     vs.append(v)
    
15. lon = array(lons)
    
16. lat = array(lats)  
    
17. v = array(vs)
    
18. 
    
19. #Get grid lighting number
    
20. glon = arange(90, 131, 1)
    
21. glat = arange(20, 51, 1)
    
22. gv = griddata((lon, lat), v, xi=(glon, glat), method='inside_count')[0]
    
23. 
    
24. #Plot
    
25. axesm()
    
26. mlayer = shaperead('D:/Temp/map/bou2_4p.shp')
    
27. geoshow(mlayer)
    
28. levs = [1,2,5,10,20,50,100,200,300,400,500,600,700,800,900,1000]
    
29. #layer = contourfm(glon, glat, gv, levs, cmap='WhBlGrYeRe')
    
30. layer = imshowm(glon, glat, gv, levs, cmap='WhBlGrYeRe')
    
31. colorbar(layer)
    
32. title('Lighting number')

#Read lighting data<br /> fn = 'D:/Temp/ascii/lighting/2009_06_06.txt'<br /> tf = open(fn)<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 />     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 /> <br /> #Get grid lighting number<br /> glon = arange(90, 131, 1)<br /> glat = arange(20, 51, 1)<br /> gv = griddata((lon, lat), v, xi=(glon, glat), method='inside_count')[0]<br /> <br /> #Plot<br /> axesm()<br /> mlayer = shaperead('D:/Temp/map/bou2_4p.shp')<br /> geoshow(mlayer)<br /> levs = [1,2,5,10,20,50,100,200,300,400,500,600,700,800,900,1000]<br /> #layer = contourfm(glon, glat, gv, levs, cmap='WhBlGrYeRe')<br /> layer = imshowm(glon, glat, gv, levs, cmap='WhBlGrYeRe')<br /> colorbar(layer)<br /> title('Lighting number')

统计一日内每小时闪电次数并绘制逐时变化图：

1. #Read lighting data
   
2. fn = 'D:/Temp/ascii/lighting/2009_06_06.txt'
   
3. tf = open(fn)
   
4. times = []
   
5. nums = []
   
6. for h in range(0, 24):
   
7.     times.append(datetime.datetime(2009,6,6,h))
   
8.     nums.append(0)
   
9. for aline in tf:   
   
10.     datalist = aline.split()
    
11.     h = int(datalist[2].split(':')[0])
    
12.     nums[h] = nums[h] + 1
    
13. 
    
14. #Plot
    
15. plot(times, nums, '-ro')
    
16. xaxis(axistype='time')
    
17. ylabel('Lighting number')
    
18. xlim(times[0], times[-1])
    
19. title('Lighting number hourly variation in ' + times[0].strftime('%Y/%m/%d'))

#Read lighting data<br /> fn = 'D:/Temp/ascii/lighting/2009_06_06.txt'<br /> tf = open(fn)<br /> times = []<br /> nums = []<br /> for h in range(0, 24):<br />     times.append(datetime.datetime(2009,6,6,h))<br />     nums.append(0)<br /> for aline in tf:    <br />     datalist = aline.split()<br />     h = int(datalist[2].split(':')[0])<br />     nums[h] = nums[h] + 1<br /> <br /> #Plot<br /> plot(times, nums, '-ro')<br /> xaxis(axistype='time')<br /> ylabel('Lighting number')<br /> xlim(times[0], times[-1])<br /> title('Lighting number hourly variation in ' + times[0].strftime('%Y/%m/%d'))

MeteoInfo

lu19320 发表于 2017-4-7 16:21
王老师，你好！由于编程基础比较弱，你编写的脚本我只能勉强看懂，自己动手写脚本还比较难。现在部分省市已 ...

1. fn = 'D:/Temp/ascii/lighting/1.csv'
   
2. table = readtable(fn, delimiter=',', format='%i%{yyyy/M/d H:mm}D%s%4f%3s%4f', \
   
3.     encoding='gb2312')
   
4. lon = table[u'经度']
   
5. lat = table[u'纬度']
   
6. v = table[u'强度(千安)']
   
7. ltype = table[u'类型']
   
8. #Plot
   
9. axesm()
   
10. mlayer = shaperead('D:/Temp/map/bou2_4p.shp')
    
11. geoshow(mlayer)
    
12. ss = makesymbolspec('point', {'value':(-10000,0), 'color':'b', 'marker':'m', 'size':6, 'caption':'Negative'}, \
    
13.     {'value':(0,10000), 'color':'r', 'marker':'+', 'size':6, 'caption':'Positive'})
    
14. layer = scatterm(lon, lat, v, symbolspec=ss)
    
15. legend(legend=layer.legend(), loc='lower left')
    
16. xlim(100, 120)
    
17. ylim(20, 30)
    
18. title('Lighting locations')

fn = 'D:/Temp/ascii/lighting/1.csv'<br /> table = readtable(fn, delimiter=',', format='%i%{yyyy/M/d H:mm}D%s%4f%3s%4f', \<br />     encoding='gb2312')<br /> lon = table[u'经度']<br /> lat = table[u'纬度']<br /> v = table[u'强度(千安)']<br /> ltype = table[u'类型']<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(100, 120)<br /> ylim(20, 30) <br /> title('Lighting locations')

好多鱼

很好的帖子，留下来了。

wuwei2163

王老师，您好！这个帖子真心点赞，不过有个问题想请您帮忙，也算是对这个帖子的一个建议:能否统计显示正闪、负闪频次数，这样更有利于业务上使用；从代码上来看：ss = makesymbolspec('point', {'value':(-10000,0), 'color':'b', 'marker':'m', 'size':6}, \{'value':(0,10000), 'color':'r', 'marker':'+', 'size':6})
是画闪电的命令，您看能不能补充一下，把闪频次数统计出来显示在图例上呢？请老师指导指导

MeteoInfo

wuwei2163 发表于 2016-10-13 09:50
王老师，您好！这个帖子真心点赞，不过有个问题想请您帮忙，也算是对这个帖子的一个建议:能否统计显示正闪 ...

你想要统计格点的闪电次数吗？比如1*1度格点？

wuwei2163

MeteoInfo 发表于 2016-10-13 10:53
你想要统计格点的闪电次数吗？比如1*1度格点？

嗯，是的，王老师。比如某一区域内的闪电次数？在这个代码的基础上，以图例的方式显示统计次数，不知道简单不简单？

MeteoInfo

wuwei2163 发表于 2016-10-13 10:57
嗯，是的，王老师。比如某一区域内的闪电次数？在这个代码的基础上，以图例的方式显示统计次数，不知道简 ...

已在一楼增加了相应的脚本。

wuwei2163

MeteoInfo 发表于 2016-10-13 12:04
已在一楼增加了相应的脚本。

王老师，您效率真高，{:eb502:}{:eb502:}这两个脚本放一起或者结合起来，我个人感觉更贴近业务需求了。王老师，以后要多多向您请教~~~

wuwei2163

MeteoInfo 发表于 2016-10-13 12:04
已在一楼增加了相应的脚本。

王老师，您好！按照你的角码，我检验了下数据，发现有错误，请您看一下：什么原因造成的？

看了看，似乎是

                               
登录/注册后可看大图

gv = griddata((lon, lat), v, xi=(glon, glat), method='inside_count')[0]这个有误？

MeteoInfo

wuwei2163 发表于 2016-10-13 16:49
王老师，您好！按照你的角码，我检验了下数据，发现有错误，请您看一下：什么原因造成的？

看了看，似 ...

你下载MeteoInfo最新版本试试。

wuwei2163

MeteoInfo 发表于 2016-10-13 19:32
你下载MeteoInfo最新版本试试。

嗯，王老师，下个新版本问题就解决了，