Isamar_Sergio_workshop.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Depositional mechanics with python\n",
    "\n",
    "\n",
    "Today we will be learning how to use python to simulate a simple 2D sediment deposition/erosion model in a river.\n",
    "\n",
    "The libraries we will be using are `numpy, maplotlib, seaborn, pandas`.\n",
    "\n",
    "The [slides for the presentation are here](https://docs.google.com/presentation/d/12Kx96Auxgq_vJuZCQK21aLqQSgJL3AkN/edit?usp=sharing&ouid=118190230180154301607&rtpof=true&sd=true)\n",
    "\n",
    "The tools we will be needing are:\n",
    "1. A heatmap to show the differences in elevation on the surface \n",
    "2. A vector field representation of the flux on the surface"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "***\n",
    "\n",
    "### Creating a heatmap with seaborn\n",
    "A heatmap is a 2D graph that uses color as its third dimension. An example looks like this"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "heat=[\n",
    "    [0,1,2,3,4],\n",
    "    [5,6,7,8,9],\n",
    "    [10,11,12,13,14]]\n",
    "\n",
    "sns.heatmap(heat)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Question 1 (Class Discussion):\n",
    "\n",
    "What am I defining in the variable \"heat\"?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Using online documentation to edit our heatplot\n",
    "We can play around with the style using seaborn documentation and online examples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.heatmap(\n",
    "    data=heat,\n",
    "    annot=True,  #show the numbers\n",
    "    cmap=\"viridis\",  #use color map viridis\n",
    "    linewidths=1,  #linewidth between blocks\n",
    "    square=True,  # make them squares (not rectangles)\n",
    "     )\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Problem 1:\n",
    "Make a heatmap of size 3 by 3 with a \"hot center\" and a 4 by 4 with a the color map \"cubehelix\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Do it here!\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "***\n",
    "\n",
    "### Formating data for input on heatmap\n",
    "2D plots use one or 3 matrices as input. We have seen what happens if we give it one matrix. If we instead give it 3, it would require: \n",
    "- the $X$ values, \n",
    "- the $Y$ values and \n",
    "- the $heat(X,Y)$ values.\n",
    "\n",
    "Before we do that let's start by making some bigger matrices for $X$ and $Y$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We first define the range of the x and y values\n",
    "xrange=np.arange(0,3,0.5) #What does np.aranage do?\n",
    "yrange=np.arange(0,10,5)\n",
    "\n",
    "#then we make them into matrices\n",
    "X,Y=np.meshgrid(xrange,yrange)\n",
    "\n",
    "print(X)\n",
    "sns.heatmap(X)\n",
    "plt.show()\n",
    "\n",
    "print(Y)\n",
    "sns.heatmap(Y)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have made the coordinates, we will make a function for the heatmap    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#remember we define functions with def\n",
    "def James(xmatrix,ymatrix):\n",
    "    return xmatrix+ymatrix #adds the matrices\n",
    "\n",
    "#we make the heat matrix with the matrices X, Y  and the function\n",
    "Z=James(X,Y)\n",
    "\n",
    "print(Z)\n",
    "sns.heatmap(Z)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "now that good and all, but those axis aren't $X$ and $Y$. To solve this we will use a pandas dataframe and a pivot table."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#we make a dictionary with the X Y and Z values and give it to pandas df\n",
    "Zdf=pd.DataFrame( {\"x\": X.flatten() ,\"y\":Y.flatten(),\"z\": Z.flatten()} );\n",
    "display(Zdf)\n",
    "\n",
    "#make a pivot table because... that's how it is\n",
    "heat=Zdf.pivot_table(columns=\"x\",index=\"y\",values=\"z\")\n",
    "\n",
    "#show\n",
    "sns.heatmap(heat)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Problem 2\n",
    "Make a heat map with the function \n",
    "$$h=elevation(x,y)=100+2x-y$$ \n",
    "on the range from $20<x<30$ and $200<y<300$. (use the grid size you think appropiate)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Here are the steps if you're feeling overwhelm follow this steps\n",
    "#1 make a mesh/grid\n",
    "\n",
    "\n",
    "\n",
    "#2 define your functions and save your heat values\n",
    "\n",
    "\n",
    "\n",
    "#3 make the df and the pivot table\n",
    "\n",
    "\n",
    "\n",
    "#4 display your results! \n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "### Vector fields with quiver\n",
    "So, what if instead you had to put an arrow at every point? a.k.a. a vector field. You need the $(x,y)$ location of the arrow and the $(u,v)$ components of the vectors. \n",
    "\n",
    "Lets start by making the $X$ and $Y$ components as before"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xrange=np.arange(0,10,0.5)\n",
    "yrange=np.arange(0,5,0.5)\n",
    "\n",
    "X,Y=np.meshgrid(xrange,yrange)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "now we define the $U$ and $V$ components of the vector field"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#define functions of the vector at every (x,y) location\n",
    "def Ucomp(xmatrix,ymatrix):\n",
    "    return xmatrix\n",
    "\n",
    "def Vcomp(xmatrix,ymatrix):\n",
    "    return ymatrix*0\n",
    "\n",
    "#assing them\n",
    "U=Ucomp(X,Y)\n",
    "V=Vcomp(X,Y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "finally, make a vector field with the `quiver` function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "axis=plt.quiver(X,Y,U,V)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Problem 3\n",
    "Make a vector field of arrows poiting up. (feel free to pick any range)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Create vector field here\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "### Exner equation: erosion and deposition of sediment\n",
    "\n",
    "Acording to wikipedia: \n",
    "\n",
    "&emsp; \"The Exner equation describes conservation of mass between sediment in the bed of a channel and sediment that is being transported\"\n",
    "\n",
    "$$ \\frac{\\partial \\eta}{\\partial t}=-\\frac{1}{\\varepsilon_0}\\nabla\\cdot q_s$$\n",
    "\n",
    "where $\\eta(x,y)=\\eta$ is the elevation and $q_s(x,y)=(U_q,V_q)$ is the vector field of sediment flux."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will try to undestand deposition and erosion on a flat surface where the sediment flux is changing. We will use our old range and make a new sediment flux"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#define functions of the vector at every (x,y) location\n",
    "def Uflux(xmatrix,ymatrix):\n",
    "    return np.sin(xmatrix)*2\n",
    "    \n",
    "def Vflux(xmatrix,ymatrix):\n",
    "    return ymatrix/3\n",
    "\n",
    "#assing them\n",
    "U=Uflux(X,Y)\n",
    "V=Vflux(X,Y)\n",
    "\n",
    "plt.quiver(X,Y,U,V,scale=30)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's a function that we will be using to apply exner's equation to the flux to find the places where deposition and erosion happen. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Exner function\n",
    "def ExnerCalculator(X,Y,Uq,Vq,eps=1,dt=0.1):\n",
    "    \"\"\"\n",
    "    ExnerCalculator takes in the X and Y locations and the vector field defined by U and V and returns the\n",
    "    deposition/erosion happening at the current time.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    X: matrix \n",
    "        x locations of the flux\n",
    "    Y: matrix \n",
    "        y locations of the flux\n",
    "    Uq: matrix \n",
    "        first components of the flux\n",
    "    Vq: matrix \n",
    "        second components of the flux\n",
    "    eps: float\n",
    "        1-porosity of the bed\n",
    "    dt: float\n",
    "        time delta for the method\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    Xs: matrix \n",
    "        x locations with n-1 locations corresponding to dh\n",
    "    Ys: matrix \n",
    "        y locations with n-1 locations corresponding to dh\n",
    "    dh: matrix\n",
    "        The difference in elevation in dt time units based on the flux\n",
    "    \"\"\"\n",
    "    m,n=X.shape\n",
    "    dx=X[0][1]-X[0][0]\n",
    "    dy=Y[1][0]-Y[0][0]\n",
    "\n",
    "    #gradient\n",
    "    dUq=[ [Uq[j][i+1]-Uq[j][i] for i in range(n-1) ] for j in range(m-1) ]\n",
    "    dVq=[ [Vq[j+1][i]-Vq[j][i] for i in range(n-1) ] for j in range(m-1) ]\n",
    "    \n",
    "    #exner\n",
    "    dh=(-dt/eps)*(dUq/dx+dVq/dy)\n",
    "\n",
    "    #new Xs and Ys\n",
    "    Xs=np.array([ [X[j][i] for i in range(n-1) ] for j in range(m-1) ])\n",
    "    Ys=np.array([ [Y[j][i] for i in range(n-1) ] for j in range(m-1) ])\n",
    "\n",
    "    return Xs,Ys,dh"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And now we will pass our vector field to this function and plot the erosion/deposition heatmap $\\Delta \\eta$ "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#use the function on the flux\n",
    "newx,newy,deltah=ExnerCalculator(X,Y,U,V)\n",
    "\n",
    "#plot the heatmap\n",
    "erosion=pd.DataFrame( {\"x\": newx.flatten() ,\"y\":newy.flatten(),\"erosion\": deltah.flatten()} )\n",
    "erosion_table=erosion.pivot_table(columns=\"x\",index=\"y\",values=\"erosion\")\n",
    "sns.heatmap(erosion_table,cmap=\"RdYlBu\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's see it next to the quiver again"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.quiver(X,Y,U,V)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "## Homework\n",
    "1. Plot the vector field $q=(U,V)$ in the domain from $-1<x<1$ and $-1<y<1$ given by\n",
    "     $$U=-X$$\n",
    "     $$V=-Y$$\n",
    "\n",
    "2. Use the `ExnerCalculator` function to display an erosional/deposition heatmap for the change in elevation $\\Delta \\eta$\n",
    "\n",
    "3. Given the data of the three variables $weight, height$ and $age$ below, find the correlation matrix between them and make a heatmap of the correlation matrix. (check the correlation example in __[here](https://python-graph-gallery.com/90-heatmaps-with-various-input-format)__  )\n",
    "\n",
    "4. Create a repository and put it up on your github."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "weight=[162.22, 156.12, 169.31, 179.76, 173.12, 165.66, 166.61, 170.91, 156.67, 170.94]\n",
    "height=[7.8954, 7.7507, 7.7312, 7.9815, 7.8738, 7.6802, 7.7422, 7.8297, 7.8992, 7.9033]\n",
    "age=[32, 33, 22, 45, 33, 49, 36, 20, 41, 25]"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "135e24ee289c60ba69d5004678bf472a320b24bde767457479ce4e07bf7a1633"
  },
  "kernelspec": {
   "display_name": "Python 3",
   "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.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}