nsw-2016-storm-impact/notebooks/04b_profile_picker.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Beach profile classifier\n",
    "This notebook is a set of functions that tries to pick a dune crest and dune toe given a beach profile."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup notebook"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Enable autoreloading of our modules. \n",
    "# Most of the code will be located in the /src/ folder, \n",
    "# and then called from the notebook.\n",
    "%matplotlib inline\n",
    "%reload_ext autoreload\n",
    "%autoreload"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.core.debugger import set_trace\n",
    "\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import os\n",
    "import decimal\n",
    "import plotly\n",
    "import plotly.graph_objs as go\n",
    "import plotly.plotly as py\n",
    "import plotly.tools as tls\n",
    "import plotly.figure_factory as ff\n",
    "from plotly import tools\n",
    "import plotly.io as pio\n",
    "from scipy import stats\n",
    "import math\n",
    "import matplotlib\n",
    "from matplotlib import cm\n",
    "import colorlover as cl\n",
    "import numpy.ma as ma\n",
    "\n",
    "from ipywidgets import widgets, Output\n",
    "from IPython.display import display, clear_output, Image, HTML\n",
    "\n",
    "from sklearn.metrics import confusion_matrix\n",
    "\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "from sklearn import linear_model, datasets\n",
    "\n",
    "from scipy.interpolate import UnivariateSpline\n",
    "from scipy.interpolate import interp1d\n",
    "from scipy.interpolate import splrep, splev\n",
    "from scipy.integrate import simps\n",
    "from scipy.stats import linregress\n",
    "from scipy.signal import find_peaks\n",
    "import json"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Matplot lib default settings\n",
    "plt.rcParams[\"figure.figsize\"] = (10,3)\n",
    "plt.rcParams['axes.grid']=True\n",
    "plt.rcParams['grid.alpha'] = 0.5\n",
    "plt.rcParams['grid.color'] = \"grey\"\n",
    "plt.rcParams['grid.linestyle'] = \"--\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Import data\n",
    "Let's first import data from our pre-processed interim data folder."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def df_from_csv(csv, index_col, data_folder='../data/interim'):\n",
    "    print('Importing {}'.format(csv))\n",
    "    return pd.read_csv(os.path.join(data_folder,csv), index_col=index_col)\n",
    "\n",
    "df_waves = df_from_csv('waves.csv', index_col=[0, 1])\n",
    "df_tides = df_from_csv('tides.csv', index_col=[0, 1])\n",
    "df_profiles = df_from_csv('profiles.csv', index_col=[0, 1, 2])\n",
    "df_sites = df_from_csv('sites.csv', index_col=[0])\n",
    "df_profile_features_crest_toes = df_from_csv('profile_features_crest_toes.csv', index_col=[0,1])\n",
    "\n",
    "# Note that the forecasted data sets should be in the same order for impacts and twls\n",
    "impacts = {\n",
    "    'forecasted': {\n",
    "        'foreshore_slope_sto06': df_from_csv('impacts_forecasted_foreshore_slope_sto06.csv', index_col=[0]),\n",
    "        'mean_slope_sto06': df_from_csv('impacts_forecasted_mean_slope_sto06.csv', index_col=[0]),\n",
    "        },\n",
    "    'observed': df_from_csv('impacts_observed.csv', index_col=[0])\n",
    "    }\n",
    "\n",
    "\n",
    "twls = {\n",
    "    'forecasted': {\n",
    "        'foreshore_slope_sto06': df_from_csv('twl_foreshore_slope_sto06.csv', index_col=[0, 1]),\n",
    "        'mean_slope_sto06':df_from_csv('twl_mean_slope_sto06.csv', index_col=[0, 1]),\n",
    "    }\n",
    "}\n",
    "print('Done!')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Implement algorithm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Get profile data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_profile_data(site, profile_type, df_profiles):\n",
    "    \"\"\"\n",
    "    Returns a list of x and z coordinates for a given site and profile_type from our list of profiles\n",
    "    \"\"\"\n",
    "\n",
    "    # Get x and z coorindates of profile\n",
    "    x = df_profiles.loc[(site, profile_type)].index.values\n",
    "    z = df_profiles.loc[(site, profile_type)].z.values\n",
    "\n",
    "    # Remove nan values\n",
    "    m = ma.masked_invalid(z)\n",
    "    x = x[~m.mask]\n",
    "    z = z[~m.mask]\n",
    "\n",
    "    #     # Remove landwards of landlim\n",
    "    #     land_lim = df_profiles.loc[(site, 'poststorm')].land_lim.loc[0]\n",
    "    #     mask = ma.masked_where(x < land_lim, x)\n",
    "    #     x = x[~mask.mask]\n",
    "    #     z = z[~mask.mask]\n",
    "\n",
    "    return x, z\n",
    "\n",
    "\n",
    "site = 'NARRA0001'\n",
    "profile_type = 'prestorm'\n",
    "x, z = get_profile_data(site, profile_type, df_profiles)\n",
    "\n",
    "print('x = {} ...'.format(x[0:5]))\n",
    "print('z = {} ...'.format(z[0:5]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fit linear interpolated spline to profile\n",
    "- Linear interpolated spline used to simplify the profile. Probably not necessary for picking dune crests/toes, but handy to simplify the profile into a list of linear lines. This helps pick up berms (for later).\n",
    "- Need to do a bit of optimization to calculate the best smoothing parameter for the profile."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0,
     7,
     16,
     21,
     26,
     38,
     88,
     110
    ]
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "def pretty_print(l):\n",
    "    print(json.dumps(l, indent=2, sort_keys=True))\n",
    "\n",
    "def pairwise(iterable):\n",
    "    \"s -> (s0,s1), (s1,s2), (s2, s3), ...\"\n",
    "    a, b = itertools.tee(iterable)\n",
    "    next(b, None)\n",
    "    return zip(a, b)\n",
    "\n",
    "\n",
    "def xz(x_start, x_end, x, z):\n",
    "    \"\"\"\n",
    "    Returns a list of x coordinates and z coorindates which lie between x_start and x_end\n",
    "    \"\"\"\n",
    "    x_seg = [val for val in x if x_start <= val <= x_end]\n",
    "    z_seg = [z_val for x_val, z_val in zip(x, z) if x_start <= x_val <= x_end]\n",
    "    return x_seg, z_seg\n",
    "\n",
    "\n",
    "def get_slope(x, z):\n",
    "    slope, intercept, r_value, p_value, std_err = linregress(x, z)\n",
    "    return slope\n",
    "\n",
    "\n",
    "def get_r_squared(x, z):\n",
    "    slope, intercept, r_value, p_value, std_err = linregress(x, z)\n",
    "    return r_value**2\n",
    "\n",
    "\n",
    "def rolling_std(a, w):\n",
    "    \"\"\"\n",
    "    Calculates a rolling window standard deviation\n",
    "    https://stackoverflow.com/a/40773366\n",
    "    \"\"\"\n",
    "    a = np.array(a)\n",
    "    nrows = a.size - w + 1\n",
    "    n = a.strides[0]\n",
    "    a2D = np.lib.stride_tricks.as_strided(a, shape=(nrows, w), strides=(n, n))\n",
    "    return np.std(a2D, 1)\n",
    "\n",
    "\n",
    "def iterate_over_smoothing_params(x, z):\n",
    "    \"\"\"\n",
    "    Finds the best smoothing parameter for a linear interpolated spline. \n",
    "    Checks a number of smoothing factors, and then tries to minimize the \n",
    "    number of knots (changepoints) while retaining decent match to the \n",
    "    original profile.\n",
    "    \"\"\"\n",
    "\n",
    "    # Store the results of our iterations in a list of dictionaries\n",
    "    results = []\n",
    "\n",
    "    # Iterate through different value of smoothing parameters\n",
    "    s_vals = [x for x in np.arange(0, 10, 0.01)]\n",
    "\n",
    "    r2_vals = []\n",
    "    n_knots = []\n",
    "    s_vals_to_keep = []\n",
    "\n",
    "    # Iterate backwards just to make further calcultions more easier.\n",
    "    for s_val in s_vals[::-1]:\n",
    "\n",
    "        # Fit spline to profile using the s_val for this iteration\n",
    "        # k=1 is used to force linear line segments between knots\n",
    "        s = UnivariateSpline(x, z, s=s_val, k=1)\n",
    "        xs = np.linspace(x[0], x[-1], 1000)\n",
    "        zs = s(xs)\n",
    "\n",
    "        # Find knots\n",
    "        x_knots = s.get_knots()\n",
    "        z_knots = s(x_knots)\n",
    "\n",
    "        # Get number of knots\n",
    "        n = len(s.get_knots())\n",
    "\n",
    "        # If we've already recorded data about this number of knots, just skip.\n",
    "        # Need to also, only have increasing results for spline interpolation\n",
    "        if n in [x['n'] for x in results\n",
    "                 ] or n < max([x['n'] for x in results], default=0):\n",
    "            continue\n",
    "\n",
    "        # Get r2, how well does the simplification fit the original profile?\n",
    "        # Create interp model from the knots at our original points.\n",
    "        f = interp1d(x_knots, z_knots)\n",
    "        z_lin = f(x)\n",
    "        r2 = get_r_squared(z, z_lin)\n",
    "\n",
    "        results.append({'n': n, 's': s_val, 'r2': r2})\n",
    "    return results\n",
    "\n",
    "\n",
    "def find_best_smoothing_param(results, min_std=0.0002):\n",
    "    \"\"\"\n",
    "    Given a list of smoothing parameters and their fits, determine what the best smoothing parameter is.\n",
    "    Larger min_std = more smoothing\n",
    "    \"\"\"\n",
    "\n",
    "    # Get the 2nd derivate of the n_knots vs r2_vals.\n",
    "    y_spl = UnivariateSpline([x['n'] for x in results],\n",
    "                             [x['r2'] for x in results],\n",
    "                             s=0)\n",
    "    y_spl_2d = y_spl.derivative(n=2)\n",
    "    y_spl_2d_vals = y_spl_2d([x['n'] for x in results])\n",
    "\n",
    "    # Get a rolling standard deviation of the derivative\n",
    "    std = rolling_std(y_spl_2d_vals, w=3)\n",
    "\n",
    "    # Best smoothing parameter will be when the standard deviation stops changing\n",
    "    best_i = np.argmax(std < min_std)\n",
    "    best_s = results[best_i]['s']\n",
    "    return best_s\n",
    "\n",
    "\n",
    "def simplify_profile(x, z, site):\n",
    "    results = iterate_over_smoothing_params(x, z)\n",
    "    best_s = find_best_smoothing_param(results)\n",
    "\n",
    "    # Fit spline to profile\n",
    "    s = UnivariateSpline(x, z, s=best_s, k=1)\n",
    "    xs = np.linspace(x[0], x[-1], 1000)\n",
    "    zs = s(xs)\n",
    "\n",
    "    # Find knots\n",
    "    x_knots = s.get_knots()\n",
    "    z_knots = s(x_knots)\n",
    "\n",
    "#     # Plot for checking\n",
    "#     plt.title('{}: Linear spline simplification'.format(site))\n",
    "#     plt.xlabel('Distance (m)')\n",
    "#     plt.ylabel('Elevation (m AHD)')\n",
    "#     plt.plot(x, z, 'k.-', label='Raw profile')\n",
    "#     plt.plot(\n",
    "#         xs, zs, 'r-', label='Interpolated profile (s={:.2f})'.format(best_s))\n",
    "#     plt.plot(\n",
    "#         x_knots,\n",
    "#         z_knots,\n",
    "#         'ro',\n",
    "#         markersize=8,\n",
    "#         markeredgewidth=1.5,\n",
    "#         markeredgecolor='orange',\n",
    "#         label='Interpolated knots (n={})'.format(len(x_knots)))\n",
    "#     plt.legend(loc='best')\n",
    "#     plt.show()\n",
    "\n",
    "    return x_knots, z_knots, best_s\n",
    "\n",
    "\n",
    "x_knots, z_knots, best_s = simplify_profile(x, z, site)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": []
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "def get_x_z_from_knots(x_knots, z_knots):\n",
    "    # Fit spline to profile\n",
    "    s = UnivariateSpline(x_knots, z_knots, s=0, k=1)\n",
    "    xs = np.linspace(x_knots[0], x_knots[-1], 1000)\n",
    "    zs = s(xs)\n",
    "    return xs, zs\n",
    "\n",
    "def simplify_segments_with_same_slopes(x_knots, x, z, site, best_s):\n",
    "    \"\"\"\n",
    "    Removes knots which have similar slopes on either side. \n",
    "    This is an extra simplifcation as the linear splines may still include\n",
    "    some nodes that aren't required.\n",
    "    \"\"\"\n",
    "    removed_x_knots = []\n",
    "    while True:\n",
    "        for x_knot1, x_knot2, x_knot3 in zip(x_knots[0:-2], x_knots[1:-1], x_knots[2:]):\n",
    "\n",
    "            # Get slope of segment behind point at x_knot2\n",
    "            x_seg, z_seg = xz(x_knot1, x_knot2, x, z)\n",
    "            slope1 = get_slope(x_seg, z_seg)\n",
    "\n",
    "            # Get slope of line in front of point at x_knot2\n",
    "            x_seg, z_seg = xz(x_knot2, x_knot3, x, z)\n",
    "            slope2 = get_slope(x_seg, z_seg)\n",
    "\n",
    "            slope_diff = abs(slope1 - slope2)\n",
    "\n",
    "            # Also check if points are too close together\n",
    "            if x_knot3 - x_knot1 < 3: # m\n",
    "                too_close = True\n",
    "            else:\n",
    "                too_close = False\n",
    "            \n",
    "            # Good, the slopes are different on each side, keep checking\n",
    "            if slope_diff > 0.01 and too_close == False:\n",
    "                continue\n",
    "            # Else bad, slopes are the same so we have to remove this point\n",
    "            else:\n",
    "#                 print('Knot at x={} removed (slope_diff={:.3f})'.format(\n",
    "#                     x_knot2, slope_diff))\n",
    "                x_knots = np.delete(x_knots, np.where(x_knots == x_knot2), axis=0)\n",
    "                removed_x_knots.append(x_knot2)\n",
    "                break\n",
    "\n",
    "        else:\n",
    "#             print(\"All segments have different slopes\")\n",
    "            x_knots_simplified = x_knots\n",
    "            z_knots_simplified = [\n",
    "                z_val for x_val, z_val in zip(x, z) if x_val in x_knots_simplified\n",
    "            ]\n",
    "            break\n",
    "    \n",
    "    # Find z location of our removed knots\n",
    "    removed_z_knots = [\n",
    "                z_val for x_val, z_val in zip(x, z) if x_val in removed_x_knots\n",
    "            ]\n",
    "    \n",
    "    # Get the new interpolated spline from our simplified knots\n",
    "    xs, zs = get_x_z_from_knots(x_knots_simplified, z_knots_simplified)\n",
    "    \n",
    "    # Plot for checking\n",
    "#     plt.title('{}: Linear spline simplification by comparing slopes'.format(site))\n",
    "#     plt.xlabel('Distance (m)')\n",
    "#     plt.ylabel('Elevation (m AHD)')\n",
    "#     plt.plot(x, z, 'k.-', label='Raw profile')\n",
    "#     plt.plot(\n",
    "#         xs, zs, 'r-', label='Interpolated profile (s={:.2f})'.format(best_s))\n",
    "#     plt.plot(\n",
    "#         x_knots_simplified,\n",
    "#         z_knots_simplified,\n",
    "#         'ro',\n",
    "#         markersize=8,\n",
    "#         markeredgewidth=1.5,\n",
    "#         markeredgecolor='orange',\n",
    "#         label='Interpolated knots (n={})'.format(len(x_knots)))\n",
    "#     plt.plot(\n",
    "#         removed_x_knots,\n",
    "#         removed_z_knots,\n",
    "#         'ro',\n",
    "#         markersize=10,\n",
    "#         markeredgewidth=2.5,\n",
    "#         markerfacecolor=\"None\",\n",
    "#         markeredgecolor='orange',\n",
    "#         label='Removed knots (n={})'.format(len(removed_x_knots)))\n",
    "#     plt.legend(loc='best')\n",
    "#     plt.show()\n",
    "\n",
    "    return x_knots_simplified,z_knots_simplified\n",
    "\n",
    "x_knots, z_knots = simplify_segments_with_same_slopes(x_knots, x, z, site, best_s)\n",
    "xs, zs = get_x_z_from_knots(x_knots, z_knots)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0,
     91
    ]
   },
   "outputs": [],
   "source": [
    "%matplotlib inline \n",
    "\n",
    "def dunes_from_d1(xx,zz,zz_d1,\n",
    "                       min_dune_face_slope=0.05,\n",
    "                       min_dune_slope_prominence=0.00):\n",
    "    \"\"\"\n",
    "    Get peaks by looking at the derivates (slope)\n",
    "    \"\"\"\n",
    "\n",
    "    # Find the peaks in the first derivatives (i.e. peaks in slopes)\n",
    "    peaks, props = find_peaks(\n",
    "        zz_d1 * -1,\n",
    "#         height=min_dune_face_slope,\n",
    "        height=(None, None),\n",
    "        width=(None, None),\n",
    "        prominence=(None, None)\n",
    "    )\n",
    "    \n",
    "#     print(peaks)\n",
    "#     plt.figure(figsize=(10, 3))\n",
    "#     plt.plot(xx, zz_d1*-1, label='raw')\n",
    "#     plt.show()\n",
    "    \n",
    "    # If no peaks found, return None\n",
    "    if len(peaks) == 0:\n",
    "        print('No dunes from d1')\n",
    "        return []\n",
    "\n",
    "    # The peaks we find are going to be near the center of any face.\n",
    "    # We need to get the left and right bases to get the peaks and crests.\n",
    "    x_crests = [xx[val] for val in props['left_bases']]\n",
    "    x_toes = [xx[val] for val in props['right_bases']]\n",
    "    heights = [val for val in props['peak_heights']]\n",
    "    x_centers = xx[peaks]\n",
    "    z_crests = []\n",
    "    z_toes = []\n",
    "    dune_face_slopes = []\n",
    "    dune_heights = []\n",
    "\n",
    "    # Let's look at each peak we identified.\n",
    "    for n, (x_crest, x_center, x_toe, height) in enumerate(zip(x_crests, x_centers, \n",
    "                                                x_toes,heights)):\n",
    "        \n",
    "        # Once the slope reduces to a certain threshold, consider this the crest/toe\n",
    "        min_slope = height / -4\n",
    "        \n",
    "        # Check the land side of the dune face\n",
    "        # Find locations where slope is greater than our minimum slope\n",
    "        land_mask = ((x_crest <= xx) & (xx <= x_center))\n",
    "        land_slope_okay = zz_d1[land_mask] < min_slope\n",
    "        x_land_slope_okay = xx[land_mask][~land_slope_okay]\n",
    "        \n",
    "        # Get last location where slope satisfies the criteria\n",
    "        if len(x_land_slope_okay) != 0:\n",
    "            x_crest = x_land_slope_okay[-1]\n",
    "            x_crests[n] = x_crest\n",
    "            \n",
    "        # Get corresponding z_value for crest\n",
    "        idx_crest = np.where(x_crest == xx)[0]\n",
    "        z_crest = zz[idx_crest][0]\n",
    "        z_crests.append(z_crest)\n",
    "\n",
    "        # Check the sea side of the dune face\n",
    "        # Find locations where slope is greater than our minimum slope\n",
    "        sea_mask = ((x_center <= xx) & (xx <= x_toe))\n",
    "        sea_slope_okay = zz_d1[sea_mask] < min_slope\n",
    "        x_sea_slope_okay = xx[sea_mask][~sea_slope_okay]\n",
    "        \n",
    "        # Get first location where slope satisfies the criteria\n",
    "        if len(x_sea_slope_okay) != 0:\n",
    "            x_toe = x_sea_slope_okay[0]\n",
    "            x_toes[n] = x_toe\n",
    "            \n",
    "        # Get corresponding z_value\n",
    "        idx_toe = np.where(x_toe == xx)[0]\n",
    "        z_toe = zz[idx_toe][0]\n",
    "        z_toes.append(z_toe)\n",
    "\n",
    "        # Recalculate new slope and dune height\n",
    "        dune_face_slopes.append((z_crest - z_toe) / (x_crest - x_toe))\n",
    "        dune_heights.append(z_crest - z_toe)\n",
    "\n",
    "    dunes = [\n",
    "        {'x_crest': x_crest,\n",
    "         'x_toe': x_toe,\n",
    "         'z_crest': z_crest,\n",
    "         'z_toe': z_toe,\n",
    "         'dune_face_slope': dune_face_slope,\n",
    "         'dune_height': dune_height} for x_crest, x_toe,z_crest,z_toe,dune_face_slope,dune_height in zip(x_crests, x_toes,z_crests,z_toes,dune_face_slopes,dune_heights)\n",
    "    ]\n",
    "\n",
    "    return dunes\n",
    "    \n",
    "def dunes_from_z_max(xx, zz, zz_d1):\n",
    "    \"\"\"\n",
    "    Estimate dune crest by from the maximum elevation of the profile.\n",
    "    Toe is estimated by based on d1 (slope)\n",
    "    Works when there is no clearly defined dune, and always gives us a dune\n",
    "    \"\"\"\n",
    "    \n",
    "    # Get location of maximum elevation. This is the crest\n",
    "    i_crest = np.argmax(zz)\n",
    "    x_crest = xx[i_crest]\n",
    "    z_crest = zz[i_crest]\n",
    "    \n",
    "    # Find peaks in slope seaward of the slope\n",
    "    peaks, props = find_peaks(\n",
    "        zz_d1 * -1,\n",
    "        height=0,\n",
    "        width=0,\n",
    "        prominence=0)\n",
    "    \n",
    "    # Check which peaks are landward of the crest.\n",
    "    # Use right base of peak prominence since we want the toe (seaward side of the dune face)\n",
    "    mask = [True if x > i_crest else False for x in peaks]\n",
    "    i_seaward_slope_peaks = peaks[mask]\n",
    "    i_right_bases = props['right_bases'][mask]\n",
    "    heights = props['peak_heights'][mask]\n",
    "                             \n",
    "    if len(i_seaward_slope_peaks) != 0:\n",
    "        i_seaward_slope_peak = i_seaward_slope_peaks[0]\n",
    "        i_right_base = i_right_bases[0]\n",
    "        height = heights[0]\n",
    "        \n",
    "        x_center = xx[i_seaward_slope_peak]\n",
    "        x_toe = xx[i_right_base]\n",
    "        z_toe = zz[i_right_base]\n",
    "        \n",
    "        # Once the slope reduces to a certain threshold, consider this the /toe\n",
    "        min_slope = height / -4\n",
    "        \n",
    "        # Check the sea side of the dune face\n",
    "        # Find locations where slope is greater than our minimum slope\n",
    "        sea_mask = ((x_center <= xx) & (xx <= x_toe))\n",
    "        sea_slope_okay = zz_d1[sea_mask] < min_slope\n",
    "        x_sea_slope_okay = xx[sea_mask][~sea_slope_okay]\n",
    "        \n",
    "        # Get first location where slope satisfies the criteria\n",
    "        if len(x_sea_slope_okay) != 0:\n",
    "            x_toe = x_sea_slope_okay[0]\n",
    "            \n",
    "        # Get corresponding z_value\n",
    "        idx_toe = np.where(x_toe == xx)[0]\n",
    "        z_toe = zz[idx_toe][0]\n",
    "                \n",
    "        dune_face_slope = (z_crest - z_toe) / (x_crest - x_toe)\n",
    "        dune_height = z_crest - z_toe           \n",
    "\n",
    "    else:\n",
    "        # No peaks were found...\n",
    "        \n",
    "        # Find position of maximum slope seaward of crest\n",
    "        min_seaward_slope = np.min(zz_d1[i_crest:])\n",
    "        i_min_seaward_slope = np.argmin(zz_d1[i_crest:])\n",
    "\n",
    "        # Check if this slope reduces from a certain threshold\n",
    "        min_threshold_slope = min_seaward_slope / 6\n",
    "\n",
    "        i_possible_toes = np.where(zz_d1[i_min_seaward_slope:] < min_threshold_slope)\n",
    "\n",
    "        if len(i_possible_toes) != 0:\n",
    "            x_toe = xx[i_min_seaward_slope:][i_possible_toes][0]\n",
    "            z_toe = zz[i_min_seaward_slope:][i_possible_toes][0]\n",
    "            dune_face_slope = (z_crest - z_toe) / (x_crest - x_toe)\n",
    "            dune_height = z_crest - z_toe\n",
    "        else:\n",
    "            x_toe = np.nan\n",
    "            z_toe = np.nan\n",
    "            dune_face_slope = np.nan\n",
    "            dune_height = np.nan\n",
    "    \n",
    "    return [{'x_crest': x_crest,\n",
    "             'x_toe': x_toe,\n",
    "             'z_crest': z_crest,\n",
    "             'z_toe': z_toe,\n",
    "             'dune_face_slope': dune_face_slope,\n",
    "             'dune_height': dune_height}]\n",
    "    \n",
    "def find_crest_and_toes(x, z, x_knots, best_s):\n",
    "\n",
    "    # Get a cubic spline so we can have continuous second order deriviates for our profile\n",
    "    spl = splrep(x, z, k=3, s=best_s)\n",
    "\n",
    "    # Sample interpolated spline and calcualte derivaives.\n",
    "    xx = np.linspace(x[0], x[-1], 1000)\n",
    "    zz = splev(xx, spl)\n",
    "    zz_d1 = splev(xx, spl, der=1)\n",
    "    zz_d2 = splev(xx, spl, der=2)\n",
    "    zz_d3 = splev(xx, spl, der=3)\n",
    "\n",
    "    # Get potential dunes from different methods\n",
    "    dunes_d1 = dunes_from_d1(xx,zz,zz_d1,\n",
    "                                         min_dune_face_slope=0.15,\n",
    "                                         min_dune_slope_prominence=0.15)\n",
    "    # Dunes from maximum elevation\n",
    "    dunes_z_max = dunes_from_z_max(xx, zz, zz_d1)\n",
    "\n",
    "    # Join list of all potential dunes\n",
    "    dunes = dunes_d1 + dunes_z_max\n",
    "    \n",
    "    # Remove dunes which have the same crest and toe x location\n",
    "    dunes = [x  for x in dunes if x['x_crest']!=x['x_toe']]\n",
    "    \n",
    "    # Calculate area under d1 curve. This'll give an indication of how important the slope is\n",
    "    for dune in dunes:\n",
    "        mask = [True if (dune['x_crest'] < x_val) and (x_val < dune['x_toe']) else False for x_val in xx]\n",
    "        x_seg = xx[mask]\n",
    "        zz_d1_seg = zz_d1[mask]\n",
    "        d1_area = simps(zz_d1_seg, x_seg)\n",
    "        dune['d1_area'] = d1_area\n",
    "\n",
    "    # Create a filtered flag for all our dunes\n",
    "    dunes = [dict(item, **{'filtered':False}) for item in dunes]\n",
    "\n",
    "    # Iterate through each potential dune and decide whether to filter or not\n",
    "    for dune in dunes:\n",
    "        \n",
    "        if dune['z_crest'] < np.median(zz)/2:\n",
    "            dune['filtered'] = True\n",
    "            \n",
    "        elif dune['dune_height'] < 0.5: # m\n",
    "            dune['filtered'] = True\n",
    "    \n",
    "        # Filter out dunes which occur landward of the highest elevation point\n",
    "    \n",
    "    \n",
    "        # TODO Decide to remove the dune toe\n",
    "        # If mean slope and std between dune toe and dune crests is similar to between dune crest and end of profile,\n",
    "        # We can remove the toe\n",
    "        \n",
    "        \n",
    "        \n",
    "#     pretty_print(dunes)\n",
    "    return xx, zz, zz_d1, zz_d2, zz_d3, dunes\n",
    "\n",
    "\n",
    "\n",
    "def plot_profile(site, x, z, xx, zz, zz_d1, zz_d2, zz_d3, dunes):\n",
    "    \n",
    "    plt.figure(figsize=(10, 12))\n",
    "    plt.subplot(411)\n",
    "    plt.plot(x, z, label='raw')\n",
    "    plt.plot(xx, zz, label='interpolated')\n",
    "    plt.plot(\n",
    "        [x['x_crest'] for x in dunes if x['filtered']==True],\n",
    "        [x['z_crest'] for x in dunes if x['filtered']==True],\n",
    "        'rv',\n",
    "        markersize=10,\n",
    "        fillstyle='none',\n",
    "        label='Dune crest (potential)')\n",
    "    plt.plot(\n",
    "        [x['x_toe'] for x in dunes if x['filtered']==True],\n",
    "        [x['z_toe'] for x in dunes if x['filtered']==True],\n",
    "        'r^',\n",
    "        markersize=10,\n",
    "        fillstyle='none',\n",
    "        label='Dune toe (potential)')\n",
    "    plt.plot(\n",
    "        [x['x_crest'] for x in dunes if x['filtered']==False],\n",
    "        [x['z_crest'] for x in dunes if x['filtered']==False],\n",
    "        'rv',\n",
    "        markersize=10,\n",
    "#         label='Dune crest (x={:.2f} m, z={:.2f} m)'.format(x_crest, z_crest)\n",
    "    )\n",
    "    plt.plot(\n",
    "        [x['x_toe'] for x in dunes if x['filtered']==False],\n",
    "        [x['z_toe'] for x in dunes if x['filtered']==False],\n",
    "        'r^',\n",
    "        markersize=10,\n",
    "#         label='Dune toe (x={:.2f} m, z={:.2f} m)'.format(x_toe, z_toe)\n",
    "    )\n",
    "    plt.title('{}: elevations'.format(site))\n",
    "    plt.legend(loc='best')\n",
    "\n",
    "    plt.subplot(412)\n",
    "    d1 = np.diff(zz) / np.diff(xx)\n",
    "    plt.plot(xx, zz_d1)\n",
    "#     plt.axvline(x_crest, color='red', label='Dune crest')\n",
    "#     plt.axvline(x_toe, color='red', label='Dune toe')\n",
    "    # plt.ylim([-0.2,0.2])\n",
    "    plt.title('first derivative - slope')\n",
    "\n",
    "    plt.subplot(413)\n",
    "    d2 = np.diff(d1) / np.diff(xx[1:])\n",
    "    plt.plot(xx, zz_d2)\n",
    "#     plt.axvline(x_crest, color='red', label='Dune crest')\n",
    "#     plt.axvline(x_toe, color='red', label='Dune toe')\n",
    "    plt.title('second derivative - change in slope')\n",
    "    #     plt.ylim([-0.025,0.025])\n",
    "\n",
    "    plt.subplot(414)\n",
    "    plt.plot(xx, zz_d3)\n",
    "#     plt.axvline(x_crest, color='red', label='Dune crest')\n",
    "#     plt.axvline(x_toe, color='red', label='Dune toe')\n",
    "    plt.title('third derivative')\n",
    "    #     plt.ylim([-0.0025,0.0025])\n",
    "\n",
    "    plt.tight_layout()\n",
    "    plt.show()\n",
    "\n",
    "\n",
    "xx, zz, zz_d1, zz_d2, zz_d3, dunes = find_crest_and_toes(x,z,x_knots,best_s)\n",
    "\n",
    "plot_profile(site, x,z, xx, zz, zz_d1, zz_d2, zz_d3, dunes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": []
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "def get_crests_and_toes(site, profile_type, df_profiles, debug_plot=False):\n",
    "    x, z = get_profile_data(site, profile_type, df_profiles)\n",
    "    x_knots, z_knots, best_s = simplify_profile(x, z, site)\n",
    "    x_knots, z_knots = simplify_segments_with_same_slopes(x_knots, x, z, site, best_s)\n",
    "    xx, zz, zz_d1, zz_d2, zz_d3, dunes = find_crest_and_toes(x,z,x_knots,best_s)\n",
    "\n",
    "    if debug_plot:\n",
    "        plot_profile(site, x,z, xx, zz, zz_d1, zz_d2, zz_d3, dunes)\n",
    "    \n",
    "    profile_data = {\n",
    "        'site': site,\n",
    "        'profile_type': profile_type,\n",
    "        'xx': xx,\n",
    "        'zz': zz,\n",
    "        'zz_d1': zz_d1,\n",
    "        'zz_d2': zz_d2,\n",
    "        'zz_d3': zz_d3,\n",
    "        'dunes': dunes,\n",
    "        'x': x,\n",
    "        'z': z\n",
    "    }\n",
    "    \n",
    "    return profile_data\n",
    "    \n",
    "    \n",
    "import random\n",
    "site = random.sample(df_profiles.index.get_level_values('site_id').unique().tolist(), 1)[0]\n",
    "\n",
    "print(site)\n",
    "# site='NINEMs0004'\n",
    "profile_data = get_crests_and_toes(site,'prestorm', df_profiles, debug_plot=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sites = df_profiles.index.get_level_values('site_id').unique().tolist()[:10]\n",
    "profile_data = [get_crests_and_toes(site, 'prestorm', df_profiles) for site in sites]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     2,
     11
    ]
   },
   "outputs": [],
   "source": [
    "%matplotlib notebook\n",
    "\n",
    "def closest_node(node, nodes):\n",
    "    \"\"\"\n",
    "    https://codereview.stackexchange.com/a/28210\n",
    "    \"\"\"\n",
    "    nodes = np.asarray(nodes)\n",
    "    deltas = nodes - node\n",
    "    dist_2 = np.einsum('ij,ij->i', deltas, deltas)\n",
    "    return np.argmin(dist_2)\n",
    "\n",
    "def get_data(i, profile_data):\n",
    "    site = profile_data[i]['site']\n",
    "    x = profile_data[i]['x']\n",
    "    z = profile_data[i]['z']\n",
    "    xx = profile_data[i]['xx']\n",
    "    zz = profile_data[i]['zz']\n",
    "    zz_d1 = profile_data[i]['zz_d1']\n",
    "    zz_d2 = profile_data[i]['zz_d2']\n",
    "    zz_d3 = profile_data[i]['zz_d3']\n",
    "    dunes = profile_data[i]['dunes']\n",
    "    \n",
    "    return site, x, z, xx, zz,zz_d1,zz_d2,zz_d3, dunes\n",
    "\n",
    "i = 0\n",
    "site, x, z, xx, zz,zz_d1,zz_d2,zz_d3, dunes = get_data(i, profile_data)\n",
    "\n",
    "fig =plt.figure(figsize=(10, 9))\n",
    "ax1 = plt.subplot(311)\n",
    "line_xz, = plt.plot(x, z, label='raw')\n",
    "line_xxzz, = plt.plot(xx, zz, label='interpolated')\n",
    "pts_filtered_crests, = plt.plot(\n",
    "    [x['x_crest'] for x in dunes if x['filtered']==True],\n",
    "    [x['z_crest'] for x in dunes if x['filtered']==True],\n",
    "    'rv',\n",
    "    markersize=10,\n",
    "    fillstyle='none',\n",
    "    label='Dune crest (filtered)')\n",
    "pts_filtered_toes, = plt.plot(\n",
    "    [x['x_toe'] for x in dunes if x['filtered']==True],\n",
    "    [x['z_toe'] for x in dunes if x['filtered']==True],\n",
    "    'r^',\n",
    "    markersize=10,\n",
    "    fillstyle='none',\n",
    "    label='Dune toe (filtered)')\n",
    "pts_potential_crests, = plt.plot(\n",
    "    [x['x_crest'] for x in dunes if x['filtered']==False],\n",
    "    [x['z_crest'] for x in dunes if x['filtered']==False],\n",
    "    'rv',\n",
    "    markersize=10,\n",
    "    label='Dune crests (potential)',\n",
    ")\n",
    "pts_potential_toes, = plt.plot(\n",
    "    [x['x_toe'] for x in dunes if x['filtered']==False],\n",
    "    [x['z_toe'] for x in dunes if x['filtered']==False],\n",
    "    'r^',\n",
    "    markersize=10,\n",
    "    label='Dune toes (potential)',\n",
    ")\n",
    "plt.title('{}: elevations'.format(site))\n",
    "plt.legend(loc='best')\n",
    "\n",
    "plt.subplot(312)\n",
    "d1 = np.diff(zz) / np.diff(xx)\n",
    "plt.plot(xx, zz_d1)\n",
    "plt.title('first derivative - slope')\n",
    "\n",
    "plt.subplot(313)\n",
    "d2 = np.diff(d1) / np.diff(xx[1:])\n",
    "plt.plot(xx, zz_d2)\n",
    "plt.title('second derivative - change in slope')\n",
    "plt.tight_layout()\n",
    "text=ax1.text(0,0, \"test\", va=\"bottom\", ha=\"left\")\n",
    "\n",
    "# Get list of potential crests and toes\n",
    "potential_crests_x = [x['x_crest'] for x in dunes if x['filtered']==False]\n",
    "potential_crests_z = [x['z_crest'] for x in dunes if x['filtered']==False]\n",
    "potential_crests = np.array([(x,z) for x,z in zip(potential_crests_x, potential_crests_z)])\n",
    "\n",
    "potential_toes_x = [x['x_toe'] for x in dunes if x['filtered']==False]\n",
    "potential_toes_z = [x['z_toe'] for x in dunes if x['filtered']==False]\n",
    "potential_toes = np.array([(x,z) for x,z in zip(potential_toes_x, potential_toes_z)])\n",
    "\n",
    "def on_click(event):\n",
    "    tx = 'button=%d, x=%d, y=%d, xdata=%f, ydata=%f' % (event.button, event.x, event.y, event.xdata, event.ydata)\n",
    "    text.set_text(tx)\n",
    "\n",
    "    # Left click = dune crest\n",
    "    if event.button == 1:\n",
    "        i_selected_crest = closest_node(np.array((event.xdata, event.ydata)), potential_crests)\n",
    "        selected_crest = potential_crests[i_selected_crest]\n",
    "        tx = 'Selected crest at {}'.format(selected_crest)\n",
    "        text.set_text(tx)\n",
    "    \n",
    "        # Save to data\n",
    "        \n",
    "    \n",
    "    # Right click = dune toe\n",
    "    if event.button == 3:\n",
    "        \n",
    "        if event.xdata > xx[-1]:\n",
    "            selected_toe = np.nan\n",
    "            tx = 'Cleared selected toe'\n",
    "        else:\n",
    "            i_selected_toe = closest_node(np.array((event.xdata, event.ydata)), potential_crests)\n",
    "            selected_toe = potential_toes[i_selected_toe]\n",
    "            tx = 'Selected toe at {}'.format(selected_toe)\n",
    "        text.set_text(tx)\n",
    "    \n",
    "def on_press(event):\n",
    "    global i\n",
    "    global potential_crests\n",
    "    global potential_toes\n",
    "\n",
    "    if event.key == 'left':\n",
    "        i = max(0, i-1)\n",
    "        \n",
    "    if event.key == 'right':\n",
    "        i = min(len(profile_data), i+1)\n",
    "        \n",
    "#     text.set_text(i)\n",
    "    site, x, z, xx, zz,zz_d1,zz_d2,zz_d3, dunes = get_data(i, profile_data)\n",
    "\n",
    "    line_xz.set_xdata(x)\n",
    "    line_xz.set_ydata(z)\n",
    "    line_xxzz.set_xdata(xx)\n",
    "    line_xxzz.set_ydata(zz)\n",
    "    pts_filtered_crests.set_xdata([x['x_crest'] for x in dunes if x['filtered']==True])\n",
    "    pts_filtered_crests.set_ydata([x['z_crest'] for x in dunes if x['filtered']==True])\n",
    "    pts_filtered_toes.set_xdata([x['x_toe'] for x in dunes if x['filtered']==True])\n",
    "    pts_filtered_toes.set_ydata([x['z_toe'] for x in dunes if x['filtered']==True])\n",
    "    pts_potential_crests.set_xdata([x['x_crest'] for x in dunes if x['filtered']==False])\n",
    "    pts_potential_crests.set_ydata([x['z_crest'] for x in dunes if x['filtered']==False])\n",
    "    pts_potential_toes.set_xdata([x['x_toe'] for x in dunes if x['filtered']==False])\n",
    "    pts_potential_toes.set_ydata([x['z_toe'] for x in dunes if x['filtered']==False])\n",
    "    \n",
    "    \n",
    "    # Get list of potential crests and toes\n",
    "    potential_crests_x = [x['x_crest'] for x in dunes if x['filtered']==False]\n",
    "    potential_crests_z = [x['z_crest'] for x in dunes if x['filtered']==False]\n",
    "    potential_crests = np.array([(x,z) for x,z in zip(potential_crests_x, potential_crests_z)])\n",
    "\n",
    "    potential_toes_x = [x['x_toe'] for x in dunes if x['filtered']==False]\n",
    "    potential_toes_z = [x['z_toe'] for x in dunes if x['filtered']==False]\n",
    "    potential_toes = np.array([(x,z) for x,z in zip(potential_toes_x, potential_toes_z)])\n",
    "    \n",
    "    ax1.set_title('{}: elevations'.format(site))\n",
    "    \n",
    "    ax1.relim()\n",
    "    ax1.autoscale_view()\n",
    "    \n",
    "fig.canvas.mpl_connect('button_press_event', on_click)\n",
    "fig.canvas.mpl_connect('key_press_event', on_press)\n",
    "    \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Find berms \n",
    "for x_knot1, x_knot2, in zip(x_knots[0:-2], x_knots[1:-1]):\n",
    "    \n",
    "    # Berms are only located seaward of the dune toe (but sometimes dune toe is undefined, so use dune crest)\n",
    "    if x_knot1 < x_crest:\n",
    "        continue\n",
    "    \n",
    "    # Get slope of segment\n",
    "    x_seg, z_seg = xz(x_knot1, x_knot2, x, z)\n",
    "    slope = get_slope(x_seg, z_seg)\n",
    "    z_mean = np.mean(z_seg)\n",
    "    width = x_knot2 - x_knot1\n",
    "    \n",
    "    if slope > -0.01:\n",
    "        print('{:.1f}m wide berm at x={:.1f}m to x={:.1f}m at z={:.1f}m'.format(width, x_knot1, x_knot2, z_mean))\n",
    "\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true
   },
   "outputs": [],
   "source": [
    "# def reject_outliers(data, m = 5):\n",
    "#     d = np.abs(data - np.median(data))\n",
    "#     mdev = np.median(d)\n",
    "#     s = d/mdev if mdev else 0.\n",
    "#     return ma.masked_where(s > m, data)\n",
    "\n",
    "# def hampel(vals_orig, k=7, t0=3):\n",
    "#     '''\n",
    "#     vals: pandas series of values from which to remove outliers\n",
    "#     k: size of window (including the sample; 7 is equal to 3 on either side of value)\n",
    "#     '''\n",
    "#     #Make copy so original not edited\n",
    "#     vals_orig = pd.Series(vals_orig)\n",
    "#     vals=vals_orig.copy()    \n",
    "#     #Hampel Filter\n",
    "#     L= 1.4826\n",
    "#     rolling_median=vals.rolling(k).median()\n",
    "#     difference=np.abs(rolling_median-vals)\n",
    "#     median_abs_deviation=difference.rolling(k).median()\n",
    "#     threshold= t0 *L * median_abs_deviation\n",
    "#     outlier_idx=difference>threshold\n",
    "#     vals[outlier_idx]=np.nan\n",
    "#     return(vals.tolist())\n",
    "\n",
    "# # Remove segements where slopes are too high, indicative of a structure\n",
    "# segments = [xz(x1, x2,x,z) for (x1, x2) in pairwise(x_knots_simplified)]\n",
    "# slopes = np.array([get_slope(x1,x2) for x1,x2 in segments])\n",
    "# mask = reject_outliers(np.array(slopes))\n",
    "# x_knots_filtered = x_knots_simplified[np.append(True,~mask.mask)]\n",
    "# z_knots_filtered = [z_val for x_val, z_val in zip(x,z) if x_val in x_knots_filtered]\n",
    "\n",
    "\n",
    "\n",
    "# # Hampel filter based on z values\n",
    "# mask = ma.masked_invalid(hampel(z_knots_filtered))\n",
    "# x_knots_filtered = x_knots_filtered[~mask.mask]\n",
    "# z_knots_filtered = [z_val for x_val, z_val in zip(x,z) if x_val in x_knots_filtered]\n",
    "\n",
    "\n",
    "\n",
    "# # plt.figure(figsize=(10,5))\n",
    "# # plt.plot(x_knots_simplified, np.append(slopes[0],slopes), '.-')\n",
    "# # plt.plot(x_knots_simplified[np.append(True,~mask.mask)], np.append(slopes[~mask.mask][0], slopes[~mask.mask]), '.-')\n",
    "# # # plt.plot(xs, zs)\n",
    "# # # plt.plot(x_knots_filtered, z_knots_filtered,'.',markersize=15)\n",
    "# # plt.show()\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "# plt.figure(figsize=(10,5))\n",
    "# plt.plot(x, z, '.-')\n",
    "# plt.plot(xs, zs)\n",
    "# plt.plot(x_knots_filtered, z_knots_filtered,'.',markersize=15)\n",
    "# plt.show()\n",
    "# z_knots_filtered"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Find dune crests"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     18
    ]
   },
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "\n",
    "def crest_by_deriv_2(z):\n",
    "    \"\"\"\n",
    "    Try get crest elevation by taking the difference in 2nd derivative of elevation.\n",
    "    \"\"\"\n",
    "\n",
    "    z = savgol_filter(z, 15, 3, mode='nearest')\n",
    "    deriv_2_diff = np.diff(np.gradient(z)) * -1\n",
    "\n",
    "    # Ensure same size arrays for plotting\n",
    "    deriv_2_diff = np.insert(deriv_2_diff, 0, deriv_2_diff[0])\n",
    "    med = np.median(deriv_2_diff)\n",
    "    std = np.std(deriv_2_diff)\n",
    "    peaks, peak_props = find_peaks(deriv_2_diff, height=med+std*2, distance=10)\n",
    "    return peaks, peak_props, deriv_2_diff\n",
    "\n",
    "\n",
    "def crest_by_elevation_peaks(z, min_dune_prominence=1.0):\n",
    "    # Try find peaks in elevation\n",
    "    peaks, peak_props = find_peaks(z, distance=10, prominence=min_dune_prominence)\n",
    "    return peaks, peak_props\n",
    "\n",
    "\n",
    "def get_dune_crest(x, z, xs,zs, x_knots, z_knots, site):\n",
    "   \n",
    "    peaks, peak_props = crest_by_elevation_peaks(z)\n",
    "    \n",
    "    if peaks.size != 0:\n",
    "        # Choose crest by most landward peak\n",
    "        x_crest = x[peaks[-1]]\n",
    "        z_crest = z[peaks[-1]]\n",
    "    \n",
    "    # If no peaks in elevation are found, find by maximum second derivative\n",
    "    else:\n",
    "        peaks, peak_props, deriv_2_diff = crest_by_deriv_2(z)\n",
    "        \n",
    "        plt.title(\n",
    "            '{}: Difference of 2nd derivative in elevation'.format(site))\n",
    "        plt.xlabel('Distance (m)')\n",
    "        plt.ylabel('Diff(2nd Derivative Elevation)')\n",
    "        plt.plot(x, deriv_2_diff)\n",
    "\n",
    "               \n",
    "        if peaks.size!= 0:\n",
    "\n",
    "            # Choose crest by highest peak in diff of 2nd derivative\n",
    "            i_best = np.argmax(peak_props['peak_heights'])\n",
    "            x_crest = x[peaks[i_best]]\n",
    "            z_crest = z[peaks[i_best]]\n",
    "            plt.plot(\n",
    "                x[peaks],\n",
    "                deriv_2_diff[peaks],\n",
    "                'o',\n",
    "                markersize=20,\n",
    "                fillstyle='none')\n",
    "            plt.show() \n",
    "            \n",
    "        else:\n",
    "            plt.show() \n",
    "            # If no peaks in diff of 2nd derivative are found, take landward most point\n",
    "            if peaks.size == 0:\n",
    "                x_crest = x_knots[0]\n",
    "                z_crest = z_knots[0]\n",
    "\n",
    "    # Move crest to closest knot\n",
    "    idx_crest = min(range(len(x_knots)), key=lambda i: abs(x_knots[i]-x_crest))\n",
    "    x_crest = x_knots[idx_crest]\n",
    "    z_crest = z_knots[idx_crest]\n",
    "    \n",
    "    # Plot for checking\n",
    "    plt.title('{}: Find dune crest by peak in diff 2nd deriv'.format(site))\n",
    "    plt.xlabel('Distance (m)')\n",
    "    plt.ylabel('Elevation (m AHD)')\n",
    "    plt.plot(x, z, 'k.-', label='Raw profile')\n",
    "    plt.plot(\n",
    "        xs,\n",
    "        zs,\n",
    "        'r-',\n",
    "        label='Interpolated profile (s={:.2f})'.format(best_s))\n",
    "    plt.plot(\n",
    "        x_knots,\n",
    "        z_knots,\n",
    "        'ro',\n",
    "        markersize=8,\n",
    "        markeredgewidth=1.5,\n",
    "        markeredgecolor='orange',\n",
    "        label='Interpolated knots (n={})'.format(len(x_knots)))\n",
    "    plt.plot(\n",
    "        x_crest,\n",
    "        z_crest,\n",
    "        'bv',\n",
    "        markersize=10,\n",
    "        label='Dune crest (selected) (x={:.1f}m,z={:.1f}m)'.format(x_crest,z_crest))\n",
    "    plt.plot(\n",
    "        x[peaks],\n",
    "        z[peaks],\n",
    "        'o',\n",
    "        markersize=10,\n",
    "        fillstyle='none',\n",
    "        label='Dune crest (potential)')\n",
    "\n",
    "    plt.legend(loc='best')\n",
    "    plt.show()\n",
    "        \n",
    "    return x_crest, z_crest\n",
    "        \n",
    "x_crest, z_crest = get_dune_crest(x, z, xs,zs, x_knots, z_knots, site)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get toe by second derivative\n",
    "def toe_by_deriv_2(z):\n",
    "    \"\"\"\n",
    "    Try get crest elevation by taking the difference in 2nd derivative of elevation.\n",
    "    \"\"\"\n",
    "\n",
    "    z = savgol_filter(z,11, 3, mode='nearest')\n",
    "    deriv_2_diff = np.diff(np.gradient(z))\n",
    "\n",
    "    # Ensure same size arrays for plotting\n",
    "    deriv_2_diff = np.insert(deriv_2_diff, 0, deriv_2_diff[0])\n",
    "    med = np.median(deriv_2_diff)\n",
    "    std = np.std(deriv_2_diff)\n",
    "    peaks, peak_props = find_peaks(deriv_2_diff, height=max(0.01,med+std*3), distance=10)\n",
    "    return peaks, peak_props, deriv_2_diff\n",
    "\n",
    "def get_dune_toe(x, z, xs,zs, x_crest,z_crest,x_knots, z_knots, site):\n",
    "    peaks, peak_props, deriv_2_diff = toe_by_deriv_2(z)\n",
    "\n",
    "    plt.title('{}: Difference of 2nd derivative in elevation'.format(site))\n",
    "    plt.xlabel('Distance (m)')\n",
    "    plt.ylabel('Diff(2nd Derivative Elevation)')\n",
    "    plt.plot(x, deriv_2_diff)\n",
    "    \n",
    "    if peaks.size!= 0:\n",
    "        \n",
    "        # Remove peaks which are behind the crest\n",
    "        # TODO What happens if all peaks are removed?\n",
    "        mask = ma.masked_where(x[peaks] <= x_crest, x[peaks])\n",
    "        peaks = peaks[~mask.mask]\n",
    "\n",
    "    if peaks.size!= 0:\n",
    "        # Choose toe by highest peak in diff of 2nd derivative\n",
    "        i_best = np.argmax(peak_props['peak_heights'])\n",
    "        x_toe = x[peaks[i_best]][0]\n",
    "        z_toe = z[peaks[i_best]][0]\n",
    "#         # Move crest to closest knot\n",
    "#         idx_toe = min(range(len(x_knots)), key=lambda i: abs(x_knots[i]-x_toe))\n",
    "#         x_toe = x_knots[idx_toe]\n",
    "#         z_toe = z_knots[idx_toe]\n",
    "\n",
    "        plt.plot(\n",
    "            x[peaks],\n",
    "            deriv_2_diff[peaks],\n",
    "            'o',\n",
    "            markersize=20,\n",
    "            fillstyle='none')\n",
    "        plt.show()\n",
    "    else:\n",
    "        x_toe = np.nan\n",
    "        z_toe = np.nan\n",
    "        plt.show()\n",
    "\n",
    "    # Plot for checking\n",
    "    plt.title('{}: Find dune crest by peak in diff 2nd deriv'.format(site))\n",
    "    plt.xlabel('Distance (m)')\n",
    "    plt.ylabel('Elevation (m AHD)')\n",
    "    plt.plot(x, z, 'k.-', label='Raw profile')\n",
    "    plt.plot(\n",
    "        xs,\n",
    "        zs,\n",
    "        'r-',\n",
    "        label='Interpolated profile (s={:.2f})'.format(best_s))\n",
    "    plt.plot(\n",
    "        x_knots,\n",
    "        z_knots,\n",
    "        'ro',\n",
    "        markersize=8,\n",
    "        markeredgewidth=1.5,\n",
    "        markeredgecolor='orange',\n",
    "        label='Interpolated knots (n={})'.format(len(x_knots)))\n",
    "    plt.plot(\n",
    "        x_crest,\n",
    "        z_crest,\n",
    "        'bv',\n",
    "        markersize=10,\n",
    "        label='Dune crest(x={:.1f} m,z={:.1f} m)'.format(x_crest,z_crest))\n",
    "    plt.plot(\n",
    "        x_toe,\n",
    "        z_toe,\n",
    "        'b^',\n",
    "        markersize=10,\n",
    "        label='Dune toe (x={:.1f} m,z={:.1f} m)'.format(x_toe,z_toe))\n",
    "\n",
    "    plt.legend(loc='best')\n",
    "    plt.show()\n",
    "    return x_toe, z_toe\n",
    "x_toe, z_toe = get_dune_toe(x, z, xs,zs, x_crest,z_crest,x_knots, z_knots, site)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def find_dune_crest_and_toe(site, profile_type, df_profiles):\n",
    "    x, z = get_profile_data(site, profile_type, df_profiles)\n",
    "    x_knots, z_knots = simplify_profile(x, z, site)\n",
    "    x_knots, z_knots = simplify_segments_with_same_slopes(x_knots, x, z, site)\n",
    "    xs, zs = get_x_z_from_knots(x_knots, z_knots)\n",
    "    x_crest, z_crest = get_dune_crest(x, z, xs,zs, x_knots, z_knots, site)\n",
    "    x_toe, z_toe = get_dune_toe(x, z, xs,zs, x_crest,z_crest,x_knots, z_knots, site)\n",
    "    \n",
    "find_dune_crest_and_toe('ENTRA0004', 'prestorm', df_profiles)\n"
   ]
  }
 ],
 "metadata": {
  "hide_input": false,
  "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.6.6"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}