examples/using_pybarsim.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "98f09d53-5177-4944-b5e6-d787fa320009",
   "metadata": {},
   "source": [
    "<div style=\"display: flex;\">\n",
    "    <span style=\"margin-left: auto; margin-right: 5px;  margin-top: 1px;\"><a href=\"https://creativecommons.org/licenses/by/4.0/\">\n",
    "        <img src=\"https://licensebuttons.net/l/by/4.0/80x15.png\" />\n",
    "    </a></span>\n",
    "    <span style=\"margin-right: auto; margin-left: 5px;\"><a href=\"https://opensource.org/licenses/MIT\">\n",
    "        <img src=\"https://img.shields.io/badge/License-MIT-green.svg\" />\n",
    "    </a></span>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c4c147f1",
   "metadata": {},
   "source": [
    "# How to use pyBarSim?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f89108b0-8cf3-4cc3-b87a-5e1a7ec62555",
   "metadata": {},
   "source": [
    "**Author:** Guillaume Rongier"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a70c651a-9517-4414-b8cd-2a012d9e9ea0",
   "metadata": {},
   "source": [
    "In this notebook, we will look at different functionalities of pyBarSim to simulate the deposits of wave-dominated shallow-marine environments in 2D."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "df6f1038-8d98-4e3a-a003-978e2d3de693",
   "metadata": {},
   "source": [
    "### Imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f76292dc",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.colors as mcolors\n",
    "import matplotlib.animation as animation\n",
    "from IPython.display import HTML\n",
    "\n",
    "from pybarsim import BarSim2D"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f48f6a34",
   "metadata": {},
   "source": [
    "## 1. Setup and run"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "72f0c691",
   "metadata": {},
   "source": [
    "Define the initial elevation and cell size (in m):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5d767177",
   "metadata": {},
   "outputs": [],
   "source": [
    "initial_elevation = np.linspace(1000., 900., 200)\n",
    "\n",
    "spacing = 100."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d337983b-c841-4afa-8a8e-7e79ea841f05",
   "metadata": {},
   "source": [
    "Define the run time (in yr) and the inflection points for the variations of sea level (in m):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ce1c3da3-2665-4eb9-a28c-944f8a169ec4",
   "metadata": {},
   "outputs": [],
   "source": [
    "run_time = 25000.\n",
    "\n",
    "sea_level_curve = np.array([\n",
    "    (0., 998.),\n",
    "    (0.25*run_time, 985.),\n",
    "    (0.5*run_time, 975.),\n",
    "    (0.75*run_time, 985.),\n",
    "    (run_time, 998.)\n",
    "])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4e7b20bf-5e5b-4ef5-adb1-6f93d5dc67c2",
   "metadata": {},
   "source": [
    "Define the inflection points for the variations of sediment supply (in m$^2$/yr):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "082bf111-eda7-4cb2-bd04-413f1a9116ac",
   "metadata": {},
   "outputs": [],
   "source": [
    "sediment_supply_curve = np.array([\n",
    "    (0., 25.),\n",
    "    (0.25*run_time, 25.),\n",
    "    (0.5*run_time, 25.),\n",
    "    (0.75*run_time, 5.),\n",
    "    (run_time, 5.)\n",
    "])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "639806d0",
   "metadata": {},
   "source": [
    "Initialize a `BarSim2D` object and run the simulation:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "60c0e28e-268d-4c7f-9c55-ac354637839d",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "<b>&#9888;</b> This takes more time to run the first time because Numba needs to compile the Python code (around 15 s against less than 1 s for the following runs).\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e96562cb",
   "metadata": {},
   "outputs": [],
   "source": [
    "barsim = BarSim2D(initial_elevation,\n",
    "                  sea_level_curve,\n",
    "                  sediment_supply_curve,\n",
    "                  spacing=spacing,\n",
    "                  max_wave_height_fair_weather=1.5,\n",
    "                  allow_storms=True,\n",
    "                  start_with_storm=False,\n",
    "                  max_wave_height_storm=6.,\n",
    "                  tidal_amplitude=2.,\n",
    "                  min_tidal_area_for_transport=100.,\n",
    "                  sediment_size=(5., 50., 125., 250.),\n",
    "                  sediment_fraction=(0.25, 0.25, 0.25, 0.25),\n",
    "                  initial_substratum=(100., (0.25, 0.25, 0.25, 0.25)),\n",
    "                  erodibility=0.1,\n",
    "                  washover_fraction=0.5,\n",
    "                  tide_sand_fraction=0.3,\n",
    "                  depth_factor_backbarrier=5.,\n",
    "                  depth_factor_shoreface=10.,\n",
    "                  local_factor_shoreface=1.5,\n",
    "                  local_factor_backbarrier=1.,\n",
    "                  fallout_rate_backbarrier=0.,\n",
    "                  fallout_rate_shoreface=0.0002,\n",
    "                  max_width_backbarrier=500.,\n",
    "                  curve_preinterpolation=None,\n",
    "                  seed=42)\n",
    "barsim.run(run_time, dt_fair_weather=15., dt_storm=1.)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2fd9cd7e",
   "metadata": {},
   "source": [
    "`run` creates `sequence_`, a [xarray](https://docs.xarray.dev/en/stable/) dataset containing the values of sea level, sediment supply, elevation, stratigraphy, and facies in time:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "84ee5e33",
   "metadata": {},
   "outputs": [],
   "source": [
    "barsim.sequence_"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d670785a-dbe2-463f-bdad-7df8cac36407",
   "metadata": {},
   "source": [
    "We can visualize all the variables and their variation through time using [xarray's plotting functions](https://docs.xarray.dev/en/stable/user-guide/plotting.html?highlight=plotting):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "19aaf204-5b33-4b9f-85c3-225e73481392",
   "metadata": {},
   "outputs": [],
   "source": [
    "barsim.sequence_['Sea level'].plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a6ce1b30-ee4c-47f7-8826-5f70e1507d6b",
   "metadata": {},
   "outputs": [],
   "source": [
    "barsim.sequence_['Elevation'].plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "22be0053-a252-4993-aadb-94a1477e0a46",
   "metadata": {},
   "source": [
    "## 2. Stratigraphy visualization"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5019d061-f1f1-49c7-b4ba-2c87e01fb62c",
   "metadata": {},
   "source": [
    "In `sequence_`, the stratigraphy is directly the final stratigraphy (i.e., it stores the remaining deposits after erosion), while the elevation corresponds to the true evolution of elevation through time. To visualize the final stratigraphy, `finalize` will update the elevation to account for erosion (new variable `Horizons`), and compute the mean grain size and the sorting term:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a8596160-4284-4d74-b47e-6fc53c52682b",
   "metadata": {},
   "outputs": [],
   "source": [
    "barsim.finalize()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "49346076-9221-4b4f-97ab-6535d099f280",
   "metadata": {},
   "outputs": [],
   "source": [
    "barsim.sequence_"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7fb5985d-a0cc-4dfe-b30e-8419fde30644",
   "metadata": {},
   "source": [
    "When the number of time steps gets too high, plotting takes longer and can become distorted. `subsample` creates `subsequence_`, a [xarray](https://docs.xarray.dev/en/stable/) dataset with a given number of time steps (here 20) merged together:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bd51217c-3d06-4661-913f-feced2c8d6a7",
   "metadata": {},
   "outputs": [],
   "source": [
    "barsim.subsample(20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "987e15c4-beef-4d69-a551-97de2435b2f8",
   "metadata": {},
   "outputs": [],
   "source": [
    "barsim.subsequence_"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ba7b1ca6-f44d-4395-982d-9877262894c6",
   "metadata": {},
   "source": [
    "Calling `finalize` again specifically on `subsequence_` computes the mean grain size, the sorting term, and the major facies:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dcba5f34-6f79-4e73-8074-dde1647f8a23",
   "metadata": {},
   "outputs": [],
   "source": [
    "barsim.finalize(on='subsequence')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0e34efad-f520-4027-b94e-917db78ce5ad",
   "metadata": {},
   "outputs": [],
   "source": [
    "barsim.subsequence_"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "141fe1e2-8943-43c5-8913-a176ff80bb81",
   "metadata": {},
   "source": [
    "Similarly to `sequence_`, we can visualize the values in time using [xarray's plotting functions](https://docs.xarray.dev/en/stable/user-guide/plotting.html?highlight=plotting):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f5d199c8-68bc-4164-899d-c979bb56ce94",
   "metadata": {},
   "outputs": [],
   "source": [
    "barsim.subsequence_['Mean grain size'].plot();"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cea595d2-18e5-4b88-8050-174985225f72",
   "metadata": {},
   "source": [
    "Or we can use the function `plot_subsequence` to plot the final stratigraphy in space:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dce69733-d879-4124-b349-8ac08d6c782c",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(12, 4))\n",
    "ax.fill_between(barsim.subsequence_['X'][:-1],\n",
    "                barsim.subsequence_['Horizons'][0, :-1],\n",
    "                barsim.subsequence_['Horizons'][0, :-1].min(),\n",
    "                color='#d9d9d9')\n",
    "c = barsim.plot_subsequence(ax, var='Mean grain size')\n",
    "fig.colorbar(c[0], ax=ax, label=r'Mean grain size ($\\mu$m)')\n",
    "ax.set(xlabel='x (m)', ylabel='z (m)');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bcc39fa4-cedd-4a10-a36d-97f1b3ff09e6",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(12, 4))\n",
    "ax.fill_between(barsim.subsequence_['X'][:-1],\n",
    "                barsim.subsequence_['Horizons'][0, :-1],\n",
    "                barsim.subsequence_['Horizons'][0, :-1].min(),\n",
    "                color='#d9d9d9')\n",
    "c = barsim.plot_subsequence(ax,\n",
    "                            var='Major facies',\n",
    "                            cmap='Set2',\n",
    "                            norm=mcolors.BoundaryNorm([0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5], 6))\n",
    "cbar = fig.colorbar(c[0], ax=ax, label=r'Major facies')\n",
    "cbar.set_ticks(ticks=[1, 2, 3, 4, 5, 6],\n",
    "               labels=barsim.subsequence_['Environment'].values[1:])\n",
    "ax.set(xlabel='x (m)', ylabel='z (m)');"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "229390f1-1258-418c-93e9-7b2e1e98dcc8",
   "metadata": {},
   "source": [
    "We can also plot specific grain sizes or facies using the `idx` parameter:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7b554b9b-c621-4a6a-ac81-a40b840d27aa",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(12, 4))\n",
    "ax.fill_between(barsim.subsequence_['X'][:-1],\n",
    "                barsim.subsequence_['Horizons'][0, :-1],\n",
    "                barsim.subsequence_['Horizons'][0, :-1].min(),\n",
    "                color='#d9d9d9')\n",
    "c = barsim.plot_subsequence(ax, var='Facies', idx=5, mask_zeros=False)\n",
    "fig.colorbar(c[0], ax=ax, label=r'Fraction of ' + str(barsim.subsequence_['Environment'][5].values))\n",
    "ax.set(xlabel='x (m)', ylabel='z (m)');"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "731cd9d4-0033-41ec-b404-6425fca4ef99",
   "metadata": {},
   "source": [
    "## 3. Stratigraphy regridding"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5af54008-d2cf-44cf-a8bc-34b494dc2b63",
   "metadata": {},
   "source": [
    "`subsequence_`'s grid is irregular and can be difficult to use in subsequent simulations. `regrid` reinterpolates BarSim's outputs on a regular grid:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "70707832-6a7a-46fb-9b3c-21be9b97aeda",
   "metadata": {},
   "outputs": [],
   "source": [
    "barsim.regrid(900., 1000., 0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f8bf4913",
   "metadata": {},
   "source": [
    "`regrid` creates `record_`, a [xarray](https://docs.xarray.dev/en/stable/) dataset containing the stratigraphy and facies in space:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "036bc2c1",
   "metadata": {},
   "outputs": [],
   "source": [
    "barsim.record_"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8c817b98-6694-4ebc-8ae0-04186a822951",
   "metadata": {},
   "source": [
    "`finalize` computes once again the mean grain size, the sorting term, and the major facies:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3cfe8090-c2ae-4046-bc41-61e71ab780b8",
   "metadata": {},
   "outputs": [],
   "source": [
    "barsim.finalize(on='record')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8a297b20-d25e-417c-a666-d4b64f0f9daf",
   "metadata": {},
   "outputs": [],
   "source": [
    "barsim.record_"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "20ac0928",
   "metadata": {},
   "source": [
    "We can visualize the resulting grid using [xarray's plotting functions](https://docs.xarray.dev/en/stable/user-guide/plotting.html?highlight=plotting):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ffb1ca54",
   "metadata": {},
   "outputs": [],
   "source": [
    "barsim.record_['Mean grain size'].plot(figsize=(12, 4))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "88507792",
   "metadata": {},
   "outputs": [],
   "source": [
    "barsim.record_['Sorting term'].plot(figsize=(12, 4))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "775c8aee",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(12, 4))\n",
    "im = barsim.record_['Major facies'].where(barsim.record_['Major facies'] > 0, np.nan).plot(ax=ax,\n",
    "                                                                                           cmap='Set2',\n",
    "                                                                                           norm=mcolors.BoundaryNorm([0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5], 7),\n",
    "                                                                                           add_colorbar=False)\n",
    "cbar = fig.colorbar(im, ax=ax, label=r'Major facies')\n",
    "cbar.set_ticks(ticks=[1, 2, 3, 4, 5, 6, 7],\n",
    "               labels=barsim.subsequence_['Environment'].values);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "128cc3c6-218c-49c7-84ed-f52166e20f9d",
   "metadata": {},
   "source": [
    "## 4. Process visualization"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f6cf32e6",
   "metadata": {},
   "source": [
    "Let's create a simple animation of the variations of sea level and elevation through time using [matplotlib](https://matplotlib.org/stable/api/animation_api.html):"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b686fe79-cb83-46d1-922f-1f7ce9f57f32",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "<b>&#9888;</b> This takes some time to run.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a7da9ab1",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(12, 4))\n",
    "\n",
    "step = 15\n",
    "time = barsim.sequence_['Time'][::step]\n",
    "sea_level = barsim.sequence_['Sea level'][::step]\n",
    "elevation = barsim.sequence_['Elevation'][::step, :-1].copy()\n",
    "x = barsim.sequence_['X'].values[:-1]\n",
    "\n",
    "def update(i):\n",
    "    label_time.set_text(str(round(int(time[i]), -min(2, int(np.log10(time[i] + 1e-8))))) + ' yr')\n",
    "    path = fill_sea.get_paths()[0]\n",
    "    path.vertices[len(elevation[i]) + 2:-1, 1] = sea_level[i]\n",
    "    path = fill_subsurface.get_paths()[0]\n",
    "    path.vertices[len(elevation[i]) + 2:-1, 1] = elevation[i][::-1]\n",
    "    k = 0\n",
    "    for j in range(i):\n",
    "        if j%15 == 0:\n",
    "            elevation_prev = elevation[j].to_numpy()\n",
    "            elevation_prev[elevation_prev > elevation[i]] = elevation[i][elevation_prev > elevation[i]]\n",
    "            lines_elevation_prev[k].set_data((x, elevation_prev))\n",
    "            k += 1\n",
    "    line_elevation.set_ydata(elevation[i])\n",
    "    return label_time, fill_sea, fill_subsurface, line_elevation, lines_elevation_prev\n",
    "\n",
    "ax.annotate('Time:', (0.85, 0.92), xycoords='axes fraction')\n",
    "label_time = ax.annotate(str(round(int(time[0]), -min(2, int(np.log10(time[0] + 1e-8))))) + ' yr', (0.965, 0.92), ha='right', xycoords='axes fraction')\n",
    "fill_sea = ax.fill_between(x, elevation.min(), sea_level[0], edgecolor='#6baed6', facecolor='#c6dbef', zorder=0)\n",
    "fill_subsurface = ax.fill_between(x, elevation.min(), elevation[0], color='#fff7bc', zorder=1)\n",
    "lines_elevation_prev = [ax.plot([], [], c='0.5', lw=0.5, zorder=2)[0] for i in range(len(elevation[::15]))]\n",
    "line_elevation, = ax.plot(x, elevation[0], c='k', lw=1.5, zorder=3)\n",
    "ax.set_xlim(x[0], x[-1])\n",
    "ax.set_ylim(elevation.min(), elevation.max() + 20.)\n",
    "ax.set_xlabel('x (m)')\n",
    "ax.set_ylabel('z (m)')\n",
    "\n",
    "plt.close()\n",
    "\n",
    "ani = animation.FuncAnimation(fig, update, len(sea_level), interval=100)\n",
    "HTML(ani.to_jshtml())"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}