Coastal Impact Beach States#

First import some necessary packages#

import logging
import os
import pathlib
import warnings

import bokeh.io
import colorcet as cc
import geopandas as gpd
import holoviews as hv
import hvplot.pandas  # noqa: API import
import ipyleaflet
import IPython
import ipywidgets as widgets
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import panel as pn
from bokeh.models import HoverTool, PanTool, WheelZoomTool
from bokeh.resources import INLINE
from IPython.display import HTML, display

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 from cloud storage
questions = cd.read_questions(
    "az://coastal-dynamics/questions/5_coastal_impact_beach_states_hashed.json",
    storage_options={"account_name": "coclico"},
)

In this notebook you might encounter some ‘FutureWarning’ or ‘PanelDeprecationWarning’. This is the result of using an old version of the holoviews package. These messages can safely be ignored. If you want to get rid of these messages, you can update the holoviews package using the code cell below (uncommented!).

# !pip install --upgrade holoviews

Morphodynamics of the upper shoreface#

Welcome to the notebook of week 5! The main topic of this notebook is beach states (section 7.3 in the book). These are discussed also in the context of coastal classification (chapter 2 in the book). This notebook consists of the following sections:

  • Coastal classification

  • Beach states

Each section contains questions for you to practice with (cross-shore) sediment transport. Let’s get started!

Coastal classification#

As you know from chapter 2 of the book and the week 1 notebook, coasts can be classified based on tectonics and processes. In this notebook we will look at different Brazilian coastal systems, and try to classify them. Let’s start with loading the dataset with data for the Brazilian coastline. This dataset contains the significant wave height (Hs), the peak period (T), the mean tidal range (MTR), the shoreface slope, and other useful values for different coordinates along the Brazilian coast. This dataset is adapted from:

Specifically, the first source is used to get values for Hs, T and MTR. The second source is used to get the RTR, D and the beachface slope for the Santa Catarina, Maranhao, and Rio de Janeiro provinces. The third source is used to get the RTR, D and the beachface slope for the Rio de Janeiro province.

Since this dataset is aggregated from multiple sources, the accuracy of its content at a local scale should be critically interpreted. Athanasiou et al. (2023) for example use a numerical model to compute offshore wave data, whereas Klein et al. (2009) use direct measurements. All of this is to say that this data is not presumed to be locally accurate. However, in this context, we use it to compare the Brazilian coastline on a larger regional scale and draw broad conclusions to illustrate the concepts introduced in this course.

Load data#

The cell below loads the data from the cloud. The first 5 entries are shown using the .head() method, so you can see what the dataframe looks like!

import pooch

fp = pooch.retrieve(
    "https://coclico.blob.core.windows.net/coastal-dynamics/5_coastal_impact_beach_states/5_data.gpkg",
    known_hash="661ddc9ad6a396dd6fe9a9cf2126a32b1134fef92075fa19f5ee1ee445125934",
)

# We load this file as a GeoDataFrame, which comes with a column containing the geometry of each entry. For this dataset, these are points (longitude, latitude)
gdf = gpd.read_file(fp)


def convert_to_numeric(value):
    if pd.isnull(value):
        return np.nan  # Convert None to NaN
    if value.startswith(">"):
        epsilon = 0.1  # Define how much greater
        return float(value[1:])
    elif value.startswith("~"):
        return float(value[1:])
    else:
        return float(value)


gdf["Omega [-]"] = gdf["Omega [-]"].apply(convert_to_numeric)

gdf.head()
Downloading data from 'https://coclico.blob.core.windows.net/coastal-dynamics/5_coastal_impact_beach_states/5_data.gpkg' to file '/home/runner/.cache/pooch/81a5cf68f6bc9db5da738afb54aeb75d-5_data.gpkg'.
Label Province Region Longitude Latitude Hs [m] T [s] MTR [m] RTR [-] D [mm] w_s [m/s] Beachface slope [degrees] Omega [-] geometry
0 Campo Bom Santa Catarina 5 -49.0572 -28.7212 1.59 8.65 0.581 0.7 0.21 0.02188 1.8 5.0 POINT (-49.0572 -28.7212)
1 Laguna Santa Catarina 4 -48.7619 -28.4884 1.45 8.55 0.583 0.9 0.19 0.01852 3.0 3.0 POINT (-48.7619 -28.4884)
2 Enseada de Pinheira Santa Catarina 4 -48.5967 -27.8547 1.34 8.51 0.647 0.9 0.19 0.01852 3.0 3.0 POINT (-48.5967 -27.8547)
3 Praia do Moçambique Santa Catarina 3 -48.4173 -27.5260 1.39 8.71 0.741 1.1 0.30 0.03720 4.1 2.0 POINT (-48.4173 -27.526)
4 Tijucas Santa Catarina 2 -48.6134 -27.2402 1.16 8.07 0.942 1.9 0.32 0.04048 5.0 1.0 POINT (-48.6134 -27.2402)

In this dataframe the fall velocity w_s was calculated using other values from the dataset. It was calculated using Soulsby (1997). You don’t have to remember this equation, but it is included here for completeness: $\( w_s = (\nu_{kin} / D) * \left( \sqrt{10.36^2 + 1.049 * (D_{*}^3)} - 10.36 \right) \)\( with \) D_{*} = \left(\frac{g * (\rho_s / \rho_w - 1)}{\nu_{kin}^2}\right) ^{1/3} * D \(, \)g=9.81\( m/s\)^2\(, \)\nu_{kin}=1.36e-6\( m\)^2\(/s, \)\rho_w=1027\( kg/m\)^3\(, \)\rho_s=2650\( kg/m\)^3$

The dimensionless fall velocity can be calculated using equation 7.8 from the book: $\( \Omega = \frac{H_b}{w_s T} \)$ We did some preliminary calculations, and the Omega values as presented in the dataset should be sufficient for qualitatively comparing different coastal systems.

Again, we would like to stress that the values should not be interpreted as highly accurate for specific beaches. They are just useful to get a feel for the surrounding area and to get first estimates of the order of magnitude.

Define plot 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 plot_brazilian_coast(plot_where="pop-out"):
    """
    change value of 'plot_where' to:
    'inline' if you would like the plot to show in the notebook
    'pop-out' if you would like the plot to show in a new tab (i.e. seperate window)
    """

    # below we build the widget
    title_bar = pn.Row(
        pn.pane.Markdown(
            "## Brazilian Coast",
            styles={"color": "black"},
            sizing_mode="fixed",
            margin=(10, 5, 10, 15),
        ),
        align="center",
    )

    # dropdown menu of coasts
    options = {
        "Pará": ["Princesa", "Atalaia", "Ajuruteua"],
        "Maranhão": ["São Luís"],
        "Rio de Janeiro": ["Marambaia"],
        "Santa Catarina": [
            "Campo Bom",
            "Laguna",
            "Enseada de Pinheira",
            "Praia do Moçambique",
            "Tijucas",
            "Balneário Camboriú",
            "Do Ubatuba",
            "Barra Velha",
        ],
    }
    coasts_dropdown = pn.widgets.Select(
        name="Coast select (grouped by province)", groups=options, value="Campo Bom"
    )

    @pn.depends(coasts_dropdown.param.value)
    def plot_coast(name, plot_size=0.04):
        beach = gdf[gdf["Label"] == name].copy()
        beach = beach.astype({"Omega [-]": str})
        lat, lon = beach[["Latitude", "Longitude"]].values.flatten()
        lat, lon = np.float64(lat), np.float64(lon)

        points = gdf.hvplot.points(
            geo=True,
            tiles="ESRI",
            ylabel="Latitude [deg]",
            xlabel="Longitude [deg]",
            xlim=(lon - plot_size / 2, lon + plot_size / 2),
            ylim=(lat - plot_size / 2, lat + plot_size / 2),
            tools=["tap"],
            hover_cols=[
                "Label",
                "Province",
                "Longitude",
                "Latitude",
                "Hs [m]",
                "T [s]",
                "MTR [m]",
                "RTR [-]",
                "D [mm]",
                "w_s [m/s]",
                "Beachface slope [degrees]",
                "Omega [-]",
            ],
            c="Province",
            cmap="Accent",
            line_color="black",
            size=300,
        )

        plot = (points).opts(width=1200, height=800, tools=["wheel_zoom"])

        return plot

    app = pn.Column(
        pn.Row(title_bar, align="center"),
        pn.Row(coasts_dropdown, align="center"),
        pn.Row(plot_coast, align="center"),
    )

    if plot_where == "inline":
        return app
    elif plot_where == "pop-out":
        app.show()
    else:
        print("please use either inline or pop-out for the plot_where variable")

Now plot the coastal data#

Execute the cell below to generate the plot by using the function we defined above. In this plot, coastal systems from 4 different Brazilian provinces are shown. These provinces (from north to south) are Pará, Maranhão, Rio de Janeiro, and Santa Catarina. 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.

logging.getLogger().setLevel(logging.ERROR)

plot_brazilian_coast(plot_where="pop-out")
Launching server at http://localhost:44473

Note: Remember that we are working with an aggregated dataset. Eeach data point should be interpreted as representative of the region around it, and not as specific to that beach.

You can freely zoom and move around in the generated plot. You can see all of the values from the dataframe when hovering over a data point with the cursor. Some of these values will be needed to answer the questions.

Using this plot, try to answer the questions below.

q1 = cd.QuestionFactory(questions["Q5-1"]).serve()
q2 = cd.QuestionFactory(questions["Q5-2"]).serve()
q3 = cd.QuestionFactory(questions["Q5-3"]).serve()
q4 = cd.QuestionFactory(questions["Q5-4"]).serve()

pn.Column(q1, q2, q3, q4)

Have a look at section 4.4.2 from the book, and figure 4.13. Let’s try to classify some of the coastal sections using the mean tidal range (MTR) and mean wave height! We plot some of the Brazilian coasts and see how they are classified using the figure. Use the Brazilian coast plot to determine relevant values for the mean tidal range and the mean wave height, and add these values to the data variable below.

# Define points here
data_413 = (
    # ("label", mean wave height [m], mean tidal range [m], relative tidal range [-]),
    ("Campo Bom", 1.590, 0.581),  # example; you can add more tuples yourself!
)

# Or uncomment the line below to get all the points
data_413 = zip(gdf.Label, gdf["Hs [m]"], gdf["MTR [m]"], gdf["RTR [-]"])

df_413 = pd.DataFrame(data_413, columns=("Label", "Hs [m]", "MTR [m]", "RTR [-]"))

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!

Again, you can hover over points to get relevant values. Don’t forget to have a look at the relative tidal ranges (RTR) for the different coastal systems! The RTR for this specific dataset is actually quite small, but in reality it can be much bigger (>15). For RTR < 3 we have the wave-dominated beaches as described in Section 4.3.3. For RTR > 15 the beaches gradually approach the pure tidal flat situation. The RTR values for the selected Para and Sao Luis beaches hint towards wave-shaped beaches with significant tidal influence and are called tide-dominated beaches in this notebook.

warnings.simplefilter("ignore", category=FutureWarning)

bokeh.io.output_notebook(INLINE)
hv.extension("bokeh")

# Load background
fp = pooch.retrieve(
    "https://coclico.blob.core.windows.net/coastal-dynamics/5_coastal_impact_beach_states/5_fig413_bg.jpg",
    known_hash="f71b11a7f30fdf49e99379a328f6d018f2181dd906b6b2c4495014336c5e1161",
)
bg = hv.RGB.load_image(fp, bounds=(0, 0, 2.5, 6)).opts(alpha=0.5)

# # create the points
points = df_413.hvplot.points(
    x="Hs [m]",
    y="MTR [m]",
    by="Label",
    size=100,
    cmap="Accent",
    line_color="black",
    hover_cols=["Label", "Hs [m]", "MTR [m]", "RTR [-]"],
)

fig413 = (bg * points).opts(
    width=700,
    height=600,
    show_grid=True,
    active_tools=[],
    toolbar=None,
    xlabel="mean wave height [m]",
    ylabel="mean tidal range [m]",
    xlim=(0, 2.5),
    ylim=(0, 6),
    show_legend=True,
)

fig413
Loading BokehJS ...
Downloading data from 'https://coclico.blob.core.windows.net/coastal-dynamics/5_coastal_impact_beach_states/5_fig413_bg.jpg' to file '/home/runner/.cache/pooch/299528ee74431ed2652539105b902e33-5_fig413_bg.jpg'.

Using this figure, and the Brazilian coast plot, answer the questions below.

q5 = cd.QuestionFactory(questions["Q5-5"]).serve()
q6 = cd.QuestionFactory(questions["Q5-6"]).serve()
q7 = cd.QuestionFactory(questions["Q5-7"]).serve()
q8 = cd.QuestionFactory(questions["Q5-8"]).serve()
q9 = cd.QuestionFactory(questions["Q5-9"]).serve()

pn.Column(q5, q6, q7, q8, q9)

This is the end of the first part of this notebook. You can continue with the next part on beach states!

Beach states#

For this section, we will focus on the wave-dominated coastal sections (i.e. the coastal sections from the Rio de Janeiro and Santa Catarina coastline). From section 7.3 of the book, you know that a wave-dominated system may be classified as reflective or dissipative, or an intermediate state in between. Note that there is a broad spectrum of intermediate beach states. In this notebook we will not look at this in detail and make some broad generalizations. For instance, we call higher intermediate states “more dissipative” and lower intermediate states “more reflective”.

A common classifier used besides the Iribarren number is the dimensionless fall velocity, which is calculated as follows:
$\( \Omega = \frac{H_b}{w_s T} \)\( where \)H_b\( is the wave height at breaking, \)T\( is the wave period and \)w_s$ is the fall velocity.

Note that the coastal sections considered as wave-dominated here all have a relative tidal range of lower than 3.

Klein et al. (2005) use beach slope and sediment size as a proxy for classifying wave-dominated systems, see also Figure 7.10 in the book. Let’s try this for ourselves! We should note that we have looked at regional values for the sediment size, and the beachface slope. Therefore some of the coastal systems from the Brazilian coast plot have equal values for the sediment size and beachface slope. This is also the likely reason that we don’t really see any reflective beaches, at least with respect to the data. Note that different coastal sections might still have varying dimensionless fall velocity, even if they are within the same region. Think about why this is the case.

The regions are displayed in the legend of the plot below, so you can see which coastal systems belong to which regions. Alternatively, a table is provided below.

Region (in Santa Catarina province)

Coastal section

D [mm]

Beachface slope [degrees]

1

Barra Velha,
Do Ubatuba

0.30

5.0

2

Balnario Camboriu,
Tijucas

0.32

5.0

3

Praia do Macambique

0.30

4.1

4

Enseada de Pinheira,
Laguna

0.19

3.0

5

Campo Bom

0.21

1.8

Finish the code below with some of the coastal sections to show where they would lie in figure 7.10.

# Define points here
data_710 = (
    # ("label", slope [degrees], mean grain size [mm], Omega [-]),
    ("Campo Bom", 2.987, 0.21, 4.385),  # example; you can add more tuples yourself!
)

# Or uncomment the lines below to get all the points
data_710 = zip(
    gdf[gdf.Province == "Santa Catarina"].Label,
    gdf[gdf.Province == "Santa Catarina"]["Beachface slope [degrees]"],
    gdf[gdf.Province == "Santa Catarina"]["D [mm]"],
    gdf[gdf.Province == "Santa Catarina"]["Omega [-]"],
)
data_710 = list(data_710)
data_710.append(
    list(
        zip(
            gdf[gdf.Province == "Rio de Janeiro"].Label,
            gdf[gdf.Province == "Rio de Janeiro"]["Beachface slope [degrees]"],
            gdf[gdf.Province == "Rio de Janeiro"]["D [mm]"],
            gdf[gdf.Province == "Rio de Janeiro"]["Omega [-]"],
        )
    )[0]
)

df_710 = pd.DataFrame(
    data=data_710, columns=("Label", "Beachface slope [degrees]", "D [mm]", "Omega [-]")
)

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!

warnings.simplefilter("ignore", category=UserWarning)

r = 1
ps = hv.Points([])

for slope, D in [(5, 0.3), (5, 0.32), (4.1, 0.3), (3, 0.19), (1.8, 0.21), (7.4, 0.6)]:
    df = df_710[df_710["Beachface slope [degrees]"] == slope][df_710["D [mm]"] == D]

    # plot the points
    if r == 6:
        point = (
            df.hvplot.points(
                x="Beachface slope [degrees]",
                y="D [mm]",
                size=100,
                c="Label",
                cmap="Accent",
                hover_cols=[
                    "Label",
                    "Beachface slope [degrees]",
                    "D [mm]",
                    "Omega [-]",
                ],
                line_color="black",
            )
            .opts(xlabel="Slope [degrees]", ylabel="Mean grain size [mm]")
            .relabel("Rio de Janeiro, Marambaia")
        )
    else:
        point = df.hvplot.points(
            x="Beachface slope [degrees]",
            y="D [mm]",
            size=100,
            cmap="Accent",
            hover_cols=["Label", "Beachface slope [degrees]", "D [mm]", "Omega [-]"],
            line_color="black",
            label=f"Santa Catarina, region {r}",
        ).opts(xlabel="Slope [degrees]", ylabel="Mean grain size [mm]")

    # add points to ps variable
    ps *= point

    # used for labelling of the regions
    r += 1

# plot horizontal and vertical lines
hlines = hv.HLine(0.25).opts(color="lightgrey") * hv.HLine(0.5).opts(color="lightgrey")
vlines = (
    hv.VLine(3.5).opts(color="lightgrey") * hv.VLine(8.5).opts(color="lightgrey")
).opts(border_line_color="lightgrey")

# plot labels
classify_labels = (
    hv.Text(1.75, 0.95, "dissipative", fontsize=10)
    * hv.Text(6, 0.95, "intermediate", fontsize=10)
    * hv.Text(13, 0.95, "reflective", fontsize=10)
)

grain_labels = (
    hv.Text(15.8, 0.75, "coarse sand", fontsize=10, halign="right")
    * hv.Text(15.8, 0.375, "medium sand", fontsize=10, halign="right")
    * hv.Text(15.8, 0.125, "fine sand", fontsize=10, halign="right")
)


fig710 = (ps * hlines * vlines * classify_labels * grain_labels).opts(
    ylim=(0, 1), xlim=(0, 16), width=1200, height=600, active_tools=[], toolbar=None
)

fig710

Using the plot above, as well as the Brazilian coast plot, try to answer the questions below.

q10 = cd.QuestionFactory(questions["Q5-10"]).serve()
q11 = cd.QuestionFactory(questions["Q5-11"]).serve()
q12 = cd.QuestionFactory(questions["Q5-12"]).serve()
q13 = cd.QuestionFactory(questions["Q5-13"]).serve()
q14 = cd.QuestionFactory(questions["Q5-14"]).serve()
q15 = cd.QuestionFactory(questions["Q5-15"]).serve()
q16 = cd.QuestionFactory(questions["Q5-16"]).serve()
q17 = cd.QuestionFactory(questions["Q5-17"]).serve()
q18 = cd.QuestionFactory(questions["Q5-18"]).serve()
q19 = cd.QuestionFactory(questions["Q5-19"]).serve()
q20 = cd.QuestionFactory(questions["Q5-20"]).serve()
q21 = cd.QuestionFactory(questions["Q5-21"]).serve()

pn.Column(q10, q11, q12, q13, q14, q15, q16, q17, q18, q19, q20, q21)

This is the end of the notebook for week 5.