Last active
September 18, 2015 16:52
-
-
Save privong/497490fa2a1f49b490dc to your computer and use it in GitHub Desktop.
sismologia.cl scraper for recent earthquakes
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
#!/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') |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment