chrisl
/
nsw-2016-storm-impact
2
1
Fork 0
Code Issues Pull Requests Releases Activity
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.
develop
master
Branches Tags
${ item.name }
Create tag ${ searchTerm }
Create branch ${ searchTerm }
from 'bb17d6c786'
${ noResults }
nsw-2016-storm-impact/notebooks/04b_profile_picker.ipynb

1375 lines
52 KiB
Plaintext
Raw Normal View History Unescape Escape

Update notebooks
6 years ago
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Beach profile classifier\n",
"Using RANSAC algorithm"
]
},
{
"cell_type": "markdown",
"metadata": {
"heading_collapsed": true
},
"source": [
"## Setup notebook"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"hidden": true
},
"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": {
"hidden": true
},
"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"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"hidden": true
},
"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": {
"heading_collapsed": true
},
"source": [
"## Import data"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"hidden": true
},
"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 profile.\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",
" # Get land limit\n",
" land_lim = df_profiles.loc[(site, 'poststorm')].land_lim.loc[0]\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",
"# 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",
"# Load profile data\n",
"# site = 'WAMBE0010'\n",
"# site = 'NINEMs0048'\n",
"# site = 'NAMB0013'\n",
"# site = 'AVOCAn0009'\n",
"# site = 'GRANTSs0014'\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": "code",
"execution_count": null,
"metadata": {
"hide_input": true
},
"outputs": [],
"source": [
"\n",
"\n",
"\n",
"# # n_samples = 1000\n",
"# # n_outliers = 50\n",
"\n",
"\n",
"# # X, y, coef = datasets.make_regression(n_samples=n_samples, n_features=1,\n",
"# # n_informative=1, noise=10,\n",
"# # coef=True, random_state=0)\n",
"\n",
"# # # Add outlier data\n",
"# # np.random.seed(0)\n",
"# # X[:n_outliers] = 3 + 0.5 * np.random.normal(size=(n_outliers, 1))\n",
"# # y[:n_outliers] = -3 + 10 * np.random.normal(size=n_outliers)\n",
"\n",
"# # # Fit line using all data\n",
"# # lr = linear_model.LinearRegression()\n",
"# # lr.fit(X, y)\n",
"\n",
"# # Robustly fit linear model with RANSAC algorithm\n",
"# ransac = linear_model.RANSACRegressor()\n",
"# ransac.fit(x, z)\n",
"# inlier_mask = ransac.inlier_mask_\n",
"# outlier_mask = np.logical_not(inlier_mask)\n",
"\n",
"# # Predict data of estimated models\n",
"# line_x = np.arange(x.min(), x.max())[:, np.newaxis]\n",
"# line_y_ransac = ransac.predict(line_x)\n",
"\n",
"# lw = 2\n",
"# plt.scatter(x[inlier_mask], z[inlier_mask], color='yellowgreen', marker='.',\n",
"# label='Inliers')\n",
"# plt.scatter(x[outlier_mask], z[outlier_mask], color='gold', marker='.',\n",
"# label='Outliers')\n",
"# plt.plot(line_x, line_y_ransac, color='cornflowerblue', linewidth=lw,\n",
"# label='RANSAC regressor')\n",
"# plt.legend(loc='lower right')\n",
"# plt.xlabel(\"Input\")\n",
"# plt.ylabel(\"Response\")\n",
"# plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Fit linear interpolated spline to profile\n",
"- Linear interpolated spline used to simplify the profile.\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": [
4,
11,
20,
25,
30,
42,
92
]
},
"outputs": [],
"source": [
"from scipy.interpolate import UnivariateSpline\n",
"from scipy.interpolate import interp1d\n",
"\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": [
0
]
},
"outputs": [],
"source": [
"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": [
1
]
},
"outputs": [],
"source": [
"## BACKUP\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 = interpolate.splrep(x, z, k=4,t=x_knots[1:-1])\n",
" spl = interpolate.splrep(x, z, k=3,s=best_s)\n",
"\n",
" xx = np.linspace(x[0], x[-1], 1000)\n",
" zz = interpolate.splev(xx, spl)\n",
" zz_d1 = interpolate.splev(xx, spl, der=1)\n",
" zz_d2 = interpolate.splev(xx, spl, der=2)\n",
" zz_d3 = interpolate.splev(xx, spl, der=3)\n",
"\n",
" try:\n",
" # Find dune crest by looking at elevation\n",
" peaks_crest, props_crest = find_peaks(zz, height=3, prominence=0.5)\n",
" \n",
" # Save potential crests\n",
" x_potential_crests = xx[peaks_crest]\n",
" z_potential_crests = zz[peaks_crest]\n",
"\n",
" # Take most seaward peak as dune crest\n",
" i_crest = peaks_crest[-1]\n",
" x_crest = xx[i_crest]\n",
" z_crest = zz[i_crest]\n",
"\n",
" print('Found {} potential dune crests from elevation'.format(len(peaks_crest)))\n",
"\n",
" except:\n",
" print('No crests found from elevation')\n",
" x_crest = np.nan\n",
" z_crest = np.nan\n",
" pass\n",
"\n",
" try:\n",
" if np.isnan(x_crest) and np.isnan(x_crest):\n",
" # Find dune crest by looking at 2nd deriviate\n",
" # Multiply by -1 since crest will have a negative 2nd derivative and we want to find peaks (i.e. positive)\n",
" peaks_crest, props_crest = find_peaks(zz_d2*-1, height=0.02,width=0)\n",
" \n",
" \n",
" # Save potential crests\n",
" x_potential_crests = xx[peaks_crest]\n",
" z_potential_crests = zz[peaks_crest]\n",
" print('Found {} potential dune crests from 2nd derivitive'.format(len(peaks_crest)))\n",
"\n",
" # Take peak with biggest height\n",
" i_biggest_crest = np.argmax(props_crest['peak_heights'])\n",
" i_crest = peaks_crest[i_biggest_crest]\n",
"# x_crest = xx[i_crest]\n",
"# z_crest = zz[i_crest]\n",
"\n",
" # Then take the left base of that peak\n",
" x_crest = xx[props_crest['left_bases'][i_biggest_crest]]\n",
" z_crest= zz[props_crest['left_bases'][i_biggest_crest]]\n",
"\n",
"\n",
" except:\n",
" # Take crest as location with maximum elevation\n",
" i_crest = np.argmax(zz)\n",
" x_crest = xx[i_crest]\n",
" z_crest = zz[i_crest]\n",
" print('No crest found, taking maximum elevation as crest')\n",
"\n",
" try:\n",
" # Find due toe\n",
" peaks_toe, props_toe = find_peaks(zz_d2, height=0, prominence=0.02)\n",
"# peaks_toe, props_toe = find_peaks(zz_d2, height=-3, prominence=0.00)\n",
"\n",
" # Save potential crests\n",
" x_potential_toes = xx[peaks_toe]\n",
" z_potential_toes = zz[peaks_toe]\n",
"\n",
" # Remove toes which are behind dune crest\n",
" heights_toe = props_toe['peak_heights'][peaks_toe > i_crest]\n",
" peaks_toe = peaks_toe[peaks_toe > i_crest]\n",
"\n",
" # Remove toes that are less than 0.5m in height from dune crest\n",
"\n",
" mask = z_crest-0.5>zz[peaks_toe]\n",
" heights_toe = heights_toe[mask]\n",
" peaks_toe = peaks_toe[mask]\n",
"\n",
" # Take toe with biggest 2nd derivative height\n",
" i_toe = peaks_toe[np.argmax(heights_toe)]\n",
" x_toe = xx[i_toe]\n",
" z_toe = zz[i_toe]\n",
" print('Found {} potential dune toes from 2nd derivitive'.format(len(peaks_toe)))\n",
"\n",
"# # The above method usually results in a toe value which is slightly too high up.\n",
"# # We can improve this by moving the toe seaward to the next negative peak in the third\n",
"# # derivative of elevation\n",
"# peaks_toe_d3, props_toe_d3 = find_peaks(zz_d3*-1, height=0.03, prominence=0.02)\n",
"\n",
"# # Filter everything landward of the previous toe position, then take the first (most landward) peak\n",
"# mask = xx[peaks_toe_d3] > x_toe\n",
"# peaks_toe_d3 = peaks_toe_d3[mask]\n",
"\n",
"# if peaks_toe_d3.size != 0:\n",
"# print('Adjusting toe position based on d3')\n",
"# x_toe = xx[peaks_toe_d3[0]]\n",
"# z_toe = zz[peaks_toe_d3[0]]\n",
"\n",
"# # Move toe forward to based on the third derivative of elevation\n",
"# print('Adjusting toe position based on d3')\n",
"# mask = (zz_d3 > -0.0008) & (xx > x_toe) & (np.r_[np.diff(zz_d3)[0],np.diff(zz_d3)] > 0)\n",
"# x_toe = xx[mask][1]\n",
"# z_toe = zz[mask][1]\n",
"\n",
" except:\n",
" x_toe = np.nan\n",
" z_toe = np.nan\n",
" print('No toe found')\n",
" \n",
" return xx,zz, zz_d1,zz_d2,zz_d3,x_crest,z_crest, x_toe,z_toe, x_potential_crests,z_potential_crests, x_potential_toes,z_potential_toes"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"code_folding": [
126
]
},
"outputs": [],
"source": [
"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 = interpolate.splrep(x, z, k=4,t=x_knots[1:-1])\n",
" spl = interpolate.splrep(x, z, k=3,s=best_s)\n",
"\n",
" xx = np.linspace(x[0], x[-1], 1000)\n",
" zz = interpolate.splev(xx, spl)\n",
" zz_d1 = interpolate.splev(xx, spl, der=1)\n",
" zz_d2 = interpolate.splev(xx, spl, der=2)\n",
" zz_d3 = interpolate.splev(xx, spl, der=3)\n",
"\n",
" # Find the biggest slopes\n",
" peaks, props = find_peaks(zz_d1*-1, height=0.15, width=0,prominence=0.1)\n",
"\n",
" if len(peaks) != 0:\n",
"\n",
" x_potential_crests =[xx[val] for val in props['left_bases']]\n",
" x_potential_toes =[xx[val] for val in props['right_bases']]\n",
" x_potential_heights = [val for val in props['peak_heights']]\n",
" x_potential_face_center = xx[peaks]\n",
"\n",
" # Put a limit on how when a slope ends.\n",
" # Dune faces has to be steeper than this value\n",
"# min_slope = -0.10\n",
"\n",
" x_new_potential_crests = []\n",
" x_new_potential_toes = []\n",
" z_new_potential_crests = []\n",
" z_new_potential_toes = []\n",
" dune_heights = []\n",
" slopes = []\n",
"\n",
" for x_crest, x_center, x_toe,height in zip(x_potential_crests, x_potential_face_center, x_potential_toes,x_potential_heights):\n",
" print('Checking x_crest={:.2f}m,x_center={:.2f}m,x_toe={:.2f}m'.format(x_crest, x_center, x_toe))\n",
"\n",
" min_slope = height/-6\n",
" # Check the land side of the dune face\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",
" if len(x_land_slope_okay) == 0:\n",
" x_new_crest = x_crest\n",
" else:\n",
" x_new_crest = x_land_slope_okay[-1]\n",
" idx_new_crest = np.where(x_new_crest==xx)[0]\n",
" z_new_crest = zz[idx_new_crest][0]\n",
"\n",
" # # TODO Check, for now, keep crest as detected level\n",
" # # Maybe we don't need to change\n",
" # x_new_crest = x_crest\n",
" # idx_new_crest = np.where(x_new_crest==xx)[0]\n",
" # z_new_crest = zz[idx_new_crest][0]\n",
" # ###\n",
"\n",
" # Check the sea side of the dune face\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",
" if len(x_sea_slope_okay) == 0:\n",
" x_new_toe = x_toe\n",
" else:\n",
" x_new_toe = x_sea_slope_okay[0]\n",
" idx_new_toe = np.where(x_new_toe==xx)[0]\n",
" z_new_toe = zz[idx_new_toe][0]\n",
"\n",
" slopes.append((z_new_crest - z_new_toe)/(x_new_crest - x_new_toe))\n",
" dune_heights.append(z_new_crest - z_new_toe)\n",
" x_new_potential_crests.append(x_new_crest)\n",
" x_new_potential_toes.append(x_new_toe)\n",
" z_new_potential_crests.append(z_new_crest)\n",
" z_new_potential_toes.append(z_new_toe)\n",
"\n",
" # Specifiy minimum dune crest elevation and minimum dune height\n",
" min_dune_ele = np.median(zz)/3 # m AHD\n",
" print('min_dune_ele: {:.2f}m'.format(min_dune_ele))\n",
" min_dune_height = 0.5 # m\n",
" mask = (np.array(z_new_potential_crests) > min_dune_ele) & (np.array(dune_heights) > min_dune_height )\n",
"\n",
" slopes = np.array(slopes)[mask]\n",
" dune_heights = np.array(dune_heights)[mask] \n",
" x_potential_crests = np.array(x_new_potential_crests)[mask]\n",
" x_potential_toes = np.array(x_new_potential_toes)[mask]\n",
" z_potential_crests = np.array(z_new_potential_crests)[mask]\n",
" z_potential_toes = np.array(z_new_potential_toes)[mask]\n",
"\n",
"# # Select dune by largest dune height\n",
"# i_dune = np.argmax(dune_heights)\n",
"\n",
" ## TODO Throws error if x_potential_crests is empty\n",
"\n",
" # Select dune by most seaward\n",
" i_dune = np.argmax(x_potential_crests)\n",
"\n",
"\n",
" \n",
"# # Select dune by biggest slope\n",
"# i_dune =np.argmax(slopes*-1)\n",
" x_crest = x_potential_crests[i_dune]\n",
" x_toe = x_potential_toes[i_dune]\n",
" z_crest = z_potential_crests[i_dune]\n",
" z_toe = z_potential_toes[i_dune]\n",
"\n",
" # try:\n",
" # # Find dune crest by looking at elevation\n",
" # peaks_crest, props_crest = find_peaks(zz, height=3, prominence=0.5)\n",
"\n",
" # # Save potential crests\n",
" # x_potential_crests = xx[peaks_crest]\n",
" # z_potential_crests = zz[peaks_crest]\n",
"\n",
" # # Take most seaward peak as dune crest\n",
" # i_crest = peaks_crest[-1]\n",
" # x_crest = xx[i_crest]\n",
" # z_crest = zz[i_crest]\n",
"\n",
" # print('Found {} potential dune crests from elevation'.format(len(peaks_crest)))\n",
"\n",
" # except:\n",
" # print('No crests found from elevation')\n",
" # x_crest = np.nan\n",
" # z_crest = np.nan\n",
" # pass\n",
"\n",
" # try:\n",
" # if np.isnan(x_crest) and np.isnan(x_crest):\n",
" # # Find dune crest by looking at 2nd deriviate\n",
" # # Multiply by -1 since crest will have a negative 2nd derivative and we want to find peaks (i.e. positive)\n",
" # peaks_crest, props_crest = find_peaks(zz_d2*-1, height=0.02,width=0)\n",
"\n",
"\n",
" # # Save potential crests\n",
" # x_potential_crests = xx[peaks_crest]\n",
" # z_potential_crests = zz[peaks_crest]\n",
" # print('Found {} potential dune crests from 2nd derivitive'.format(len(peaks_crest)))\n",
"\n",
" # # Take peak with biggest height\n",
" # i_biggest_crest = np.argmax(props_crest['peak_heights'])\n",
" # i_crest = peaks_crest[i_biggest_crest]\n",
" # # x_crest = xx[i_crest]\n",
" # # z_crest = zz[i_crest]\n",
"\n",
" # # Then take the left base of that peak\n",
" # x_crest = xx[props_crest['left_bases'][i_biggest_crest]]\n",
" # z_crest= zz[props_crest['left_bases'][i_biggest_crest]]\n",
"\n",
"\n",
" # except:\n",
" # # Take crest as location with maximum elevation\n",
" # i_crest = np.argmax(zz)\n",
" # x_crest = xx[i_crest]\n",
" # z_crest = zz[i_crest]\n",
" # print('No crest found, taking maximum elevation as crest')\n",
"\n",
" # try:\n",
" # # Find due toe\n",
" # peaks_toe, props_toe = find_peaks(zz_d2, height=0, prominence=0.02)\n",
" # # peaks_toe, props_toe = find_peaks(zz_d2, height=-3, prominence=0.00)\n",
"\n",
" # # Save potential crests\n",
" # x_potential_toes = xx[peaks_toe]\n",
" # z_potential_toes = zz[peaks_toe]\n",
"\n",
" # # Remove toes which are behind dune crest\n",
" # heights_toe = props_toe['peak_heights'][peaks_toe > i_crest]\n",
" # peaks_toe = peaks_toe[peaks_toe > i_crest]\n",
"\n",
" # # Remove toes that are less than 0.5m in height from dune crest\n",
"\n",
" # mask = z_crest-0.5>zz[peaks_toe]\n",
" # heights_toe = heights_toe[mask]\n",
" # peaks_toe = peaks_toe[mask]\n",
"\n",
" # # Take toe with biggest 2nd derivative height\n",
" # i_toe = peaks_toe[np.argmax(heights_toe)]\n",
" # x_toe = xx[i_toe]\n",
" # z_toe = zz[i_toe]\n",
" # print('Found {} potential dune toes from 2nd derivitive'.format(len(peaks_toe)))\n",
"\n",
" # # # The above method usually results in a toe value which is slightly too high up.\n",
" # # # We can improve this by moving the toe seaward to the next negative peak in the third\n",
" # # # derivative of elevation\n",
" # # peaks_toe_d3, props_toe_d3 = find_peaks(zz_d3*-1, height=0.03, prominence=0.02)\n",
"\n",
" # # # Filter everything landward of the previous toe position, then take the first (most landward) peak\n",
" # # mask = xx[peaks_toe_d3] > x_toe\n",
" # # peaks_toe_d3 = peaks_toe_d3[mask]\n",
"\n",
" # # if peaks_toe_d3.size != 0:\n",
" # # print('Adjusting toe position based on d3')\n",
" # # x_toe = xx[peaks_toe_d3[0]]\n",
" # # z_toe = zz[peaks_toe_d3[0]]\n",
"\n",
" # # # Move toe forward to based on the third derivative of elevation\n",
" # # print('Adjusting toe position based on d3')\n",
" # # mask = (zz_d3 > -0.0008) & (xx > x_toe) & (np.r_[np.diff(zz_d3)[0],np.diff(zz_d3)] > 0)\n",
" # # x_toe = xx[mask][1]\n",
" # # z_toe = zz[mask][1]\n",
"\n",
" # except:\n",
" # x_toe = np.nan\n",
" # z_toe = np.nan\n",
" # print('No toe found')\n",
"\n",
" \n",
" else:\n",
" # Take crest as location with maximum elevation\n",
" i_crest = np.argmax(zz)\n",
" x_crest = xx[i_crest]\n",
" z_crest = zz[i_crest]\n",
" x_potential_crests = np.nan\n",
" z_potential_crests = np.nan\n",
" print('No crest found, taking maximum elevation as crest')\n",
" \n",
" x_toe = np.nan\n",
" z_toe = np.nan\n",
" x_potential_toes=np.nan\n",
" z_potential_toes=np.nan\n",
" \n",
" # Check if toe is necessary. If slope of the dune face is similar to the median slope seaward of the dune\n",
" # toe, let's remove the dune toe\n",
" try:\n",
" if np.isnan(x_toe) == False:\n",
" dune_slope =(z_crest - z_toe) / (x_crest - x_toe)\n",
" median_foreshore_slope = np.median(zz_d1[xx > x_toe])\n",
" std_foreshore_slope = np.std(zz_d1[xx>x_toe])\n",
" print('foreshore_std:{:.4f}'.format(std_foreshore_slope))\n",
" print('foreshore_med_slope:{:.4f}'.format(median_foreshore_slope))\n",
" print('dune_slope:{:.4f}'.format(dune_slope))\n",
" if abs(dune_slope - median_foreshore_slope) < 0.075 and std_foreshore_slope <0.025:\n",
" x_toe = np.nan\n",
" z_toe = np.nan\n",
" except:\n",
" pass\n",
" \n",
" return xx,zz, zz_d1,zz_d2,zz_d3,x_crest,z_crest, x_toe,z_toe, x_potential_crests,z_potential_crests, x_potential_toes,z_potential_toes \n",
" \n",
" \n",
"def plot_profile(site, x,z, xx,zz, zz_d1,zz_d2,zz_d3,x_crest,z_crest, x_toe,z_toe, x_potential_crests,z_potential_crests, x_potential_toes,z_potential_toes):\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(x_potential_crests,z_potential_crests,'rv',markersize=10,fillstyle='none',label='Dune crest (potential)')\n",
" plt.plot(x_potential_toes,z_potential_toes,'r^',markersize=10,fillstyle='none',label='Dune toe (potential)')\n",
" plt.plot(x_crest,z_crest,'rv',markersize=10,label='Dune crest (x={:.2f} m, z={:.2f} m)'.format(x_crest,z_crest))\n",
" plt.plot(x_toe,z_toe,'r^',markersize=10,label='Dune toe (x={:.2f} m, z={:.2f} m)'.format(x_toe,z_toe))\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",
"# xx,zz, x_crest,z_crest, x_toe,z_toe, x_potential_crests,z_potential_crests, x_potential_toes,z_potential_toes = find_crest_and_toes(x,z,x_knots)\n",
"\n",
"# plot_profile(site, x,z, xx,zz, x_crest,z_crest, x_toe,z_toe, x_potential_crests,z_potential_crests, x_potential_toes,z_potential_toes)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"code_folding": []
},
"outputs": [],
"source": [
"def get_crests_and_toes(site, profile_type, df_profiles):\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",
"# xs, zs = get_x_z_from_knots(x_knots, z_knots)\n",
" xx,zz, zz_d1,zz_d2,zz_d3,x_crest,z_crest, x_toe,z_toe, x_potential_crests,z_potential_crests, x_potential_toes,z_potential_toes = 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,x_crest,z_crest, x_toe,z_toe, x_potential_crests,z_potential_crests, x_potential_toes,z_potential_toes)\n",
"\n",
" \n",
"# site = 'WAMBE0010'\n",
"# site = 'NINEMs0048'\n",
"# site = 'NAMB0013'\n",
"# site = 'AVOCAn0009'\n",
"# site = 'GRANTSs0014'\n",
"site = 'NARRA0001'\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='NSHORE_s0014'\n",
"get_crests_and_toes(site,'prestorm', df_profiles)"
]
},
{
"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": [
"from scipy.signal import find_peaks\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
}