Plotting Seismograms on the Station MapΒΆ
It may be useful to compare seismograms interactively on a map. Although gmprocess does not have this capability, it is relatively straightforward to do using folium for generating the map and altair for creating the figures.
Note that the altair package is needed for the interactive plotting, but is not used by gmprocess and will need to be installed separately. The package can be installed via conda or mamba.
We first take care of the imports.
import numpy as np
import pandas as pd
import folium
import altair as alt
import json
from gmprocess.io.asdf.stream_workspace import StreamWorkspace
from gmprocess.utils.constants import DATA_DIR
We are using a ground motion dataset of earthquakes from the 2014 Napa earthquake. The processed waveforms in the ASDF file will be used for this tutorial.
# Open ASDF file containing Nisqually dataset
data_path = DATA_DIR / 'demo' / 'nc72282711' / 'workspace.h5'
workspace = StreamWorkspace.open(data_path)
Lets first list the stations stored in the ASDF file
ds = workspace.dataset
station_list = ds.waveforms.list()
print(station_list)
['CE.68150']
As an illustration, we will plot the vertical component on a map for a few stations such that clicking on the station will popup a figure of the corresponding vertical component time-series.
event = workspace.get_event('nc72282711')
pstreams = workspace.get_streams('nc72282711', labels=['default'])
In order to handle lengthy time-series, we need disable the maximum row count of 5000 that altair uses to keep performance in check. There are alternative options that can be used to improve performance if we were working with a larger dataset. See the following documentation for more information:
alt.data_transformers.disable_max_rows()
DataTransformerRegistry.enable('default')
The code below will generate the map and could be customized as the user desires.
station_map = folium.Map(
location=[event.latitude, event.longitude], zoom_start=7, control_scale=True, tiles='cartodb positron'
)
lats = np.array([stream[0].stats.coordinates["latitude"] for stream in pstreams])
lons = np.array([stream[0].stats.coordinates["longitude"] for stream in pstreams])
stnames = np.array([stream[0].stats.station for stream in pstreams])
networks = np.array([stream[0].stats.network for stream in pstreams])
stations = pd.DataFrame(
{
"stnames": stnames,
"network": networks,
"coords": zip(lats, lons),
}
)
for i, r in stations.iterrows():
st = ds.waveforms[r["network"]+'.'+r["stnames"]]['nc72282711_default']
dd = pd.DataFrame({'times':st[2].times(),'Acceleration':st[2].data})
chart = alt.Chart(dd).mark_line().encode(
x='times',
y='Acceleration',
color=alt.value("#000000"))
popup = folium.Popup()
folium.features.VegaLite(chart, height=250, width=250).add_to(popup)
folium.CircleMarker(
location=r["coords"],
tooltip=r["stnames"],
popup=popup,
color='green',
fill=True,
radius=6,
).add_to(station_map)
Finally, we display the map. When a station icon is clicked on, the corresponding vertical time-series will popup in a window.
station_map