绘制闪电密度图-2

qq200602116

辟如（72.1,10.1）、（72.3,10.1）这两个数据都该在第一个网格内的，按您的算法后面数据就不在第一网格了

1099221723

支持楼主。

koala54

楼主考虑的很严谨，把闪电密度定义成单位面积的闪电频次更加合理，用经纬网格统计得到的闪电频次除以球面多边形的面积，这样计算也很容易实现。我之前的思路是：先选择一个中心经纬度，建立2km x 2km 的水平网格坐标，再将闪电记录的经纬度坐标转换成该水平网格坐标，最后统计网格中的闪电数。这样涉及坐标的转换，实现起来也更复杂。好贴顶起~

helloqingfeng

学习了，谢谢楼主。

hqf2008jy

感谢樓主分享！没用过MeteoInfo，需要画闪电密度图，要慢慢研究啦！非常感謝！

hqf2008jy

前面下载脚本，扣了三次金币都没下载下来，金币不够了

呼小喵喵

王老师，我用的数据格式稍微有些不同，我把脚本里的格式修改运行后，提示错误如下：

提示错误

0514-ic.txt (9.5 MB, 下载次数: 0)

2017-5-15 16:54 上传

点击文件名下载附件
原始数据
希望老师能提示指导一下。

MeteoInfo

呼小喵喵 发表于 2017-5-15 16:55
王老师，我用的数据格式稍微有些不同，我把脚本里的格式修改运行后，提示错误如下：
希望老师能提示指导一 ...

参考这里：http://bbs.06climate.com/forum.p ... &extra=page%3D1

呼小喵喵

王老师，我用1楼的程序来画闪电密度分布，程序在计算格点面积那提示错误，麻烦您有空帮我看看是什么问题，谢谢。
#--------------------------------------------------------        
# Author: Yaqiang Wang                                          
# Date: 2013-5-5                                                
# Purpose: Read and plot lighting frequency data  
# Note: Sample                                                   
#-----------------------------------------------------------  
#---- Import clr and classes
import clr
clr.AddReferenceByPartialName("System")
clr.AddReferenceByPartialName("System.Drawing")
from System import *
from System.Collections.Generic import *
clr.AddReference("MeteoInfoC.dll")
from MeteoInfoC import *
from MeteoInfoC.Data import *
from MeteoInfoC.Data.MeteoData import *
from MeteoInfoC.Legend import *
from MeteoInfoC.Shape import *
from MeteoInfoC.Geoprocess import *
import os.path

#---- Set directory
dataDir = "F:\\Meteinfo\\MeteoInfo\\"

#---- Define grid data
print 'Define grid data...'
x0 = 108
y0 = 29
xdelt = 0.2
ydelt = 0.2
xnum = 88
ynum = 71
gData = GridData(x0,xdelt,xnum,y0,ydelt,ynum)

#---- Get area of each grid
print 'Get area of each grid...'
areaData = GridData(x0,xdelt,xnum,y0,ydelt,ynum)
for i in range(0,ynum):
        for j in range(0,xnum):
                pList = List[PointD]()
                x = x0 - xdelt / 2 + j * xdelt
                y = y0 - ydelt / 2 + i * ydelt
                pList.Add(PointD(x, y))
                y = y + ydelt
                pList.Add(PointD(x, y))
                x = x + xdelt
                pList.Add(PointD(x, y))
                y = y - ydelt
                pList.Add(PointD(x, y))
                pList.Add(pList[0])
                area = GeoComputation.SphericalPolygonArea(pList)
                areaData.Data[i, j] = area / 1000000   #Square meters to square kilometers

#---- Read lighting data file
print 'Read lighting data file...'
fn = dataDir + "0514CG.txt"
tf = open(fn)
for aline in tf:        
        datalist = aline.split()
        lat = float(datalist[3])
        lon = float(datalist[4])               
        xIdx = int((lon - x0 + xdelt / 2) / xdelt)
&#160; &#160;&#160; &#160;&#160;&#160;if xIdx < 0 or xIdx > xnum -1:
&#160; &#160;&#160; &#160;&#160; &#160;&#160; &#160;&#160; &#160; continue
&#160; &#160;&#160; &#160;&#160;&#160;yIdx = int((lat - y0 + ydelt / 2) / ydelt)
&#160; &#160;&#160; &#160;&#160;&#160;if yIdx < 0 or yIdx > ynum - 1:
&#160; &#160;&#160; &#160;&#160; &#160;&#160; &#160;&#160; &#160; continue
&#160; &#160;&#160; &#160;&#160;&#160;
&#160; &#160;&#160; &#160;&#160;&#160;gData.Data[yIdx, xIdx] += 1

#---- Get lighting times per square kilometers
for i in range(0,ynum):
&#160; &#160;&#160; &#160;&#160;&#160;for j in range(0,xnum):
&#160; &#160;&#160; &#160;&#160; &#160;&#160; &#160;&#160; &#160; gData.Data[i, j] = gData.Data[i, j] / areaData.Data[i, j]

#---- Ouput grid data as Surfer ASCII data file
outfn = dataDir + "0514CG.grd"
gData.SaveAsSurferASCIIFile(outfn)

#---- Create shaded layer from the grid data
#aLS = LegendScheme(ShapeTypes.Polygon)
#aLS.ImportFromXMLFile(dataDir + "test1.lgs")
aLS = LegendManage.CreateLegendSchemeFromGridData(gData, LegendType.GraduatedColor, ShapeTypes.Polygon)
aLayer = DrawMeteoData.CreateShadedLayer(gData, aLS, "Lighting_Shaded", "Times")
bLayer = DrawMeteoData.CreateGridFillLayer(gData, aLS, "Lighting_Grid", "Times")
aLS = LegendManage.CreateLegendSchemeFromGridData(areaData, LegendType.GraduatedColor, ShapeTypes.Polygon)
areaLayer = DrawMeteoData.CreateGridFillLayer(areaData, aLS, "Area_Grid", "Area")
mipy.MapDocument.ActiveMapFrame.AddLayer(aLayer)
mipy.MapDocument.ActiveMapFrame.MoveLayer(aLayer, 0)
mipy.MapDocument.ActiveMapFrame.AddLayer(bLayer)
mipy.MapDocument.ActiveMapFrame.MoveLayer(bLayer, 0)
mipy.MapDocument.ActiveMapFrame.AddLayer(areaLayer)
mipy.MapDocument.ActiveMapFrame.MoveLayer(areaLayer, 0)

print 'Finished!'

提出错误：
File "<iostream>", line 39
    ? ?? ???for j in range(0,xnum):
    ^
SyntaxError: no viable alternative at character '?'
我用的数据格式是这样的，根据格式，我对脚本里读闪电数据修改了。
ID        日期        时间        lat        lon        强度        陡度
1        2015/5/14        0:05:20        29.289        117.7208        55.2        12.6
2        2015/5/14        0:08:40        29.2667        117.7798        110.9        17.5
3        2015/5/14        0:09:41        29.3316        117.8805        74.7        18.8
4        2015/5/14        0:10:36        28.8085        117.7493        -54.3        -25.7
5        2015/5/14        0:12:22        29.2603        117.7879        92.9        15.7
6        2015/5/14        0:14:33        28.7751        117.8876        -56.1        -26.6
7        2015/5/14        0:14:33        29.0603        117.0943        -46.7        -15.3
8        2015/5/14        0:18:34        29.2807        117.7976        102.2        17
9        2015/5/14        0:19:22        27.7117        114.0031        -46.1        -13.7
10        2015/5/14        0:21:39        29.2238        117.8531        131.8        19.4

阿凡达

谢谢楼主了