用Python分析过去10年地震分布~

回复“书籍”即可获赠Python从入门到进阶共10本电子书

大家好，我是才哥。

前段时间在黄海海域发生了5.2级地震，当时朋友圈华东一带的朋友们都“震动”了。

今天，就带大家一起看看过去10年里地震的分布情况！

目录：

• 1. 数据获取

• 2. 数据处理

• 2.1. 经纬度转换

• 2.2. 发震时刻数据处理

• 2.3. 地震所在地（省市）

• 3. 统计与可视化

• 3.1. 历年地震数

• 3.2. 各省地震数

• 3.3. 各城市地震数

• 3.4. 地震分布散点图

• 3.5. 地址分布热力图

• 3.6. 近10年每月地震分布图（动态）

• 4. 补充知识

1. 数据获取

直接在中国地震台网就可以下载到过去10年的地震数据。

中国地震台网：

https://news.ceic.ac.cn/index.html?time=1637318776

大家直接在历史查询里选择时间范围进行查询后，点击保存到本地即可获取数据。

存在本地的数据预览如下：

2. 数据处理

在数据获取的时候， 我们发现获取的数据包含时间、震级、经纬度、深度以及参考位置。而经纬度是百度地图下的，考虑到不同地图的经纬度会存在差异，而我们后续绘图采用的是高德地图，所以这里需要对经纬度进行转换。

2.1. 经纬度转换

高德开放平台有对应的经纬度转换的API，我们可以直接调用（免费额度够）。

import requests
import pandas as pd
# 读取数据
df = pd.read_excel(r'eqList.xlsx')

longitude_list = []
latitude_list = []
# 百度地图经纬度转换为高德地图经纬度
for i , location in enumerate(df[['经度(°)','纬度(°)']].values):
    location = str(location[0])+','+str(location[1])
    url = 'https://restapi.amap.com/v3/assistant/coordinate/convert?'
    
    parames = {
        'locations':location,
        'coordsys':'baidu',
        'key':'你的key',
        }
    
    r = eval(requests.get(url, params=parames).json()['locations'])
    # 经度
    longitude_list.append(r[0])
    # 纬度
    latitude_list.append(r[1])
    print(f'\r{i+1}',end='')

df['经度(°)'] = longitude_list
df['纬度(°)'] = latitude_list

如此，我们就成功将百度地图系的经纬度坐标转换为了高德地图系的经纬度坐标。

2.2. 发震时刻数据处理

对于发震时刻来说，我希望的是精确到月份以及小时，用于后续统计分析。

# 转换为时间格式
df['发震时刻'] = pd.to_datetime(df['发震时刻'])
# 获取年月
df['月份'] = df['发震时刻'].apply(lambda x: str(x)[:7])
# 获取小时
df['小时'] = df['发震时刻'].dt.hour

这里我感觉我的获取年月的操作有点low，有小伙伴有更好的方案还请留言哈，我要学习。

2.3. 地震所在地（省市）

由于在原始数据中参考位置无法便捷的解析出省和市，我这边打算用经纬度信息通过高德的API来进行获取。参考此前《你知道海底捞在全国有多少家店吗？》。

citys = []
provinces = []
for i , location in enumerate(df[['经度(°)','纬度(°)']].values):
    location = str(location[0])+','+str(location[1])
    url = 'https://restapi.amap.com/v3/geocode/regeo?'
    params = {
        'location':location,
        'key':'你的key',
        'extensions':'base',
        'batch':'false',
        'roadlevel':0,
        }
    
    r = requests.get(url, params=params)
    data = r.json()['regeocode']
    city = data['addressComponent']['city']
    province = data['addressComponent']['province']
    if len(city)==0:
        city = province
    citys.append(city)
    provinces.append(province)
    print(f'\r{i+1}',end='')
df['城市'] = citys
df['省'] = provinces

经过以上的处理，我们最终得到了以下数据：

3. 统计与可视化

这部分我们只进行简单的统计展示，不做类似预测或者别的什么深度分析，毕竟基于现有原始数据也得不出太大的有价值结论。

3.1. 历年地震数

中小地震次数

按照[2, 4.6]作为中小地震等级区间，近10年一共6188次，年均600余次！

# 中小地震（[2,4.6]）
df['年份'] = df['发震时刻'].dt.year
df_cn = df[df['省']!='[]']
df_xiao = df_cn.query('2<=`震级(M)`<=4.6')
df_xiao.groupby('年份')['发震时刻'].count().to_frame('次数').reset_index()

破坏性地震次数

按照[4.7, ∞]作为中小地震等级区间，近10年一共505次，年均50余次！

3.2. 各省地震数

考虑到一些地震发生在海域，这里统一化为中国海域，另外对于发生在国外的咱们这里不参与。

另外，我们只看破坏性地震分布，可以发现在我国新疆、西藏、云南、台湾和四川是地震高发带！

# 各省 破坏性地震
df_province = df_cn.query('`震级(M)`>=4.7').groupby('省')['发震时刻'].count().to_frame('次数').sort_values(by='次数',ascending=False).reset_index()

3.3. 各城市地震数

在我国233个城市发生过地震，发生地震次数最多的城市大多集中在新疆、四川。

注：这里统计的是全部地震数据

像新疆的克孜勒苏柯尔克孜自治州、和田地区

四川的宜宾和雅安等

3.4. 地震分布散点图

我们只看破坏性地震分布，根据经纬度坐标信息绘制

import folium.plugins as plugins
import folium

df = df[df['震级(M)']>=4.7]
ss1 = [[latitude,longitude] for latitude,longitude in df[['纬度(°)', '经度(°)']].values.tolist()]
m = folium.Map([39.904989, 116.405285],
              tiles='https://webrd01.is.autonavi.com/appmaptile?&x={x}&y={y}&z={z}&lang=zh_cn&size=1&scale=1&style=7',
               zoom_start=6,
               attr='高德')
groups = folium.FeatureGroup('')
for l,label in zip(ss1,df['参考位置'].to_list()):
    groups.add_child(
        folium.CircleMarker(
            location=l,
            radius=1,
            color='red',
            fill=True,
            fill_color='red',
            fill_opacity=0.4,
            popup=folium.Popup(html=label,max_width=100),
        )
    )
m.add_child(groups)
m.add_child(folium.LatLngPopup()) 
m.save('地震分布.html')               

3.5. 地址分布热力图

可以看到新疆西南部，和台湾省东部海域属于地震多发地段。

data_all = df[['纬度(°)', '经度(°)','震级(M)']].values.tolist()
m = folium.Map([39.904989, 116.405285],
               tiles='https://webrd01.is.autonavi.com/appmaptile?&x={x}&y={y}&z={z}&lang=zh_cn&size=1&scale=1&style=7',
               zoom_start=6,
               attr='高德')
hm = plugins.HeatMap(data_all, radius=10, 
                     gradient={.1: 'green', .6: 'yellow', 1: 'red'},
                     )
hm.add_to(m)
hm.save('地震热力图.html')

3.6. 近10年每月地震分布图（动态）

注：这里统计的是全部地震数据

绘制代码

data_move = []
date_list = df['月份'].sort_values().unique()
for month in date_list:
    data_move.append(df[df['月份'] == month][['纬度(°)', '经度(°)','震级(M)']].values.tolist())
m = folium.Map([39.904989, 116.405285],
               tiles='https://webrd01.is.autonavi.com/appmaptile?&x={x}&y={y}&z={z}&lang=zh_cn&size=1&scale=1&style=7',
               zoom_start=6,
               attr='高德')
time_index = df['月份'].sort_values().unique().tolist()
hm = plugins.HeatMapWithTime(data_move, 
                             index=time_index,
                             radius=10)

hm.add_to(m)
hm.save('地震动态热力图.html')

4. 补充知识

地震带分布图

我们筛选 5级以上地震绘制分布图如下：

基本吻合~

    小伙伴们，快快用实践一下吧！如果在学习过程中，有遇到任何问题，欢迎加我好友，我拉你进Python学习交流群共同探讨学习。

------------------- End -------------------

往期精彩文章推荐：

• Python中encode和encoding的区别是什么？

• 手把手教你使用openpyxl库从Excel文件中提取指定的数据并生成新的文件（附源码）

• 手把手教你用Python来模拟绘制自由落体运动过程中的抛物线（附源码）

• 盘点3种Python网络爬虫过程中的中文乱码的处理方法
  

欢迎大家点赞，留言，转发，转载，感谢大家的相伴与支持

想加入Python学习群请在后台回复【入群】

万水千山总是情，点个【在看】行不行

/今日留言主题/

随便说一两句吧~~