<jupyter_start><jupyter_text>**QuakeMAP**The following python script is used to process and visualize seismic event and seismic station data in the Kiskatinaw Seismic Monitoring and Mitigation Area (KSMMA), north-east British Columbia, Canada from 2014-2018.It is also intended to highlight the functionality of Pandas, GeoPandas and Bokeh for processing and visualizing open source georeferenced data using open source software, in particular for producing an interactive time series map for explaratory analysis.This python script will be updated with seismic event data from 2019-2020 once it is publicly available. **Data** KSMMA boundary data from the British Columbia Oil and Gas Commission (BCOGC) (open source)- ftp://ftp.bcogc.ca/outgoing/Induced%20Seismicity%20in%20KSMMA_Report%20Appendices/Appendix%20D%20-%20Mapping%20and%20Pressure%20Data%20Files/Shapefiles/General/ Seismic event data from the Geological Survey of Canada (GSC) (open source)- https://geoscan.nrcan.gc.ca/starweb/geoscan/servlet.starweb?path=geoscan/fulle.web&search1=R=306292 (2014-2016)- https://geoscan.nrcan.gc.ca/starweb/geoscan/servlet.starweb?path=geoscan/fulle.web&search1=R=326015 (2017-2018) Seismic station data from the Incoporated Research Institutions for Seismology (IRIS) (open source)- https://ds.iris.edu/gmap/network=*&starttime=2014-01-01&maxlat=56.4811&maxlon=-119.9268&minlat=55.6528&minlon=-121.2608&drawingmode=box&planet=earth (note to only select stations where seismic data is available (A), or partially available (P), and not restricted (R), view 'more information' on the station icon. This also highlights which stations are co-located (broadband and accelerometer rather than just broadband) by the presence of a Titan instrument. For this data, this includes stations: BCH1B, BCH2B, MONT1, MONT2, MONT3, MONT6, MONT8, MONT9, MONTA, KSM01, KSM10) **Packages**- jupyterlab 2.2.9- pandas 1.1.3- xlrd 1.2.0- geopandas 0.8.1- bokeh 2.2.3<jupyter_code>import os
import pandas as pd
import geopandas as gpd
import numpy as np
from bokeh.models import ColumnDataSource, CDSView, GroupFilter, Title, Range1d, WMTSTileSource, Div, Slider, CustomJS, HoverTool
from bokeh.plotting import output_file, figure, show
from bokeh.layouts import gridplot<jupyter_output><empty_output><jupyter_text>**Data pre-processing**<jupyter_code>#############################################################################################################
## pre-process data
#############################################################################################################

## change cwd to folder containing 'QuakeMAP.ipynb' and 'data'
os.chdir("/home/tom/python/QuakeMAP")
#print(os.getcwd())

## create tuple for cleaning datetime columns for consistency
non_numeric = ("`", "¬", "!", "\"", "£", "$", "%", "^", "&", "*", "(", ")", "-", "_", "=", "+",
               "[", "{", "]", "}",
               ";", ":", "'", "@", "#", "~",
               "\\", "|", ",", "<", ".", ">", "/", "?",
               " ",
               "a", "b", "c", "d", "e", "f", "g", "h", "i", "j", "k", "l", "m", "n", "o", "p", "q", "r", "s",
               "t", "u", "v", "w", "x", "y", "z",
               "A", "B", "C", "D", "E", "F", "G", "H", "I", "J", "K", "L", "M", "N", "O", "P", "Q", "R", "S",
               "T", "U", "V", "W", "X", "Y", "Z")

##################################################
## seismic events
##################################################

## read seismic event data into pd's
pd1 = pd.read_excel("./data/of_8335_Appendix1.xlsx").astype(str)
pd2 = pd.read_excel("./data/of_8718_Appendix_1.xlsx", header = 1).astype(str)
#print(pd1.head())
#print(pd2.head())

## format data of interest
pd1["datetime"] = pd1["Y-M-D"] + " " + pd1["H:M:S"].str[:8]
pd1["lat"] = pd1["Lat. (N)"]
pd1["lon"] = pd1["Lon. (E)"]
pd1["mag"] = pd1["ML"]
pd1["depth_km"] = pd1["Depth"]
pd1["origin"] = "GSC"
pd2["datetime"] = pd2["datetime"].str[:19]
pd2["depth_km"] = pd2["depth"]
pd2["origin"] = "GSC"

## copy columns from each pd to a list and concatenate list to new pd
l12 = [pd1[["datetime", "lat", "lon", "mag", "depth_km", "origin"]].copy(),
       pd2[["datetime", "lat", "lon", "mag", "depth_km", "origin"]].copy()]
pd_SE = pd.concat(l12)

## format datetime and convert back to string
for char in non_numeric:
    pd_SE["datetime"] = pd_SE["datetime"].str.replace(char, "")
pd_SE["datetime"] = pd.to_datetime(pd_SE["datetime"], format = "%Y%m%d%H%M%S")
pd_SE.sort_values(by = ["datetime"], inplace = True, ascending = True)
pd_SE["datetime"] = pd_SE["datetime"].astype(str)

## format numeric data for consistency
pd_SE[["lat","lon","mag","depth_km"]] = pd_SE[["lat","lon","mag","depth_km"]].apply(pd.to_numeric)
pd_SE["lat"] = [f"{lat:.6f}" for lat in pd_SE["lat"]]
pd_SE["lon"] = [f"{lon:.6f}" for lon in pd_SE["lon"]]
pd_SE["mag"] = [f"{mag:.2f}" for mag in pd_SE["mag"]]
pd_SE["depth_km"] = [f"{depth_km:.2f}" for depth_km in pd_SE["depth_km"]]

## reset index and print key info
pd_SE.reset_index(drop = True, inplace = True)
print(pd_SE["datetime"].count(), "seismic events")
print(pd_SE["datetime"].min())
print(pd_SE["datetime"].max())
print(pd_SE.head())

## [OPTIONAL] write to csv
pd_SE2014 = pd_SE[pd_SE["datetime"].str.match("2014")]
pd_SE2015 = pd_SE[pd_SE["datetime"].str.match("2015")]
pd_SE2016 = pd_SE[pd_SE["datetime"].str.match("2016")]
pd_SE2017 = pd_SE[pd_SE["datetime"].str.match("2017")]
pd_SE2018 = pd_SE[pd_SE["datetime"].str.match("2018")]
l_name = ["SE2014_2018", "SE2014", "SE2015", "SE2016", "SE2017", "SE2018"]
l_year = [pd_SE, pd_SE2014, pd_SE2015, pd_SE2016, pd_SE2017, pd_SE2018]
for i in range(0, len(l_year)):
    l_year[i].to_csv(l_name[i] + ".csv", index = False)
    
##################################################
## seismic stations
##################################################

## read seismic station data into pd
pd_SS = pd.read_csv("./data/gmap-stations.txt", sep = "|", skiprows = [2]).astype(str)
#print(pd_SS.head())

## format data of interest
pd_SS.columns = pd_SS.columns.str.strip()
pd_SS = pd_SS.rename(columns = {"#Network":"network",
                                "Station":"station",
                                "Latitude":"lat",
                                "Longitude":"lon",
                                "Elevation":"elev_m",
                                "StartTime":"starttime",
                                "EndTime":"endtime"})
## only include stations with available (A) or partially available (P) seismic data (see website)
pd_SS = pd_SS[pd_SS.network != "XL"]
## add column for co-located stations (see website)
pd_SS.set_index(pd_SS.station, inplace = True)
pd_SS.loc[["BCH1B","BCH2B", "MONT1", "MONT2", "MONT3", "MONT6", "MONT8", "MONT9", "MONTA", "KSM01", "KSM10"],
          "co_loc"] = 1
pd_SS.co_loc.fillna(value = 0, inplace = True)
pd_SS["origin"] = "IRIS"

## format start time and convert back to string
for char in non_numeric:
    pd_SS["starttime"] = pd_SS["starttime"].str.replace(char, "")
    pd_SS["endtime"] = pd_SS["endtime"].str.replace(char, "")
pd_SS["starttime"] = pd.to_datetime(pd_SS["starttime"], format = "%Y%m%d%H%M%S")
pd_SS.sort_values(by = ["starttime"], inplace = True, ascending = True)
pd_SS["starttime"] = pd_SS["starttime"].astype(str)
## end time has dates out of bounds so we must split and concatenate to parse desired format
pd_SS["year"] = pd_SS["endtime"].str[:4]
pd_SS["month"] = pd_SS["endtime"].str[4:]
pd_SS["month"] = pd.to_datetime(pd_SS["month"], format = "%m%d%H%M%S")
pd_SS["endtime"] = pd_SS["year"] + pd_SS["month"].astype(str).str[4:]
pd_SS = pd_SS.drop(["Sitename", "year", "month"], axis = 1)

## format numeric data for consistency
pd_SS[["lat","lon","elev_m","co_loc"]] = pd_SS[["lat", "lon", "elev_m", "co_loc"]].apply(pd.to_numeric)
pd_SS["lat"] = [f"{lat:.6f}" for lat in pd_SS["lat"]]
pd_SS["lon"] = [f"{lon:.6f}" for lon in pd_SS["lon"]]
pd_SS["elev_m"] = [f"{elev_m:.0f}" for elev_m in pd_SS["elev_m"]]
pd_SS["co_loc"] = [f"{co_loc:.0f}" for co_loc in pd_SS["co_loc"]]

## reset index and print key info
pd_SS.reset_index(drop = True, inplace = True)
print(pd_SS["station"].count(), "seismic stations")
print(pd_SS["starttime"].min())
print(pd_SS["starttime"].max())
print(pd_SS.head())

## [OPTIONAL] write to csv
pd_SS.to_csv("SS2014_2020.csv", index = False)<jupyter_output><empty_output><jupyter_text>**Data visualisation**<jupyter_code>#############################################################################################################
## prepare data
#############################################################################################################

##################################################
## KSMMA boundary
##################################################

## load shapefile as gpd, simplify columns to geometry only, change crs to Web Mercator for mapping in Bokeh
gpd_KSMMA = gpd.read_file("./data/KSMMA Outline_plygn.shp")
gpd_KSMMA.drop(gpd_KSMMA.iloc[:, 0:8], inplace = True, axis = 1)
gpd_KSMMA_WM = gpd_KSMMA.to_crs("epsg:3857")

## extract max extent for plot
pd_bounds = gpd_KSMMA_WM.bounds.round(-3)
minx = pd_bounds.iloc[0]["minx"]
miny = pd_bounds.iloc[0]["miny"]
maxx = pd_bounds.iloc[0]["maxx"]
maxy = pd_bounds.iloc[0]["maxy"]

## extract x and y to new columns, remove geometry column (note this does not change gpd to pd)
coordinates_array = np.asarray(gpd_KSMMA_WM.geometry.iloc[0].exterior.coords.xy).tolist()
gpd_KSMMA_WM["x"] = [coordinates_array[0]]
gpd_KSMMA_WM["y"] = [coordinates_array[1]]
gpd_KSMMA_WM = gpd_KSMMA_WM.drop("geometry", axis = 1).copy()

## view gpd
#print(gpd_KSMMA_WM.head())

## create cds for KSMMA boundary data
cds_KSMMA_WM = ColumnDataSource(gpd_KSMMA_WM)

##################################################
## seismic events
##################################################

## transform pd to gpd (set current crs in order to change it to Web Mercator for mapping in Bokeh)
gpd_SE = gpd.GeoDataFrame(pd_SE, crs = "epsg:4326", geometry = gpd.points_from_xy(pd_SE.lon, pd_SE.lat))
gpd_SE_WM = gpd_SE.to_crs("epsg:3857")

## [OPTIONAL] restrict seismic event data to KSMMA boundary (switch off to include all seismic event data)
gpd_SE_WM = gpd_SE_WM.cx[minx:maxx, miny:maxy]

## extract x and y to new columns, remove geometry column (note this does not change gpd to pd)
gpd_SE_WM["x"] = gpd_SE_WM.centroid.x
gpd_SE_WM["y"] = gpd_SE_WM.centroid.y
gpd_SE_WM = gpd_SE_WM.drop("geometry", axis = 1).copy()

## create new columns to use slider
gpd_SE_WM["month"] = gpd_SE_WM["datetime"].str[:7]
gpd_SE_WM["month_count"] = pd.factorize(gpd_SE_WM["datetime"].str[:7])[0]
## group by BCOGC operating thresholds
gpd_SE_WM["mag"] = gpd_SE_WM["mag"].apply(pd.to_numeric)
conditions = [(gpd_SE_WM["mag"] < 1.5),
              (gpd_SE_WM["mag"] >= 1.5) & (gpd_SE_WM["mag"] < 3.0),
              (gpd_SE_WM["mag"] > 3.0)]
values = ["0", "1", "2"]
gpd_SE_WM["threshold"] = np.select(conditions, values)
## visualise events by magnitude
gpd_SE_WM["mag_size"] = gpd_SE_WM["mag"] / 0.003
gpd_SE_WM["mag"] = [f"{mag:.2f}" for mag in gpd_SE_WM["mag"]]
## plot depth
gpd_SE_WM["depth_km"] = gpd_SE_WM["depth_km"].apply(pd.to_numeric)

## view gpd
#print(gpd_SE_WM.head())

## create an array of unique date values
dates = gpd_SE_WM.month.unique()

## create empty lists for each unique date value
xs = [[] for _ in range(len(dates))]
ys = [[] for _ in range(len(dates))]
datetimes = [[] for _ in range(len(dates))]
latitudes = [[] for _ in range(len(dates))]
longitudes = [[] for _ in range(len(dates))]
magnitudes = [[] for _ in range(len(dates))]
depth_kms = [[] for _ in range(len(dates))]
origins = [[] for _ in range(len(dates))]
thresholds = [[] for _ in range(len(dates))]
mag_sizes = [[] for _ in range(len(dates))]

## create a list of the number of unique date values
date = [i for i in range(len(dates))]

## append monthly data to each empty list
for _, row in gpd_SE_WM.iterrows():
    i = int(row.month_count)
    xi = row.x
    yi = row.y
    datetimei = row.datetime
    latitudei = row.lat
    longitudei = row.lon
    magnitudei = row.mag
    depth_kmi = row.depth_km
    origini = row.origin
    thresholdi = row.threshold
    mag_sizei = row.mag_size
    xs[i].append(xi)
    ys[i].append(yi)
    datetimes[i].append(datetimei)
    latitudes[i].append(latitudei)
    longitudes[i].append(longitudei)
    magnitudes[i].append(magnitudei)
    depth_kms[i].append(depth_kmi)
    origins[i].append(origini)
    thresholds[i].append(thresholdi)
    mag_sizes[i].append(mag_sizei)

## create variables for data in the first month
x = xs[0]
y = ys[0]
datetime = datetimes[0]
latitude = latitudes[0]
longitude = longitudes[0]
magnitude = magnitudes[0]
depth_km = depth_kms[0]
origin = origins[0]
threshold = thresholds[0]
mag_size = mag_sizes[0]

## cds columns must be same length, so we select the month with maximum number of events
month_max = gpd_SE_WM.groupby(["month_count"]).size().reset_index(name = "counts")
## to determine coefficients a and b
if month_max.counts.max() >= len(dates):
    a = month_max.counts.max() - month_max.counts.values[0]
else:
    a = len(dates) - month_max.counts.values[0]
b = month_max.counts.max() - len(dates)
x += [[]] * a
y += [[]] * a
datetime += [[]] * a
latitude += [[]] * a
longitude += [[]] * a
magnitude += [[]] * a
depth_km += [[]] * a
origin += [[]] * a
threshold += [[]] * a
mag_size += [[]] * a
xs += [[]] * b
ys += [[]] * b
datetimes += [[]] * b
latitudes += [[]] * b
longitudes += [[]] * b
magnitudes += [[]] * b
depth_kms += [[]] * b
origins += [[]] * b
thresholds += [[]] * b
mag_sizes += [[]] * b

## create cds for seismic event data
cds_SE_WM = ColumnDataSource(data = dict(x = x, y = y, xs = xs, ys = ys,
                                         datetime = datetime, datetimes = datetimes,
                                         latitude = latitude, latitudes = latitudes,
                                         longitude = longitude, longitudes = longitudes,
                                         magnitude = magnitude, magnitudes = magnitudes,
                                         depth_km = depth_km, depth_kms = depth_kms,
                                         origin = origin, origins = origins,
                                         threshold = threshold, thresholds = thresholds,
                                         mag_size = mag_size, mag_sizes = mag_sizes))

## filter cds by BCOGC operating thresholds
vcds_SE_WM_s = CDSView(source = cds_SE_WM, filters = [GroupFilter(column_name = "threshold", group = "0")])
vcds_SE_WM_m = CDSView(source = cds_SE_WM, filters = [GroupFilter(column_name = "threshold", group = "1")])
vcds_SE_WM_l = CDSView(source = cds_SE_WM, filters = [GroupFilter(column_name = "threshold", group = "2")])

##################################################
## seismic stations
##################################################

## transform pd to gpd (set current crs in order to change it to Web Mercator for mapping in Bokeh)
gpd_SS = gpd.GeoDataFrame(pd_SS, crs = "epsg:4326", geometry = gpd.points_from_xy(pd_SS.lon, pd_SS.lat))
gpd_SS_WM = gpd_SS.to_crs("epsg:3857")

## extract x and y to new columns, remove geometry column (note this does not change gpd to pd)
gpd_SS_WM["x"] = gpd_SS_WM.centroid.x
gpd_SS_WM["y"] = gpd_SS_WM.centroid.y
gpd_SS_WM = gpd_SS_WM.drop("geometry", axis = 1).copy()

## create new columns to use with the same slider we use for seismic events
gpd_SS_WM["month"] = gpd_SS_WM["starttime"].str[:7]

## month count to use with slider must match month count from seismic event data
gpd_SS_WM = gpd_SS_WM.merge(gpd_SE_WM[["month", "month_count"]].drop_duplicates(),
                            left_on="month",
                            right_on="month",
                            how="left")
## if seismic stations were installed prior to the first seismic event we set missing values to zero
gpd_SS_WM.loc[:gpd_SS_WM["month_count"].first_valid_index()-1,"month_count"] = 0
## if missing values exist in the middle of the dataset we backfill
gpd_SS_WM["month_count"].fillna(method="bfill", inplace=True)
## if seismic stations were installed after the last seismic event we remove data with missing values
gpd_SS_WM = gpd_SS_WM.dropna()
## set month count values as integer
gpd_SS_WM["month_count"] = gpd_SS_WM["month_count"].astype(int)
## create column with values indicating the number of months each seismic station was present
gpd_SS_WM["n"] = len(dates) - gpd_SS_WM["month_count"]
## replicate each row by n
gpd_SS_WM = gpd_SS_WM.loc[gpd_SS_WM.index.repeat(gpd_SS_WM.n)].copy()
## create column with values indicating the number of times each row has been repeated
gpd_SS_WM["repeat_id"] = gpd_SS_WM.groupby(level=0).cumcount()
## reset index to avoid duplicate values
gpd_SS_WM = gpd_SS_WM.reset_index(drop=True)
## replace monthly count column with values indicating the month each seismic station was present
gpd_SS_WM["month_count"] = gpd_SS_WM["month_count"] + gpd_SS_WM["repeat_id"]
## drop unecessary columns
gpd_SS_WM = gpd_SS_WM.drop(["n", "repeat_id"], axis=1).copy()

## view gpd
#print(gpd_SS_WM.head())

## create empty lists for each unique date value (seismic event date values)
xs2 = [[] for _ in range(len(dates))]
ys2 = [[] for _ in range(len(dates))]
networks2 = [[] for _ in range(len(dates))]
stations2 = [[] for _ in range(len(dates))]
latitudes2 = [[] for _ in range(len(dates))]
longitudes2 = [[] for _ in range(len(dates))]
elevation_ms2 = [[] for _ in range(len(dates))]
starttimes2 = [[] for _ in range(len(dates))]
endtimes2 = [[] for _ in range(len(dates))]
co_locs2 = [[] for _ in range(len(dates))]
origins2 = [[] for _ in range(len(dates))]

## append monthly data to each empty list
for _, row in gpd_SS_WM.iterrows():
    i = int(row.month_count)
    x2i = row.x
    y2i = row.y
    network2i = row.network
    station2i = row.station
    latitude2i = row.lat
    longitude2i = row.lon
    elevation_m2i = row.elev_m
    starttime2i = row.starttime
    endtime2i = row.endtime
    co_loc2i = row.co_loc
    origin2i = row.origin
    xs2[i].append(x2i)
    ys2[i].append(y2i)
    networks2[i].append(network2i)
    stations2[i].append(station2i)
    latitudes2[i].append(latitude2i)
    longitudes2[i].append(longitude2i)
    elevation_ms2[i].append(elevation_m2i)
    starttimes2[i].append(starttime2i)
    endtimes2[i].append(endtime2i)
    co_locs2[i].append(co_loc2i)
    origins2[i].append(origin2i)

## create variables for data in the first month
x2 = xs2[0]
y2 = ys2[0]
network2 = networks2[0]
station2 = stations2[0]
latitude2 = latitudes2[0]
longitude2 = longitudes2[0]
elevation_m2 = elevation_ms2[0]
starttime2 = starttimes2[0]
endtime2 = endtimes2[0]
co_loc2 = co_locs2[0]
origin2 = origins2[0]

## cds columns must be same length, so we select the month with maximum number of events
month_max = gpd_SS_WM.groupby(["month_count"]).size().reset_index(name = "counts")
## to determine coefficients a2 and b2
if month_max.counts.max() >= len(dates):
    a2 = month_max.counts.max() - month_max.counts.values[0]
else:
    a2 = len(dates) - month_max.counts.values[0]
b2 = month_max.counts.max() - len(dates)
x2 += [[]] * a2
y2 += [[]] * a2
network2 += [[]] * a2
station2 += [[]] * a2
latitude2 += [[]] * a2
longitude2 += [[]] * a2
elevation_m2 += [[]] * a2
starttime2 += [[]] * a2
endtime2 += [[]] * a2
co_loc2 += [[]] * a2
origin2 += [[]] * a2
xs2 += [[]] * b2
ys2 += [[]] * b2
networks2 += [[]] * b2
stations2 += [[]] * b2
latitudes2 += [[]] * b2
longitudes2 += [[]] * b2
elevation_ms2 += [[]] * b2
starttimes2 += [[]] * b2
endtimes2 += [[]] * b2
co_locs2 += [[]] * b2
origins2 += [[]] * b2

## create cds for seismic event data
cds_SS_WM = ColumnDataSource(data = dict(x2 = x2, y2 = y2, xs2 = xs2, ys2 = ys2,
                                         network2 = network2, networks2 = networks2,
                                         station2 = station2, stations2 = stations2,
                                         latitude2 = latitude2, latitudes2 = latitudes2,
                                         longitude2 = longitude2, longitudes2 = longitudes2,
                                         elevation_m2 = elevation_m2, elevation_ms2 = elevation_ms2,
                                         starttime2 = starttime2, starttimes2 = starttimes2,
                                         endtime2 = endtime2, endtimes2 = endtimes2,
                                         co_loc2 = co_loc2, co_locs2 = co_locs2,
                                         origin2 = origin2, origins2 = origins2))

## filter cds by co-located stations
vcds_SS_WM_0 = CDSView(source = cds_SS_WM, filters = [GroupFilter(column_name = "co_loc2", group = "0")])
vcds_SS_WM_1 = CDSView(source = cds_SS_WM, filters = [GroupFilter(column_name = "co_loc2", group = "1")])

#############################################################################################################
## output to standalone HTML file
#############################################################################################################

output_file("QuakeMAP.html")

#############################################################################################################
## create plots
#############################################################################################################

## select tools for Bokeh toolbar
tools = "pan,wheel_zoom,reset"

plot1 = figure(title = "QUAKE MAP",
               tools = tools,
               toolbar_location = "right",
               active_scroll = "wheel_zoom",
               plot_width = 480,
               plot_height = 480,
               x_axis_type = None,
               ## show Web Mercator as latitude
               y_axis_type = "mercator",
               y_axis_label = "Latitude",
               ## confine range to KSMMA boundary plus or minus buffer
               x_range = (minx - 6000, maxx + 6000),
               y_range = (miny - 8000, maxy + 29000),
               background_fill_color = "black",
               border_fill_color = "black",
               outline_line_color = "black",
               ## adjust margin to avoid white lines within gridplot
               margin = (0, - 1, - 1, 0))
plot1.ygrid.grid_line_color = None
plot1.yaxis.axis_label_text_color = "white"
plot1.yaxis.axis_label_text_font_style = "italic"
plot1.yaxis.axis_label_text_alpha = 0.5
plot1.yaxis.axis_label_text_font_size = "13px"
plot1.yaxis.major_label_text_color = "white"
plot1.yaxis.major_label_text_font_style = "normal"
plot1.yaxis.major_label_text_alpha = 0.5
plot1.yaxis.major_label_text_font_size = "11px"
plot1.title.text_color = "white"
plot1.title.text_font_style = "bold"
plot1.title.text_alpha = 0.5
plot1.title.text_font_size = "40px"
## add additional title
plot1.add_layout(
    Title(
        text = "Kiskatinaw Seismic Monitoring and Mitigation Area",
        text_color = "white",
        text_font_style = "bold",
        text_alpha = 0.5,
        text_font_size = "11px",
        align = "center",
        text_line_height = 0),
    'above')

plot2 = figure(toolbar_location = None,
               plot_width = 455,
               plot_height = 124,
               ## show Web Mercator as longitude
               x_axis_type = "mercator",
               x_axis_label = "Longitude",
               y_axis_label = "Depth (km)",
               ## confine range to match plot 1
               x_range = plot1.x_range,
               background_fill_color = "black",
               border_fill_color = "black",
               outline_line_color = "black",
               ## increase border between plot 1 and 2
               min_border_top = 30)
plot2.xgrid.grid_line_color = None
plot2.ygrid.grid_line_color = "white"
plot2.ygrid.grid_line_alpha = 0.5
plot2.ygrid.grid_line_dash = [6, 4]
plot2.axis.axis_label_text_color = "white"
plot2.axis.axis_label_text_font_style = "italic"
plot2.axis.axis_label_text_alpha = 0.5
plot2.axis.axis_label_text_font_size = "13px"
plot2.yaxis.axis_label_standoff = 13
plot2.axis.major_label_text_color = "white"
plot2.axis.major_label_text_font_style = "normal"
plot2.axis.major_label_text_alpha = 0.5
plot2.axis.major_label_text_font_size = "11px"
## flip direction of y axis
plot2.y_range.flipped = True
plot2.y_range = Range1d(15, 0)
plot2.yaxis.ticker = [0, 5, 10, 15]

plot3 = figure(toolbar_location = None,
               plot_width = 85,
               plot_height = 479,
               y_axis_type = None,
               x_axis_label = "Depth (km)",
               x_axis_location = "above",
               ## confine range to match plot 1
               y_range = plot1.y_range,
               background_fill_color = "black",
               border_fill_color = "black",
               outline_line_color = "black",
               ## increase border between plot 1 and 3
               min_border_left = 10,
               min_border_right = 20)
plot3.xgrid.grid_line_color = "white"
plot3.xgrid.grid_line_alpha = 0.5
plot3.xgrid.grid_line_dash = [6, 4]
plot3.xaxis.axis_label_text_color = "white"
plot3.xaxis.axis_label_text_font_style = "italic"
plot3.xaxis.axis_label_text_alpha = 0.5
plot3.xaxis.axis_label_text_font_size = "13px"
plot3.xaxis.major_label_text_color = "white"
plot3.xaxis.major_label_text_font_style = "normal"
plot3.xaxis.major_label_text_alpha = 0.5
plot3.xaxis.major_label_text_font_size = "11px"
plot3.x_range = Range1d(0, 15)
plot3.xaxis.ticker = [0, 5, 10, 15]

#############################################################################################################
## add glyphs
#############################################################################################################

##################################################
## background map
##################################################

plot1.add_tile(
    WMTSTileSource(
        url = "https://abcd.basemaps.cartocdn.com/dark_all/{z}/{x}/{y}{@2x}.png",
        attribution =
        """
        <span style="color:#808080;background-color:black;padding:10px;margin:-2px">
        &copy <a style="color:#808080;text-decoration:none" href="http://www.openstreetmap.org/copyright">
        OpenStreetMap contributors </a/span> | <span>
        &copy <a style="color:#808080;text-decoration:none" href="https://carto.com/attributions">
        CARTO </a/span> | <span>
        &copy <a style="color:#808080;text-decoration:none" href="https://github.com/TomSwinscoe">
        created by <NAME> </a/span>
        """
    )
)

##################################################
## KSMMA boundary
##################################################

plot1.patches("x", "y", source = cds_KSMMA_WM,
              fill_color = None, line_color = "white", line_width = 1, line_alpha = 0.6, legend_label = "KSMMA")

##################################################
## seismic events
##################################################

plot1.circle("x", "y", source = cds_SE_WM,
             fill_color = "green", fill_alpha = 0.5, line_color = "green", line_width = 0.5, line_alpha = 1,
             radius = "mag_size", legend_label = "<1.5 ML", name = "s", view = vcds_SE_WM_s)
plot1.circle("x", "y", source = cds_SE_WM,
             fill_color = "orange", fill_alpha = 0.5, line_color = "orange", line_width = 0.5, line_alpha = 1,
             radius = "mag_size", legend_label = ">=1.5 - <3.0 ML", name = "m", view = vcds_SE_WM_m)
plot1.circle("x", "y", source = cds_SE_WM,
             fill_color = "red", fill_alpha = 0.5, line_color = "red", line_width = 0.5, line_alpha = 1,
             radius = "mag_size", legend_label = ">=3.0 ML", name = "l", view = vcds_SE_WM_l)

plot2.circle("x", "depth_km", source = cds_SE_WM,
             fill_color = "green", fill_alpha = 0.5, line_color = "green", line_width = 0.5, line_alpha = 1,
             radius = "mag_size", name = "s", view = vcds_SE_WM_s)
plot2.circle("x", "depth_km", source = cds_SE_WM,
             fill_color = "orange", fill_alpha = 0.5, line_color = "orange", line_width = 0.5, line_alpha = 1,
             radius = "mag_size", name = "m", view = vcds_SE_WM_m)
plot2.circle("x", "depth_km", source = cds_SE_WM,
             fill_color = "red", fill_alpha = 0.5, line_color = "red", line_width = 0.5, line_alpha = 1,
             radius = "mag_size", name = "l", view = vcds_SE_WM_l)

plot3.circle("depth_km", "y", source = cds_SE_WM,
             fill_color = "green", fill_alpha = 0.5, line_color = "green", line_width = 0.5, line_alpha = 1,
             radius = "mag_size", name = "s", view = vcds_SE_WM_s, radius_dimension='y')
plot3.circle("depth_km", "y", source = cds_SE_WM,
             fill_color = "orange", fill_alpha = 0.5, line_color = "orange", line_width = 0.5, line_alpha = 1,
             radius = "mag_size", name = "m", view = vcds_SE_WM_m, radius_dimension='y')
plot3.circle("depth_km", "y", source = cds_SE_WM,
             fill_color = "red", fill_alpha = 0.5, line_color = "red", line_width = 0.5, line_alpha = 1,
             radius = "mag_size", name = "l", view = vcds_SE_WM_l, radius_dimension='y')

##################################################
## seismic stations
##################################################

plot1.triangle("x2", "y2", source = cds_SS_WM,
               fill_color = "white", fill_alpha = 0.6, line_color = "black", line_width = 0.5, line_alpha = 0.6,
               size = 10, legend_label = "Seismic Station (SS)", name = "n", view = vcds_SS_WM_0)
plot1.hex("x2", "y2", source = cds_SS_WM,
          fill_color = None, fill_alpha = 0.6, line_color = "white", line_width = 2, line_alpha = 0.6,
          size = 20, legend_label = "SS + Accelerometer", name = "y", view = vcds_SS_WM_1)

#############################################################################################################
## add tools
#############################################################################################################

## switch legend items on/off
plot1.legend.click_policy = "hide"
plot1.legend.location = "bottom_left"
plot1.legend.label_text_font = "times"
plot1.legend.label_text_font_style = "italic"
plot1.legend.background_fill_color = None
plot1.legend.border_line_color = None
plot1.legend.inactive_fill_alpha = 0.3
plot1.legend.inactive_fill_color = "black"

## format slider text and create list of unique months
div = Div(text = "Jan 2014",
          style = {"font-size":"32.5px", "color":"white", "background":"black", "padding":"16px",
                   "text-align":"center", "color":"#808080", "font-family":"times"},
          align = ("end"),
          margin = (-1, 0, 0, -25), height_policy = "auto", width = 110, height = 135, sizing_mode = "fixed") 
str_list = gpd_SE_WM["month"].unique().tolist()
str_list = pd.to_datetime(str_list, format = "%Y-%m").strftime("%b %Y")

## add monthly slider
slider = Slider(start = 0, end = max(date), value = 0, step = 1, title = None, bar_color = "black",
                background = "black", tooltips = False, margin = (-130, 0, 0, 80), width_policy = "fixed",
                width = 350, height = 20)

## we parse five arguments to our CustomJS callback
callback = CustomJS(args = dict(s1 = cds_SE_WM, s2 = cds_SS_WM, slider = slider, div = div, str_list = str_list),
                    code =
"""
    const data = s1.data
    const data2 = s2.data
    const date = slider.value
    div.text = str_list[date]
    
    const old_x = data['x']
    const old_y = data['y']
    const old_datetime = data['datetime']
    const old_latitude = data['latitude']
    const old_longitude = data['longitude']
    const old_magnitude = data['magnitude']
    const old_depth_km = data['depth_km']
    const old_origin = data['origin']
    const old_threshold = data['threshold']
    const old_mag_size = data['mag_size']
    
    const old_x2 = data2['x2']
    const old_y2 = data2['y2']
    const old_network2 = data2['network2']
    const old_station2 = data2['station2']
    const old_latitude2 = data2['latitude2']
    const old_longitude2 = data2['longitude2']
    const old_elevation_m2 = data2['elevation_m2']
    const old_starttime2 = data2['starttime2']
    const old_endtime2 = data2['endtime2']
    const old_co_loc2 = data2['co_loc2']
    const old_origin2 = data2['origin2']
    
    const x = data['xs'][date]
    const y = data['ys'][date]
    const datetime = data['datetimes'][date]
    const latitude = data['latitudes'][date]
    const longitude = data['longitudes'][date]
    const magnitude = data['magnitudes'][date]
    const depth_km = data['depth_kms'][date]
    const origin = data['origins'][date]
    const threshold = data['thresholds'][date]
    const mag_size = data['mag_sizes'][date]
    
    const x2 = data2['xs2'][date]
    const y2 = data2['ys2'][date]
    const network2 = data2['networks2'][date]
    const station2 = data2['stations2'][date]
    const latitude2 = data2['latitudes2'][date]
    const longitude2 = data2['longitudes2'][date]
    const elevation_m2 = data2['elevation_ms2'][date]
    const starttime2 = data2['starttimes2'][date]
    const endtime2 = data2['endtimes2'][date]
    const co_loc2 = data2['co_locs2'][date]
    const origin2 = data2['origins2'][date]
    
    for (var i = 0; i < old_x.length; i++) {old_x[i] = x[i]}
    for (var i = 0; i < old_y.length; i++) {old_y[i] = y[i]}
    for (var i = 0; i < old_datetime.length; i++) {old_datetime[i] = datetime[i]}
    for (var i = 0; i < old_latitude.length; i++) {old_latitude[i] = latitude[i]}
    for (var i = 0; i < old_longitude.length; i++) {old_longitude[i] = longitude[i]}
    for (var i = 0; i < old_magnitude.length; i++) {old_magnitude[i] = magnitude[i]}
    for (var i = 0; i < old_depth_km.length; i++) {old_depth_km[i] = depth_km[i]}
    for (var i = 0; i < old_origin.length; i++) {old_origin[i] = origin[i]}
    for (var i = 0; i < old_threshold.length; i++) {old_threshold[i] = threshold[i]}
    for (var i = 0; i < old_mag_size.length; i++) {old_mag_size[i] = mag_size[i]}
    
    for (var i = 0; i < old_x2.length; i++) {old_x2[i] = x2[i]}
    for (var i = 0; i < old_y2.length; i++) {old_y2[i] = y2[i]}
    for (var i = 0; i < old_network2.length; i++) {old_network2[i] = network2[i]}
    for (var i = 0; i < old_station2.length; i++) {old_station2[i] = station2[i]}
    for (var i = 0; i < old_latitude2.length; i++) {old_latitude2[i] = latitude2[i]}
    for (var i = 0; i < old_longitude2.length; i++) {old_longitude2[i] = longitude2[i]}
    for (var i = 0; i < old_elevation_m2.length; i++) {old_elevation_m2[i] = elevation_m2[i]}
    for (var i = 0; i < old_starttime2.length; i++) {old_starttime2[i] = starttime2[i]}
    for (var i = 0; i < old_endtime2.length; i++) {old_endtime2[i] = endtime2[i]}
    for (var i = 0; i < old_co_loc2.length; i++) {old_co_loc2[i] = co_loc2[i]}
    for (var i = 0; i < old_origin2.length; i++) {old_origin2[i] = origin2[i]}
    
    s1.change.emit()
    s2.change.emit()
"""
                   )
slider.js_on_change("value", callback)

## add hover tools (seismic events and seismic stations)
hover_se = HoverTool(names = ["s", "m", "l"],
                     tooltips =
"""
<div style = "font-family:times;font-size:12px;color:#808080;background-color:black;padding:10px;margin:-2px">
    <table>
        <tr>
            <td style = "text-align:right;color:#cccccc">DateTime:</td>
            <td>@datetime</td>
        </tr>
        <tr>
            <td style = "text-align:right;color:#cccccc">Latitude:</td>
            <td>@latitude</td>
        </tr>
        <tr>
            <td style = "text-align:right;color:#cccccc">Longitude:</td>
            <td>@longitude</td>
        </tr>
        <tr>
            <td style = "text-align:right;color:#cccccc">Magnitude (Ml):</td>
            <td>@magnitude</td>
        </tr>
        <tr>
            <td style = "text-align:right;color:#cccccc">Depth (km):</td>
            <td>@depth_km</td>
        </tr>
        <tr>
            <td style = "text-align:right;color:#cccccc">Origin:</td>
            <td>@origin</td>
        </tr>
    </table>
</div>
"""
                    )

hover_ss = HoverTool(names = ["n", "y"],
                     toggleable = False,
                     tooltips =
"""
<div style = "font-family:times;font-size:12px;color:#808080;background-color:black;padding:10px;margin:-2px">
    <table>
        <tr>
            <td style = "text-align:right;color:#cccccc">Network:</td>
            <td>@network2</td>
        </tr>
        <tr>
            <td style = "text-align:right;color:#cccccc">Station:</td>
            <td>@station2</td>
        </tr>
        <tr>
            <td style = "text-align:right;color:#cccccc">Latitude:</td>
            <td>@latitude2</td>
        </tr>
        <tr>
            <td style = "text-align:right;color:#cccccc">Longitude:</td>
            <td>@longitude2</td>
        </tr>
        <tr>
            <td style = "text-align:right;color:#cccccc">Elevation (m):</td>
            <td>@elevation_m2</td>
        </tr>
        <tr>
            <td style = "text-align:right;color:#cccccc">Start Time:</td>
            <td>@starttime2</td>
        </tr>
        <tr>
            <td style = "text-align:right;color:#cccccc">End Time:</td>
            <td>@endtime2</td>
        </tr>
        <tr>
            <td style = "text-align:right;color:#cccccc">Origin:</td>
            <td>@origin2</td>
        </tr>
    </table>
</div>
"""
                    )
plot1.add_tools(hover_se, hover_ss)

## set autohide to true to only show the toolbar when mouse is over plot
plot1.toolbar.autohide = True

#############################################################################################################
## show results
#############################################################################################################

layout = gridplot([[plot1, plot3],[plot2, div],[slider]], merge_tools = False)
show(layout)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.2117355022710468, 'fraction_numerical': 0.027411006654695258, 'mean_word_length': 3.603574033552152, 'pattern_counts': {'":': 15, '<': 141, '<?xml version=': 0, '>': 122, 'https://': 6, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 41, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '122802', 'license': 'MIT', 'n_tokens_mistral': 13622, 'n_tokens_neox': 12937, 'n_words': 2890, 'path': 'QuakeMAP.ipynb', 'repo_name': 'TomSwinscoe/QuakeMAP', 'size': '37877'}