闪电数据的统计处理

chengxf

前几天有网友提到闪电数据的统计处理问题，因此有此帖。

闪电频次统计程序(LightningStatistics.exe)，是本人设计制作的一个Surfer软件的应用程序。所用程序语言是Pascal，编程环境为Embarcadero 的 RAD Studio XE3。如果仅仅根据闪电数据来统计计算闪电频次，生成一个自定义的网格(.grd)文件和一个典型的X(中心经度)、Y(中心纬度)、Z(频次)的数据(.dat)文件，并不需要Surfer软件。如果还要运行本程序的作图部分，需要在你的计算机上安装Surfer11软件。

下面就程序的运行使用方法做一些说明。

1、运行环境，Windows7系统，屏幕分辨率至少1024x768，当程序界面最大化时，如果有1280x800的分辨率则是更好的选择。根据运行功能的多少，安装或不安装Surfer软件，我是在Surfer11的版本下运行的，其他版本的Windows和Surfer是否适用，还须网友来检验。

2、程序界面纵向分为三个部分。最左边一部分是关于闪电数据的导入部分，这个数据文件就是原始的闪电数据，论坛上网友曾提供过一个这样的文件，共有8万8千多条记录，每一条记录就是一个经度和纬度指定的一个地理位置点。这里需要说明的是，原始数据是以多个空格作为分隔符，每条记录的前面有几个空格，然后是经度数据，再是几个空格，再是纬度数据。很遗憾也很抱歉，这样的格式本程序无法读入。目前本程序只能处理以逗号为分隔符的数据集。如下图所示，左边是原始格式，右边是处理规整以后的数据格式。
DataFormat0.jpg (22.76 KB, 下载次数: 21)

2013-5-10 14:56 上传

点击文件名下载附件
阅读权限: 10          DataFormat1.jpg (26.11 KB, 下载次数: 5)

2013-5-10 14:56 上传

点击文件名下载附件
阅读权限: 10

这样的规整工作可以在Surfer中实现，在Surfer中用“文件｜打开”命令，在“工作表”中打开原始文件，然后用保存命令保存成dat文件，有关参数取默认值即可。这样就得到了本程序需要的数据文件。

这一部分有两个按钮，“选择闪电记录数据文件”，仅仅是选择文件，用来判断是否是合格的数据文件。“导入并分析数据”按钮才进行真正的数据导入和分析工作。

3、程序界面的中间部分，是统计闪电频次的，你要自己确定“经度方向格点数”，例如，确定为50，就是说将经度方向分成50个等分。如果你将“自动确定纬度方向参数”复选框选中，则程序自动为你确定纬度方向的格点数和格点间隔。使得纬度间隔和经度间隔在地理长度上基本相等。我们知道，纬度间隔的地理长度在任何经度上都是一样的，但是某个经度间隔在不同的纬度上有不同的地理长度，越是接近两极，相同的经度间隔距离越小。具体算法这里不说了。“自动确定”就会尽量做到经纬度间隔的实际地理长度相等。使得统计频次时尽量在一个正方形的地理位置中。

当然，你也可以分别单独确定经纬度格点的多寡，从而单独确定网格的间隔。

4、按一下“统计闪电频次”，88000个数据的统计结果立刻出现在下面的表格中(大约耗时1秒)。除了表格数据外，还有频次的极大极小值、格点总数和频次总数的统计。这里的频次总数是各个中心位置附近的闪电次数合计而来，应该和左边部分的“记录数”一致，也以另一种方式验证了程序运行的正确性。

利用得到的频次统计结果，可以将表格数据“保存为GRD网格文件”，注意，这里保存的.grd网格文件是ASCII码格式的，并非二进制格式。它是已经网格化好了的文件了，当然在Surfer中照样可用，立刻就可以用它作等值线图了。但是注意这里的网格格点数是有限的，它就是经度方向格点数x纬度方向格点数，例子中只有50x41=2050个，这样做出来的图必然是粗糙的。

还可以将表格数据“保存为频次数据文件”。这个按钮保存的是标准的XYZ数据文件，后缀名为.dat，你可以用它在Surfer中显示成散点图(意义不大)，也可以用它选各种插值算法再做插值运算，这样再得出一个.grd网格文件，这里就可以将格点数调整大一些，例如200x164，这样得出的图形就要细致一些。

5、程序界面的最右边部分，是利用Surfer程序作图的部分。可以选择是直接利用统计时保存确定的网格文件作图，还是“利用频次数据重新网格化”来作图。选择前者，就没有“插值算法”可选；选择后者，除了可以选择插值算法以外，还可以重新确定X、Y方向的格点数，例如调整为200和164，或者更高。要作图，还需要底图文件和白化文件，也需要你确定。程序的例子中这两个文件都是chongqing.bln，既作底图，也是白化文件。

按了“调用Surfer作图”按钮，等待的时间要看你用的什么算法，以及节点数的多少。克里金算法较其他的算法耗时较长，节点数越多耗时越多。一般十几秒到一分钟，你可以看那个鼠标光标，小圆圈转啊转的，就是绘图没有结束。

随主程序LightningStatistics.exe一起打包的还包括：

1、关于程序的说明.pdf：即本文；

2、Chongqing.bln：例子中用到的重庆的底图文件，也是白化用到的文件，和主程序放在一起；

3、GrayScale.clr；Bathymetry.clr；Geology.clr；RainbowPastel.clr；RainbowLight.clr：5个等值线着色方案文件，这5个文件需要复制到Surfer11安装目录的ColorScales文件夹中，我的机器上是如下路径，参照自己实际的路径作修改：

C:\Program Files\Golden Software\Surfer 11\ColorScales

如果你要自定义一个彩色方案文件，请将该文件命名为LS_Custom.clr，就可以在本程序中使用了。

考虑到论坛上原始数据的提供者的意见，原始数据文件没有打包进来。(如果哪位网友可以提供测试用的原始闪电数据文件，我这里代大家表示感谢。)

建议你将主程序复制到某个文件夹中，自己再做一个快捷方式到桌面，用起来就方便了。

看到气象论坛有闪电统计的帖子，自己给自己下了个课题，用了几天时间编了这个小应用。气象界的同行们，如果这个小程序能对你在今后的部分工作中有所帮助，我是非常高兴的。这个程序你可以免费使用，你也可以copy介绍给其他朋友，但是你不能修改它和反编译它。当然最重要的是如果发现了什么错误，还麻烦你及时在本气象论坛上通知我。

论坛的MeteoInfo板块中有MeteoInfo大侠的“绘制闪电密度图”的帖子，绘出的图形很漂亮，值得学习，特此推荐。

闪电统计程序.rar (1.29 MB, 下载次数: 36)

2013-5-10 15:11 上传

点击文件名下载附件
阅读权限: 10
发现主程序中一个四舍五入的错误，对不起，下载这个吧。--
闪电统计程序(20130511).rar (1.29 MB, 下载次数: 93)

2013-5-11 14:48 上传

点击文件名下载附件
阅读权限: 10
程序补充和修改已完成，下面提供新的程序下载，上面的两个下载文件都不需要。有关说明见下文。
闪电统计程序(20130517).rar (1.55 MB, 下载次数: 204)

2013-5-17 14:31 上传

点击文件名下载附件
阅读权限: 10
程序界面.jpg (206.51 KB, 下载次数: 5)

2013-5-10 14:57 上传

点击文件名下载附件
阅读权限: 10                                                             

补充:

关于闪电统计软件中的计算方法，应一些网友的要求，在此做一个补充说明。

1、关于格点间距的计算

我们知道，地球是椭球体，球心到赤道的长半轴为 6378245米，而到两极的短半轴只有6356863米，扁率为298.3。以上数据根据Surfer软件中以北京-1954为大地基准、Krasovsky-1940的参考椭球体得到的。但在绝大多数情况下，我们还是可以将地球作为一个球体来看待，忽略它的椭球特性。因而我们取长短半轴的平均值6367554为地球的半径来进行下面的计算。

为了获得某一地区的闪电频次，我们要将该地区划分成很多矩形区域，以便统计该区域中发生的闪电次数。一般情况下，我们将经度方向的最大经度，减去最小经度，然后拿你确定的网格点数来除，就得到了经度方向的间隔值。这里我们只修改网格点数，而由程序来计算出间隔的经度值。也同样是网友的要求，我觉得如果能将这个经度间隔换算成具体的距离显示出来，可能更有意义，所以程序在这次添加了间隔单位为长度“米”的显示框。

那么，每个经度代表的长度是多少呢？如果在赤道，这个问题很简单：地球赤道的周长/360。而赤道的周长=Pi, 这里D是地球赤道的直径，6367554 x 2。特别说明的是pi值，由于采用程序的定义值，小数点后面大约是20几位，并非自己输入的3.14159等等。这样算出来的地球赤道周长是40008.522千米，大约为4万公里(题外话：这就是毛主席诗词“坐地日行八万里”的科学依据)。将此数值除以360，就得到每一经度在赤道的距离是111.135千米。由于我们假设地球是正球体，因此，穿越两极的任何一条经线，周长都是这么多，也就是说每一纬度的距离也是111.135千米。

但是如果我们不在赤道，例如我们在重庆，每一个经度的距离又是多少呢？每一个纬度的距离又是多少呢？每一纬度的距离还是不变，还是那么多。然而每一经度的距离就不一样了，因为在不同的纬度上，纬度线的周长不一样了。看下面的示意图。

闪电说明.jpg (35.27 KB, 下载次数: 4)

2013-5-17 14:31 上传

点击文件名下载附件
阅读权限: 10

图中的Q点假定是重庆。通过Q作一条垂直于赤道半径OA的垂线QB。可见CQ=OB，OQ=OA，通过Q点的绿色的纬度线周长，应该等于2 x CQ x pi，而红色的赤道的周长等于2 x OA x pi，两个周长之比就是CQ：OA，应该等于OB：OQ，那就是cos(∠AOQ)，就是Q点纬度值的余弦。

这是一个很重要的结论。经纬度间隔数值相同时，实际的地理距离数，南北方向上要大于东西方向，愈是接近两极这种差别就愈大。我们为了得到一个正方形的地理区域，经纬度间隔就不能取相同的值，当纬度间隔数值=经度间隔数值x cos(当地纬度)时，我们可以得到一个近似的正方形。

程序中，改变经度方向的格点数时，经度方向的格点间距=(经度极大值-经度极小值)/经度方向的格点数。同样，纬度方向的格点间距=(纬度极大值-纬度极小值)/纬度方向的格点数。当你将“自动确定纬度方向参数”打勾选中时，纬度方向参数计算方法不同，首先计算格点间距，再反过来计算格点数。纬度方向的格点间距=经度方向的格点间距x cos(图形纬度范围的中间值)，这个中间值=(纬度极大值-纬度极小值)/2。取中间值也是为了计算方便，对于一个东西方向跨度较大而南北方向跨度不大的地区，例如重庆市、广东湖北山东等省这样的地区，应该问题不大，对于南北方向跨度较大的地区，我还没有想到怎样处理，程序将会变得更加复杂。纬度方向的格点间距计算好了以后，再确定纬度方向的格点数，它=(纬度极大值-纬度极小值)/纬度方向的格点间距，得出的结果要整数，小数部分只进不舍。原来程序这里做了四舍五入，结果不行，有错误，已经改正了。

纬度方向的格点数和格点间距计算方法自动和手动时不相同，因此你可能发现相同的格点数在自动和手动计算时会有不同的纬度格点间距，这是正常的。

这次程序添加的格点间距(单位：米)的显示框，计算方法是：

纬度方向格点间距(单位:米)=地球半径 x pi x 纬度格点间距(单位:度)/180；

经度方向格点间距(单位:米)=地球半径 x pi x 经度格点间距(单位:度) x cos(图形纬度范围的中间值)/180 (注，这个中间值单位是弧度，要换算一下) 。

有了这个实际距离显示框，可以比较直观地了解你确定的网格的大小。这里要附带说明一下，网格面积不要选择得太小，有网友说要确定为1个平方千米的范围。我觉得对于闪电数据统计来说是太小了，这样反而不能得到比较好的等值线。不过这是个人意见，不算数的。

2、关于闪电频次的统计

这个都是用的落点统计办法，很多程序很多网友都做过，任何语言都可以实现，统计结果差别不大，不多说了。

以上是这次程序修订的一些补充说明，这次程序修改以后，在“自动确定纬度方向参数”的计算中，更加尽量地将经纬度间隔的矩形近似为正方形。因此在格点数和格点间隔上的计算结果与以前的程序不同。

谢谢使用者。更谢谢你们提出宝贵意见。

程贤辅 2013/5/17

感谢网友的关心，程序也作了一些修改，具体的修改内容请见附件内的说明文件，这里就不多说了。chengxf 2013/7/6
闪电统计程序(20130706).rar (1.67 MB, 下载次数: 303)

2013-7-6 16:56 上传

点击文件名下载附件
阅读权限: 10

这个闪电统计程序根据网友提出的意见作了进一步修改，最新版的如下(你如果没有耐心看完帖子，就不会看到这个版本——嘻嘻)：
闪电统计程序(20140928).rar (1.67 MB, 下载次数: 76)

2014-9-30 12:50 上传

点击文件名下载附件
阅读权限: 10

根据网友提出来安徽省常用的DL-II型闪电定位仪有两种格式的数据，一种是：
00000001 001001534 +114.72248 +31.90601 -102.887 0 1 +035.478 000015.864 000001.853 024.185 004 003 004 008 010
    这里，从第20个字符起，依次是经度值，纬度值，以及强度值。
另一种是：
00000001 2005 09 19 00 30 16 844 +32.80877 +118.26498 -031.232 +002.5 -002.6 00000.48 1
    这里，从第34个字符起，依次是纬度值，经度值，以及强度值。
因此，修改了程序，使得这两种数据都能够读取。
新的程序如下：
LightningStatistics20160516.rar (1.05 MB, 下载次数: 92)

2016-5-19 09:41 上传

点击文件名下载附件
阅读权限: 10
chengxf 2016.5.19
-------------------------------------------------
该程序的源代码已经公开，请参阅“ 闪电统计DelphiXE3源程序”的帖子。    程贤辅2014.4.2

有关该程序的使用实例以及该程序的最新版本，请参考： “闪电数据的统计处理”程序的使用一例     程贤辅2014.11.3

chengxf

宋秀鹏 发表于 2013-6-13 12:12

                               
登录/注册后可看大图

很清楚，明白。请问楼上的统计软件是自己开发的？还是有现成的？

最初是为了解决气象家园论坛上一些网友的具体问题而做的软件，也没有注明软件作者。后来功能更全面一些，也就注明了作者姓名，那个作者就是本人。5月10日放到论坛，至今才1个月零几天。程序肯定有不少纰漏，还望您指正。当前能够处理的闪电数据就是两种，一个是国家气象局行业标准，另一个是安徽省很多台站用的DL-II型闪电定位仪的数据，程序中大家可以看得到。闪电定位仪有很多种，各自产生的数据格式不尽相同。如果仔细看了软件中小小的帮助文档，可以知道我欢迎网友提供各种闪电数据，然后我来改进程序，使得它也能适合你的具体情况，在工作中助你一臂之力。
最后说一下，这一切都是免费的，我不收任何费用，图的是快乐——能够帮助你，我就快乐。

宋秀鹏

chengxf 发表于 2013-6-9 10:52

                               
登录/注册后可看大图

你要什么方面的详细说明?

这是用surfer处理的吗？

MeteoInfo

楼主够有心的，侠义情怀，赞一个！

mbt125

程老师真是爱钻研，热心啊

MeteoInfo

不过我比较好奇程老师怎么会有这么多时间，哈哈

chengxf

MeteoInfo 发表于 2013-5-10 17:10

                               
登录/注册后可看大图

不过我比较好奇程老师怎么会有这么多时间，哈哈

这不奇怪，您到了我这年龄同样有用不完的时间。哈哈。

kongfeng0824

谢谢程老师如此详细的帖子！！！

mofangbao

才看到，好帖被一堆提问帖子挡住了，得推起来才行，谢谢程老师

▍|斯|_▍

呀呀呀！ 太棒了太棒了！ 谢谢程老师！ 通过一个问题学习了好多东西！

chengxf

▍|斯|_▍ 发表于 2013-5-11 10:40

                               
登录/注册后可看大图

呀呀呀！ 太棒了太棒了！ 谢谢程老师！ 通过一个问题学习了好多东西！

我就盼望着你能试用一下(你有原始数据)，看看和你的期望还差多少。也欢迎提出意见。

▍|斯|_▍

chengxf 发表于 2013-5-11 10:56

                               
登录/注册后可看大图

我就盼望着你能试用一下(你有原始数据)，看看和你的期望还差多少。也欢迎提出意见。

         
程老师，我用了你做的闪电频次统计程序，选择的网格数是50x41，间距为0.1度x0.1度得到的闪电频次，和我用同样的网格设置用fortran程序求的闪电频次不一样，但是都刚好是两倍。