首页 > NumPy 阅读数:65

NumPy实例:地震数据的统计分析

利用 NumPy 可以对大数据进行简单的统计分析,包括数据文件的读取、数据列的提取、数据类型的转换和各种统计计算,还有数据的排序、搜索和计数等。现以中国历史地震数据分析为例来介绍简单的统计分析思路。

【例 1】将中国历史地震数据文件命名为 earthquakes.csv 文件。该文件中包含的数据列有日期、时间、地震发生纬度与经度、深度(km)、震级、烈度和精确度等。数据包括 1969 年 12 月 22 日以前中国历史上所发生的地震数据,其中部分数据列表显示如下。

日期,时间(北京时),纬度(度-分),经度(度-分),深度(km),震级,烈度,精确度
408/00/00,::.,"39°00""","100°30""",,4.75,VI,9,
416/00/00,::.,"34°18""","105°30""",,5,VI,9,
421/00/00,::.,"41°36""","120°24""",,5,VI,9,
462/08/16,::.,"35°00""","116°48""",,6,VIII,9,
495/03/31,::.,"37°30""","121°30""",,5.5,VII,9,
…………………………………………………………………………………………………………………
1969/12/7,52:39.0,"23°42""","121°54""",3.3,5.1,,9
1969/12/17,00:04.0,"18°30""","110°36""",,5.1,,3
1969/12/19,51:05.0,"31°12""","88°12""",,4,,2
1969/12/20,09:15.0,"18°12""","110°18""",2.4,5.2,,2
1969/12/21,38:57.0,"39°54""","76°54""",3.3,4,,9


请完成下列数据分析任务:

(1)计算历史地震数据中震级和烈度的最大值、最小值、算术平均值和中位数。

(2)统计不同震级发现次数。

(3)统计发生6级以上震级次数的总和。

(4)列出7.5级以上地震的发生日期、时间、经度、纬度和震级信息。

中国历史地震数据统计分析思路如下

1. 文件数据列的读取

利用 NumPy 中的 loadtxt 文件读取命令,将 earthquakes.csv 文件中的震级和烈度这两列数据分别读取并存储到 magnitude 数组和 lido 数组中。

在读取文件时,loadtxt 命令需要使用 delimiter、usecols、dtype 和 skiprows 参数,其中,delimiter 表示文件数据列的分割符为逗号,usecols 是确定读取文件中的列数位置,dtype 是确定数组的数据类型,skiprows 表示跳过第一行的行数。

由于地震的烈度分为 12 级,分别用罗马数字 I,II,…,XII 来表示。但是,在统计计算时需要将罗马数字转换成对应的 1~12 数字表示,因此,在读取烈度列时,loadtxt 命令还需要使用 converters 创建烈度列与转换函数之间进行映射的字典。同时,创建将烈度列中的罗马数字转换成对应的 1~12 数字的自定义函数 lido_conv(x)。

2. 自定义函数 lido_conv(x)

在自定义函数 lido_conv(x) 中,先创建数组 keys 用于存储空格和 1~12 的罗马数字,并将数组 keys 转换成列表 keys_list,而列表中每个元素的索引正是该元素的罗马数字对应的阿拉伯数字。由于从 CSV 文件中读取的数据是二进制字符串,因此需要将其转换成字符串。然后,使用 NumPy 的 where() 函数,将罗马数字转换成对应的阿拉伯数字。

3. 计算历史地震数据中震级和烈度的最大值、最小值、算术平均值和中位数

计算最大值、最小值、算术平均值和中位数的统计函数分别为 max()、min()、mean() 和 median() 等。由于烈度列数据中有许多空格,在计算烈度的平均值时要将烈度为 0 的项排除,因此,计算烈度平均值不能使用 mean() 函数,而是要使用 sum() 函数求出烈度数据总和,再使用 count_nonzero() 函数求出烈度数据不为 0 的总个数,然后将总和与总个数相除。

4. 统计不同震级发现的次数

利用 NumPy 中的 unique() 函数来实现。

5. 统计发生 6 级以上震级次数的总和

利用 NumPy 中的 where() 和 count_nonzero() 函数来实现。

6. 列出 7.5 级以上地震的发生日期、时间、经度、纬度和震级信息

利用 NumPy 中的 where() 函数来实现。

其程序代码 test1 如下。
# -*- coding: utf-8 -*-
import numpy as np
def lido_conv(x):
    keys = np.array(['','I','II','III','IIII','V','VI','VII','VIII',
                         'IX','X','XI','XII'])
    keys_list = keys.tolist()       #将数组转换成列表

    a = str(x, encoding="utf-8")    #将二进制字符转换成字符
    # 使用NumPy的where函数,将罗马数字转换成对应的阿拉伯数字
    result = np.where(a in keys,keys_list.index(a),0)
    return result

#取出震级存入magnitude数组和烈度存入lido数组
magnitude= np.loadtxt('earthquakes.csv',delimiter = ",",
                      usecols=(5,),dtype=float,skiprows=1)
lido = np.loadtxt('earthquakes.csv',delimiter = ",",
                  usecols=(6,),dtype=np.int,skiprows=1,
                  converters={6:lido_conv})
print(magnitude,lido)
#计算震级和烈度的最大值、最小值、算术平均值和中位数
print ("最大震级 =", np.max(magnitude))
print ("最大震级 =", np.amax(magnitude))
print ("最大烈度 =", np.max(lido))
print ("最小震级 =", np.min(magnitude))
print ("最小震级 =", np.amin(magnitude))
print ("最小烈度 =", np.min(lido))
print ("震级平均值 =", np.mean(magnitude))
print ("烈度平均值 =", np.mean(lido))
print ("震级中位数 =", np.median(magnitude))
print ("烈度平均值 =", np.sum(lido)/np.count_nonzero(lido))

#统计不同震级发生的次数
print ('不同震级震级数组:',np.unique(magnitude))
print('不同震级发现次数:')
arr_u,u_inverse = np.unique(magnitude, return_counts=True)
print(u_inverse)

#统计发生6级以上的震级次数总和
m_index = np.where(magnitude>6)      #震级大于6的元素的索引
m_sum = np.count_nonzero(m_index)
print(m_sum)

#列出7.5级以上地震的发生日期、时间、经度、纬度和震级信息
arr = np.loadtxt('earthquakes.csv',delimiter = ',',
                 dtype=np.string_,usecols=(0,1,2,3),skiprows=1)
print('读取的地震数据:',arr)
cond = np.where(magnitude>7.5)
print ('7.5级以上地震的元素索引值:',cond)
print ('使用元素索引值提取元素:')
print ('7.5级以上地震数据:',arr[cond].astype(np.string_))