国家地震信息中心 (NEIC) 是美国地质调查局 (USGS) 的一部分,监测全球地震并记录其位置和大小。该地震活动数据库可供公众用于研究和教育目的。用户可以通过 API 访问这些数据并获取不同格式的数据,例如 CSV、GeoJSON、KML 等...

这个博客显示,

• 如何在python中使用USGS API key获取地震数据。

• 了解 GeoJson 格式。

• 将数据导出为 GeoPackage 并在 QGIS 中可视化。

从 USGS 获取地震数据

USGS记录全球地震活动并将数据公开,我们可以从地震.usgs.gov访问这些数据。我们还可以获取_Past Hour_、Past Day、Past 7 Days 和_Past 30 Days_ 的数据。由于这些实时数据每分钟更新一次,如果我们需要定期获取数据进行分析,使用 API 获取数据效率更高。编写完 Python 脚本后,我们使用 CronJob/windows Scheduler 运行程序并获得空间数据层,无需任何人工干预。

现在,Lawrence Hayes 博士无需任何人工操作即可获取每小时的地震活动数据作为空间层。 圣安地列斯的电影场景

让我们编码。

python 模块是我们可以在脚本中安装和使用的可重用代码块。以下模块用于此项目。

os模块用于连接目录和文件名的路径。所以我们不必担心这段代码运行的是哪个操作系统。 IT 为我们做了正确的事。

requests模块用于获取特定的URL值,这里是USGS地震活动数据。

json模块用于读取请求返回的对象并将其转换为 JSON 格式。

time模块用于获取当前时间,因为我们打算最终运行一个 cronJob,命名文件。更重要的是,它用于将地震活动发生时间转换为人类可读的格式。

geopandas模块是流行的pandas库的空间版本。我们可以创建 GeoDataFrame,读写矢量图层。按照惯例,它被导入为gpd

time,geopandas许多在默认安装中不可用。我发现Anaconda对解决安装包依赖问题非常有帮助。谢谢蟒蛇! :)

import os
import json
import time 
import requests
import geopandas as gpd

访问地震.usgs.gov网页以获取 GeoJSON 摘要 API URL。在这里您可以选择不同的格式,例如 KML 准备使用矢量文件,我们可以将其加载到Google 地球以进行 3D 可视化。 USGS 已设置此 KML 文件的样式。但是对于这个项目,让我们坚持 GeoJSON 总结。它将以结构化格式向我们返回一个包含所有详细信息的文本文件。

API - 应用程序接口允许我们与远程系统或软件通信并获取信息。令人惊讶的是,在最流行的 API 中,有两个与地理空间世界有关。谷歌地图和Accuweather。使用 requests 模块中的.get()方法可以让我们与 USGS 交互并查询我们需要的数据。由于API是一个特定的结构,我们不必担心操作系统的操作。

response =requests.get('https://earthquake.usgs.gov/earthquakes/feed/v1.0/summary/all_hour.geojson')

任何request调用的实际输出格式都是 JSON。当我们添加geometry -> type, and coordinates等空间属性时,就称为GeoJSON。这些格式是人类可读的,并且可以使用索引和键在 python 中快速完成数据抽象。

response = response.json()
response['features'][0]

# Displaying first feature from the GeoJSON data.
{'type': 'Feature',
 'properties': {'mag': 1.8,
  'place': '30 km SW of Skwentna, Alaska',
  'time': 1632888597630,
  'updated': 1632888797590,
  'tz': None,
  'url': 'https://earthquake.usgs.gov/earthquakes/eventpage/ak021chwm6mr',
  'detail': 'https://earthquake.usgs.gov/earthquakes/feed/v1.0/detail/ak021chwm6mr.geojson',
  'felt': None,
  'cdi': None,
  'mmi': None,
  'alert': None,
  'status': 'automatic',
  'tsunami': 0,
  'sig': 50,
  'net': 'ak',
  'code': '021chwm6mr',
  'ids': ',ak021chwm6mr,',
  'sources': ',ak,',
  'types': ',origin,',
  'nst': None,
  'dmin': None,
  'rms': 0.69,
  'gap': None,
  'magType': 'ml',
  'type': 'earthquake',
  'title': 'M 1.8 - 30 km SW of Skwentna, Alaska'},
 'geometry': {'type': 'Point', 'coordinates': [-151.807, 61.7972, 99.3]},
 'id': 'ak021chwm6mr'}

GeoPandas 改变了空间数据与 python 的工作方式,因为它建立在强大的模块pandas之上。操作与 pandas 模块非常相似。

让我们创建一个空的 GeoDataFrame 并将每个地震记录添加到其中。

exportDf = gpd.GeoDataFrame()

这个循环将为每个特征创建一个新的 DataFrame,并且只存储震级、地点、时间、警报、海啸状态和几何形状。我们将所有记录保存为带有纬度和经度信息的点要素。由于都是GPS坐标,所以CRS是EPSG:4326。 USGS 给出的time以毫秒为单位。所以我们可以通过将其除以 1000 将其转换为纪元时间并转换为人类可以理解的格式。

for i, data in enumerate(response['features']):  
    gdf = gpd.GeoDataFrame()
    coord = data['geometry']['coordinates']
    geometry = gpd.points_from_xy([coord[0]], [coord[1]])
    gdf = gpd.GeoDataFrame(data['properties'], index = [i], geometry = geometry, crs='EPSG:4326')
    gdf = gdf[['mag', 'place', 'time', 'alert','status','tsunami', 'geometry']]
    gdf['time'] = time.strftime('%Y-%m-%d', time.gmtime(gdf['time'][i]/1000))
    exportDf = exportDf.append(gdf)

由于主 GeoDataFrame 没有定义 CRS,让我们为其定义 CRS。

if exportDf.crs == None:
    exportDf.set_crs('epsg:4326', inplace=True)

在导出过程中,我们可以将其保存在特定目录中。另外,由于我们每小时都在获取数据,如果每次都生成一个新的 GeoPackage,那将是一团糟,因此通过创建一个一天的 GeoPackage,我们可以将每小时的数据存储为一个层。

对于 2021 年 9 月 30 日,文件名为20210930,图层名为00_13、01_12、02_06....24_10。小时后的数字指示器是该小时内发生的地震总数。

directory = 'earthquakes'
name = strftime("%Y%m%d", time.localtime())
filename = os.path.join(directory, name + '.gpkg')
layername = "_".join([strftime("%H", time.localtime()), str(len(response['features']))])

exportDf.to_file(driver='GPKG', 
                 filename=filename, 
                 layer=layername, 
                 encoding='utf-8')
print(f'Export {name} {layername} sucessfull')

下载笔记本

阅读我的下一篇博客,了解如何设置 cron 作业并使用空间连接过滤到特定区域。

快乐学习:)