Revised Earthquake Monitoring

Recently, there was a M8.3 earthquake North of Santiago, which caused a substantial amount of damage near the epicenter and a few deaths. The earthquake was felt, with only moderate shaking and no damage, down here in Concepción.

I had been running my earthquake checking script periodically, but the large earthquake did not show up in the list because I had only been searching within a 500 km radius of Concepción. My first thought was to simply increase the search radius to, say, 1000 km. The large earthquake then appeared in the list, but so did numerous other small earthquakes that I would not have felt here. Next I thought about adding a magnitude cutoff, so the script would only return information on the stronger earthquakes. But I can sometimes feel M3.0 earthquakes if they are close to me, but not M6.0 that are up by Santiago.

How do I select only those earthquakes which I might have felt?

While thinking about this, I remembered the USGS “did you feel it” information, collected for various earthquakes. If you go to that link and select the Intensity Vs. Distance tab, it shows a plot of the reported intensity as a function of distance from the epicenter, along with some fits to the data. It occurred to me that, with the formula, I could use the reported magnitude and calculated distance to estimate the intensity felt in Concepción and only show those earthquakes above a certain threshold.

I opted to use the California shake model on the Did You Feel It? page, as it seemed to better approximate the intensity–distance relation of the Chilean earthquakes. So, with a modification to the script, you can now select only earthquakes based on maximum distance, minimum magnitude, minimum estimated intensity at your position, or some combination of those.

I have been running this updated script for a few days now, and I haven’t seen an earthquake reported which I have not felt (I set no distance or magnitude limits, but did set a minimum Modified Mercalli Intensity of > 3.5). Of course, I have also not felt any earthquakes since then, so I have not been able to see if I set the minimum intensity too high. I suspect I’ll have a datapoint in a month or two…

Edit: Literally minutes after I finished writing and posting this, I felt an earthquake. A quick check suggested it was this M6.2 earthquake, but setting a minimum MMI of 3.5 did not pull it up. Neither did 3, or 2.5. Setting the MMI minimum to 2 resulted in the earthquake being displayed. So, I guess up high in my apartment building I can feel earthquakes with an MMI > 2. For what it’s worth, the M6.2 earthquake was nearly 700km away from me.

Here’s the (latest) code:

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