sismologia.cl scraper for recent earthquakes

sismologia.py

#!/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')