lessons/permamodel/PFdepth-dependent-diffusion.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import rasterio as rio\n",
    "import xarray as xr\n",
    "from landlab import RasterModelGrid\n",
    "from landlab.components import DepthDependentDiffuser\n",
    "from matplotlib.colors import LightSource\n",
    "from netCDF4 import Dataset\n",
    "from scipy import ndimage"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Coupling with Landlab components"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook gives one *very simplistic* example for how Ku can be used alongside landscape geomorphology models. We will use the active layer depth results from Ku to drive a depth-dependent hillslope diffusion model over the Eight Mile Lake study site."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Eight Mile Lake DEM is a subset of the 2m ArcticDEM from the Polar Geospatial Center (see: https://www.pgc.umn.edu/data/arcticdem/). To run the following code cells, put the DEM of your choice in the data/ directory, and change the file names and paths accordingly. Alternatively, explore idealized examples by ignoring this section and altering the grid's topographic elevation field directly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load a DEM\n",
    "with rio.open(\"/data/eightmile-lake-dem-cropped.tif\") as f:\n",
    "    data = f.read(1)\n",
    "\n",
    "    # ArcticDEM tiles are rotated, so let's rotate it back\n",
    "    fine_dem = ndimage.rotate(data[:3280, :4530], 90)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From this point onwards, the computational requirements of almost every step will scale with the size of this DEM. To save time, we can rescale it with typical image processing tools. I'll use xarray's coarsen() method here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Coarsen the DEM\n",
    "xr_dem = xr.DataArray(fine_dem)\n",
    "scale = 10\n",
    "dem = xr_dem.coarsen(dim_0=scale, dim_1=scale).mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When visualizing DEMs, it is often helpful to use matplotlib's LightSource.hillshade(). For a bit of fun, mess with directional lighting using the 'azdeg' argument."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Show the DEM\n",
    "fig, ax = plt.subplots(figsize=(15, 15))\n",
    "\n",
    "ls = LightSource(azdeg=90, altdeg=35)\n",
    "shaded_dem = ls.hillshade(np.array(dem))\n",
    "\n",
    "ax.imshow(shaded_dem, cmap=\"Greys\")\n",
    "plt.axis(\"off\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we need to load the output files from Ku. If you haven't already completed a Ku model run, do that before running the following cells. Here, we're going to focus on depth-dependent diffusion, so we will only need the active layer thickness output."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# First, load the netCDF4 file\n",
    "alt_nc = Dataset(\"./outputs/eightmile-ALT.nc\")\n",
    "\n",
    "# Then, extract a time series of the active layer thickness data\n",
    "alt = alt_nc.variables[\"data\"][:, 0, 0]\n",
    "\n",
    "# Set any negative values to NaN\n",
    "alt[np.where(alt < 0)] = np.nan\n",
    "\n",
    "# Interpolate through the NaN values\n",
    "mask = np.isnan(alt)\n",
    "alt[mask] = np.interp(np.flatnonzero(mask), np.flatnonzero(~mask), alt[~mask])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that, you can ignore the warning that pops up here, code is still working properly below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the active layer thickness dataset\n",
    "plt.plot(alt)\n",
    "plt.xlabel(\"Simulation Year\")\n",
    "plt.ylabel(\"Active layer thickness (m)\")\n",
    "plt.title(\"Ku model output: Active layer thickness\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Landlab components"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is the simplest possible example, only using the Depth-Dependent Diffuser and no other inputs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# First, make the model grid\n",
    "grid = RasterModelGrid(dem.shape)\n",
    "\n",
    "# Add fields for elevation and soil depth\n",
    "z = grid.add_field(\"node\", \"topographic__elevation\", dem.astype(np.float64), units=\"m\")\n",
    "soilH = grid.add_field(\n",
    "    \"node\", \"soil__depth\", np.full(grid.shape, alt[0], dtype=np.float64), units=\"m\"\n",
    ")\n",
    "production = grid.add_zeros(\"soil_production__rate\", at=\"node\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set parameter values\n",
    "linear_diffusivity = 1e-2  # m^2 / yr\n",
    "transport_decay_depth = 1  # m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initialize the Landlab component\n",
    "ddiff = DepthDependentDiffuser(grid, linear_diffusivity, transport_decay_depth)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run the diffusion model\n",
    "\n",
    "# Set the total runtime and time step size\n",
    "nt = 1000  # years\n",
    "dt = 1  # years\n",
    "\n",
    "# Prepare an array for outputs\n",
    "out = np.ndarray([nt, grid.shape[0], grid.shape[1]])\n",
    "\n",
    "# The main loop\n",
    "for i in range(nt):\n",
    "    # If array lengths match, index alt directly with the loop index\n",
    "    if nt == len(alt):\n",
    "        soilH[:] = alt[i]\n",
    "\n",
    "    # Else, rescale the indices to match the length of the ALT array\n",
    "    else:\n",
    "        idx = int((i / nt) * (len(alt)))\n",
    "        soilH[:] = alt[idx]\n",
    "\n",
    "    # Step the diffusion component\n",
    "    ddiff.run_one_step(dt)\n",
    "\n",
    "    # Save out the elevation grid\n",
    "    out[i] = grid.field_values(\"node\", \"topographic__elevation\").reshape(grid.shape)\n",
    "\n",
    "    if i % 100 == 0:\n",
    "        print(\"Completed year \" + str(i))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Visualize the final landscape, after the model completed its run\n",
    "# Note that there will be oddities around the edges, as we did not specify boundary conditions\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(10, 10))\n",
    "\n",
    "ls = LightSource(azdeg=90, altdeg=35)\n",
    "shaded_dem = ls.hillshade(out[-1])\n",
    "\n",
    "ax.imshow(shaded_dem, cmap=\"Greys\")\n",
    "plt.axis(\"off\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Crop away the edges of the DEM\n",
    "cropped = out[0, :-25, 10:] - out[-1, :-25, 10:]\n",
    "\n",
    "# Then, visualize the difference from beginning to end at each grid cell\n",
    "fig, ax = plt.subplots(figsize=(10, 10))\n",
    "# im = ax.imshow(cropped, cmap = 'coolwarm')\n",
    "im = ax.imshow(cropped, cmap=\"jet\")\n",
    "plt.axis(\"off\")\n",
    "fig.colorbar(im, ax=ax, label=\"Elevation change (m)\")\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Ivy",
   "language": "python",
   "name": "ivy"
  },
  "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.9.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}