It goes without saying that Chile has a fair number of earthquakes.
Here in Concepción we seem to feel a reasonably sized earthquake (M4.5–5) every 1–2 months.
So, I wanted a quick and easy way to check for information about earthquakes from publicly available data sources.
USGS Database
Initially I worked with the USGS Earthquake page and their associated RSS feed.
I wrote a short script to grab the list of earthquakes from the past 24 hours and filter them by distance from a specified latitude/longitude pair:
#!/usr/bin/env python
# Load the daily USGS atom feed and see if there are any earthquakes within
# a specified distance of a specified location
import feedparser
import argparse
import math
import sys
# compute the great circle distance using the haversine formula
def GreatCircDist(loc1, loc2):
dlat = loc2[0] - loc1[0]
dlong = loc2[1] - loc1[1]
rdlat = math.radians(dlat)
rdlong = math.radians(dlong)
ha = math.sin(rdlat/2) * math.sin(rdlat/2) + \
math.cos(math.radians(loc1[0])) * \
math.cos(math.radians(loc2[0]))*math.sin(rdlong/2) * \
math.sin(rdlong/2)
hc = 2 * math.atan2(math.sqrt(ha), math.sqrt(1-ha))
return 6378.1*hc
parser = argparse.ArgumentParser(description="Load the daily USGS atom feed \
and see if there are any earthquakes within a specified distance of \
a specified location.")
parser.add_argument('-l', type=float, default=False,
help='Latitude of interest.')
parser.add_argument('-g', type=float, default=False,
help='Longitude of interest.')
parser.add_argument('-d', type=float, default=False,
help='Radius of interest (in km).')
args = parser.parse_args()
if not(args.l and args.g and args.d):
sys.stderr.write("Error. Incorrect arguments. Exiting.\n")
sys.exit(1)
usgs = feedparser.parse('http://earthquake.usgs.gov/earthquakes/feed/v1.0/summary/2.5_day.atom')
number = 0
for entry in usgs['entries']:
pos = [float(x) for x in entry['georss_point'].split()]
dist = GreatCircDist([args.l, args.g], pos)
if args.d > dist:
number += 1
sys.stdout.write(
entry['title'].split(' -')[0] + ' earthquake ' +
'{0:0.0f}'.format(dist) + ' km away,' +
entry['title'].split('-')[1] + ' (' +
str(pos[0]) + ', ' + str(pos[1]) + ') at an elevation of ' +
'{0:0.1f}'.format(float(entry['georss_elev'])/1000.) +
'km.\n')
if not(number):
sys.stdout.write('No earthquakes within ' + str(args.d) + ' km of (' +
str(args.l) + ', ' + str(args.g) + ') in the past 24h.\n')
else:
sys.stdout.write('\n' + str(number) + ' earthquake')
if number != 1:
sys.stdout.write('s')
sys.stdout.write(' within ' +
str(args.d) + ' km of (' + str(args.l) + ', ' +
str(args.g) + ') in the past 24h.\n')
For example, to find earthquakes within 1000 km of Concepción, you can invoke the script as so:
$ usgs_daily_earthquake_check.py -l -36.82 -g -73.04 -d 1000
No earthquakes within 1000.0 km of (-36.82, -73.04) in the past 24h.
Sismologia.cl
I was using that for a while, but eventually discovered that it was missing some earthquakes.
In a discussion with folks here, I learned of the Chilean government’s site, sismologia.cl , which has more earthquakes, but limited to Chile.
So, I wrote a similar script for that website:
#!/usr/bin/env python
#
# sismologia.py
#
# Scrape the sismologia.cl website list of most recent earthquakes > 3.0
# (http://www.sismologia.cl/links/ultimos_sismos.html) and localize them.
#
# Copyright (C) 2015 George C. Privon
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
import pycurl
import math
import numpy as np
import re
import sys
import argparse
from io import BytesIO
def GreatCircDist(loc1, loc2):
"""
compute the great circle distance between two lat/long pairs using the
haversine formula
Arguments:
- loc1, loc1: tuple containing (lat,long) pairs
Returns:
- great circle distance in km
"""
dlat = loc2[0] - loc1[0]
dlong = loc2[1] - loc1[1]
rdlat = math.radians(dlat)
rdlong = math.radians(dlong)
ha = math.sin(rdlat/2) * math.sin(rdlat/2) + \
math.cos(math.radians(loc1[0])) * \
math.cos(math.radians(loc2[0]))*math.sin(rdlong/2) * \
math.sin(rdlong/2)
hc = 2 * math.atan2(math.sqrt(ha), math.sqrt(1-ha))
return 6378.1*hc
def estimateIntensity(mag, dist):
"""
Estimate earthquake intensity (Modified Mercalli Scale) based on the
distance and magnitude of an earthquake. This intensity is estimated
using the "Did you feel it" technique from the USGS:
http://earthquake.usgs.gov/research/dyfi/
The Western USA equation seems to better match the data:
CDI = 1.15 + 1.01*Magnitude - 0.00054*Distance - 1.72*(log Distance)/(log 10)
Arguments:
- mag: earthquake magnitude
- dist: distance to the epicenter
Returns:
- value of the estimated intensity
"""
return 1.15 + 1.01 * mag - 0.00054 * dist - 1.72 * np.log(dist) / np.log(10)
parser = argparse.ArgumentParser(description="Load the sismologia.cl \
recent earthquakes page and see if there are any earthquakes within \
specified distance of a specified location.")
parser.add_argument('-l', '--latitude', type=float, default=False,
required=True,
help='Latitude of interest.')
parser.add_argument('-g', '--longitude', type=float, default=False,
required=True,
help='Longitude of interest.')
parser.add_argument('-d', '--distance', type=float, default=None,
help='Radius of interest (in km).')
parser.add_argument('-m', '--magnitude', type=float, default=None,
help='Minimum earthquake magnitude.')
parser.add_argument('-i', '--intensity', type=int, default=None,
help='Show only earthquakes whose estimated intensity \
(Modified Mercalli Intensity Scale) at the provided position is greater than \
this value.')
args = parser.parse_args()
if not(args.distance) and not(args.magnitude) and not(args.intensity):
parser.print_help()
sys.stderr.write("\nYou must specify either a distance, magnitude, or \
intensity.\n")
sys.exit(1)
sys.stdout.write('Searching sismologia.cl for recent earthquakes with: \n')
if args.distance:
sys.stdout.write('\tdistance from ({0:1.2f}, {1:1.2f}) < {2:1.0f} km\n'.
format(args.latitude, args.longitude, args.distance))
if args.magnitude:
sys.stdout.write('\tmagnitude > {0:1.1f}\n'.format(args.magnitude))
if args.intensity:
sys.stdout.write('\testimated Mod. Mercalli Intensity at ({0:1.2f}, \
{1:1.2f}) > {2:1.1f}\n'.format(args.latitude, args.longitude, args.intensity))
sys.stdout.write('\n')
ultimo = 'http://www.sismologia.cl/links/ultimos_sismos.html'
buf = BytesIO()
curl = pycurl.Curl()
curl.setopt(pycurl.URL, ultimo)
curl.setopt(pycurl.WRITEDATA, buf)
curl.perform()
curl.close()
body = buf.getvalue().decode('iso-8859-1').splitlines()
grab = False
number = 0
for line in body:
if re.search('<tbody>', line):
grab = True
continue
elif re.search('</tbody>', line):
grab = False
break
if grab:
match = True
if re.search('_self', line):
line = line.split('</td><td>')
eqpos = [float(line[2]), float(line[3])]
dist = GreatCircDist([args.latitude, args.longitude],
eqpos)
mag = float(line[5].split()[0])
intensity = estimateIntensity(mag, dist)
if args.distance and dist > args.distance:
match *= False
if args.magnitude and mag < args.magnitude:
match *= False
if args.intensity and estimateIntensity(mag, dist) < args.intensity:
match *= False
if match:
number += 1
sys.stdout.write(line[0].split('>')[3].split('</a')[0] + ': ' +
line[5] + ' earthquake ' +
'{0:0.0f}'.format(dist) + ' km away' +
' (' + line[2] + ', ' + line[3] + '; ' +
line[7].split('<')[0] +
') at a depth of ' + line[4] + ' km.\n')
if not(number):
sys.stdout.write('No recent earthquakes meet the search criteria.\n')
else:
sys.stdout.write('\n{0:1d} recent earthquake'.format(number))
if number != 1:
sys.stdout.write('s')
sys.stdout.write(' meet')
if number == 1:
sys.stdout.write('s')
sys.stdout.write(' the search criteria.\n')
The usage is very similar:
$ sismologia.py -l -36.82 -g -73.04 -d 500
...
3 earthquakes within 500.0 km of (-36.82, -73.04) in the past 24h.
Both scripts are available for use, and it should be fairly easy to deconstruct/adapt them to other uses.