绘制闪电密度图

MeteoInfo

已有所有闪电发生的位置来计算和绘制闪电密度图。数据见此贴：http://bbs.06climate.com/forum.php?mod=viewthread&tid=14099&extra=&page=1

用脚本程序来实现，首先定义一个格点数据（GridData），主要参数是格点的X/Y方向的起始位置、格点间距、格点数。然后读取闪电位置数据文件，计算出每个闪电位于X/Y方向上第几个格点，将该格点的值加1，遍历所有闪电后就可以算出每个格点的闪电数目。可以将此格点数据存为一个文本格点数据文件（比如Surfer ASCII 格点数据），格点数据可以依据图例（LegendScheme）设置来生成等值线填充图层（Shaded）。


可以验证一下上面计算出的闪电密度数据是否正确。先用脚本程序生成一个闪电位置点图层：


再用上面保存的格点数据文件生成一个Grid_Fill图层：


将该图层每个图元进行标注（也就是每个格点中闪电发生的次数），图形放大后可以数数一些格点内的闪电个数和标注数字是否相同（相同当然说明计算正确了）


脚本程序如下：

1. #--------------------------------------------------------        
   
2. # Author: Yaqiang Wang                                          
   
3. # Date: 2013-5-5                                                
   
4. # Purpose: Read and plot lighting frequency data  
   
5. # Note: Sample                                                   
   
6. #-----------------------------------------------------------  
   
7. #---- Import clr and classes
   
8. import clr
   
9. clr.AddReferenceByPartialName("System")
   
10. clr.AddReferenceByPartialName("System.Drawing")
    
11. from System import *
    
12. clr.AddReference("MeteoInfoC.dll")
    
13. from MeteoInfoC import *
    
14. from MeteoInfoC.Data import *
    
15. from MeteoInfoC.Data.MeteoData import *
    
16. from MeteoInfoC.Legend import *
    
17. from MeteoInfoC.Shape import *
    
18. import os.path
    
19. 
    
20. #---- Set directory
    
21. dataDir = "D:\\Temp\\ascii\\"
    
22. 
    
23. #---- Define grid data
    
24. print 'Define grid data...'
    
25. x0 = 105
    
26. y0 = 28
    
27. xdelt = 0.1
    
28. ydelt = 0.1
    
29. xnum = 55
    
30. ynum = 45
    
31. gData = GridData(x0,xdelt,xnum,y0,ydelt,ynum)
    
32. 
    
33. #---- Read lighting data file
    
34. print 'Read lighting data file...'
    
35. fn = dataDir + "cunshuju.txt"
    
36. tf = open(fn)
    
37. for aline in tf:       
    
38.         datalist = aline.split()
    
39.         lat = float(datalist[1])
    
40.         lon = float(datalist[0])       
    
41.         xIdx = int((lon - x0 + xdelt / 2) / xdelt)
    
42.         if xIdx < 0 or xIdx > xnum -1:
    
43.                 continue
    
44.         yIdx = int((lat - y0 + ydelt / 2) / ydelt)
    
45.         if yIdx < 0 or yIdx > ynum - 1:
    
46.                 continue
    
47.        
    
48.         gData.Data[yIdx, xIdx] += 1
    
49. 
    
50. #---- Ouput grid data as Surfer ASCII data file
    
51. outfn = dataDir + "test.grd"
    
52. gData.SaveAsSurferASCIIFile(outfn)
    
53. 
    
54. #---- Create shaded layer from the grid data
    
55. aLS = LegendScheme(ShapeTypes.Polygon)
    
56. aLS.ImportFromXMLFile(dataDir + "test.lgs")
    
57. aLayer = DrawMeteoData.CreateShadedLayer(gData, aLS, "Lighting_Shaded", "Frequency")
    
58. aLayer.IsMaskout = True
    
59. mipy.MapDocument.ActiveMapFrame.AddLayer(aLayer)
    
60. mipy.MapDocument.ActiveMapFrame.MoveLayer(aLayer, 0)
    
61. 
    
62. print 'Finished!'

#--------------------------------------------------------        <br /> # Author: Yaqiang Wang                                           <br /> # Date: 2013-5-5                                                 <br /> # Purpose: Read and plot lighting frequency data  <br /> # Note: Sample                                                   <br /> #-----------------------------------------------------------  <br /> #---- Import clr and classes<br /> import clr<br /> clr.AddReferenceByPartialName("System")<br /> clr.AddReferenceByPartialName("System.Drawing")<br /> from System import *<br /> clr.AddReference("MeteoInfoC.dll")<br /> from MeteoInfoC import *<br /> from MeteoInfoC.Data import *<br /> from MeteoInfoC.Data.MeteoData import *<br /> from MeteoInfoC.Legend import *<br /> from MeteoInfoC.Shape import *<br /> import os.path<br /> <br /> #---- Set directory<br /> dataDir = "D:\\Temp\\ascii\\"<br /> <br /> #---- Define grid data<br /> print 'Define grid data...'<br /> x0 = 105<br /> y0 = 28<br /> xdelt = 0.1<br /> ydelt = 0.1<br /> xnum = 55<br /> ynum = 45<br /> gData = GridData(x0,xdelt,xnum,y0,ydelt,ynum)<br /> <br /> #---- Read lighting data file<br /> print 'Read lighting data file...'<br /> fn = dataDir + "cunshuju.txt"<br /> tf = open(fn)<br /> for aline in tf:        <br />         datalist = aline.split()<br />         lat = float(datalist[1])<br />         lon = float(datalist[0])        <br />         xIdx = int((lon - x0 + xdelt / 2) / xdelt)<br />         if xIdx < 0 or xIdx > xnum -1:<br />                 continue<br />         yIdx = int((lat - y0 + ydelt / 2) / ydelt)<br />         if yIdx < 0 or yIdx > ynum - 1:<br />                 continue<br />         <br />         gData.Data[yIdx, xIdx] += 1<br /> <br /> #---- Ouput grid data as Surfer ASCII data file<br /> outfn = dataDir + "test.grd"<br /> gData.SaveAsSurferASCIIFile(outfn)<br /> <br /> #---- Create shaded layer from the grid data<br /> aLS = LegendScheme(ShapeTypes.Polygon)<br /> aLS.ImportFromXMLFile(dataDir + "test.lgs")<br /> aLayer = DrawMeteoData.CreateShadedLayer(gData, aLS, "Lighting_Shaded", "Frequency")<br /> aLayer.IsMaskout = True<br /> mipy.MapDocument.ActiveMapFrame.AddLayer(aLayer)<br /> mipy.MapDocument.ActiveMapFrame.MoveLayer(aLayer, 0)<br /> <br /> print 'Finished!'

1. #--------------------------------------------------------        
   
2. # Author: Yaqiang Wang                                          
   
3. # Date: 2013-5-5                                                
   
4. # Purpose: Read and plot lighting data  
   
5. # Note: Sample                                                   
   
6. #-----------------------------------------------------------  
   
7. #---- Import clr and classes
   
8. import clr
   
9. clr.AddReferenceByPartialName("System")
   
10. clr.AddReferenceByPartialName("System.Drawing")
    
11. from System import *
    
12. from System.Drawing import *
    
13. from System.Collections.Generic import *
    
14. clr.AddReference("MeteoInfoC.dll")
    
15. from MeteoInfoC import *
    
16. from MeteoInfoC.Data.MapData import *
    
17. from MeteoInfoC.Geoprocess import *
    
18. from MeteoInfoC.Layer import *
    
19. from MeteoInfoC.Shape import *
    
20. from MeteoInfoC.Legend import *
    
21. from MeteoInfoC.Drawing import *
    
22. import os.path
    
23. 
    
24. #---- Create a point layer to plot lighting data
    
25. print 'Create lighting point layer...'
    
26. pLayer = VectorLayer(ShapeTypes.Point)
    
27. pLayer.LayerName = "Lighting_Point"
    
28. pLayer.Visible = True
    
29. pLayer.EditAddField("ID", String)
    
30. 
    
31. #---- Read lighting data file
    
32. print 'Read lighting data file...'
    
33. fn = "D:\\Temp\\ascii\\cunshuju.txt"
    
34. tf = open(fn)
    
35. id = 1
    
36. for aline in tf:       
    
37.         datalist = aline.split()
    
38.         lat = float(datalist[1])
    
39.         lon = float(datalist[0])       
    
40.         aPS = PointShape()
    
41.         aPoint = PointD(lon, lat)
    
42.         aPS.Point = aPoint       
    
43.         shapeNum = pLayer.ShapeNum
    
44.         if pLayer.EditInsertShape(aPS, shapeNum):
    
45.                 pLayer.EditCellValue("ID", shapeNum, str(id))
    
46.        
    
47.         id += 1
    
48. 
    
49. pLayer.UpdateLegendScheme(LegendType.SingleSymbol, "ID")
    
50. for legend in pLayer.LegendScheme.LegendBreaks:
    
51.         legend.Style = PointStyle.Plus
    
52.         legend.DrawOutline = False
    
53.         legend.Color = Color.LightBlue
    
54.         legend.Size = 8
    
55. 
    
56. mipy.MapDocument.ActiveMapFrame.AddLayer(pLayer)   
    
57. print 'Finished!'

#--------------------------------------------------------        <br /> # Author: Yaqiang Wang                                           <br /> # Date: 2013-5-5                                                 <br /> # Purpose: Read and plot lighting data  <br /> # Note: Sample                                                   <br /> #-----------------------------------------------------------  <br /> #---- Import clr and classes<br /> import clr<br /> clr.AddReferenceByPartialName("System")<br /> clr.AddReferenceByPartialName("System.Drawing")<br /> from System import *<br /> from System.Drawing import *<br /> from System.Collections.Generic import *<br /> clr.AddReference("MeteoInfoC.dll")<br /> from MeteoInfoC import *<br /> from MeteoInfoC.Data.MapData import *<br /> from MeteoInfoC.Geoprocess import *<br /> from MeteoInfoC.Layer import *<br /> from MeteoInfoC.Shape import *<br /> from MeteoInfoC.Legend import *<br /> from MeteoInfoC.Drawing import *<br /> import os.path<br /> <br /> #---- Create a point layer to plot lighting data<br /> print 'Create lighting point layer...'<br /> pLayer = VectorLayer(ShapeTypes.Point)<br /> pLayer.LayerName = "Lighting_Point"<br /> pLayer.Visible = True<br /> pLayer.EditAddField("ID", String)<br /> <br /> #---- Read lighting data file<br /> print 'Read lighting data file...'<br /> fn = "D:\\Temp\\ascii\\cunshuju.txt"<br /> tf = open(fn)<br /> id = 1<br /> for aline in tf:        <br />         datalist = aline.split()<br />         lat = float(datalist[1])<br />         lon = float(datalist[0])        <br />         aPS = PointShape()<br />         aPoint = PointD(lon, lat)<br />         aPS.Point = aPoint        <br />         shapeNum = pLayer.ShapeNum<br />         if pLayer.EditInsertShape(aPS, shapeNum):<br />                 pLayer.EditCellValue("ID", shapeNum, str(id))<br />         <br />         id += 1<br /> <br /> pLayer.UpdateLegendScheme(LegendType.SingleSymbol, "ID")<br /> for legend in pLayer.LegendScheme.LegendBreaks:<br />         legend.Style = PointStyle.Plus<br />         legend.DrawOutline = False<br />         legend.Color = Color.LightBlue<br />         legend.Size = 8<br /> <br /> mipy.MapDocument.ActiveMapFrame.AddLayer(pLayer)   <br /> print 'Finished!'

chengxf

学习了.谢谢版主.

▍|斯|_▍

努力学习，谢谢版主

koikaze

版主的程序关键地方可以加上注释吗？想研究一下，如何写脚本划分网格，谢谢。

MeteoInfo

koikaze 发表于 2013-5-6 15:26

                               
登录/注册后可看大图

版主的程序关键地方可以加上注释吗？想研究一下，如何写脚本划分网格，谢谢。

不是已经有注释了吗

koikaze

MeteoInfo 发表于 2013-5-6 16:10

                               
登录/注册后可看大图

不是已经有注释了吗

请教版主，代码中：cunshuju.txt
是怎样一种格式的？

原帖中的数据已经给lz删除了，想研究一下如何用代码处理网格和数据统计，感激不尽。

MeteoInfo

koikaze 发表于 2013-5-7 16:46

                               
登录/注册后可看大图

请教版主，代码中：cunshuju.txt
是怎样一种格式的？

原帖34楼已经讲了格式。

MeteoInfo

koikaze 发表于 2013-5-7 16:46

                               
登录/注册后可看大图

请教版主，代码中：cunshuju.txt
是怎样一种格式的？

这个帖子里有类似的数据（求助：有闪电定位资料，如何在grads中绘制闪电的分布情况？
http://bbs.06climate.com/forum.p ... 187&fromuid=106
），不过是纬度在前，经度在后，稍微改改脚本程序就可以了。

koikaze

MeteoInfo 发表于 2013-5-7 16:54

                               
登录/注册后可看大图

原帖34楼已经讲了格式。

了解，就是经纬度坐标数组。
还有一点疑问需要请教：
代码中：test.grd
是否为经纬度坐标下的地闪次数，数组：经度,纬度,地闪？

如果地闪数据是在某一时刻出现的次数是大于1次（如：113.435837,23.916666,4），如果只根据坐标出现的次数进行统计，会漏项，脚本应该如何修改比较完善呢？

MeteoInfo

koikaze 发表于 2013-5-7 17:02

                               
登录/注册后可看大图

了解，就是经纬度坐标数组。
还有一点疑问需要请教：
代码中：test.grd

test.grd是已经计算了闪电次数的格点数据。

把次数也读成一个变量，在给格点数据加闪电次数时用这个变量就可以了。