From 5850871c14dc0b6e4bf439a7e24bc089cb2507b7 Mon Sep 17 00:00:00 2001 From: Chris Leaman Date: Fri, 8 Feb 2019 14:49:18 +1100 Subject: [PATCH] Update notebooks --- .gitattributes | 2 + notebooks/01_exploration.ipynb | 45 +- notebooks/04b_profile_picker.ipynb | 1025 +++++++++-------- notebooks/04c_profile_picker_manual.ipynb | 242 ++++ notebooks/05_twl_exceedence.ipynb | 866 +++++++++++++- notebooks/06_change_in_slope.ipynb | 270 ++++- notebooks/07-longshore-plotting.ipynb | 348 ++++++ notebooks/08-narr-topo-bathy-slope-test.ipynb | 767 ++++++++++++ 8 files changed, 3027 insertions(+), 538 deletions(-) create mode 100644 notebooks/04c_profile_picker_manual.ipynb create mode 100644 notebooks/07-longshore-plotting.ipynb create mode 100644 notebooks/08-narr-topo-bathy-slope-test.ipynb diff --git a/.gitattributes b/.gitattributes index d9d6885..ff7e90b 100644 --- a/.gitattributes +++ b/.gitattributes @@ -1,3 +1,5 @@ + + *.ipynb filter=nbstripout *.ipynb diff=ipynb diff --git a/notebooks/01_exploration.ipynb b/notebooks/01_exploration.ipynb index 0430e6f..c034e09 100644 --- a/notebooks/01_exploration.ipynb +++ b/notebooks/01_exploration.ipynb @@ -611,7 +611,46 @@ "metadata": {}, "outputs": [], "source": [ - "df_sites.site_no.to_csv('temp.csv')" + "# Get twl information at maximum R_high for a site\n", + "# site_id = 'NARRA0008'\n", + "site_id = 'NARRA0012'\n", + "\n", + "print('TWLs:')\n", + "t = twls['forecasted']['mean_slope_sto06'].xs(site_id,level='site_id').R_high.idxmax()\n", + "print(twls['forecasted']['mean_slope_sto06'].loc[(site_id, t)])\n", + "\n", + "print('\\nforecast regime:')\n", + "print(impacts['forecasted']['mean_slope_sto06'].loc[site_id])\n", + "\n", + "print('\\nobserved regime:')\n", + "print(impacts['observed'].loc[site_id])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "site_id = 'NARRA0010'\n", + "\n", + "print(site_id)\n", + "# Storm duration\n", + "Hs0_storm = 3.0\n", + "df_twls = twls['forecasted']['mean_slope_sto06'].xs(site_id,level='site_id')\n", + "t_start = np.argmax(df_twls.Hs0 > Hs0_storm)\n", + "t_end = np.argmax((df_twls.Hs0 > Hs0_storm)[::-1])\n", + "i_start = df_twls.index.get_loc(t_start)\n", + "i_end = df_twls.index.get_loc(t_end)\n", + "\n", + "df_storm = df_twls.iloc[i_start:i_end]\n", + "print('Storm length: {} hrs'.format(len(df_storm)))\n", + "\n", + "# Get hours above a certain elevation\n", + "z_critical = 2.4\n", + "\n", + "n_impact_hrs = np.sum(df_storm.R_high >z_critical)\n", + "print('Number of hours before peak water level with R_high > {}m: {}hrs'.format(z_critical,n_impact_hrs))\n" ] }, { @@ -879,10 +918,10 @@ "height": "656px", "left": "508px", "top": "90px", - "width": "218.797px" + "width": "282.797px" }, "toc_section_display": true, - "toc_window_display": true + "toc_window_display": false }, "varInspector": { "cols": { diff --git a/notebooks/04b_profile_picker.ipynb b/notebooks/04b_profile_picker.ipynb index 160904f..a78e620 100644 --- a/notebooks/04b_profile_picker.ipynb +++ b/notebooks/04b_profile_picker.ipynb @@ -5,14 +5,12 @@ "metadata": {}, "source": [ "# Beach profile classifier\n", - "Using RANSAC algorithm" + "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": { - "heading_collapsed": true - }, + "metadata": {}, "source": [ "## Setup notebook" ] @@ -20,9 +18,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "hidden": true - }, + "metadata": {}, "outputs": [], "source": [ "# Enable autoreloading of our modules. \n", @@ -36,9 +32,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "hidden": true - }, + "metadata": {}, "outputs": [], "source": [ "from IPython.core.debugger import set_trace\n", @@ -69,15 +63,21 @@ "import numpy as np\n", "from matplotlib import pyplot as plt\n", "\n", - "from sklearn import linear_model, datasets" + "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": { - "hidden": true - }, + "metadata": {}, "outputs": [], "source": [ "# Matplot lib default settings\n", @@ -90,19 +90,16 @@ }, { "cell_type": "markdown", - "metadata": { - "heading_collapsed": true - }, + "metadata": {}, "source": [ - "## Import data" + "## Import data\n", + "Let's first import data from our pre-processed interim data folder." ] }, { "cell_type": "code", "execution_count": null, - "metadata": { - "hidden": true - }, + "metadata": {}, "outputs": [], "source": [ "def df_from_csv(csv, index_col, data_folder='../data/interim'):\n", @@ -156,35 +153,27 @@ "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", + " 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", - " # 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", + " # # 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", - "# 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", @@ -193,63 +182,12 @@ "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", + "- 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." ] }, @@ -258,20 +196,22 @@ "execution_count": null, "metadata": { "code_folding": [ - 4, - 11, - 20, - 25, - 30, - 42, - 92 + 0, + 7, + 16, + 21, + 26, + 38, + 88, + 110 ] }, "outputs": [], "source": [ - "from scipy.interpolate import UnivariateSpline\n", - "from scipy.interpolate import interp1d\n", + "%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", @@ -396,22 +336,22 @@ " 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", + "# # 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", @@ -424,12 +364,12 @@ "cell_type": "code", "execution_count": null, "metadata": { - "code_folding": [ - 0 - ] + "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", @@ -468,14 +408,14 @@ " 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", + "# 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", + "# 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", @@ -528,417 +468,316 @@ "execution_count": null, "metadata": { "code_folding": [ - 1 + 0, + 91 ] }, "outputs": [], "source": [ - "## BACKUP\n", - "def find_crest_and_toes(x,z,x_knots,best_s):\n", + "%matplotlib inline \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", + "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", - " 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", + " # 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", - " # 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", + " # 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", - " # 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", + " # 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", - " 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_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", + "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", - "# # Select dune by largest dune height\n", - "# i_dune = np.argmax(dune_heights)\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", - " ## TODO Throws error if x_potential_crests is empty\n", + " # Check if this slope reduces from a certain threshold\n", + " min_threshold_slope = min_seaward_slope / 6\n", "\n", - " # Select dune by most seaward\n", - " i_dune = np.argmax(x_potential_crests)\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", - " \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", + " # 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", - " 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", + " # 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", - " # 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", + " 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", - " 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", + "# 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(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.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.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.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", + " # 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.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", + " # 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)" + "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)" ] }, { @@ -949,29 +788,217 @@ }, "outputs": [], "source": [ - "def get_crests_and_toes(site, profile_type, df_profiles):\n", + "%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", - "# 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", + " 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", - "# 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)" + "# 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" ] }, { @@ -1082,7 +1109,7 @@ }, "outputs": [], "source": [ - "from scipy.signal import find_peaks\n", + "\n", "\n", "\n", "def crest_by_deriv_2(z):\n", diff --git a/notebooks/04c_profile_picker_manual.ipynb b/notebooks/04c_profile_picker_manual.ipynb new file mode 100644 index 0000000..22b922d --- /dev/null +++ b/notebooks/04c_profile_picker_manual.ipynb @@ -0,0 +1,242 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Profile picker" + ] + }, + { + "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\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": { + "hidden": true + }, + "outputs": [], + "source": [ + "# Matplot lib default settings\n", + "plt.rcParams[\"figure.figsize\"] = (10,6)\n", + "plt.rcParams['axes.grid']=True\n", + "plt.rcParams['grid.alpha'] = 0.5\n", + "plt.rcParams['grid.color'] = \"grey\"\n", + "plt.rcParams['grid.linestyle'] = \"--\"\n", + "plt.rcParams['axes.grid']=True\n", + "\n", + "# https://stackoverflow.com/a/20709149\n", + "matplotlib.rcParams['text.usetex'] = True\n", + "\n", + "matplotlib.rcParams['text.latex.preamble'] = [\n", + " r'\\usepackage{siunitx}', # i need upright \\micro symbols, but you need...\n", + " r'\\sisetup{detect-all}', # ...this to force siunitx to actually use your fonts\n", + " r'\\usepackage{helvet}', # set the normal font here\n", + " r'\\usepackage{amsmath}',\n", + " r'\\usepackage{sansmath}', # load up the sansmath so that math -> helvet\n", + " r'\\sansmath', # <- tricky! -- gotta actually tell tex to use!\n", + "] " + ] + }, + { + "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_profiles = df_from_csv('profiles.csv', index_col=[0, 1, 2])\n", + "df_profile_features_crest_toes = df_from_csv('profile_features_crest_toes.csv', index_col=[0,1])\n", + "\n", + "print('Done!')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Manually pick features" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib notebook\n", + "\n", + "sites = df_profiles.index.get_level_values('site_id').unique()\n", + "\n", + "\n", + "fig =plt.figure(figsize=(10, 3))\n", + "\n", + "df_prestorm = df_profiles.xs((sites[0],'prestorm'),level=('site_id','profile_type'))\n", + "df_poststorm = df_profiles.xs((sites[0],'poststorm'),level=('site_id','profile_type'))\n", + "line_prestorm, = plt.plot(df_prestorm.index, df_prestorm.z, label='prestorm')\n", + "line_poststorm, = plt.plot(df_prestorm.index, df_prestorm.z, label='poststorm')\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# df_profiles.xs((sites[0],'prestorm'),level=('site_id','profile_type'))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "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 +} diff --git a/notebooks/05_twl_exceedence.ipynb b/notebooks/05_twl_exceedence.ipynb index 848a527..2055570 100644 --- a/notebooks/05_twl_exceedence.ipynb +++ b/notebooks/05_twl_exceedence.ipynb @@ -52,11 +52,41 @@ "import matplotlib\n", "from matplotlib import cm\n", "import colorlover as cl\n", - "\n", + "from tqdm import tqdm_notebook\n", "from ipywidgets import widgets, Output\n", "from IPython.display import display, clear_output, Image, HTML\n", + "from scipy import stats\n", + "from sklearn.metrics import confusion_matrix\n", + "import matplotlib.pyplot as plt\n", + "from scipy.interpolate import interp1d" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Matplot lib default settings\n", + "plt.rcParams[\"figure.figsize\"] = (10,6)\n", + "plt.rcParams['axes.grid']=True\n", + "plt.rcParams['grid.alpha'] = 0.5\n", + "plt.rcParams['grid.color'] = \"grey\"\n", + "plt.rcParams['grid.linestyle'] = \"--\"\n", + "\n", "\n", - "from sklearn.metrics import confusion_matrix" + "# https://stackoverflow.com/a/20709149\n", + "matplotlib.rcParams['text.usetex'] = True\n", + "matplotlib.rcParams['font.family'] = 'sans-serif'\n", + "matplotlib.rcParams['font.sans-serif'] = 'Helvetica'\n", + "\n", + "matplotlib.rcParams['text.latex.preamble'] = [\n", + " r'\\usepackage{siunitx}', # i need upright \\micro symbols, but you need...\n", + " r'\\sisetup{detect-all}', # ...this to force siunitx to actually use your fonts\n", + " r'\\usepackage{helvet}', # set the normal font here\n", + " r'\\usepackage{sansmath}', # load up the sansmath so that math -> helvet\n", + " r'\\sansmath' # <- tricky! -- gotta actually tell tex to use!\n", + "] " ] }, { @@ -105,8 +135,8 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "## Calculate vertical distribution of wave count\n", - "For each site, calculate how many waves reached a certain elevation (store as a binned histogram)." + "## Calculate impact hours\n", + "- For each site_id, determine the R2 elevation." ] }, { @@ -115,11 +145,306 @@ "metadata": {}, "outputs": [], "source": [ - "# Helper functions\n", - "def find_nearest(array, value):\n", - " array = np.asarray(array)\n", - " idx = np.nanargmin((np.abs(array - value)))\n", - " return array[idx], idx" + "# Create figure to plot results\n", + "fig = tools.make_subplots(\n", + " rows=2,\n", + " cols=2,\n", + " specs=[[{}, {}], [{}, {}]],\n", + " subplot_titles=('Swash/Swash', 'Swash/Collision', 'Collision/Swash',\n", + " 'Collision/Collision'),\n", + " shared_xaxes=True,\n", + " shared_yaxes=True,\n", + " horizontal_spacing=0.05,\n", + " vertical_spacing=0.1,\n", + " print_grid=False)\n", + "\n", + "# Iterate through each site\n", + "print('Calculating cumulate frequency of R_high for each site:')\n", + "site_ids = twls['forecasted']['mean_slope_sto06'].index.get_level_values(\n", + " 'site_id').unique().values\n", + "for site_id in tqdm_notebook(site_ids):\n", + "\n", + " # Put data into a temporary dataframe, shorter name is easier to work with\n", + " df_impacts = impacts['forecasted']['mean_slope_sto06'].loc[site_id]\n", + " df_twls = twls['forecasted']['mean_slope_sto06'].loc[site_id]\n", + "\n", + " D_low = df_impacts.dune_toe_z\n", + " if np.isnan(D_low):\n", + " continue\n", + "\n", + " # Get R_high elevations minus dune toe\n", + " R_high_ts = df_twls.R_high.dropna().values\n", + " R_high_D_low_ts = R_high_ts - D_low\n", + "\n", + " # Get SWL minus dune toe\n", + " SWL_D_low_ts = df_twls['tide'].dropna().values - D_low\n", + " DSWL_D_low_ts = (df_twls['tide'] + df_twls['setup']).dropna().values - D_low\n", + "\n", + " # Get cumulative freq\n", + " cumfreq = stats.cumfreq(R_high_D_low_ts, numbins=100)\n", + "# cumfreq = stats.cumfreq(DSWL_D_low_ts, numbins=100)\n", + "\n", + " # Calculate space of values for x\n", + " bin_vals = cumfreq.lowerlimit + np.linspace(\n", + " 0, cumfreq.binsize * cumfreq.cumcount.size, cumfreq.cumcount.size)\n", + "\n", + " # Check which subplot we should put this site on\n", + " forecasted_regime = impacts['forecasted']['mean_slope_sto06'].loc[\n", + " site_id].storm_regime\n", + " observed_regime = impacts['observed'].loc[site_id].storm_regime\n", + "\n", + " if forecasted_regime == 'swash' and observed_regime == 'swash':\n", + " x_col = 1\n", + " y_col = 1\n", + " elif forecasted_regime == 'collision' and observed_regime == 'collision':\n", + " x_col = 2\n", + " y_col = 2\n", + " elif forecasted_regime == 'swash' and observed_regime == 'collision':\n", + " x_col = 2\n", + " y_col = 1\n", + " elif forecasted_regime == 'collision' and observed_regime == 'swash':\n", + " x_col = 1\n", + " y_col = 2\n", + " else:\n", + " continue\n", + "\n", + " fig.append_trace(\n", + " go.Scattergl(\n", + " x=bin_vals,\n", + " y=[max(cumfreq.cumcount) - x for x in cumfreq.cumcount],\n", + " name=site_id,\n", + " line=dict(\n", + " color=('rgba(22, 22, 22, 0.2)'),\n", + " width=0.5,\n", + " )), x_col, y_col)\n", + "\n", + "print('Finalizing plot:')\n", + "# Change some formatting for the plot\n", + "layout = go.Layout(\n", + " xaxis=dict(domain=[0, 0.45]),\n", + " yaxis=dict(\n", + " domain=[0, 0.45],\n", + " type='log',\n", + " ),\n", + " xaxis2=dict(domain=[0.55, 1]),\n", + " xaxis4=dict(domain=[0.55, 1], anchor='y4'),\n", + " yaxis3=dict(\n", + " domain=[0.55, 1],\n", + " type='log',\n", + " ),\n", + " yaxis4=dict(\n", + " domain=[0.55, 1],\n", + " anchor='x4',\n", + " type='log',\n", + " ))\n", + "\n", + "fig['layout'].update(\n", + " showlegend=False,\n", + " title='Impact hours',\n", + " height=800,\n", + ")\n", + "\n", + "for ax in ['yaxis', 'yaxis2']:\n", + "# fig['layout'][ax]['range'] = [0, 400]\n", + " fig['layout'][ax]['range'] = [0, 15]\n", + "\n", + "for ax in ['xaxis', 'xaxis2']:\n", + "# fig['layout'][ax]['range'] = [-2.5, 2.5]\n", + " fig['layout'][ax]['range'] = [-1, 1]\n", + "\n", + "fig['layout']['xaxis'].update(title='R_high - D_low')\n", + "fig['layout']['xaxis2'].update(title='R_high - D_low')\n", + "fig['layout']['yaxis'].update(title='No. of Hours')\n", + "fig['layout']['yaxis2'].update(title='No. of Hours')\n", + "\n", + "# pio.write_image(fig, 'fig2.png')\n", + "\n", + "go.FigureWidget(fig)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This gives an overview of the distribution of impact hours. Try to calculate the confidence interval bounds for each prediction/observed combination.\n", + "\n", + "The following cell looks at combining all the observations from each CDF into one large CDF and calculating a confidence interval from it, but I'm not sure if this is a valid method." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "code_folding": [ + 3 + ] + }, + "outputs": [], + "source": [ + "from statsmodels.distributions.empirical_distribution import ECDF\n", + "from statsmodels.distributions.empirical_distribution import _conf_set\n", + "\n", + "df_twls = twls['forecasted']['mean_slope_sto06']\n", + "df_forecasted_impacts = impacts['forecasted']['mean_slope_sto06']\n", + "df_observed_impacts = impacts['observed']\n", + "\n", + "plt.figure(figsize=(6,8))\n", + "\n", + "# Do some data rearranging and joining to make it easier\n", + "df_joined = df_twls.reset_index()\n", + "df_joined = df_joined.set_index('site_id')\n", + "df_joined = df_joined.merge(\n", + " df_observed_impacts[['storm_regime']],\n", + " left_on='site_id',\n", + " right_on='site_id').rename({\n", + " 'storm_regime': 'observed_regime'\n", + " },\n", + " axis='columns')\n", + "df_joined = df_joined.merge(\n", + " df_forecasted_impacts[['storm_regime', 'dune_toe_z']],\n", + " left_on='site_id',\n", + " right_on='site_id').rename({\n", + " 'storm_regime': 'forecasted_regime'\n", + " },\n", + " axis='columns')\n", + "\n", + "regime_combinations = [\n", + " ('swash', 'swash', '#2b83ba'),\n", + " ('collision', 'swash', '#abdda4'),\n", + " ('swash', 'collision', '#fdae61'),\n", + " ('collision', 'collision', '#d7191c'),\n", + "]\n", + "\n", + "for comb in regime_combinations:\n", + "\n", + " observed_regime = comb[0]\n", + " forecasted_regime = comb[1]\n", + " color = comb[2]\n", + "\n", + " # Get values of elevation difference to plot\n", + " query = '(observed_regime==\"{}\") & (forecasted_regime==\"{}\")'.format(\n", + " observed_regime, forecasted_regime)\n", + " df = df_joined.query(query)\n", + " R_high_D_low = (df.R_high - df.dune_toe_z).values\n", + " R_high_D_low = R_high_D_low[~np.isnan(R_high_D_low)]\n", + "\n", + " ecdf = ECDF(R_high_D_low)\n", + "\n", + " y = ecdf.y\n", + " lower, upper = _conf_set(y, alpha=0.05)\n", + " x = ecdf.x\n", + "\n", + " avg_hrs = df.groupby('site_id').count().R_high.mean()\n", + " y = [avg_hrs - v * avg_hrs for v in y]\n", + " lower = [avg_hrs - v * avg_hrs for v in lower]\n", + " upper = [avg_hrs - v * avg_hrs for v in upper]\n", + "\n", + " plt.step(\n", + " x,\n", + " y,\n", + " color=color,\n", + " label='Pred={}, Obs={}'.format(forecasted_regime, observed_regime))\n", + " plt.fill_between(\n", + " x, y, upper, color='grey', alpha=0.2, interpolate=False, step='pre')\n", + " plt.fill_between(\n", + " x, y, lower, color='grey', alpha=0.2, interpolate=False, step='pre')\n", + "\n", + "# # Plot for checking\n", + "\n", + "plt.title('Empirical CDF with 95\\% confidence intervals')\n", + "plt.xlabel('$R_{high} - D_{low} (m)$')\n", + "plt.ylabel('Hours of Elevation Exceedence')\n", + "plt.xlim([-1, 1])\n", + "plt.ylim([0, 25])\n", + "plt.legend(loc='best')\n", + "\n", + "# Print to figure\n", + "plt.savefig('05-empirical-cdf.png', dpi=600, bbox_inches='tight') \n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The plot above shows:\n", + "- collision if R_high - D_low > 0.25 m for 6 hours\n", + "- swash if R_high - D_low < -0.8m for 7 hours\n", + "\n", + "additionaly:\n", + "- collision if R_high > D_low for more than 10 hours\n", + " \n", + "Let's test how these new critera would perform." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Calculate elevation exceedence for each hours we are interested in\n", + "ele_exceedence_6hr = twls['forecasted']['mean_slope_sto06'].sort_values(['R_high'],ascending=False).groupby('site_id').R_high.nth(6-1).rename('ele_exceedence_6hr')\n", + "\n", + "ele_exceedence_7hr = twls['forecasted']['mean_slope_sto06'].sort_values(['R_high'],ascending=False).groupby('site_id').R_high.nth(7-1).rename('ele_exceedence_7hr')\n", + "\n", + "\n", + "ele_exceedence_2hr = twls['forecasted']['mean_slope_sto06'].sort_values(['R_high'],ascending=False).groupby('site_id').R_high.nth(2-1).rename('ele_exceedence_2hr')\n", + "\n", + "ele_exceedence_1hr = twls['forecasted']['mean_slope_sto06'].sort_values(['R_high'],ascending=False).groupby('site_id').R_high.nth(0).rename('ele_exceedence_1hr')\n", + "\n", + "\n", + "# Get our dune toes\n", + "dune_toes = df_profile_features_crest_toes.xs('prestorm',level='profile_type')['dune_toe_z']\n", + "\n", + "# Get our observed regimes\n", + "observed_regime = impacts['observed'].storm_regime.rename('observed_regime')\n", + "\n", + "# Concat into one data frame\n", + "df = pd.concat([dune_toes, ele_exceedence_6hr, ele_exceedence_7hr, ele_exceedence_1hr, ele_exceedence_2hr, observed_regime],axis=1)\n", + "\n", + "# Get predicted regime based on old criteria\n", + "df.loc[df.ele_exceedence_1hr < df.dune_toe_z, 'forecasted_regime'] = 'swash'\n", + "df.loc[df.ele_exceedence_1hr > df.dune_toe_z , 'forecasted_regime'] = 'collision'\n", + "\n", + "\n", + "regime_combinations = [\n", + " ('swash','swash'),\n", + " ('collision','swash'),\n", + " ('swash','collision'),\n", + " ('collision','collision'),\n", + "]\n", + "\n", + "print('Original')\n", + "for comb in regime_combinations:\n", + " query = 'forecasted_regime==\"{}\" & observed_regime==\"{}\"'.format(comb[0], comb[1])\n", + " print('Forecasted: {}, Observed: {}, Count: {}'.format(comb[0], comb[1], len(df.query(query))))\n", + "\n", + "\n", + "# Get predicted regime based on our new criteria\n", + "\n", + "adjust_swash_criteria = (df.forecasted_regime == 'swash') & (df.ele_exceedence_7hr - df.dune_toe_z > -0.8)\n", + "adjust_collision_criteria = (df.forecasted_regime == 'collision') & (df.ele_exceedence_6hr - df.dune_toe_z < 0.25)\n", + "df.loc[adjust_swash_criteria, 'forecasted_regime'] = 'collision'\n", + "df.loc[adjust_collision_criteria, 'forecasted_regime'] = 'swash'\n", + "\n", + "# df.loc[(df.ele_exceedence_1hr - df.dune_toe_z <= -0.15 ),'forecasted_regime'] = 'swash'\n", + "# df.loc[(df.ele_exceedence_1hr - df.dune_toe_z > -0.15 ),'forecasted_regime'] = 'collision'\n", + "\n", + "\n", + "print('\\nAfter adjustment')\n", + "for comb in regime_combinations:\n", + " query = 'forecasted_regime==\"{}\" & observed_regime==\"{}\"'.format(comb[0], comb[1])\n", + " print('Forecasted: {}, Observed: {}, Count: {}'.format(comb[0], comb[1], len(df.query(query))))\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Looking at the adjusted values, we can see these criteria actually make it worse. There must be something wrong with the technique - maybe the way of calculating the confidence intervals is wrong? Let's try calculate confidence intervals for each regime combination.\n", + "\n", + "*This cell I don't think is used...*\n" ] }, { @@ -128,7 +453,65 @@ "metadata": {}, "outputs": [], "source": [ - "df_profile_features_crest_toes.loc[(site_id,'prestorm'),'dune_toe_z']" + "def mean_confidence_interval(data, confidence=0.95):\n", + " a = 1.0 * np.array(data)\n", + " n = len(a)\n", + " m, se = np.mean(a), stats.sem(a)\n", + " h = se * stats.t.ppf((1 + confidence) / 2., n-1)\n", + " return m, m-h, m+h\n", + "\n", + "# Add columns indicating how many n hrs was this the largest record\n", + "df = twls['forecasted']['mean_slope_sto06'].sort_values(['R_high'],ascending=False)\n", + "df['n_hrs_largest']= df.groupby('site_id').cumcount()+1\n", + "\n", + "# Join observed and forecast impacts and dune toe elevation\n", + "observed_regime = impacts['observed'].storm_regime.rename('observed_regime').to_frame()\n", + "forecasted_regime = impacts['forecasted']['mean_slope_sto06'].storm_regime.rename('forecasted_regime').to_frame()\n", + "dune_info = df_profile_features_crest_toes.xs('prestorm', level='profile_type')\n", + "\n", + "df['datetime'] = df.index.get_level_values('datetime')\n", + "df = df.merge(observed_regime,left_on=['site_id'],right_on='site_id')\n", + "df = df.merge(forecasted_regime,left_on=['site_id'],right_on='site_id')\n", + "df = df.merge(dune_info,left_on=['site_id'],right_on='site_id')\n", + "\n", + "# Make new column for R_high minus D_low\n", + "df['R_high_D_low_diff'] = df.R_high - df.dune_toe_z\n", + "\n", + "\n", + "regime_combinations = [\n", + " ('swash','swash'),\n", + " ('swash','collision'),\n", + " ('collision','swash'),\n", + " ('collision','collision'),\n", + "]\n", + "\n", + "print('Calculating hr exceedence elevations for each combination:')\n", + "exceedence_data = []\n", + "for hr in tqdm_notebook([x for x in range(1,101)]):\n", + " \n", + " for comb in regime_combinations:\n", + " \n", + " vals = df.loc[(df.n_hrs_largest==hr) & (df.observed_regime==comb[0]) & (df.forecasted_regime==comb[1])].R_high_D_low_diff.dropna().values\n", + " \n", + " ci = mean_confidence_interval(vals)\n", + "\n", + " exceedence_data.append({\n", + " 'observed_regime': comb[0],\n", + " 'forecasted_regime': comb[1],\n", + " 'exceedence_hr': hr,\n", + " 'ci_mean': ci[0],\n", + " 'ci_lower': ci[1],\n", + " 'ci_upper': ci[2],\n", + " })\n", + " \n", + " " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Let's try a different apporach and try split the observed swash and collision regimes at each impact duration hour. " ] }, { @@ -136,6 +519,443 @@ "execution_count": null, "metadata": {}, "outputs": [], + "source": [ + "from scipy.stats import norm\n", + "\n", + "best_split = []\n", + "exceedence_hrs = []\n", + "swash_mean = []\n", + "swash_95_lower = []\n", + "swash_95_upper = []\n", + "collision_mean = []\n", + "collision_95_lower = []\n", + "collision_95_upper = []\n", + "swash_median = []\n", + "swash_q1 = []\n", + "swash_q3 = []\n", + "collision_median = []\n", + "collision_q1 = []\n", + "collision_q3 = []\n", + "\n", + "for hr in tqdm_notebook([x for x in range(1,101)]):\n", + " \n", + " dists = []\n", + " plt.figure(figsize=(10,2))\n", + " for observed_regime in ['swash','collision']:\n", + " \n", + " vals = df.loc[(df.n_hrs_largest==hr) &\n", + " (df.observed_regime==observed_regime)].R_high_D_low_diff.dropna().values\n", + " \n", + " if observed_regime =='collision':\n", + " color = 'red'\n", + " label='collision'\n", + " else:\n", + " color = 'blue'\n", + " label='swash'\n", + " \n", + " plt.hist(vals, bins='auto',color=color, alpha=0.5,label=label) \n", + " plt.title(\"{} hour exceedence TWL\".format(hr))\n", + " plt.xlim([-2.5,2.5])\n", + " \n", + " dists.append(norm.fit(vals))\n", + " \n", + " # Find which elevation best splits swash and collision\n", + "# eles = [x for x in np.linspace(-2,2,1000)]\n", + "# total_cdfs = []\n", + "# for ele in eles:\n", + "# swash_cdf = norm.cdf(ele,*dists[0])\n", + "# collision_cdf = 1 - norm.cdf(ele,*dists[1])\n", + "# total_cdfs.append(swash_cdf + collision_cdf)\n", + "\n", + "# i_max = np.argmax(total_cdfs)\n", + "# best_ele = eles[i_max]\n", + "\n", + "# exceedence_hrs.append(hr)\n", + "# best_split.append(best_ele)\n", + "\n", + " # Find which elevation best splits swash and collision\n", + " eles = [x for x in np.linspace(-2,2,100)]\n", + " total_cdfs = []\n", + " swash_vals = df.loc[(df.n_hrs_largest==hr) &\n", + " (df.observed_regime=='swash')].R_high_D_low_diff.dropna().values\n", + " collision_vals = df.loc[(df.n_hrs_largest==hr) &\n", + " (df.observed_regime=='collision')].R_high_D_low_diff.dropna().values\n", + " for ele in eles:\n", + " swash_samples = np.sum( swash_vals < ele) / len(swash_vals)\n", + " collision_samples = np.sum( collision_vals > ele) / len(collision_vals) \n", + " total_cdfs.append(swash_samples + collision_samples)\n", + " \n", + " i_max = np.argmax(total_cdfs)\n", + " best_ele = eles[i_max]\n", + "\n", + " exceedence_hrs.append(hr)\n", + " best_split.append(best_ele) \n", + " \n", + " \n", + " # Store stastistics\n", + " swash_mean.append(dists[0][0])\n", + " swash_95_lower.append(norm.interval(0.5, *dists[0])[0])\n", + " swash_95_upper.append(norm.interval(0.5, *dists[0])[1])\n", + " collision_mean.append(dists[1][0])\n", + " collision_95_lower.append(norm.interval(0.5, *dists[1])[0])\n", + " collision_95_upper.append(norm.interval(0.5, *dists[1])[1])\n", + " \n", + " swash_median.append(np.percentile(swash_vals, 50))\n", + " swash_q1.append(np.percentile(swash_vals, 25))\n", + " swash_q3.append(np.percentile(swash_vals, 75))\n", + " collision_median.append(np.percentile(collision_vals, 50))\n", + " collision_q1.append(np.percentile(collision_vals, 25))\n", + " collision_q3.append(np.percentile(collision_vals, 75)) \n", + " \n", + " plt.axvline(best_ele, label='Best split (x={:.2f} m)'.format(best_ele))\n", + " plt.legend(loc='upper right', prop={'size': 10} )\n", + " plt.xlabel('$R_{high} - D_{low}$')\n", + " plt.ylabel('No. of sites')\n", + " plt.xlim([-2,2])\n", + " if hr == 80 or hr < 5 or hr==90:\n", + " plt.show()\n", + " \n", + " plt.close()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now, let's plot our distributions for swash/collision and the best seperation between them." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plt.figure(figsize=(5,5))\n", + "plt.plot(best_split,exceedence_hrs, label='Best split', color='#000000', linestyle='--')\n", + "\n", + "# plt.plot(swash_mean, exceedence_hrs, label='Swash', color='#2b83ba')\n", + "# plt.fill_betweenx(exceedence_hrs,swash_95_lower,swash_95_upper, color='#2b83ba', alpha=0.2, interpolate=False)\n", + "\n", + "# plt.plot(collision_mean, exceedence_hrs, label='Collision', color='#d7191c')\n", + "# plt.fill_betweenx(exceedence_hrs,collision_95_lower,collision_95_upper, color='#d7191c', alpha=0.2, interpolate=False,label='plus 50')\n", + "\n", + "\n", + "plt.plot(swash_median, exceedence_hrs, label='Swash', color='#2b83ba')\n", + "plt.fill_betweenx(exceedence_hrs,swash_q1,swash_q3, color='#2b83ba', alpha=0.2, interpolate=False,label='Swash IQR')\n", + "\n", + "plt.plot(collision_median, exceedence_hrs, label='Collision', color='#d7191c')\n", + "plt.fill_betweenx(exceedence_hrs,collision_q1,collision_q3, color='#d7191c', alpha=0.2, interpolate=False,label='Collision IQR')\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "#===\n", + "# # Let's plot one site as well, just to check\n", + "# import random\n", + "# site_ids = list(impacts['observed'].index.unique().values)\n", + "# site_id = random.choice(site_ids)\n", + "\n", + "\n", + "# site_id = 'TREACH0011'\n", + "# site_predicted_regime = impacts['forecasted']['mean_slope_sto06'].loc[site_id].storm_regime\n", + "# site_observed_regime = impacts['observed'].loc[site_id].storm_regime\n", + "# df_site = df.loc[site_id]\n", + "# plt.plot(df_site.R_high_D_low_diff, df_site.n_hrs_largest,label='site_id={}\\n(pred={},obs={})'.format(site_id,site_predicted_regime, site_observed_regime),color='#ffffff', linestyle='--')\n", + "\n", + "\n", + "plt.title('Observed Swash/Collision - Best Split')\n", + "plt.xlabel('$R_{high} - D_{low}$ (m)')\n", + "plt.ylabel('Exceedance hours')\n", + "plt.ylim([0,100])\n", + "plt.xlim([-2,2])\n", + "plt.legend()\n", + "\n", + "# Print to figure\n", + "plt.savefig('05-best-split.png', dpi=600, bbox_inches='tight') \n", + "\n", + "plt.show()\n", + "plt.close()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Plot above shows that if Rhigh = Dlow plus/minus 0.25m, we should say the storm regime is uncertain, rather than trying to make an incorrect prediction." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "df_for = impacts['forecasted']['mean_slope_sto06']\n", + "df_obs = impacts['observed']\n", + "\n", + "# Join forecasted and observed impacts into same dataframe\n", + "df_for = df_for.rename(columns={'storm_regime': 'forecasted_regime'})\n", + "df_obs = df_obs.rename(columns={'storm_regime': 'observed_regime'})\n", + "df_for = df_for.merge(\n", + " df_obs.observed_regime.to_frame(), left_index=True, right_index=True)\n", + "\n", + "# Get wrong forecasts\n", + "incorrect_for = df_for.forecasted_regime != df_for.observed_regime\n", + "\n", + "# How many wrong/correct forecasts\n", + "print('There were {} correct forecasts'.format(len(df_for[~incorrect_for])))\n", + "print('There were {} incorrect forecasts'.format(len(df_for[incorrect_for])))\n", + "print('')\n", + "\n", + "# How many of these forecasts were where R_high was near D_low?\n", + "close_eles = ((df.R_high > df.dune_toe_z - 0.25) &\n", + " (df.R_high < df.dune_toe_z + 0.25))\n", + "\n", + "s = 'R_high and D_low elevations were close at {} correctly forecasted sites'\n", + "print(s.format(len(df_for[~incorrect_for & close_eles])))\n", + "\n", + "s = 'R_high and D_low elevations were close at {} wrongly forecasted sites'\n", + "print(s.format(len(df_for[incorrect_for & close_eles])))\n", + "\n", + "# df[(df.R_high>df.dune_toe_z-0.25)&(df.R_high helvet\n", + " r'\\sansmath', # <- tricky! -- gotta actually tell tex to use!\n", + "] " + ] + }, { "cell_type": "markdown", "metadata": {}, @@ -140,6 +167,12 @@ "# Get x and z at poststorm dune toe for each site\n", "df_dune_toe_poststorm = df_profile_features_crest_toes.xs('poststorm', level='profile_type')[['dune_toe_x','dune_toe_z']]\n", "\n", + "# If there is no poststorm dune toe defined, use the dune crest\n", + "df_dune_crest_poststorm = df_profile_features_crest_toes.xs('poststorm', level='profile_type')[['dune_crest_x','dune_crest_z']]\n", + "df_dune_toe_poststorm.dune_toe_x = df_dune_toe_poststorm.dune_toe_x.fillna(df_dune_crest_poststorm.dune_crest_x)\n", + "df_dune_toe_poststorm.dune_toe_z = df_dune_toe_poststorm.dune_toe_z.fillna(df_dune_crest_poststorm.dune_crest_z)\n", + "\n", + "\n", "# Join df for mhw and dune toe\n", "df = df_mhw_poststorm.join(df_dune_toe_poststorm)\n", "df['beta'] = -(df['dune_toe_z'] - df['z_mhw']) / (df['dune_toe_x'] -df['x_mhw'])\n", @@ -176,26 +209,89 @@ "df" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We also should add the change in beach width between prestorm and post storm profiles" + ] + }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ - "df_data.index" + "ele = 0.7\n", + "data = []\n", + "for site_id, df_site in df_profiles.groupby('site_id'):\n", + " \n", + " # Beach width should be measured from dune toe (or crest if doesn't exist) to MHW\n", + " \n", + " dune_toe_x = np.nanmax([\n", + " df_profile_features_crest_toes.loc[(site_id,'prestorm')].dune_crest_x,\n", + " df_profile_features_crest_toes.loc[(site_id,'prestorm')].dune_toe_x\n", + " ])\n", + " \n", + " \n", + " # TODO This probably should take the closest value to ele starting from the seaward end of the profile\n", + " temp = df_site.xs('prestorm',level='profile_type').dropna(subset=['z'])\n", + " prestorm_width = temp.iloc[(temp.z - ele).abs().argsort()[0]].name[1] - dune_toe_x\n", + " \n", + " temp = df_site.xs('poststorm',level='profile_type').dropna(subset=['z'])\n", + " poststorm_width = temp.iloc[(temp.z - ele).abs().argsort()[0]].name[1] - dune_toe_x\n", + " \n", + " width_change = prestorm_width - poststorm_width\n", + " data.append(\n", + " {\n", + " 'site_id': site_id,\n", + " 'width_change': width_change,\n", + " 'prestorm_width': prestorm_width,\n", + " 'poststorm_width': poststorm_width\n", + " })\n", + " \n", + " \n", + " \n", + " \n", + "df_width_change = pd.DataFrame(data)\n", + "df_width_change = df_width_change.set_index(['site_id'])\n", + "\n", + "# Join with the data\n", + "df = df.merge(df_width_change, left_on=['site_id'], right_on=['site_id'])\n" ] }, { - "cell_type": "markdown", + "cell_type": "code", + "execution_count": null, "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "metadata": { + "heading_collapsed": true + }, "source": [ - "Plot our data" + "## Plot our data in a confusion matrix\n", + "Superseded" ] }, { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "hidden": true + }, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "hidden": true + }, "outputs": [], "source": [ "fig = tools.make_subplots(\n", @@ -279,7 +375,9 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "hidden": true + }, "source": [ "Looking at the above plot:\n", "- In general, we can see that the prestorm mean slope is flatter than the poststorm mean slope. This can be explained by the presence of prestorm berms, which increase the prestorm mean slope. During the storm, these berms get eroded and decrease the slope.\n", @@ -287,6 +385,166 @@ "- **Swash/Collision**: Where we predict collision but observe swash, we can see that the prestorm mean slopes >0.1 generate high TWLs. \n", "\n" ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Plot our data in a confusion matrix\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "df[cc_mask].loc[df[cc_mask].poststorm_beta+0.05< df[cc_mask].prestorm_beta]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "f, ([ax1, ax2], [ax3, ax4],) = plt.subplots(\n", + " 2,\n", + " 2,\n", + " sharey=True,\n", + " sharex=True,\n", + " figsize=(8, 7))\n", + "\n", + "\n", + "ss_mask = (df.observed_regime=='swash') & (df.forecasted_regime=='swash')\n", + "sc_mask = (df.observed_regime=='swash') & (df.forecasted_regime=='collision')\n", + "cs_mask = (df.observed_regime=='collision') & (df.forecasted_regime=='swash')\n", + "cc_mask = (df.observed_regime=='collision') & (df.forecasted_regime=='collision')\n", + "\n", + "# Define colormap for our observations\n", + "cm = plt.cm.get_cmap('plasma')\n", + "\n", + "params = {'edgecolors': '#999999',\n", + " 's': 12,\n", + " 'linewidth': 0.1, \n", + " 'cmap':cm,\n", + " 'vmin':0, \n", + " 'vmax':60\n", + " }\n", + "\n", + "sc=ax1.scatter(df[ss_mask].prestorm_beta, df[ss_mask].poststorm_beta, c=df[ss_mask].width_change,**params)\n", + "ax1.set_title('Swash/Swash')\n", + "ax1.set_ylabel('Observed swash')\n", + "\n", + "ax2.scatter(df[sc_mask].prestorm_beta, df[sc_mask].poststorm_beta, c=df[sc_mask].width_change,**params)\n", + "ax2.set_title('Swash/Collision')\n", + "\n", + "ax3.scatter(df[cs_mask].prestorm_beta, df[cs_mask].poststorm_beta, c=df[cs_mask].width_change,**params)\n", + "ax3.set_title('Collision/Swash')\n", + "ax3.set_ylabel('Observed collision')\n", + "ax3.set_xlabel('Predicted swash')\n", + "\n", + "ax4.scatter(df[cc_mask].prestorm_beta, df[cc_mask].poststorm_beta, c=df[cc_mask].width_change,**params)\n", + "ax4.set_title('Collision/Collision')\n", + "ax4.set_xlabel('Predicted collision')\n", + "\n", + "for ax in [ax1,ax2,ax3,ax4]:\n", + " ax.plot([0,0.2],[0,0.2], 'k--')\n", + " ax.set_xlim([0,0.2])\n", + " ax.set_ylim([0,0.2])\n", + "\n", + " \n", + "# Create a big ax so we can use common axis labels\n", + "# https://stackoverflow.com/a/36542971\n", + "f.add_subplot(111, frameon=False)\n", + "plt.tick_params(labelcolor='none', top=False, bottom=False, left=False, right=False)\n", + "plt.grid(False)\n", + "plt.xlabel(\"Prestorm mean slope (-)\", labelpad=25)\n", + "plt.ylabel(\"Poststorm mean slope (-)\", labelpad=25)\n", + " \n", + "# Layout adjustment\n", + "plt.tight_layout()\n", + "plt.subplots_adjust(hspace=0.25, bottom=0.1,right=0.9)\n", + "\n", + "# Add colorbar\n", + "cbar_ax = f.add_axes([0.95, 0.15, 0.05, 0.7])\n", + "cb = f.colorbar(sc, cax=cbar_ax)\n", + "cb.set_label(r'$\\varDelta$ beach width at MHW (m)')\n", + "\n", + "# Save and show figure\n", + "plt.savefig('06-confusion-change-in-slope.png'.format(beach), dpi=600, bbox_inches='tight') \n", + "plt.show()\n", + "plt.close()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Plot for single beach" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "beach = 'NARRA'\n", + "\n", + "df_beach = df.loc[df.index.str.contains(beach)]\n", + "\n", + "# Get index\n", + "n = [x for x in range(len(df_beach))][::-1]\n", + "n_sites = [x for x in df_beach.index][::-1]\n", + "\n", + "f, (ax1,ax2,ax3,ax4) = plt.subplots(1,4, sharey=True,figsize=(10, 8))\n", + "\n", + "ax1.plot(df_beach.prestorm_beta,n,label='Prestorm slope',color='#4d9221')\n", + "ax1.plot(df_beach.poststorm_beta,n,label='Poststorm slope',color='#c51b7d')\n", + "ax1.set_title('Mean beach slope')\n", + "ax1.legend(loc='center', bbox_to_anchor=(0.5, -0.15))\n", + "\n", + "# Replace yticks with site_ids\n", + "yticks = ax1.get_yticks().tolist()\n", + "yticks = [n_sites[int(y)] if 0 <= y <= len(n_sites) else y for y in yticks ]\n", + "ax1.set_yticklabels(yticks)\n", + "ax1.set_xlabel(r'Slope (-)')\n", + "\n", + "ax2.plot(df_beach.prestorm_width,n,label='Prestorm width',color='#4d9221')\n", + "ax2.plot(df_beach.poststorm_width,n, label='Poststorm width',color='#c51b7d')\n", + "# ax2.set_xlim([200,300])\n", + "ax2.set_xlabel(r'Beach width (m)')\n", + "ax2.set_title('Beach width\\nat MHW')\n", + "ax2.legend(loc='center', bbox_to_anchor=(0.5, -0.15))\n", + "\n", + "ax3.plot(df_beach.width_change,n,color='#999999')\n", + "ax3.set_xlim([0,75])\n", + "ax3.set_title('Change in MHW\\nbeach width')\n", + "ax3.set_xlabel(r'$\\varDelta$ Beach width (m)')\n", + "\n", + "\n", + "ax4.plot(df_beach.poststorm_beta / df_beach.prestorm_beta,n,color='#999999')\n", + "ax4.set_title('Ratio between pre and\\npost storm mean slopes')\n", + "\n", + "plt.tight_layout()\n", + "f.subplots_adjust(top=0.88)\n", + "f.suptitle(beach)\n", + "\n", + "# Print to figure\n", + "plt.savefig('06-change-in-slope-{}.png'.format(beach), dpi=600, bbox_inches='tight') \n", + "plt.show()\n", + "plt.close()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "df_beach" + ] } ], "metadata": { diff --git a/notebooks/07-longshore-plotting.ipynb b/notebooks/07-longshore-plotting.ipynb new file mode 100644 index 0000000..1f5a954 --- /dev/null +++ b/notebooks/07-longshore-plotting.ipynb @@ -0,0 +1,348 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Longshore plots of each beach\n", + "- Need to create a longshore plot of each beach to see how the variables change alongshore." + ] + }, + { + "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", + "from tqdm import tqdm_notebook\n", + "from ipywidgets import widgets, Output\n", + "from IPython.display import display, clear_output, Image, HTML\n", + "from scipy import stats\n", + "from sklearn.metrics import confusion_matrix\n", + "import matplotlib.pyplot as plt\n", + "from scipy.interpolate import interp1d\n", + "from pandas.api.types import CategoricalDtype" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Matplot lib default settings\n", + "plt.rcParams[\"figure.figsize\"] = (10,6)\n", + "plt.rcParams['axes.grid']=True\n", + "plt.rcParams['grid.alpha'] = 0.5\n", + "plt.rcParams['grid.color'] = \"grey\"\n", + "plt.rcParams['grid.linestyle'] = \"--\"\n", + "plt.rcParams['axes.grid']=True\n", + "\n", + "# https://stackoverflow.com/a/20709149\n", + "matplotlib.rcParams['text.usetex'] = True\n", + "\n", + "matplotlib.rcParams['text.latex.preamble'] = [\n", + " r'\\usepackage{siunitx}', # i need upright \\micro symbols, but you need...\n", + " r'\\sisetup{detect-all}', # ...this to force siunitx to actually use your fonts\n", + " r'\\usepackage{helvet}', # set the normal font here\n", + " r'\\usepackage{amsmath}',\n", + " r'\\usepackage{sansmath}', # load up the sansmath so that math -> helvet\n", + " r'\\sansmath', # <- tricky! -- gotta actually tell tex to use!\n", + "] " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Import data" + ] + }, + { + "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": [ + "## Generate plot for each beach" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "beach = 'NARRA'\n", + "\n", + "# Get the dataframe\n", + "df = impacts['forecasted']['mean_slope_sto06']\n", + "df = df.rename(columns={'storm_regime': 'forecasted_regime'})\n", + "\n", + "df_beach = df.loc[df.index.str.contains(beach)]\n", + "\n", + "# Add information about hydrodynamics at max(R_high) time\n", + "df_beach = df_beach.merge(\n", + " twls['forecasted']['mean_slope_sto06'].drop(columns=['R_high', 'R_low']),\n", + " left_on=['site_id', 'datetime'],\n", + " right_on=['site_id', 'datetime'])\n", + "\n", + "# Add information about observed impacts\n", + "obs_impacts = impacts['observed'].rename(columns={\n", + " 'storm_regime': 'observed_regime'\n", + "}).observed_regime.to_frame()\n", + "df_beach = df_beach.merge(obs_impacts, left_on='site_id', right_on='site_id')\n", + "\n", + "# Convert storm regimes to categorical datatype\n", + "cat_type = CategoricalDtype(\n", + " categories=['swash', 'collision', 'overwash', 'inundation'], ordered=True)\n", + "df_beach.forecasted_regime = df_beach.forecasted_regime.astype(cat_type)\n", + "df_beach.observed_regime = df_beach.observed_regime.astype(cat_type)\n", + "\n", + "# Get index\n", + "n = [x for x in range(len(df_beach))][::-1]\n", + "n_sites = [x for x in df_beach.index][::-1]\n", + "\n", + "f, (ax1, ax2, ax3, ax4, ax5, ax6, ax7, ax8) = plt.subplots(\n", + " 1,\n", + " 8,\n", + " sharey=True,\n", + " figsize=(14, 8),\n", + " gridspec_kw={'width_ratios': [4, 4, 2, 2, 2, 2,2,2]})\n", + "\n", + "# Specify colors for storm regimes\n", + "cmap = {\n", + " 'swash': '#1a9850',\n", + " 'collision': '#fee08b',\n", + " 'overwash': '#d73027'\n", + "}\n", + "\n", + "colors = [cmap.get(x) for x in df_beach.observed_regime]\n", + "colors = ['#d73027' if c is None else c for c in colors]\n", + "\n", + "# Plot forecasted and observed storm regime\n", + "ax1.scatter(\n", + " df_beach.observed_regime.cat.codes.replace(-1,np.NaN),\n", + " n,\n", + " color=colors,\n", + " marker='o',\n", + " label='Observed regime')\n", + "\n", + "ax1.scatter(\n", + " df_beach.forecasted_regime.cat.codes.replace(-1,np.NaN),\n", + " n,\n", + " color='b',\n", + " marker='o',\n", + " edgecolors='black',\n", + " facecolors='none',\n", + " label='Forecasted regime')\n", + "\n", + "ax1.set_title('Storm\\nregime')\n", + "ax1.set_xticks([0,1,2,3])\n", + "ax1.set_xticklabels(['swash','collision','overwash','inundation'])\n", + "ax1.tick_params(axis='x', rotation=45)\n", + "ax1.legend(loc='center', bbox_to_anchor=(0.5, -0.15))\n", + "\n", + "# Replace yticks with site_ids\n", + "yticks = ax1.get_yticks().tolist()\n", + "yticks = [n_sites[int(y)] if 0 <= y <= len(n_sites) else y for y in yticks ]\n", + "ax1.set_yticklabels(yticks)\n", + "\n", + "# Water levels\n", + "ax2.plot(df_beach.R_high, n, color='#2c7bb6')\n", + "ax2.plot(df_beach.R_low, n, color='#2c7bb6')\n", + "ax2.fill_betweenx(\n", + " n, df_beach.R_low, df_beach.R_high, alpha=0.2, color='#2c7bb6', label='$R_{low}$ to $R_{high}$')\n", + "\n", + "# Dune elevations\n", + "ax2.plot(df_beach.dune_crest_z, n, color='#fdae61')\n", + "ax2.plot(df_beach.dune_toe_z, n, color='#fdae61')\n", + "ax2.fill_betweenx(\n", + " n, df_beach.dune_toe_z, df_beach.dune_crest_z, alpha=0.2, color='#fdae61', label='$D_{low}$ to $D_{high}$')\n", + "\n", + "ax2.set_title('TWL \\& Dune\\nElevations')\n", + "ax2.legend(loc='center',bbox_to_anchor=(0.5,-0.15))\n", + "ax2.set_xlabel('Elevation (m AHD)')\n", + "\n", + "# Plot R_high - D_low\n", + "ax3.plot(df_beach.R_high - df_beach.dune_toe_z,n,color='#999999')\n", + "ax3.axvline(x=0,color='black',linestyle=':')\n", + "ax3.set_title('$R_{high}$ - $D_{low}$')\n", + "ax3.set_xlabel('Height (m)')\n", + "ax3.set_xlim([-2,2])\n", + "\n", + "# Wave height, wave period, beach slope\n", + "ax4.plot(df_beach.Hs0, n,color='#377eb8')\n", + "ax4.set_title('$H_{s0}$')\n", + "ax4.set_xlabel('Sig. wave height (m)')\n", + "ax4.set_xlim([3,5])\n", + "\n", + "ax5.plot(df_beach.Tp, n,color='#e41a1c')\n", + "ax5.set_title('$T_{p}$')\n", + "ax5.set_xlabel('Peak wave period (s)')\n", + "ax5.set_xlim([8,14])\n", + "\n", + "ax6.plot(df_beach.tide, n,color='#a6cee3')\n", + "ax6.set_title('Tide')\n", + "ax6.set_xlabel('Elevation (m AHD)')\n", + "ax6.set_xlim([0,2])\n", + "\n", + "ax7.plot(df_beach.beta, n,color='#4daf4a')\n", + "ax7.set_title(r'$\\beta$')\n", + "ax7.set_xlabel('Mean prestorm\\nbeach slope')\n", + "ax7.set_xlim([0,0.15])\n", + "\n", + "ax8.plot(df_beach.R2, n,color='#6a3d9a')\n", + "ax8.set_title(r'$R_{2\\%}$')\n", + "ax8.set_xlabel('Height (m)')\n", + "\n", + "plt.tight_layout()\n", + "f.subplots_adjust(top=0.88)\n", + "f.suptitle(beach)\n", + "\n", + "\n", + "# Print to figure\n", + "plt.savefig('07-{}.png'.format(beach), dpi=600, bbox_inches='tight') \n", + "\n", + "plt.show()\n", + "plt.close()" + ] + } + ], + "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 +} diff --git a/notebooks/08-narr-topo-bathy-slope-test.ipynb b/notebooks/08-narr-topo-bathy-slope-test.ipynb new file mode 100644 index 0000000..2e02e23 --- /dev/null +++ b/notebooks/08-narr-topo-bathy-slope-test.ipynb @@ -0,0 +1,767 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Narrabeen Slope Test\n", + "With full topo and bathy combined" + ] + }, + { + "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", + "from tqdm import tqdm_notebook\n", + "from ipywidgets import widgets, Output\n", + "from IPython.display import display, clear_output, Image, HTML\n", + "from scipy import stats\n", + "from sklearn.metrics import confusion_matrix\n", + "import matplotlib.pyplot as plt\n", + "from scipy.interpolate import interp1d\n", + "from pandas.api.types import CategoricalDtype\n", + "from scipy.interpolate import UnivariateSpline\n", + "from shapely.geometry import Point, LineString" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Matplot lib default settings\n", + "plt.rcParams[\"figure.figsize\"] = (10,6)\n", + "plt.rcParams['axes.grid']=True\n", + "plt.rcParams['grid.alpha'] = 0.5\n", + "plt.rcParams['grid.color'] = \"grey\"\n", + "plt.rcParams['grid.linestyle'] = \"--\"\n", + "plt.rcParams['axes.grid']=True\n", + "\n", + "# https://stackoverflow.com/a/20709149\n", + "matplotlib.rcParams['text.usetex'] = True\n", + "\n", + "matplotlib.rcParams['text.latex.preamble'] = [\n", + " r'\\usepackage{siunitx}', # i need upright \\micro symbols, but you need...\n", + " r'\\sisetup{detect-all}', # ...this to force siunitx to actually use your fonts\n", + " r'\\usepackage{helvet}', # set the normal font here\n", + " r'\\usepackage{amsmath}',\n", + " r'\\usepackage{sansmath}', # load up the sansmath so that math -> helvet\n", + " r'\\sansmath', # <- tricky! -- gotta actually tell tex to use!\n", + "] " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Import .csv data" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "data_filename = '08-narr-topo-bathy-slope-test-full-profiles.csv'\n", + "\n", + "df_profiles = pd.read_csv(data_filename).set_index(['site_id','x'])\n", + "df_profiles = df_profiles[~df_profiles.index.duplicated(keep='first')]\n", + "print('df_profiles:')\n", + "df_profiles.head()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Manually cut off the prestorm topo \n", + "cuts = {'NARRA0004': {'prestorm_topo_max_x': 330,\n", + " 'poststorm_topo_max_x': 250},\n", + " 'NARRA0008': {'prestorm_topo_max_x': 290,\n", + " 'poststorm_topo_max_x': 250},\n", + " 'NARRA0012': {'prestorm_topo_max_x': 300,\n", + " 'poststorm_topo_max_x': 250},\n", + " 'NARRA0016': {'prestorm_topo_max_x': 300,\n", + " 'poststorm_topo_max_x': 225},\n", + " 'NARRA0021': {'prestorm_topo_max_x': 280,\n", + " 'poststorm_topo_max_x': 225},\n", + " 'NARRA0023': {'prestorm_topo_max_x': 275,\n", + " 'poststorm_topo_max_x': 215},\n", + " 'NARRA0027': {'prestorm_topo_max_x': 260,\n", + " 'poststorm_topo_max_x': 225},\n", + " 'NARRA0031': {'prestorm_topo_max_x': 260,\n", + " 'poststorm_topo_max_x': 225},\n", + " }\n", + "\n", + "for site_id in cuts:\n", + " mask1 = df_profiles.index.get_level_values('site_id') == site_id\n", + " mask2 = df_profiles.index.get_level_values('x') > cuts[site_id]['prestorm_topo_max_x']\n", + " df_profiles.loc[(mask1)&(mask2), 'pre_topo'] = np.nan\n", + " \n", + " mask3 = df_profiles.index.get_level_values('x') > cuts[site_id]['poststorm_topo_max_x']\n", + " df_profiles.loc[(mask1)&(mask3), 'post_topo'] = np.nan\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# for site_id,df_site in df_profiles.groupby('site_id'):\n", + "# f, (ax1) = plt.subplots(1,1, figsize=(6, 3))\n", + "# ax1.set_title(site_id)\n", + " \n", + "# ax1.plot(df_site.index.get_level_values('x'),\n", + "# df_site.pre_topo,\n", + "# label='Pre Topo',\n", + "# color='#2c7bb6')\n", + "# ax1.plot(df_site.index.get_level_values('x'),\n", + "# df_site.pre_bathy,\n", + "# label='Pre Bathy',\n", + "# color='#abd9e9')\n", + "\n", + "# ax1.plot(df_site.index.get_level_values('x'),\n", + "# df_site.post_topo,\n", + "# label='Post Topo',\n", + "# color='#d7191c')\n", + "# ax1.plot(df_site.index.get_level_values('x'),\n", + "# df_site.post_bathy,\n", + "# label='Post Bathy',\n", + "# color='#fdae61')\n", + "\n", + "# ax1.legend()\n", + "# plt.show()\n", + "# plt.close()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "df_profiles = df_profiles.dropna(\n", + " subset=['post_topo', 'post_bathy', 'pre_bathy', 'pre_topo'], how='all')\n", + "\n", + "df_profiles = df_profiles.reset_index()\n", + "df_profiles = df_profiles.melt(id_vars=['site_id','x','lat','lon'],\n", + " value_vars=['post_topo','post_bathy','pre_bathy','pre_topo']).rename(columns={'variable':'profile_type', 'value':'z'})\n", + "\n", + "df_profiles = df_profiles.dropna(subset=['z'])\n", + "\n", + "df_profiles.loc[df_profiles.profile_type=='post_topo','profile_type']='poststorm'\n", + "df_profiles.loc[df_profiles.profile_type=='post_bathy','profile_type']='poststorm'\n", + "df_profiles.loc[df_profiles.profile_type=='pre_topo','profile_type']='prestorm'\n", + "df_profiles.loc[df_profiles.profile_type=='pre_bathy','profile_type']='prestorm'\n", + "\n", + "df_profiles = df_profiles.set_index(['site_id', 'profile_type', 'x'])\n", + "df_profiles = df_profiles[~df_profiles.index.duplicated(keep='first')]\n", + "\n", + "df_profiles = df_profiles.sort_index()\n", + "\n", + "print('df_profiles:')\n", + "df_profiles.head()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Just plots each site's x and z values\n", + "for site_id,df_site in df_profiles.groupby('site_id'):\n", + " f, (ax1) = plt.subplots(1,1, figsize=(6, 3))\n", + " ax1.set_title(site_id)\n", + " \n", + " prestorm=df_site.index.get_level_values('profile_type') == 'prestorm'\n", + " ax1.plot(df_site[prestorm].index.get_level_values('x'),\n", + " df_site[prestorm].z,\n", + " label='Pre Topo',\n", + " color='#2c7bb6')\n", + "\n", + " \n", + " poststorm=df_site.index.get_level_values('profile_type') == 'poststorm'\n", + " ax1.plot(df_site[poststorm].index.get_level_values('x'),\n", + " df_site[poststorm].z,\n", + " label='Post Topo',\n", + " color='#d7191c')\n", + "\n", + "\n", + " ax1.legend()\n", + " plt.show()\n", + " plt.close()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Get dune faces" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "code_folding": [] + }, + "outputs": [], + "source": [ + "# Manually define dune x coordinates and work out slope\n", + "\n", + "dune_data = [\n", + " {\n", + " 'site_id': 'NARRA0004',\n", + " 'dune_crest_x': 180,\n", + " 'dune_toe_x': 205\n", + " },\n", + " {\n", + " 'site_id': 'NARRA0008',\n", + " 'dune_crest_x': 180,\n", + " 'dune_toe_x': 205\n", + " },\n", + " {\n", + " 'site_id': 'NARRA0012',\n", + " 'dune_crest_x': 195,\n", + " 'dune_toe_x': 205\n", + " },\n", + " {\n", + " 'site_id': 'NARRA0016',\n", + " 'dune_crest_x': 190,\n", + " 'dune_toe_x': 200\n", + " },\n", + " {\n", + " 'site_id': 'NARRA0021',\n", + " 'dune_crest_x': 205,\n", + " 'dune_toe_x': 210\n", + " },\n", + " {\n", + " 'site_id': 'NARRA0023',\n", + " 'dune_crest_x': 205,\n", + " 'dune_toe_x': 215\n", + " },\n", + " {\n", + " 'site_id': 'NARRA0027',\n", + " 'dune_crest_x': 210,\n", + " 'dune_toe_x': 219\n", + " },\n", + " {\n", + " 'site_id': 'NARRA0031',\n", + " 'dune_crest_x': 210,\n", + " 'dune_toe_x': 218\n", + " },\n", + "]\n", + "\n", + "for site_dune in dune_data:\n", + " df_site = df_profiles.xs(site_dune['site_id'], level='site_id').xs('prestorm',level='profile_type')\n", + " \n", + " dune_crest_x = site_dune['dune_crest_x']\n", + " dune_toe_x = site_dune['dune_toe_x']\n", + " dune_crest_z = df_site.iloc[df_site.index.get_loc(site_dune['dune_crest_x'],method='nearest')].z\n", + " dune_toe_z = df_site.iloc[df_site.index.get_loc(site_dune['dune_toe_x'],method='nearest')].z\n", + "\n", + " dune_slope = (dune_crest_z - dune_toe_z)/(dune_crest_x - dune_toe_x)\n", + " \n", + " site_dune['dune_crest_z'] = dune_crest_z\n", + " site_dune['dune_toe_z'] = dune_toe_z\n", + " site_dune['dune_slope'] = dune_slope\n", + " \n", + " \n", + "# Join back into main data\n", + "df_dunes = pd.DataFrame(dune_data).set_index('site_id')\n", + "print('df_dunes:')\n", + "df_dunes.head()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# # Just plots each site's x and z values\n", + "# for site_id,df_site in df_profiles.xs('prestorm',level='profile_type').groupby('site_id'):\n", + "# f, (ax1) = plt.subplots(1,1, figsize=(6, 3))\n", + "# ax1.set_title(site_id)\n", + "# ax1.plot(df_site.index.get_level_values('x'),\n", + "# df_site.z)\n", + "# ax1.plot([df_dunes.loc[site_id].dune_crest_x, df_dunes.loc[site_id].dune_toe_x],\n", + "# [df_dunes.loc[site_id].dune_crest_z, df_dunes.loc[site_id].dune_toe_z],\n", + "# 'r.-')\n", + "# ax1.set_xlim([150,250])\n", + "# ax1.set_ylim([0,15])\n", + "# plt.show()\n", + "# plt.close()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Get prestorm slope" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "z_ele = 0.7\n", + "debug=False\n", + "\n", + "def find_nearest_idx(array, value):\n", + " array = np.asarray(array)\n", + " idx = (np.abs(array - value)).argmin()\n", + " return idx\n", + "\n", + "prestorm_slope_data =[]\n", + "for site_id, df_site in df_profiles.xs('prestorm',level='profile_type').groupby('site_id'):\n", + " \n", + " # Find index of our z_ele\n", + " idx = np.where(df_site.z.values>=z_ele)[0][-1]\n", + " \n", + " prestorm_end_x = df_site.iloc[idx].name[1]\n", + " prestorm_end_z = df_site.iloc[idx].z\n", + " \n", + " prestorm_start_x = df_dunes.loc[site_id].dune_toe_x\n", + " prestorm_start_z = df_dunes.loc[site_id].dune_toe_z\n", + " \n", + " prestorm_slope = (prestorm_end_z-prestorm_start_z)/(prestorm_end_x-prestorm_start_x)\n", + " \n", + " prestorm_slope_data.append({\n", + " 'site_id': site_id,\n", + " 'prestorm_end_x': prestorm_end_x,\n", + " 'prestorm_end_z': prestorm_end_z,\n", + " 'prestorm_start_x': prestorm_start_x,\n", + " 'prestorm_start_z': prestorm_start_z,\n", + " 'prestorm_slope': prestorm_slope\n", + " })\n", + " \n", + "df_prestorm_slope = pd.DataFrame(prestorm_slope_data).set_index(['site_id'])\n", + "print('df_prestorm_slope:')\n", + "df_prestorm_slope.head()\n", + " " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Get shelf slope\n", + "At 10 m contour" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "code_folding": [] + }, + "outputs": [], + "source": [ + "# Elevation to take shelf slope at\n", + "z_ele = -9\n", + "debug=False\n", + "\n", + "def find_nearest_idx(array, value):\n", + " array = np.asarray(array)\n", + " idx = (np.abs(array - value)).argmin()\n", + " return idx\n", + "\n", + "def slope_at_point(x, z, z_ele,debug=False):\n", + " # Smooth profile a bit\n", + " # TODO the smoothing factor will change based on the number of data points\n", + " # Need to fix\n", + " s = UnivariateSpline(x, z, s=50)\n", + " xs = np.linspace(min(x),max(x),1000)\n", + " zs = s(xs)\n", + "\n", + " # Calculate derivates of spline\n", + " dzdx = np.diff(zs)/np.diff(xs)\n", + "\n", + " # Find index of z_ele\n", + " idx = find_nearest_idx(zs, z_ele)\n", + " slope = dzdx[idx]\n", + " shelf_x = xs[idx]\n", + "\n", + "\n", + " \n", + " # For checking how much smoothing is going on\n", + " if debug:\n", + " f, (ax1) = plt.subplots(1,1, figsize=(6, 3))\n", + " ax1.plot(x,z)\n", + " ax1.plot(xs,zs)\n", + " plt.show()\n", + " plt.close()\n", + " \n", + " return slope, shelf_x, z_ele\n", + " \n", + "shelf_data = []\n", + "for site_id, df_site in df_profiles.xs('prestorm',level='profile_type').groupby('site_id'):\n", + " shelf_slope, shelf_x, shelf_z = slope_at_point(df_site.index.get_level_values('x').values,\n", + " df_site.z, \n", + " z_ele, debug=debug)\n", + " shelf_data.append({\n", + " 'site_id': site_id,\n", + " 'shelf_slope': shelf_slope,\n", + " 'shelf_x': shelf_x,\n", + " 'shelf_z': shelf_z\n", + " })\n", + " \n", + "df_shelf = pd.DataFrame(shelf_data).set_index(['site_id'])\n", + "\n", + "df_shelf.loc['NARRA0004','shelf_slope'] = -0.02\n", + "\n", + "print('df_shelf:')\n", + "df_shelf.head()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Do geometry\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "df_site" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "for site_id, df_site in df_profiles.groupby('site_id'):\n", + "\n", + " # Project the dune face outwards\n", + " dune_face_toe = Point(df_dunes.loc[site_id].dune_toe_x,\n", + " df_dunes.loc[site_id].dune_toe_z)\n", + " dune_face_sea = Point(\n", + " df_dunes.loc[site_id].dune_toe_x + 1000,\n", + " # df_dunes.loc[site_id].dune_toe_z +1000 * -1\n", + " df_dunes.loc[site_id].dune_toe_z +\n", + " 1000 * df_dunes.loc[site_id].dune_slope)\n", + " dune_line = LineString([dune_face_toe, dune_face_sea])\n", + "\n", + " # Project the shelf slope landwards\n", + " shelf_point = Point(df_shelf.loc[site_id].shelf_x,\n", + " df_shelf.loc[site_id].shelf_z)\n", + " shelf_land = Point(\n", + " df_shelf.loc[site_id].shelf_x - 1000, df_shelf.loc[site_id].shelf_z -\n", + " 1000 * df_shelf.loc[site_id].shelf_slope)\n", + " shelf_sea = Point(\n", + " df_shelf.loc[site_id].shelf_x + 1000, df_shelf.loc[site_id].shelf_z +\n", + " 1000 * df_shelf.loc[site_id].shelf_slope)\n", + " shelf_line = LineString([shelf_land, shelf_point, shelf_sea])\n", + "\n", + " # Find intersection between to lines\n", + " dune_shelf_int = dune_line.intersection(shelf_line)\n", + " dist_toe_to_int = dune_face_toe.distance(dune_shelf_int)\n", + "\n", + " # Plots\n", + " f, (ax1) = plt.subplots(1, 1, figsize=(12, 4))\n", + "\n", + " # Raw profile prestorm\n", + " ax1.plot(\n", + " df_site.xs('prestorm',\n", + " level='profile_type').index.get_level_values('x'),\n", + " df_site.xs('prestorm', level='profile_type').z,\n", + " label='Prestorm profile')\n", + "\n", + " # Raw profile poststorm\n", + " ax1.plot(\n", + " df_site.xs('poststorm',\n", + " level='profile_type').index.get_level_values('x'),\n", + " df_site.xs('poststorm', level='profile_type').z,\n", + " label='Poststorm profile')\n", + "\n", + " # Dune face\n", + " ax1.plot(\n", + " [df_dunes.loc[site_id].dune_crest_x, df_dunes.loc[site_id].dune_toe_x],\n", + " [df_dunes.loc[site_id].dune_crest_z, df_dunes.loc[site_id].dune_toe_z],\n", + " linestyle=':',\n", + " color='#999999',\n", + " label='Dune face ({:.2f})'.format(-df_dunes.loc[site_id].dune_slope))\n", + "\n", + " # Projected dune face\n", + " ax1.plot(\n", + " dune_line.xy[0],\n", + " dune_line.xy[1],\n", + " linestyle='--',\n", + " color='#999999',\n", + " label='Dune face (projected)')\n", + "\n", + " # Projected shelf slope\n", + " ax1.plot(\n", + " shelf_line.xy[0],\n", + " shelf_line.xy[1],\n", + " linestyle='--',\n", + " color='#999999',\n", + " label='Shelf slope (projected)')\n", + "\n", + " # Intersection\n", + " ax1.scatter(\n", + " dune_shelf_int.xy[0],\n", + " dune_shelf_int.xy[1],\n", + " marker='x',\n", + " color='#999999',\n", + " label='Dune/shelf projected intersection')\n", + "\n", + " # Prestorm slope\n", + " ax1.plot([\n", + " df_prestorm_slope.loc[site_id].prestorm_start_x,\n", + " df_prestorm_slope.loc[site_id].prestorm_end_x\n", + " ], [\n", + " df_prestorm_slope.loc[site_id].prestorm_start_z,\n", + " df_prestorm_slope.loc[site_id].prestorm_end_z\n", + " ],\n", + " color='violet',\n", + " label='Prestorm slope ({:.2f})'.format(\n", + " -df_prestorm_slope.loc[site_id].prestorm_slope))\n", + "\n", + " # # Find best slope based on distance form toe to intersection?\n", + " # best_slope_toe = shelf_line.interpolate(\n", + " # shelf_line.project(intersection) - 4 * dist_toe_to_int)\n", + " # best_slope = (dune_face_toe.xy[1][0] - best_slope_toe.xy[1][0]) / (\n", + " # dune_face_toe.xy[0][0] - best_slope_toe.xy[0][0])\n", + "\n", + " # # Best slope toe\n", + " # ax1.scatter(\n", + " # best_slope_toe.xy[0], best_slope_toe.xy[1], marker='o', color='g')\n", + "\n", + " # # Best slope\n", + " # ax1.plot([dune_face_toe.xy[0], best_slope_toe.xy[0]],\n", + " # [dune_face_toe.xy[1], best_slope_toe.xy[1]],\n", + " # color='g',\n", + " # label='Best slope ({:.3f})'.format(-best_slope))\n", + "\n", + " # Find best slope based on intersection of prestorm slope and surf zone slope\n", + " prestorm_slope_line = LineString([\n", + " Point(\n", + " df_prestorm_slope.loc[site_id].prestorm_start_x,\n", + " df_prestorm_slope.loc[site_id].prestorm_start_z,\n", + " ),\n", + " Point(\n", + " df_prestorm_slope.loc[site_id].prestorm_start_x + 10000,\n", + " df_prestorm_slope.loc[site_id].prestorm_start_z +\n", + " 10000 * df_prestorm_slope.loc[site_id].prestorm_slope)\n", + " ])\n", + "\n", + " # Where prestorm slope projection intersects shelf line\n", + " prestorm_slope_shelf_int = prestorm_slope_line.intersection(shelf_line)\n", + "\n", + " # Distance between dune/shelf intersection and prestorm/shelf intersection\n", + " dist_shelf_prestorm_ints = prestorm_slope_shelf_int.distance(\n", + " dune_shelf_int)\n", + "\n", + " best_slope_pt = shelf_line.interpolate(\n", + " shelf_line.project(dune_shelf_int) + 0.3 * (shelf_line.project(prestorm_slope_shelf_int) -\n", + " shelf_line.project(dune_shelf_int)))\n", + " \n", + " best_slope =(df_prestorm_slope.loc[site_id].prestorm_start_z-best_slope_pt.xy[1][0])/(df_prestorm_slope.loc[site_id].prestorm_start_x-best_slope_pt.xy[0][0])\n", + " \n", + " if not prestorm_slope_shelf_int.is_empty:\n", + " ax1.plot(\n", + " prestorm_slope_shelf_int.xy[0],\n", + " prestorm_slope_shelf_int.xy[1],\n", + " marker='x',\n", + " color='#999999',\n", + " label='Prestorm slope/shelf\\nprojected intersection')\n", + " ax1.plot(\n", + " prestorm_slope_line.xy[0],\n", + " prestorm_slope_line.xy[1],\n", + " color='#999999',\n", + " linestyle='--',\n", + " label='Prestorm slope projected line')\n", + " ax1.plot(\n", + " [df_prestorm_slope.loc[site_id].prestorm_start_x,\n", + " best_slope_pt.xy[0][0]],\n", + " [df_prestorm_slope.loc[site_id].prestorm_start_z,\n", + " best_slope_pt.xy[1][0]],\n", + " color='red',\n", + " linestyle='--',\n", + " label='Best slope ({:.3f})'.format(-best_slope))\n", + " \n", + " # TEMP Target slopes\n", + " target_slopes = {\n", + " 'NARRA0004': 0.076,\n", + " 'NARRA0008': 0.093,\n", + " 'NARRA0012': 0.060,\n", + " 'NARRA0016': 0.11,\n", + " 'NARRA0021': 0.063,\n", + " 'NARRA0023': 0.061,\n", + " 'NARRA0027': 0.060,\n", + " 'NARRA0031': 0.057,\n", + " }\n", + "\n", + " target_direction = {\n", + " 'NARRA0004': \"flatter\",\n", + " 'NARRA0008': \"steeper\",\n", + " 'NARRA0012': \"flatter\",\n", + " 'NARRA0016': \"flatter\",\n", + " 'NARRA0021': \"steeper\",\n", + " 'NARRA0023': \"steeper\",\n", + " 'NARRA0027': \"steeper\",\n", + " 'NARRA0031': \"steeper\",\n", + " }\n", + " ax1.plot([dune_face_toe.xy[0][0], dune_face_toe.xy[0][0] + 1000], [\n", + " dune_face_toe.xy[1][0],\n", + " dune_face_toe.xy[1][0] - 1000 * target_slopes[site_id]\n", + " ],\n", + " color='red',\n", + " label='Target slope\\n({} than {:.3f})'.format(\n", + " target_direction[site_id], target_slopes[site_id]))\n", + "\n", + " ax1.set_xlim([100, 800])\n", + " ax1.set_ylim([-15, 12])\n", + "# ax1.set_xlim([100, 600])\n", + "# ax1.set_ylim([-10, 12])\n", + "\n", + " # ax1.set_xlim([df_dunes.loc[site_id].dune_crest_x - 50,\n", + " # intersection.xy[0][0] + 50])\n", + " # ax1.set_ylim([intersection.xy[1][0] -3,\n", + " # df_dunes.loc[site_id].dune_crest_z + 3])\n", + "\n", + " ax1.set_title(site_id)\n", + " ax1.legend(loc='upper right', prop={'size': 10})\n", + " f.savefig('08-{}.png'.format(site_id), dpi=600)\n", + " plt.show()\n", + " plt.close()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "dune_shelf_int" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "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 +}