地震数据的可视化

25年来全球六级地震的可视化

从2000年算起，全球共发生了六级以上地震3627次（参考USGS）。

我写了一个小脚本来展示地震的发生地，详见视频：

代码思路

思路是使用matplotlib绘图，依次展示全球DEM、全球行政区划边界和地震点，地震点的颜色根据地震等级变化，

1. 导入所需的库

   包括用于数据处理的 pandas 和 geopandas，用于绘图的 matplotlib 和 rasterio，以及其他一些库。

   import geopandas as gpd
   import matplotlib.pyplot as plt
   import matplotlib.animation as animation
   from matplotlib.colors import ListedColormap, BoundaryNorm
   from matplotlib import rcParams
   from datetime import datetime
   import pandas as pd
   from tqdm import tqdm
   import numpy as np
   import rasterio
   from rasterio.plot import show
   from pandas import Period
   

2. 读取数据

   读取 DEM_TIF 文件，全球边界 Shapefile 文件和地震数据 Shapefile 文件。地震相关的矢量数据下载见参考网址2，DEM数据使用的是ETOPO 的全球DEM，全球各国的矢量数据使用的是中科院的各国矢量数据。

   # 读取Shapefile文件
   DEM_TIF=r"E:\\ruiduobao\\博客\\地震数据可视化\\矢量\\DEM.tif"
   world = gpd.read_file('E:\\ruiduobao\\博客\\地震数据可视化\\矢量\\全球边界wgs84.shp')
   earthquakes = gpd.read_file('E:\\ruiduobao\\博客\\地震数据可视化\\矢量\\地震6级以上.shp')
   # 读取 DEM_TIF 文件
   with rasterio.open(DEM_TIF) as src:
       dem = src.read(1)
       transform = src.transform
   

3. 数据预处理

   将地震数据的 ‘Mag’ 列转换为浮点数，将 ‘Year’, ‘Month’, ‘Day’, ‘Time’ 列合并成一个 datetime 对象，并按 ‘Datetime’ 列排序。

   earthquakes['Mag'] = earthquakes['Mag'].astype(float)
   earthquakes['Datetime'] = pd.to_datetime(earthquakes['Year'].astype(str) + '-' + earthquakes['Month'].astype(str) + '-' + earthquakes['Day'].astype(str) + ' ' + earthquakes['Time'])
   earthquakes = earthquakes.sort_values('Datetime')
   

4. 创建图形和子图

   创建一个新的图形，并在图形上创建两个子图。一个子图用于绘制地图，另一个子图用于绘制折线图。

   fig = plt.figure(figsize=(20, 10))
   ax1 = plt.subplot2grid((1, 3), (0, 0), colspan=2)
   

5. 绘制初始地图

   在地图子图上显示 DEM 数据，绘制世界地图，设置地图的边界，并绘制初始地震点。

   show(dem, ax=ax1, transform=transform, cmap='gray', vmin=-5000, vmax=7000)
   world.plot(edgecolor='lightblue', facecolor='none', zorder=1, ax=ax1)
   ax1.set_xlim(-180, 180)
   ax1.set_ylim(-90, 90)
   

6. 计算每个月的地震次数

   将地震数据按月份分组，并计算每个月的地震次数。

   earthquakes['YearMonth'] = earthquakes['Datetime'].dt.to_period('M')
   monthly_counts = earthquakes['YearMonth'].value_counts().sort_index()
   monthly_counts.index = monthly_counts.index.to_timestamp().date
   

7. 创建折线图

   在折线图子图上创建一个线对象，用于绘制折线图。

   ax2 = plt.subplot2grid((1, 3), (0, 2), colspan=1)
   line, = ax2.plot([], [], lw=2)
   

8. 创建动画

   定义动画的每一帧，每一帧中，获取这个月的地震数据，更新地震点，更新日期文本，重新设置地图的边界，更新折线图的数据，更新进度条。

   grouped = earthquakes.groupby(earthquakes['Datetime'].dt.to_period('M'))
   def animate(month):
       # 动画的每一帧的代码
   ani = animation.FuncAnimation(fig, animate, frames=monthly_counts.index, interval=50)
   

9. 保存动画

   将创建的动画保存为 mp4 格式的视频。

   # 保存动画为视频
   ani.save('earthquakes.mp4', writer='ffmpeg')
   
   # plt.show()
   

详细代码

import geopandas as gpd
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from matplotlib.colors import ListedColormap, BoundaryNorm
from matplotlib import rcParams
from datetime import datetime
import pandas as pd
from tqdm import tqdm
import numpy as np
import rasterio
from rasterio.plot import show
from pandas import Period

# 设置字体为 SimHei 显示中文
rcParams['font.sans-serif'] = 'SimHei'
rcParams['axes.unicode_minus'] = False


# 读取Shapefile文件
DEM_TIF=r"E:\\ruiduobao\\博客\\地震数据可视化\\矢量\\DEM.tif"
world = gpd.read_file('E:\\ruiduobao\\博客\\地震数据可视化\\矢量\\全球边界wgs84.shp')
earthquakes = gpd.read_file('E:\\ruiduobao\\博客\\地震数据可视化\\矢量\\地震6级以上.shp')
# 读取 DEM_TIF 文件
with rasterio.open(DEM_TIF) as src:
    dem = src.read(1)
    transform = src.transform
# 将 'Mag' 列转换为整数
earthquakes['Mag'] = earthquakes['Mag'].astype(float)
# 按 'Datetime' 列排序
# 将 'Year', 'Month', 'Day', 'Time' 列合并成一个 datetime 对象
earthquakes['Datetime'] = pd.to_datetime(earthquakes['Year'].astype(str) + '-' + earthquakes['Month'].astype(str) + '-' + earthquakes['Day'].astype(str) + ' ' + earthquakes['Time'])
# 按 'Datetime' 列排序
earthquakes = earthquakes.sort_values('Datetime')
# 创建自定义颜色映射
cmap = ListedColormap(['#3285BA', '#14A542', '#E6F658',  '#D91F1F'])
# 创建边界规范
norm = BoundaryNorm([6, 6.7, 7.4, 8.2, 10], cmap.N)
#处理进度条
progress_bar = tqdm(total=len(earthquakes))

# 创建一个新的图形
fig = plt.figure(figsize=(20, 10))

# 创建一个子图用于绘制地图
ax1 = plt.subplot2grid((1, 3), (0, 0), colspan=2)

# 使用 transform 参数将 DEM 数据的像素坐标转换为地理坐标
show(dem, ax=ax1, transform=transform, cmap='gray', vmin=-5000, vmax=7000)
# 绘制世界地图
world.plot(edgecolor='lightblue', facecolor='none', zorder=1, ax=ax1)

# 设置固定的边界
ax1.set_xlim(-180, 180)
ax1.set_ylim(-90, 90)


# 创建一个空的 GeoDataFrame，用于存储已经显示的地震点
displayed_earthquakes = earthquakes.iloc[[0]]
# 绘制初始地震点并创建图例
displayed_earthquakes.plot(ax=ax1, column='Mag', cmap='RdYlGn_r', markersize=5, vmin=4, vmax=6.5, legend=False, zorder=2)

# 在图形左侧添加一个文本对象，用于显示日期
date_text = ax1.text(-170, 80, '', fontsize=12)

# 计算每个月的地震次数
earthquakes['YearMonth'] = earthquakes['Datetime'].dt.to_period('M')
monthly_counts = earthquakes['YearMonth'].value_counts().sort_index()
# 将 PeriodIndex 转换为日期类型
monthly_counts.index = monthly_counts.index.to_timestamp().date


# 计算累计次数
cumulative_counts = monthly_counts.cumsum()

# 创建一个子图用于绘制折线图
ax2 = plt.subplot2grid((1, 3), (0, 2), colspan=1)

# 创建一个用于绘制折线图的线对象
line, = ax2.plot([], [], lw=2)

# 设置折线图的标题和标签
ax2.set_title('6级以上地震发生次数')
ax2.set_xlabel('年-月')
ax2.set_ylabel('次数')
# 设置 y 轴的下限为 0，去掉负数轴
ax2.set_ylim(bottom=0)
# 在开始时创建一个 scatter 对象
scatter = ax1.scatter([earthquakes.iloc[0].geometry.x], [earthquakes.iloc[0].geometry.y], c=[earthquakes.iloc[0]['Mag']], cmap=cmap, norm=norm, s=5, zorder=2)
# 添加颜色图例
cbar = plt.colorbar(scatter, ax=ax1)
cbar.set_label('地震等级')

# 将地震数据按月份分组
grouped = earthquakes.groupby(earthquakes['Datetime'].dt.to_period('M'))


# 定义动画的每一帧

def animate(month):
        # 获取这个月的地震数据
        monthly_earthquakes = grouped.get_group(Period(month, freq='M'))
        
        # 获取新的地震点的坐标和震级
        new_offsets = np.array([monthly_earthquakes.geometry.x, monthly_earthquakes.geometry.y]).T
        new_array = np.array(monthly_earthquakes['Mag'])
        
        # 将新的地震点添加到 scatter 对象中
        scatter.set_offsets(np.concatenate([scatter.get_offsets(), new_offsets]))
        scatter.set_array(np.concatenate([scatter.get_array(), new_array]))

        # 更新日期文本
        date_text.set_text(month.strftime('%Y-%m'))
        date_text.set_color('white')
        # 重新设置地图的边界
        ax1.set_xlim(-180, 180)
        ax1.set_ylim(-90, 90)
        
        # 更新折线图的数据
        line.set_data(monthly_counts.index[:monthly_counts.index.get_loc(month)+1], monthly_counts.iloc[:monthly_counts.index.get_loc(month)+1])
        ax2.set_xlim(monthly_counts.index[0], monthly_counts.index[-1])
        ax2.set_ylim(0, monthly_counts.max())
        
        # 更新 tqdm 进度条
        progress_bar.update(len(monthly_earthquakes))

# 创建动画
ani = animation.FuncAnimation(fig, animate, frames=monthly_counts.index, interval=50)
# 保存动画为视频
ani.save('earthquakes.mp4', writer='ffmpeg')

# plt.show()

参考：

https://earthexplorer.usgs.gov/

https://ds.iris.edu/ieb/index.html