python地震数据处理_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_))