Umbria EarthQuake on 30th October 2016.¶
This notebook is to complement our report. Figures from our report were produced with this notebook.
In [1]:
from obspy import UTCDateTime
from obspy.clients.fdsn import Client
In [2]:
import folium
## This cell contains functions useful for plotting events and stations
def get_depth_color(depth):
colors = [
"#FF0000",
"#FF3300",
"#FF6600",
"#FFCC00",
"#CCFF00",
"#66FF00",
"#00FF00",
]
# <=20000 <=10000 <=5000 <=1000 <=500 <=100 >100
limits = [20000, 10000, 5000, 1000, 500, 100]
for i, val in enumerate(limits):
if depth and depth >= val:
return colors[i]
return colors[-1]
def plot_events(events, origin=[0, 0], zoom=2, comments="", carte=None, outfile=None):
"""Example of use : if you need to add events to an existing map :
m = plot_events(events_1, color="red")
plot_events(events_2, color="green", outfile="map.html", carte=m)"""
if carte is None:
carte = folium.Map(location=origin, zoom_start=zoom)
for event in events:
for origin, magnitude in zip(event.origins, event.magnitudes):
lat, lon, depth, mag = (
origin.latitude,
origin.longitude,
origin.depth,
magnitude.mag,
)
infos = "(%s %s) depth=%s m mag=%s (%s)" % (lat, lon, depth, mag, comments)
folium.CircleMarker(
location=[lat, lon],
radius=50 * 2 ** (mag) / 2**10,
tooltip=infos,
color=get_depth_color(depth),
fill=True,
fill_color="#FF8C00",
).add_to(carte)
if outfile:
carte.save(outfile)
# webbrowser.open(outfile, new=2, autoraise=True)
# time.sleep(1)
return carte
def plot_event_stations(stations, events, origin=[0, 0], zoom=2, comments=""):
map = folium.Map(location=origin, zoom_start=zoom)
for event in events:
for origin, magnitude in zip(event.origins, event.magnitudes):
lat, lon, depth, mag = (
origin.latitude,
origin.longitude,
origin.depth,
magnitude.mag,
)
infos = "(%s %s) depth=%s m mag=%s (%s)" % (lat, lon, depth, mag, comments)
folium.CircleMarker(
location=[lat, lon],
radius=50 * 2 ** (mag) / 2**10,
tooltip=infos,
color=get_depth_color(depth),
fill=True,
fill_color="#FF8C00",
).add_to(map)
for net, sta, lat, lon, elev in stations:
name = ".".join([net, sta])
infos = "%s (%s, %s) %s m" % (name, lat, lon, elev)
folium.features.RegularPolygonMarker(
location=[lat, lon],
tooltip=infos,
color="blue",
fill_color="#FF8C00",
number_of_sides=3,
radius=10,
fill_opacity=0.3,
rotation=30,
).add_to(map)
return map
In [3]:
# Defining a range to narrow down our search for the event.
maxlatitude = 43
minlatitude = 41
maxlongtitude = 14
minlongtitude=12
starttime = UTCDateTime('2016-10-30')
endtime = UTCDateTime('2016-11-1')
minmagnitude = 6
umbria_event = Client('INGV').get_events(starttime= starttime, endtime=endtime, minmagnitude=6,
minlatitude=minlatitude,
maxlatitude=maxlatitude,
minlongitude=minlongtitude,
maxlongitude=maxlongtitude)
assert len(umbria_event) == 1, 'Event not found'
eqo = umbria_event[0].origins[0]
print('Event Magnitude: ', umbria_event[0].magnitudes[0].mag)
print('Event Latitude: ', eqo.latitude)
print('Event Longitude: ', eqo.longitude)
print('Event Time: ', eqo.time)
plot_events(umbria_event, origin=[43,15], zoom=7)
Event Magnitude: 6.5 Event Latitude: 42.8303 Event Longitude: 13.1092 Event Time: 2016-10-30T06:40:17.320000Z
Out[3]:
Make this Notebook Trusted to load map: File -> Trust Notebook
Getting the stations¶
We have found the seismic event at Umbria, near the city of Norcia. Using the epicenter reported by INGV, we will look for stations in a 0.2 degrees of latitude range of the epicenter. This translates to about 22 kilometers in radius.
We want channels EH?, for high frequency sampling, high gain seismometers in any direction. Any network within the radius is acceptable.
In [4]:
inventory = Client('INGV').get_stations(
level='channel',
network='*',
starttime=starttime,
endtime=endtime,
channel='EH?',
latitude=eqo.latitude,
longitude=eqo.longitude,
maxradius=0.2)
print(inventory)
Inventory created at 2024-12-13T18:14:47.223000Z Created by: INGV-ONT WEB SERVICE: fdsnws-station | version: 1.1.61 /exist/apps/fdsn-station/fdsnws/station/1/query?starttime=2016-10-... Sending institution: eXistDB (INGV-ONT) Contains: Networks (3): 3A, 8P, IV Stations (15): 3A.MZ103 (Accumoli campo sportivo - ENEA) 3A.MZ104 (Tino - ENEA) 3A.MZ105 (San Giovanni - ENEA) 8P.T1202 (Villanova (RI)) 8P.T1212 (Avendita PG) 8P.T1213 (Savelli PG) 8P.T1214 (Forca Canapine AP) 8P.T1215 (Reggiano PG) 8P.T1216 (Castelvecchio PG) 8P.T1217 (Poggiodomo PG) 8P.T1218 (Civita PG) 8P.T1221 (Spina Nuova PG) 8P.T1244 (Spelonga) IV.MC2 (Monte Cornaccione) IV.MMO1 (Montemonaco) Channels (45): 3A.MZ103..EHZ, 3A.MZ103..EHN, 3A.MZ103..EHE, 3A.MZ104..EHZ, 3A.MZ104..EHN, 3A.MZ104..EHE, 3A.MZ105..EHZ, 3A.MZ105..EHN, 3A.MZ105..EHE, 8P.T1202..EHZ, 8P.T1202..EHN, 8P.T1202..EHE, 8P.T1212..EHZ, 8P.T1212..EHN, 8P.T1212..EHE, 8P.T1213..EHZ, 8P.T1213..EHN, 8P.T1213..EHE, 8P.T1214..EHZ, 8P.T1214..EHN, 8P.T1214..EHE, 8P.T1215..EHZ, 8P.T1215..EHN, 8P.T1215..EHE, 8P.T1216..EHZ, 8P.T1216..EHN, 8P.T1216..EHE, 8P.T1217..EHZ, 8P.T1217..EHN, 8P.T1217..EHE, 8P.T1218..EHZ, 8P.T1218..EHN, 8P.T1218..EHE, 8P.T1221..EHZ, 8P.T1221..EHN, 8P.T1221..EHE, 8P.T1244..EHZ, 8P.T1244..EHN, 8P.T1244..EHE, IV.MC2..EHZ, IV.MC2..EHN, IV.MC2..EHE, IV.MMO1..EHZ, IV.MMO1..EHN, IV.MMO1..EHE
In [5]:
networks = [network for network in inventory.networks]
stations = []
for network in networks:
for station in network.stations:
stations.append(
[network.code, station.code, station.latitude, station.longitude, station.elevation]
)
print(f'We have {len(stations)} in the defined radius')
We have 15 in the defined radius
Plotting the stations around the epicenter¶
We have customised the Folium code. Our epicenter is marked as the red circle. Seismic stations as the blue triangles.
In [6]:
plot = plot_event_stations(events=umbria_event, stations=stations, zoom=11, origin=[eqo.latitude, eqo.longitude])
plot
Out[6]:
Make this Notebook Trusted to load map: File -> Trust Notebook
Using MassDownloader to download wavefronts within a radius¶
This is preferable over the get_waveforms function as the data is written to my filesystem rather than memory.
Makes it much faster to apply filters to the stream in place, and restore the raw trace whenever needed.
CircularDomain is used in this case, as we are going with a radial search for stations. We limit our stations to the inventory defined in the cells above. Without the restriction, it should still be confined to the same inventory nonetheless.
We reject channels with gaps or overlaps. There ought to be at least 95% of the stream available within the requested window.
The files are stored in a new report directory, with subdirectories being created for the trace and station metadata.
In [7]:
from obspy.clients.fdsn.mass_downloader import CircularDomain
from obspy.clients.fdsn.mass_downloader import Restrictions
from obspy.clients.fdsn.mass_downloader import MassDownloader
domain = CircularDomain(longitude=eqo.longitude,
latitude=eqo.latitude,
minradius=0,
maxradius=0.2)
restrictions = Restrictions(starttime=UTCDateTime(eqo.time),
endtime=UTCDateTime(eqo.time) + 60*3,
reject_channels_with_gaps=True,
minimum_length=0.95,
network='*',
channel='EH?',
limit_stations_to_inventory=inventory
)
mdl = MassDownloader(providers=["INGV"])
mdl.download(
domain,
restrictions,
mseed_storage="./report/waveforms",
stationxml_storage="./report/stations",
)
[2024-12-13 19:19:10,179] - obspy.clients.fdsn.mass_downloader - INFO: Initializing FDSN client(s) for INGV. [2024-12-13 19:19:10,223] - obspy.clients.fdsn.mass_downloader - INFO: Successfully initialized 1 client(s): INGV. [2024-12-13 19:19:10,224] - obspy.clients.fdsn.mass_downloader - INFO: Total acquired or preexisting stations: 0 [2024-12-13 19:19:10,225] - obspy.clients.fdsn.mass_downloader - INFO: Client 'INGV' - Requesting unreliable availability. [2024-12-13 19:19:10,562] - obspy.clients.fdsn.mass_downloader - INFO: Client 'INGV' - Successfully requested availability (0.34 seconds) [2024-12-13 19:19:10,563] - obspy.clients.fdsn.mass_downloader - INFO: Client 'INGV' - Found 14 stations (42 channels). [2024-12-13 19:19:10,564] - obspy.clients.fdsn.mass_downloader - INFO: Client 'INGV' - Will attempt to download data from 14 stations. [2024-12-13 19:19:10,565] - obspy.clients.fdsn.mass_downloader - INFO: Client 'INGV' - Status for 33 time intervals/channels before downloading: EXISTS [2024-12-13 19:19:10,566] - obspy.clients.fdsn.mass_downloader - INFO: Client 'INGV' - Status for 9 time intervals/channels before downloading: NEEDS_DOWNLOADING [2024-12-13 19:19:10,697] - obspy.clients.fdsn.mass_downloader - INFO: Client 'INGV' - No data available for request. HTTP Status code: 204 [2024-12-13 19:19:10,698] - obspy.clients.fdsn.mass_downloader - INFO: Client 'INGV' - Launching basic QC checks... [2024-12-13 19:19:10,699] - obspy.clients.fdsn.mass_downloader - INFO: Client 'INGV' - Downloaded 0.0 MB [0.00 KB/sec] of data, 0.0 MB of which were discarded afterwards. [2024-12-13 19:19:10,699] - obspy.clients.fdsn.mass_downloader - INFO: Client 'INGV' - Status for 9 time intervals/channels after downloading: DOWNLOAD_FAILED [2024-12-13 19:19:10,700] - obspy.clients.fdsn.mass_downloader - INFO: Client 'INGV' - Status for 33 time intervals/channels after downloading: EXISTS [2024-12-13 19:19:10,712] - obspy.clients.fdsn.mass_downloader - INFO: Client 'INGV' - No station information to download. [2024-12-13 19:19:10,714] - obspy.clients.fdsn.mass_downloader - INFO: ============================== Final report [2024-12-13 19:19:10,715] - obspy.clients.fdsn.mass_downloader - INFO: 33 MiniSEED files [2.6 MB] already existed. [2024-12-13 19:19:10,715] - obspy.clients.fdsn.mass_downloader - INFO: 11 StationXML files [0.4 MB] already existed. [2024-12-13 19:19:10,716] - obspy.clients.fdsn.mass_downloader - INFO: Client 'INGV' - Acquired 0 MiniSEED files [0.0 MB]. [2024-12-13 19:19:10,716] - obspy.clients.fdsn.mass_downloader - INFO: Client 'INGV' - Acquired 0 StationXML files [0.0 MB]. [2024-12-13 19:19:10,716] - obspy.clients.fdsn.mass_downloader - INFO: Downloaded 0.0 MB in total.
Out[7]:
{'INGV': <obspy.clients.fdsn.mass_downloader.download_helpers.ClientDownloadHelper at 0x11d82ae40>}
Plotting the raw traces¶
Using the same functions from our TP, we plot the power spectral density.
We use the built in plot function from obspy to plot the seismic trace.
We do this for each channel in a for loop.
For the sake of brevity, we terminate early after the 3rd channel trace is plotted.
In [8]:
import os
from obspy import read_inventory, read
from tp_obspy_utils import get_periodogram
waveforms_directory = './report/waveforms'
for index, file in enumerate(os.listdir(waveforms_directory)):
if file.endswith('.mseed'):
st = read(os.path.join(waveforms_directory, file))
if index == 0:
fs = st[0].stats.sampling_rate
get_periodogram(st[0].data, fs, semilog=True, show=True)
break
# Limit ourselves to the first three channels.
if index == 3:
break
st.plot()
Plotting Traces with Band Pass Filter¶
Removing the instrument response would require information of the instruments used for the channel itself, which is codified in the station XML metadata file, which we have retrieved with MassDownloader.
We remove the instrument response, de-mean the trace and apply a band-pass filter which attenuates frequencies outside the range to a minimum.
Once again, we terminate after the 3rd plot for the reader's convenience.
In [9]:
waveforms_directory = './report/waveforms'
inventory_directory = './report/stations'
# Loop through all the files in the directory
for index, file in enumerate(os.listdir(waveforms_directory)):
# Check if the file is an mseed file
if file.endswith('.mseed'):
seed_code = file.split('..')[0]
# Read the mseed file using obspy
st = read(os.path.join(waveforms_directory, file))
sampling_rate = st[0].stats.sampling_rate
inv = read_inventory(os.path.join(inventory_directory, seed_code + '.xml'))
st.remove_response(inv)
st.detrend('demean')
st.filter("bandpass", freqmin= 0.1, freqmax=1)
if index == 0:
get_periodogram(st[0].data, sampling_rate, semilog=True, show=True)
# Limit ourselves to the first three channels.
if index == 3:
break
st.plot()
