GSEASRX

For Earthquake Science

0%

绘制事件的时间分布图

发表于 更新于 分类于 阅读次数： Disqus：

  本文介绍利用python在时间坐标下对地震事件进行统计和可视化。包括地震频度的统计、频度图的绘制以及M-t图的绘制。
  重点需要掌握，时间轴上统计和绘图时需要完成的时间段和时间点的转换，以及时间轴样式的设置。

绘制M-t图

简介

在python中，利用matplotlib.pyplot.stem可以画茎叶图，stem的参数可以改变垂直线的类型，顶点的颜色大小等。

1
stem(x,y, linefmt=None, markerfmt=None, basefmt=None)

说明

  • x, y分别是横纵坐标；
  • linefmt:垂直线的颜色和类型；linefmt=‘r-’，代表红色的实线;
  • basefmt指y=0那条直线；
  • markerfmt设置顶点的类型和颜色，比如C3.，C(大写字母C)是默认的，后面数字应该是0-9，改变颜色，最后的.或者o（小写字母o）分别可以设置顶点为小实点或者大实点，空格表示没有顶点。

以下为线型对应的字符

字符 线型
‘-’ solid line
‘–’ dashed line
‘-.’ dash-dot line
‘:’ dotted line

绘制N-t图

简介

python中可以利用matplotlib.pyplot.bar，通过绘制bar实现地震频度图的绘制。
工具：

  • matplotlib.pyplot
  • matplotlib.dates
  • pandas
1
2
3
import matplotlib.pyplot as plt

plt.bar(x,height,width=0.8,bottom=None,align='center',data=Nonecolor='b',alpha=1...)

参数说明

xerr, yerr:
scalar or array-like of shape(N,) or shape(2, N), optional；If not None, add horizontal / vertical errorbars to the bar tips. The values are +/- sizes relative to the data:

  • scalar: symmetric +/- values for all bars
  • shape(N,): symmetric +/- values for each bar
  • shape(2, N): Separate - and + values for each bar. First row contains the lower errors, the second row contains the upper
    errors.
  • None: No errorbar. (Default)

实例

Notice!
本节最关键的问题在于地震数目的统计，和时间坐标的控制

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
#!/bin/python

import pandas as pd
import matplotlib.dates as mdate
import matplotlib.pyplot as plt

# macos下解决中文乱码
plt.rcParams['font.sans-serif'] = ['Arial Unicode MS']
plt.rcParams['axes.unicode_minus'] = False

mc = 3  # 设置最小震级
# 读入目录为DataFrame
eq = pd.read_csv('gsty.txt',sep=' ',
                 encoding='gbk',skipinitialspace=True,
                 names=['date','time','lat','lon','dep','mag','diming'])
eq['datetime'] = eq.apply(lambda x: x.date+' '+x.time, axis=1)
eq['datetime_f'] = pd.to_datetime(eq['datetime'])

# 设置图片布局
fig = plt.figure(figsize=(10,5),dpi=100)

# 利用stem绘制M-t图
ax1 = fig.add_subplot(211)
x = eq[eq.mag>=mc].datetime_f
y = eq[eq.mag>=mc].mag
ax1.stem(x,y,linefmt='k-', basefmt='r-', markerfmt='C2 ',label='地震')
plt.ylim(mc-1,eq.mag.max()+1)
plt.ylabel('震级(Ml)')
plt.gcf().autofmt_xdate()
plt.xlabel("时间")
plt.legend()

# 绘制频度图前的准备
# 将时间设置为index，便于按时间段进行数量统计
eq.set_index('datetime_f',inplace=True,drop=True)
# eq_p = eq.to_period('M') # 按月显示
# print(eq_p.index.asfreq('M'))
afreq = eq.mag.resample('A').count()  # 按年统计
mfreq = eq.mag.resample('M').count()  # 按月统计
dfreq = eq.resample('d').count().to_period('d')  # 按天统计

# 利用bar绘制频度图
ax2 = fig.add_subplot(212)
x = dfreq.index.start_time
y = dfreq.mag
#设置x轴为时间格式，这句非常重要，否则x轴显示的将是类似于‘736268’这样的转码后的数格式
ax2.xaxis.set_major_formatter(mdate.DateFormatter('%Y-%m-%d'))
plt.xticks(pd.date_range(dfreq.index[0].start_time,dfreq.index[-1]end_time,freq='10d'), rotation=45)
ax2.bar(x,y,color='c',align='edge',alpha=0.8)
plt.ylabel('频数')
plt.xlabel("时间")
plt.show()

运行以上代码生成如下图片

fig1

说明

以下为在时间坐标下绘制数据时的关键顺序

  1. 读取数据后，利用pd.set_index设置时间列为index
  2. 利用时间index对表数据进行统计eq.resample('d').count().toperiod('d')，此处d表示按天采样，获取每日地震数量。
    注意：此时生成的dataframe中index为时间段，不可直接用于时间坐标（时间坐标用时间点）
  3. 设置x轴时间格式，ax.xaxis.set_major_formatter(mdate.DateFormatter('%Y-%m-%d'))
  4. 对index时间段取start_time作为时间点，利用始末时间点生成时间序列，作为x轴坐标；pd.date_range
  5. 绘图

👉 Notice:
可选参数color， edgecolor， linewidth，xerr和yerr可以是标量或长度等于bar数目的序列。这是绘制条形图作为堆叠条形图或烛台图的基础。:wstar:
详细信息：xerr和yerr被直接传递到:meth:‘errorbar’，所以他们也可以有形状2xN,独立规定上下误差。