Pu’u O’o vent collapse and the first fissure appearance during the 2018 eruption#
In this exercise, we will plot:
a map of the Kilauea and the rift zone
locations of the eruption activities, and associated seismic events
snapshot of the seismic activity
%matplotlib inline
import cartopy.crs as ccrs
import cartopy.io.img_tiles as cimgt
import matplotlib.pyplot as plt
import numpy as np
from obspy import UTCDateTime
From the reports, populate the following details:
## Specify location of Pu'u O'o vent
puu_oo_lat = ??????
puu_oo_lon = ??????
## Specify location of first fissure eruption
puna_eruption_lat = ??????
puna_eruption_lon = ??????
## Specify start time and end time of the eruption series (as a string)
# Obspy scripts read data and time using function UTCDateTime.
# Seismogram recordings and event catalogues are always in UTC (universal) time.
# use format "YYYY-MM-DDTHH:mm:ss" e.g.: 2024-04-26T00:00:00
starttime = UTCDateTime(???????) # <--- pick time a day before Puuoo collapse
endtime = UTCDateTime(???????) # <-- pick time a day after the first fissure appears
# for plotting purposes:
# limit events occuring within the map boundary based on Pu'u O'o vent location
minlongitude = puu_oo_lon-0.2
maxlongitude = puu_oo_lon+0.3
minlatitude = puu_oo_lat-0.2
maxlatitude = puu_oo_lat+0.2
# pick the background image (either satellite image or street)
mapbox_satellite = cimgt.MapboxTiles(map_id='satellite',
access_token='pk.eyJ1IjoibG91aXNtb3Jlc2kiLCJhIjoiY2pzeG1mZzFqMG5sZDQ0czF5YzY1NmZ4cSJ9.lpsUzmLasydBlS0IOqe5JA')
mapbox_streets = cimgt.MapboxTiles(map_id='streets',
access_token='pk.eyJ1IjoibG91aXNtb3Jlc2kiLCJhIjoiY2pzeG1mZzFqMG5sZDQ0czF5YzY1NmZ4cSJ9.lpsUzmLasydBlS0IOqe5JA')
google_maps_street = cimgt.GoogleTiles(style="street")
google_maps_satellite = cimgt.GoogleTiles(style="satellite")
map_tiles = google_maps_satellite
Now run the following cell to plot the map. Can you recognize the location of the caldera and various vents along the rift zone?
# plot map
fig = plt.figure(dpi=200)
# set the projection and extent of map
projection=map_tiles.crs
ax = fig.add_subplot(111,projection=projection)
plot_extent = (minlongitude,maxlongitude,minlatitude,maxlatitude) #map boundary set by puu_oo_lat
ax.set_extent(plot_extent)
#plot satellite / street image
ax.add_image(map_tiles, 12)
# add gridlines to the plot
gl = ax.gridlines(draw_labels=True, alpha=0.2, linestyle='--')
gl.top_labels = None; gl.right_labels = None
# add title
ax.set_title('Kilauea')
# plot location of Pu'u O'o vent
x, y, _ = projection.transform_points(ccrs.Geodetic(), np.array(puu_oo_lon), np.array(puu_oo_lat)).T
ax.scatter(x, y, 30, color="#C68E17", marker="o", edgecolor=None, zorder=3)
plt.text(x, y, 'vent', ha="right",va="bottom", color="#d3d3d3", family="monospace", fontsize='small')
# plot location of 2018 Lower puna eruption
x, y, _ = projection.transform_points(ccrs.Geodetic(), np.array(puna_eruption_lon), np.array(puna_eruption_lat)).T
ax.scatter(x, y, 30, color="#E47200", marker="s", edgecolor=None, zorder=3)
plt.text(x, y, 'fissure \neruption', color="#d3d3d3", ha="right", va="bottom", family="monospace", fontsize='small')
plt.show()
Now we will get an inventory of the seismic stations operated by the Hawaiian Volcano Observatory during the time of eruption.
The specific obspy function to use is ‘get_stations’.
We will specifically plot HVO seismometers (with network code “HV”) that has a EHZ channel. EHZ is the channel name for the short period vertical component of the seismometer, which is more sensitive to the weak seismic events.
from obspy.clients.fdsn import Client
client = Client("IRIS") # we use IRIS as IRIS holds the seismic data and metadata (e.g., inventory info)
network="HV" # network code for Hawaiian Volcano Observatory
station="*" # * means include all stations
channel="EHZ" # component of the seismometer.
inv = client.get_stations(network=network,station="*",location="",
channel=channel,level="channel",
starttime=starttime, endtime=endtime)
print(inv)
Run the next cell to make a new map, which includes the seismic stations you just found.
# plot map with stations operated by HVO
fig = plt.figure(dpi=200)
# set the projection and extent of map
projection=map_tiles.crs
ax = fig.add_subplot(111,projection=projection)
plot_extent = (minlongitude,maxlongitude,minlatitude,maxlatitude) #map boundary set by puu_oo_lat
ax.set_extent(plot_extent)
#plot satellite / street image
ax.add_image(map_tiles, 12)
# add gridlines to the plot
gl = ax.gridlines(draw_labels=True, alpha=0.2, linestyle='--')
gl.top_labels = None; gl.right_labels = None
# add title
ax.set_title('Kilauea')
# plot location of Pu'u O'o vent
x, y, _ = projection.transform_points(ccrs.Geodetic(), np.array(puu_oo_lon), np.array(puu_oo_lat)).T
ax.scatter(x, y, 30, color="#C68E17", marker="o", edgecolor=None, zorder=3)
plt.text(x, y, 'vent', ha="right",va="bottom", color="#d3d3d3", family="monospace", fontsize='small')
# plot location of 2018 Lower puna eruption
x, y, _ = projection.transform_points(ccrs.Geodetic(), np.array(puna_eruption_lon), np.array(puna_eruption_lat)).T
ax.scatter(x, y, 30, color="#E47200", marker="s", edgecolor=None, zorder=3)
plt.text(x, y, 'fissure \neruption', color="#d3d3d3", ha="right", va="bottom", family="monospace", fontsize='small')
# HV stations
stations = []; stla = []; stlo = []
for sta in inv.networks[0].stations:
#print(sta.code)
stations.append(sta.code)
stla.append(sta.latitude)
stlo.append(sta.longitude)
x, y, _ = projection.transform_points(ccrs.Geodetic(), np.array(stlo), np.array(stla)).T
ax.scatter(x, y, 30, color="#D22B2B", marker="v", zorder=3)
for i in range(len(stations)):
t = plt.text(x[i], y[i], stations[i], color="#d3d3d3", va="top", family="monospace",fontsize='xx-small',clip_on=True)
t.clipbox = ax.bbox
plt.show()
We will now download a day-long time series from the station closest to the vent using the ‘get_waveforms’ function.
We will need two inputs:
station name (4 digit code)
a start time (for the time series). Pick a time ~12 hours before the first major activity at the vent (in UTC!)
If you get error messages saying FDSNNNoDataException, that means those stations are not operating (hence no data).
stationIN = ???????? # <--- input station name (as a string)
# specify time frame for downloading data
st_starttime = UTCDateTime(????????) # <---- adjust start time!
st_endtime = st_starttime + 24 * 60 * 60
stream = client.get_waveforms(network='HV',station = stationIN,
location='',channel='EHZ',
starttime=st_starttime, endtime=st_endtime)
print(stream)
There are many ways to plot the seismogram. One easy way to look a day-long time series is plotting a “day plot” using the obpys function. Run the following cell to see the plot:
stream.plot(type='dayplot',interval=60)
Combined with recordings from other stations, seismologists are able to detect many seismic activities and put together an event catalog detailing their magnitude, origin time, and location.
# We will download a catalog of events occurring near the eruption site before the vent collapse
# and shortly after first fissure appearance.
# The catalog is compiled by U.S. Geological Survey, using these HV stations.
client_catalog = Client("USGS") # switch client to USGS
cat = client_catalog.get_events(starttime = starttime, endtime = endtime,
minlatitude = minlatitude, maxlatitude = maxlatitude,
minlongitude = minlongitude, maxlongitude = maxlongitude)
print(cat.__str__(print_all=True))
We will make a final plot showing the events on the map (colored by time). Run the following cell:
from matplotlib import cm
from matplotlib.dates import DateFormatter,DayLocator
# plot map with seismic events
fig = plt.figure(dpi=200)
# set the projection and extent of map
projection=map_tiles.crs
ax = fig.add_subplot(111,projection=projection)
plot_extent = (minlongitude,maxlongitude,minlatitude,maxlatitude) #map boundary set by puu_oo_lat
ax.set_extent(plot_extent)
#plot satellite / street image
ax.add_image(map_tiles, 12)
# add gridlines to the plot
gl = ax.gridlines(draw_labels=True, alpha=0.2, linestyle='--')
gl.top_labels = None; gl.right_labels = None
# add title
ax.set_title(' Seismicity between %s and %s ' %(starttime.strftime('%b %d'), endtime.strftime('%b %d')),fontsize="medium")
# plot location of seismic events
otime_cb = []; evlat = []; evlon = []; evmag = []
for event in cat:
otime_cb.append(event.origins[0].time.matplotlib_date)
evlat.append(event.origins[0].latitude)
evlon.append(event.origins[0].longitude)
evmag.append(event.magnitudes[0].mag)
total_event = len(cat)
colormap = cm.get_cmap('magma',total_event)
x, y, _ = projection.transform_points(ccrs.Geodetic(), np.array(evlon), np.array(evlat)).T
im = ax.scatter(x, y, s=evmag, c=otime_cb, cmap=colormap, marker="o", edgecolor=None, zorder=3)
cbar = plt.colorbar(im,orientation='horizontal',shrink=0.5,
ticks=DayLocator(interval=1),
format=DateFormatter('%b %d'))
cbar.ax.tick_params(labelsize=6)
# plot location of Pu'u O'o vent
x, y, _ = projection.transform_points(ccrs.Geodetic(), np.array(puu_oo_lon), np.array(puu_oo_lat)).T
ax.scatter(x, y, 30, color="#C68E17", marker="o", edgecolor=None, zorder=3)
plt.text(x, y, 'vent', ha="right",va="bottom", color="#d3d3d3", family="monospace", fontsize='small')
# plot location of 2018 Lower puna eruption
x, y, _ = projection.transform_points(ccrs.Geodetic(), np.array(puna_eruption_lon), np.array(puna_eruption_lat)).T
ax.scatter(x, y, 30, color="#E47200", marker="s", edgecolor=None, zorder=3)
plt.text(x, y, 'fissure \neruption', color="#d3d3d3", ha="right", va="bottom", family="monospace", fontsize='small')
plt.show()
# save images for every hour:
kk = 0
for hrIN in np.arange(starttime,endtime,3*60*60):
print('Plotting ', hrIN)
# plot location of seismic events
otime_cbIN = []; evlatIN = []; evlonIN = []; evmagIN = []
for event in cat:
if event.origins[0].time >= starttime and event.origins[0].time < hrIN:
otime_cbIN.append(event.origins[0].time.matplotlib_date)
evlatIN.append(event.origins[0].latitude)
evlonIN.append(event.origins[0].longitude)
evmagIN.append(event.magnitudes[0].mag)
fig = plt.figure(dpi=200)
# set the projection and extent of map
projection=map_tiles.crs
ax = fig.add_subplot(111,projection=projection)
plot_extent = (minlongitude,maxlongitude,minlatitude,maxlatitude) #map boundary set by puu_oo_lat
ax.set_extent(plot_extent)
#plot satellite / street image
ax.add_image(map_tiles, 12)
# add gridlines to the plot
gl = ax.gridlines(draw_labels=True, alpha=0.2, linestyle='--')
gl.top_labels = None; gl.right_labels = None
# add title
ax.set_title('From Pu\'u \'O\'o vent collapse to first fissure eruption',fontsize="medium")
# plot location of Pu'u O'o vent
x, y, _ = projection.transform_points(ccrs.Geodetic(), np.array(puu_oo_lon), np.array(puu_oo_lat)).T
ax.scatter(x, y, 30, color="#C68E17", marker="o", edgecolor=None, zorder=3)
plt.text(x, y, 'vent', ha="right",va="bottom", color="#d3d3d3", family="monospace", fontsize='small')
# plot location of 2018 Lower puna eruption
x, y, _ = projection.transform_points(ccrs.Geodetic(), np.array(puna_eruption_lon), np.array(puna_eruption_lat)).T
ax.scatter(x, y, 30, color="#E47200", marker="s", edgecolor=None, zorder=3)
plt.text(x, y, 'fissure \neruption', color="#d3d3d3", ha="right", va="bottom", family="monospace", fontsize='small')
# plot events:
colormap = cm.get_cmap('magma',total_event)
x, y, _ = projection.transform_points(ccrs.Geodetic(), np.array(evlonIN), np.array(evlatIN)).T
im = ax.scatter(x, y, s=evmagIN, c=otime_cbIN, cmap=colormap, marker="o", edgecolor=None, zorder=3)
## Add time stamp to the top left corner
plt.text(0.05, 0.95, hrIN.ctime(), color="#d3d3d3",transform=ax.transAxes,fontsize=8, verticalalignment='top')
fileOUT = './Puuoo_eq_%03d.png' %(kk)
kk +=1
plt.savefig(fileOUT,dpi=250)
plt.close()
Now view the eruption as an animation!
from PIL import Image
import glob
# Create the frames
frames = []
imgs = glob.glob('./Puuoo_eq_*.png')
imgs.sort()
for i in imgs:
new_frame = Image.open(i)
frames.append(new_frame)
# Save into a GIF file that loops once!
frames[0].save('./Pu_uo_o_eruption.gif',
format='GIF',
append_images=frames[1:],
save_all=True,
duration=2e-6, loop=1)