You cannot select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
1402 lines
51 KiB
Plaintext
1402 lines
51 KiB
Plaintext
{
|
|
"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
|
|
}
|