绘制闪电密度图-2

MeteoInfo · 发表于 2013-5-30 11:29:20

登录后查看更多精彩内容~

您需要登录才可以下载或查看，没有帐号？立即注册

x

本帖最后由 MeteoInfo 于 2013-5-30 15:29 编辑

之前发了一个帖子利用MeteoInfo脚本绘制闪电密度图（绘制闪电密度图http://bbs.06climate.com/forum.php?mod=viewthread&tid=14130&fromuid=106
），看到紫夜怡静发的数据（http://bbs.06climate.com/forum.p ... tra=page%3D1&page=4）把脚本改了改来绘制闪电定位点图和0.5度*0.5度格点的闪电发生次数图，当然格点的设置可以自己在脚本里改。

先看看生成闪电定位点图脚本：

demo_lighting-1.py (2.74 KB, 下载次数: 142)

2013-5-30 11:05 上传

点击文件名下载附件
1. 创建一个点图层，加入相应的属性字段.
2. 打开闪电数据文件，读取每个闪电的相应字段信息。经纬度信息被用来生成点图元，其它字段放在属性库里。“陡度”数据里有一些'$'字符，所以做了点特殊处理，“定位方式”里的中文也要做相应的处理。
3. 设置点图层的图例。
4. 将点图层加入地图中。

如前所述，MeteoInfo中的图层是有属性信息的（GIS软件的特征），可以浏览所有闪电的属性，也可以看某一个闪电的属性信息：

还可以将属性信息标注出来：

也可以利用属性信息来编辑图例：

再来看计算格点闪电次数并生成等值线图层的脚本：

demo_lighting_f-1.py (2.14 KB, 下载次数: 69)

2013-5-30 11:22 上传

点击文件名下载附件

test1.lgs (3.77 KB, 下载次数: 57)

2013-5-30 13:42 上传

点击文件名下载附件
1. 定义格点，通过起始经纬度、经纬度间隔、格点x/y方向个数这些参数来定义格点。
2. 打开闪电数据文件，读取每个闪电的经纬度信息，判断该闪电落在哪个格点中，将所在格点的闪电次数加1。遍历所有闪电后就可以得到每个格点的闪电发生次数。
3. 利用上述步骤生成的格点数据（gData）生成等值线填充图层和Grid_Fill图层。图例可以用格点数据自动生成一个，也可以加载已经编辑好的图例，或者在脚本中详细的设置出来。

生成的等值线填充图层：

Grid_Fill图层：

可以将Grid_Fill图层的每个格点的闪电次数标注出来，和闪电定位点图进行对比验证：

脚本文件上面已经给出了，这里把代码也再贴一遍：

# This Python file uses the following encoding: GB2312
#--------------------------------------------------------
# Author: Yaqiang Wang
# Date: 2013-5-5
# Purpose: Read and plot lighting data
# Note: Sample
#-----------------------------------------------------------
#---- Import clr and classes
import clr
clr.AddReferenceByPartialName("System")
clr.AddReferenceByPartialName("System.Drawing")
from System import *
from System.Drawing import *
from System.Collections.Generic import *
clr.AddReference("MeteoInfoC.dll")
from MeteoInfoC import *
from MeteoInfoC.Data.MapData import *
from MeteoInfoC.Geoprocess import *
from MeteoInfoC.Layer import *
from MeteoInfoC.Shape import *
from MeteoInfoC.Legend import *
from MeteoInfoC.Drawing import *
import os.path
#---- Create a point layer to plot lighting data
print 'Create lighting point layer...'
pLayer = VectorLayer(ShapeTypes.Point)
pLayer.LayerName = "Lighting_Point"
pLayer.Visible = True
pLayer.EditAddField("ID", Int32)
pLayer.EditAddField("日期", String)
pLayer.EditAddField("时间", String)
pLayer.EditAddField("强度", Double)
pLayer.EditAddField("陡度", Double)
pLayer.EditAddField("误差", Double)
pLayer.EditAddField("定位方式", String)
#---- Read lighting data file
print 'Read lighting data file...'
fn = "D:\\Temp\\ascii\\2009_06_06.txt"
tf = open(fn)
for aline in tf:
datalist = aline.split()
if aline.strip() == '':
continue
id = int(datalist[0])
date = datalist[1]
time = datalist[2]
lat = float(datalist[3].split('=')[1])
lon = float(datalist[4].split('=')[1])
intensity = float(datalist[5].split('=')[1])
gradient = float(datalist[6].split('=')[1].replace('$', '0'))
error = float(datalist[7].split('=')[1])
position = unicode(datalist[8], 'GB2312')
position = position.split('：')[1]
aPS = PointShape()
aPoint = PointD(lon, lat)
aPS.Point = aPoint
shapeNum = pLayer.ShapeNum
if pLayer.EditInsertShape(aPS, shapeNum):
pLayer.EditCellValue("ID", shapeNum, id)
pLayer.EditCellValue("日期", shapeNum, date)
pLayer.EditCellValue("时间", shapeNum, time)
pLayer.EditCellValue("强度", shapeNum, intensity)
pLayer.EditCellValue("陡度", shapeNum, gradient)
pLayer.EditCellValue("误差", shapeNum, error)
pLayer.EditCellValue("定位方式", shapeNum, position)
pLayer.UpdateLegendScheme(LegendType.SingleSymbol, "ID")
for legend in pLayer.LegendScheme.LegendBreaks:
legend.Style = PointStyle.Plus
legend.DrawOutline = False
legend.Color = Color.Blue
legend.Size = 8
mipy.MapDocument.ActiveMapFrame.AddLayer(pLayer)
print 'Finished!'

# This Python file uses the following encoding: GB2312 
#--------------------------------------------------------         
# Author: Yaqiang Wang                                           
# Date: 2013-5-5                                                 
# Purpose: Read and plot lighting data   
# Note: Sample                                                    
#-----------------------------------------------------------   
#---- Import clr and classes 
import clr 
clr.AddReferenceByPartialName("System") 
clr.AddReferenceByPartialName("System.Drawing") 
from System import * 
from System.Drawing import * 
from System.Collections.Generic import * 
clr.AddReference("MeteoInfoC.dll") 
from MeteoInfoC import * 
from MeteoInfoC.Data.MapData import * 
from MeteoInfoC.Geoprocess import * 
from MeteoInfoC.Layer import * 
from MeteoInfoC.Shape import * 
from MeteoInfoC.Legend import * 
from MeteoInfoC.Drawing import * 
import os.path 
 
#---- Create a point layer to plot lighting data 
print 'Create lighting point layer...' 
pLayer = VectorLayer(ShapeTypes.Point) 
pLayer.LayerName = "Lighting_Point" 
pLayer.Visible = True 
pLayer.EditAddField("ID", Int32) 
pLayer.EditAddField("日期", String) 
pLayer.EditAddField("时间", String) 
pLayer.EditAddField("强度", Double) 
pLayer.EditAddField("陡度", Double) 
pLayer.EditAddField("误差", Double) 
pLayer.EditAddField("定位方式", String) 
 
#---- Read lighting data file 
print 'Read lighting data file...' 
fn = "D:\\Temp\\ascii\\2009_06_06.txt" 
tf = open(fn) 
for aline in tf:         
        datalist = aline.split() 
        if aline.strip() == '': 
                continue 
                
        id = int(datalist[0]) 
        date = datalist[1] 
        time = datalist[2]         
        lat = float(datalist[3].split('=')[1]) 
        lon = float(datalist[4].split('=')[1])         
        intensity = float(datalist[5].split('=')[1]) 
        gradient = float(datalist[6].split('=')[1].replace('$', '0')) 
        error = float(datalist[7].split('=')[1]) 
        position = unicode(datalist[8], 'GB2312') 
        position = position.split('：')[1] 
         
        aPS = PointShape() 
        aPoint = PointD(lon, lat) 
        aPS.Point = aPoint         
        shapeNum = pLayer.ShapeNum 
        if pLayer.EditInsertShape(aPS, shapeNum): 
                pLayer.EditCellValue("ID", shapeNum, id) 
                pLayer.EditCellValue("日期", shapeNum, date) 
                pLayer.EditCellValue("时间", shapeNum, time) 
                pLayer.EditCellValue("强度", shapeNum, intensity) 
                pLayer.EditCellValue("陡度", shapeNum, gradient) 
                pLayer.EditCellValue("误差", shapeNum, error) 
                pLayer.EditCellValue("定位方式", shapeNum, position) 
 
pLayer.UpdateLegendScheme(LegendType.SingleSymbol, "ID") 
for legend in pLayer.LegendScheme.LegendBreaks: 
        legend.Style = PointStyle.Plus 
        legend.DrawOutline = False 
        legend.Color = Color.Blue 
        legend.Size = 8 
 
mipy.MapDocument.ActiveMapFrame.AddLayer(pLayer)    
print 'Finished!'

#--------------------------------------------------------
# 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 *
clr.AddReference("MeteoInfoC.dll")
from MeteoInfoC import *
from MeteoInfoC.Data import *
from MeteoInfoC.Data.MeteoData import *
from MeteoInfoC.Legend import *
from MeteoInfoC.Shape import *
import os.path
#---- Set directory
dataDir = "D:\\Temp\\ascii\\"
#---- Define grid data
print 'Define grid data...'
x0 = 72
y0 = 10
xdelt = 0.5
ydelt = 0.5
xnum = 126
ynum = 88
gData = GridData(x0,xdelt,xnum,y0,ydelt,ynum)
#---- Read lighting data file
print 'Read lighting data file...'
fn = dataDir + "2009_06_06.txt"
tf = open(fn)
for aline in tf:
datalist = aline.split()
lat = float(datalist[3].split('=')[1])
lon = float(datalist[4].split('=')[1])
xIdx = int((lon - x0 + xdelt / 2) / xdelt)
if xIdx < 0 or xIdx > xnum -1:
continue
yIdx = int((lat - y0 + ydelt / 2) / ydelt)
if yIdx < 0 or yIdx > ynum - 1:
continue
gData.Data[yIdx, xIdx] += 1
#---- Ouput grid data as Surfer ASCII data file
outfn = dataDir + "test.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")
mipy.MapDocument.ActiveMapFrame.AddLayer(aLayer)
mipy.MapDocument.ActiveMapFrame.MoveLayer(aLayer, 0)
mipy.MapDocument.ActiveMapFrame.AddLayer(bLayer)
mipy.MapDocument.ActiveMapFrame.MoveLayer(bLayer, 0)
print 'Finished!'

#--------------------------------------------------------         
# 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 * 
clr.AddReference("MeteoInfoC.dll") 
from MeteoInfoC import * 
from MeteoInfoC.Data import * 
from MeteoInfoC.Data.MeteoData import * 
from MeteoInfoC.Legend import * 
from MeteoInfoC.Shape import * 
import os.path 
 
#---- Set directory 
dataDir = "D:\\Temp\\ascii\\" 
 
#---- Define grid data 
print 'Define grid data...' 
x0 = 72 
y0 = 10 
xdelt = 0.5 
ydelt = 0.5 
xnum = 126 
ynum = 88 
gData = GridData(x0,xdelt,xnum,y0,ydelt,ynum) 
 
#---- Read lighting data file 
print 'Read lighting data file...' 
fn = dataDir + "2009_06_06.txt" 
tf = open(fn) 
for aline in tf:         
        datalist = aline.split() 
        lat = float(datalist[3].split('=')[1]) 
        lon = float(datalist[4].split('=')[1])                
        xIdx = int((lon - x0 + xdelt / 2) / xdelt) 
        if xIdx < 0 or xIdx > xnum -1: 
                continue 
        yIdx = int((lat - y0 + ydelt / 2) / ydelt) 
        if yIdx < 0 or yIdx > ynum - 1: 
                continue 
         
        gData.Data[yIdx, xIdx] += 1 
 
#---- Ouput grid data as Surfer ASCII data file 
outfn = dataDir + "test.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") 
mipy.MapDocument.ActiveMapFrame.AddLayer(aLayer) 
mipy.MapDocument.ActiveMapFrame.MoveLayer(aLayer, 0) 
mipy.MapDocument.ActiveMapFrame.AddLayer(bLayer) 
mipy.MapDocument.ActiveMapFrame.MoveLayer(bLayer, 0) 
 
print 'Finished!'

MeteoInfo · 发表于 2013-5-30 15:16:36

本帖最后由 MeteoInfo 于 2013-5-30 15:25 编辑

之前说过用次数这种变量受到格点大小的影响太大，如果能算出单位面积的闪电发生次数会更合理些。MeteoInfo里有面积计算的功能，所以对脚本做了些扩展计算出每平方公里的闪电次数。

demo_lighting_f-2.py (3.23 KB, 下载次数: 54)

2013-5-30 15:12 上传

点击文件名下载附件

由于面积的单位是每平方公里，而格点的划分是等经纬度，其实需要计算每个等经纬度格点的地球球面面积。MeteoInfoC类库里的GeoProcess命名空间的GeoComputation类里有里SphericalPolygonArea方法来计算球面面积，参数是多边形所有顶点的列表。在原来脚本的基础上增加创建一个GridData（areaData），获得每个格点四个顶点的经纬度坐标生成顶点列表，计算出每个格点的球面面积。在计算出每个格点的闪电次数后将每个格点的次数除以面积获得单位面积（每平方公里）的闪电次数，然后在生成等值线图层。

从areaData生成的格点多边形图层，可以看出对于等经纬度网格，低纬度比高纬度格点的面积大，同纬度的格点面积一样。

具体脚本代码如下。如上所述，其实同纬度上的所有格点计算其中一个格点的面积就可以了，其它的都是一样，因此脚本中计算面积的代码还可以改得更高效一些（有兴趣的朋友可以试试）。

#--------------------------------------------------------
# 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 = "D:\\Temp\\ascii\\"
#---- Define grid data
print 'Define grid data...'
x0 = 72
y0 = 10
xdelt = 0.5
ydelt = 0.5
xnum = 126
ynum = 88
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 + "2009_06_06.txt"
tf = open(fn)
for aline in tf:
datalist = aline.split()
lat = float(datalist[3].split('=')[1])
lon = float(datalist[4].split('=')[1])
xIdx = int((lon - x0 + xdelt / 2) / xdelt)
if xIdx < 0 or xIdx > xnum -1:
continue
yIdx = int((lat - y0 + ydelt / 2) / ydelt)
if yIdx < 0 or yIdx > ynum - 1:
continue
gData.Data[yIdx, xIdx] += 1
#---- Get lighting times per square kilometers
for i in range(0,ynum):
for j in range(0,xnum):
gData.Data[i, j] = gData.Data[i, j] / areaData.Data[i, j]
#---- Ouput grid data as Surfer ASCII data file
outfn = dataDir + "test.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!'

#--------------------------------------------------------         
# 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 = "D:\\Temp\\ascii\\" 
 
#---- Define grid data 
print 'Define grid data...' 
x0 = 72 
y0 = 10 
xdelt = 0.5 
ydelt = 0.5 
xnum = 126 
ynum = 88 
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 + "2009_06_06.txt" 
tf = open(fn) 
for aline in tf:         
        datalist = aline.split() 
        lat = float(datalist[3].split('=')[1]) 
        lon = float(datalist[4].split('=')[1])                
        xIdx = int((lon - x0 + xdelt / 2) / xdelt) 
        if xIdx < 0 or xIdx > xnum -1: 
                continue 
        yIdx = int((lat - y0 + ydelt / 2) / ydelt) 
        if yIdx < 0 or yIdx > ynum - 1: 
                continue 
         
        gData.Data[yIdx, xIdx] += 1 
 
#---- Get lighting times per square kilometers 
for i in range(0,ynum): 
        for j in range(0,xnum): 
                gData.Data[i, j] = gData.Data[i, j] / areaData.Data[i, j] 
 
#---- Ouput grid data as Surfer ASCII data file 
outfn = dataDir + "test.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!'

紫夜怡静 · 发表于 2013-5-30 17:13:13

拜谢，楼主！不过MeteoInfo还未学习，收藏后慢慢研究！再次感谢！

kongfeng0824 · 发表于 2013-5-30 21:11:29

学习了谢谢程老师！

指点江山zcl · 发表于 2013-5-30 21:52:53

同谢程老师，以后应该会用的着

姜勇 · 发表于 2013-8-13 16:00:51

大开眼界啊，目前正在学习阶段

{:eb515:}{:eb515:}

rceclx · 发表于 2013-11-17 15:08:43

本帖最后由 rceclx 于 2013-11-17 15:10 编辑

王老师，您好！
“2. 打开闪电数据文件，读取每个闪电的经纬度信息，判断该闪电落在哪个格点中，将所在格点的闪电次数加1。gData.Data[yIdx, xIdx] += 1” 这部分要是改成“将所在格点的不同经纬度的闪电的强度平均”，需要怎么修改脚本？看了您的这个脚本，实在没有想到好的方法，请老师帮助一下。谢谢！

MeteoInfo · 发表于 2013-11-17 23:14:12

rceclx 发表于 2013-11-17 15:08

登录/注册后可看大图

王老师，您好！
“2. 打开闪电数据文件，读取每个闪电的经纬度信息，判断该闪电落在哪个格点中，将所在格点 ...

仔细看看此贴的沙发

rceclx · 发表于 2013-11-17 23:46:01

MeteoInfo 发表于 2013-11-17 23:14

登录/注册后可看大图

仔细看看此贴的沙发

谢谢王老师！仔细看了一下，计算每个格点的闪电次数是“gData.Data[yIdx, xIdx] += 1，但不知道如何对每个格点内不同经纬度的闪电的强度（intensity）这个属性进行平均计算？MeteoInfo里相关的方法是什么？对脚本用到的方法还不熟悉，请王老师讲一下。

MeteoInfo · 发表于 2013-11-18 09:06:55

rceclx 发表于 2013-11-17 23:46

登录/注册后可看大图

谢谢王老师！仔细看了一下，计算每个格点的闪电次数是“gData.Data[yIdx, xIdx] += 1，但不知道如何对每个 ...

在原脚本基础上再定义一个GridData来存放闪电强度的累加值，在每次判断出一个闪电所在的格点后，给该格点的闪电次数加一次同时将该闪电的闪电强度加入该格点（两个GridData，一个存闪电次数，另一个存闪电强度的累加，格点一一对应），所有闪电遍历完成后用闪电强度累加数据除以闪电次数就可以获得闪电强度平均。比如存闪电次数的格点数据为gData（参照原脚本），存闪电强度累加值的格点数据是sData，闪电强度平均值格点数据：

sData = sData / gData

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

绘制闪电密度图-2

登录后查看更多精彩内容~

评分