1a: Tectonic classification¶
Welcome to the first week of notebook exercises for Coastal Systems (TU Delft, MSc Coastal Engineering)! This is the third year that we use notebooks in Coastal Systems. With these notebooks, we hope to support the Coastal Dynamics Open Textbook with interactive material that helps you better understand the coastal processes and concepts explained in this course. Please let us know how you like the notebooks - we appreciate your feedback!
Each week of the course has at least 1 accompanying notebook. Some weeks may have multiple notebooks. This week, we will have 2 notebooks: notebooks 1a and 1b. These two notebooks are about the material in Chapter 2 of the Open Textbook, which describes the large geographical variation of coasts across the world. This chapter explains how the coasts that we have today are shaped by both present-day processes and processes millions of years ago. It distinguishes between three different orders of features, which are associated with different orders of time. In the notebooks of this first week, we will look at coastal systems at these different orders of scale.
Notebook 1a starts with the broadest (or first-order) features of the coast that cover large geographical distances (thousands of kilometres) and are linked to the long-term geological process of plate tectonics. We will do so by using earthquake data from the USGS. The dataset that we load contains a sample (10%) of observed earthquakes between Dec 2018 and Jan 2000. Why earthquake data? Earthquake data could help reveal the geophysical processes on earth, which is not only insightful to geologists, but to coastal researchers as well.
We will first prepare the visualization of the data. Next, you will be asked to answer 12 multiple choice or multiple selection questions.
Import libraries that we use for our analysis¶
Let’s first import the libraries that we use for our analysis by running the next cells.
from pathlib import Path
import warnings
from bokeh.models import PanTool, WheelZoomTool
import colorcet as cc
import geopandas as gpd
import holoviews as hv
import hvplot.pandas # noqa: API import
import numpy as np
import pandas as pd
import panel as pn
import pooch
import coastal_dynamics as cd
# Silence DeprecationWarning # Future TODO: in spring 2024 rm this silence and check status
with warnings.catch_warnings():
warnings.filterwarnings("ignore", category=DeprecationWarning)
import dask.dataframe as dd
# Activate Panel extension to make interactive visualizations
pn.extension()# read questions locally
questions = cd.read_questions(Path("../hashed_questions/1_coastal_classification_hashed.json"))
question_industry = cd.QuestionIndustry(questions)
cd.UseAnswersApp("1a").serve()Load the earthquake data and prepare the visualization¶
Load the earthquake data¶
We load the data (tabular data including geometries) and index the columns to keep only the data that we actually need in memory. In total, the dataset contains 2.2 million earthquakes, but here we use a sample (10%), so our subset contains approx. 220k earthquake entries. If you find that the interactive panel responds slowly to how you adjust the widgets, please consider taking another sample. You can do so by uncommenting the sample line in the next cell. If you set frac=0.1 you will therefore have a dataframe with approx. 22k unique earthquakes all over the world.
earthquakes_fp = Path("../database/1_coastal_classification/1_earthquakes_sample.parquet")
WEB_MERCATOR_LIMITS = (
-20037508.342789244,
20037508.342789244,
) # max polar latitudes that can be handled in World Mercator
df = (
dd.read_parquet(earthquakes_fp)
.sample(
frac=0.1
) # uncomment this line if loading the data takes too long on your computer
.set_index("time")
.compute()
.tz_localize(None)
.sort_index()
)
# To save memory we drop most of the columns. Also we drop the polar latitudes that cannot be displayed in the web mercator projection.
df = df[["mag", "depth", "latitude", "longitude", "place", "type"]][
df["northing"] < WEB_MERCATOR_LIMITS[1]
]Load the bathymetric contours¶
To link the data on earthquakes with the physical shape of the ocean, we need to load in some bathymetry (underwater topography). We take the bathymetric contours (also known as isobaths) for a water depth of -200 m from ArcGIS, which we will use as a proxy to approximate the boundary of the continental shelf (see Chapter 2.3.3 in the textbook). Because we want to maintain interactive plots, all shorter isobaths are dropped. Computing lengths is a metric operation, so the data first has to be reprojected from a geographic coordinate system in degrees latitude/longitude to a planar projection.
isobath_fp = Path("../database/1_coastal_classification/1_isobaths200.gpkg")
data200 = gpd.read_file(isobath_fp)
data200["length"] = data200.to_crs("EPSG:3857").geometry.length
data200 = data200[data200["length"] > 5 * 10**6]Define the plot function for the earthquake data¶
To explore the earthquake data, we use tools from the Holoviz project, which makes high-level tools to simplify visualization in Python. Run the cell below to generate the show_earthquakes function.
Note: Although you don’t have to understand the plot method, we include it here so you can see how these interactive plots are made!
def show_earthquakes():
"""
Creates an app showing a global map with earthquake data.
"""
# Below we build the earthquake widget
title_bar = pn.pane.Markdown(
"## Global distribution of earthquakes",
styles={"color": "black"},
width=400,
# margin=(10, 5, 10, 15),
)
# define widgets that can be used to index the data
magnitude_slider = pn.widgets.RangeSlider(
name="Earthquake magnitude [Richter]", start=0.1, end=10
)
depth_slider = pn.widgets.RangeSlider(
name="Earthquake depth [km]", start=0.1, end=650
)
date_slider = pn.widgets.DateRangeSlider(
name="Date", start=df.index[0], end=df.index[-1]
)
column_types = pn.widgets.Select(
name="Show earthquake magnitude or depth?", options=["mag", "depth"]
)
plot_isobaths = pn.widgets.Select(
name="Plot isobaths -200m?", options=["no", "yes"]
)
@pn.depends(
magnitude_slider.param.value_start,
magnitude_slider.param.value_end,
depth_slider.param.value_start,
depth_slider.param.value_end,
date_slider.param.value_start,
date_slider.param.value_end,
column_types.param.value,
plot_isobaths.param.value,
)
def plot_earthquake_panel(
magnitude_start,
magnitude_end,
depth_start,
depth_end,
date_start,
date_end,
column_type,
plot_isobath,
):
panel = df[
(df.mag > magnitude_start)
& (df.mag < magnitude_end)
& (df.depth > depth_start)
& (df.depth < depth_end)
& (df.index >= pd.Timestamp(date_start))
& (df.index <= pd.Timestamp(date_end))
]
# inverted fire colormap from colorcet
cmap = cc.CET_L4[::-1]
colorbar_labels = {
"mag": "Magnitude [Richter]",
"depth": "Earthquake depth [km]",
}
p = panel.hvplot.points(
x="longitude",
y="latitude",
geo=True,
color=column_type,
# global_extent=True,
tiles="EsriOceanBase",
# frame_width=900,
ylabel="Latitude [deg]",
xlabel="Longitude [deg]",
cmap=cmap,
tools=["tap"],
hover_cols=["place", "time"],
logz=True,
clim=(1, None),
clabel=colorbar_labels[column_type],
framewise=True,
xlim=(-180, 180),
ylim=(-65, 77)
)
if plot_isobath == "yes":
baths = data200.hvplot(
geo=True, line_width=2, line_color="white", line_dash="dashed"
)
p = p * baths
p.opts(tools=["wheel_zoom"])
p.opts(frame_width=800)
# p.opts(height=600)
p.opts()
return p
earthquake_panel = pn.Column(
pn.Row(
pn.Column(
pn.Row(title_bar, align="start"),
pn.Row(plot_isobaths, align="start"),
pn.Row(column_types, align="start"),
),
pn.Column(
pn.Row(magnitude_slider, align="start"),
pn.Row(depth_slider, align="start"),
pn.Row(date_slider, align="start"),
),
pn.Column(),
),
pn.Row(plot_earthquake_panel, align="center"),
)
return earthquake_panelPlot the earthquake data¶
Execute the cell below to generate the plot by using the function we defined above. Please note that altering the slider positions or selecting different options from the dropdown menus may trigger a warning; it can safely be ignored, and possibly silenced by adjusting the logging warning level.
After running the cells below you will have a panel with several widgets to index the earthquake data; by magnitude, depth and time, while the colours on the map show either the magnitude or the depth of the earthquakes. Here we consider the depth of an earthquake to be the distance below the earth’s surface at which the epicentre of the earthquake occurs. Earthquakes occur in the crust or upper mantle, which ranges from the earth’s surface to about 800 km deep. If you look very carefully, you might observe that according to the dataset, an unexpectedly high number of earthquakes occur at a depth of 10 km. Ten kilometres is a “fixed depth”. Sometimes data are too scarce or low-quality to compute a reliable depth for an earthquake. In such cases, the depth is assigned to be 10 km (source). Whenever you examine real-world datasets, it is always a good idea to keep an eye out for such assumptions.
If the visualization is too slow, please adjust the sliders such that less data is shown.
The earthquake data is shown on top of a bathymetric and topographic map, which is made available by ESRI. In general, the lighter blue colour in the figure shows the shallower part of the ocean (i.e., the continental shelf). Note that in shallower water (i.e., when zooming in to regions close to the coast and in estuaries), a darker shade of blue is shown again. This behaviour relates to how the map data is stored and displayed, which is much more efficient for us to load, but with the unfortunate side effect that (at least in this case) the colourmap is not constant for all scales.
For efficiency, the plots are generated without the -200m isobathymetry by default. Enable this feature if you would like to see detailed water depth contours, but note that it will increase the plotting time.
app = show_earthquakes()
cd.launch_app(app)'See the app at: http://localhost:39045'Questions¶
By running the below cells you get four blocks of in total 12 questions. You can use the earthquake panel in combination with your textbook to answer the questions.
Question block 1: Identifying the tectonic plate setting from earthquake data¶
As coastal engineers we may not primarily be interested in the earthquake data by itself. However, earthquake data show evidence of the most fundamental processes in geology: plate tectonics. Although plate tectonics is a relatively slow process that acts on the geological time scale, it has had an enormous impact on the formation of coastlines and determines the broadest features of the coast.
The first three questions will help you identify the location and type of plate boundaries by means of the earthquake data. These questions can be answered by using the earthquake panel and Figure 2.2 from the book. You can also cross-check your answers using Figure 2.4. For convenience, we suggest to open the Panel dashboard in a separate browser tab to view the questions and plot side by side. Run the code cell below to access the questions.
q = [
"Q1a-earthquakes_diverging_boundaries",
"Q1a-earthquakes_250km",
"Q1a-earthquakes_divergent_vs_convergent",
]
question_industry.serve(q)Question block 2: Tectonic influence on coasts¶
In 1971, Inman & Nordstrom used plate tectonics to classify the coast (Chapter 2.3.2). They distinguished between three main types of coasts: leading edge, trailing edge and marginal sea. In the questions below, match the mentioned coasts to each class.
q = [
"Q1a-wilmington",
"Q1a-shanghai",
"Q1a-maputo",
"Q1a-lima",
]
question_industry.serve(q)Question block 3: The influence of sediment supply¶
It’s not just the location of the coast relative to plate boundaries that matters; Inman, D. L. & Nordstrom (1971) further distinguish Afro-trailing-edge coasts and Amero-trailing-edge coasts based on differences in sediment supplies. Let’s take a look at some of the sites mentioned in Table 2.3 with this in mind.
q = [
"Q1a-large_river_deltas",
"Q1a-amazon_and_mekong",
"Q1a-redfish_pass",
]
question_industry.serve(q)Question block 4: Australian Gold Coast¶
For the upcoming two questions, consider the east Australian “Gold Coast”, which is located around the easternmost point of Australia (south of Brisbane).
q = [
"Q1a-gold_coast_shelf_width",
"Q1a-gold_coast_marginal",
]
question_industry.serve(q)The end¶
You have reached the end of this Notebook 1a. You can continue with this week’s second notebook, Notebook 1b on process-based classification.