{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "## Investigate \"collision protection volume\" concept" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "ExecuteTime": { "end_time": "2018-12-03T00:37:28.874517Z", "start_time": "2018-12-03T00:37:28.528594Z" } }, "outputs": [], "source": [ "%matplotlib inline\n", "%reload_ext autoreload\n", "%autoreload" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "ExecuteTime": { "end_time": "2018-12-03T00:37:30.229243Z", "start_time": "2018-12-03T00:37:28.875519Z" } }, "outputs": [], "source": [ "from IPython.core.debugger import set_trace\n", "\n", "import pandas as pd\n", "import numpy as np\n", "import os\n", "\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", "import plotly.io as pio" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Load data\n", "Load data from the `./data/interim/` folder and parse into `pandas` dataframes." ] }, { "cell_type": "code", "execution_count": 83, "metadata": { "ExecuteTime": { "end_time": "2018-12-03T02:20:29.132099Z", "start_time": "2018-12-03T02:20:22.217681Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Importing profiles.csv\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "C:\\Users\\z5189959\\AppData\\Local\\Continuum\\anaconda3\\lib\\site-packages\\numpy\\lib\\arraysetops.py:472: FutureWarning:\n", "\n", "elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison\n", "\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Importing profile_features.csv\n", "Importing impacts_forecasted_foreshore_slope_sto06.csv\n", "Importing impacts_forecasted_mean_slope_sto06.csv\n", "Importing impacts_observed.csv\n", "Importing twl_foreshore_slope_sto06.csv\n", "Importing twl_mean_slope_sto06.csv\n", "Done!\n" ] } ], "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 = df_from_csv('profile_features.csv', index_col=[0])\n", "\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", "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", "\n", "print('Done!')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Lets define a function to calculate the \"collision protection volume\" based on prestorm profiles." ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "ExecuteTime": { "end_time": "2018-12-03T00:37:37.244970Z", "start_time": "2018-12-03T00:37:36.924017Z" }, "code_folding": [] }, "outputs": [], "source": [ "from shapely.geometry import Point, LineString, Polygon\n", "\n", "def collision_protection_vol(x,z, d_low_x, d_low_z, lower_z, angle):\n", " # First, get the bounding line strings of our protection volume\n", " lower_line = LineString([Point(min(x), lower_z), Point(max(x), lower_z)])\n", " profile_line = LineString([Point(x_coord, z_coord) for x_coord, z_coord in zip(x,z) \n", " if all([not np.isnan(x_coord), not np.isnan(z_coord)])])\n", " slope_line = LineString([Point(d_low_x, d_low_z), \n", " Point(max(x), d_low_z - max(x) * np.sin(np.deg2rad(angle)))])\n", "\n", " # Work out where our lower line and slope line intersect\n", " lower_profile_intersection = lower_line.intersection(profile_line)\n", " d_protected_intersection = lower_line.intersection(slope_line)\n", " \n", " # Define the perimeter of the protection area\n", " profile_protected = LineString([Point(x_coord, z_coord) for x_coord, z_coord \n", " in zip(profile_line.xy[0],profile_line.xy[1]) \n", " if d_low_x < x_coord < lower_profile_intersection.xy[0][0]] \n", " + [lower_profile_intersection]\n", " + [d_protected_intersection]\n", " + [Point(d_low_x, d_low_z)])\n", " \n", " # Convert to polygon and return the area (m3/m)\n", " protection_area_poly = Polygon(profile_protected)\n", " protection_area_vol = protection_area_poly.area\n", " return protection_area_vol\n", "\n", "site_id = 'NARRA0018'\n", "profile_type = 'prestorm'\n", "query = \"site_id == '{}' and profile_type == '{}'\".format(site_id, profile_type)\n", "prestorm_profile = df_profiles.query(query)\n", "profile_features = df_profile_features.query(query)\n", "\n", "x = prestorm_profile.index.get_level_values('x')\n", "z = prestorm_profile.z\n", "d_low_x = profile_features.dune_toe_x.tolist()[0]\n", "d_low_z = profile_features.dune_toe_z.tolist()[0]\n", "angle = 60 # degrees from the horizontal\n", "lower_z = 0.5 # from mhw\n", "\n", "vol = collision_protection_vol(x,z, d_low_x, d_low_z, lower_z, angle)" ] }, { "cell_type": "code", "execution_count": 58, "metadata": { "ExecuteTime": { "end_time": "2018-12-03T01:24:18.270476Z", "start_time": "2018-12-03T01:24:18.119986Z" } }, "outputs": [], "source": [ "from datetime import timedelta\n", "\n", "def wl_time(t, wl, z_lower, z_upper):\n", " \"\"\"\n", " Returns the amount of time the water level is between two elevations.\n", " \"\"\"\n", " df_wl = pd.DataFrame.from_records([(t_val, R2_val) for t_val, R2_val in zip(t,R2)], columns=['datetime','wl'])\n", " df_wl.set_index(pd.DatetimeIndex(df_wl['datetime']),inplace=True)\n", " df_wl.drop(columns=['datetime'], inplace=True)\n", " \n", " # Assumes that each record is one hour... probably need to check this\n", " hours = len(df_wl.query('{} < wl < {}'.format(z_lower, z_upper)))\n", " return timedelta(hours=hours)\n", "\n", "def wave_power(t, wl, z_lower, z_upper, Hs0, Tp):\n", " \"\"\"\n", " Returns the cumulative wave power when the water level is between two elevations.\n", " \"\"\"\n", " df_wl = pd.DataFrame.from_records([(t_val, R2_val,Hs0_val,Tp_val) for t_val, R2_val,Hs0_val,Tp_val in zip(t,R2,Hs0,Tp)], columns=['datetime','wl', 'Hs0','Tp'])\n", " df_wl.set_index(pd.DatetimeIndex(df_wl['datetime']),inplace=True)\n", " df_wl.drop(columns=['datetime'], inplace=True)\n", " \n", " # Assumes that each record is one hour... probably need to check this\n", " rho = 1025 # kg/m3\n", " g = 9.8 # m/s2\n", " df_wl_times = df_wl.query('{} < wl < {}'.format(z_lower, z_upper))\n", " power = rho * g ** 2 / 64 / np.pi * df_wl_times.Hs0 ** 2 * df_wl_times.Tp\n", " return power.sum()\n", "\n", "df_twl = twls['forecasted']['mean_slope_sto06']\n", "df_twl_site = df_twl.query(\"site_id == '{}'\".format(site_id))\n", "\n", "R2 = df_twl_site.R2.tolist()\n", "t = df_twl_site.index.get_level_values('datetime')\n", "z_lower = 0.5\n", "z_upper = d_low_z\n", "\n", "exposed_time = wl_time(t, R2, z_lower,z_upper)\n", "\n", "Hs0 = df_twl.Hs0.tolist()\n", "Tp = df_twl.Tp.tolist()\n", "wave_p = wave_power(t, R2, z_lower,z_upper,Hs0, Tp)\n" ] }, { "cell_type": "code", "execution_count": 59, "metadata": { "ExecuteTime": { "end_time": "2018-12-03T01:24:24.215223Z", "start_time": "2018-12-03T01:24:24.210209Z" } }, "outputs": [ { "data": { "text/plain": [ "8554715.596323118" ] }, "execution_count": 59, "metadata": {}, "output_type": "execute_result" } ], "source": [ "wave_p" ] }, { "cell_type": "code", "execution_count": 21, "metadata": { "ExecuteTime": { "end_time": "2018-12-03T00:54:57.439956Z", "start_time": "2018-12-03T00:54:57.425899Z" } }, "outputs": [], "source": [ "def dune_toe_elevation_change(site_id, df_profile_features):\n", " query = \"site_id == '{}'\".format(site_id)\n", " profile_features = df_profile_features.query(query)\n", " prestorm_dune_toe_z = profile_features.query(\"profile_type=='prestorm'\").dune_toe_z.tolist()[0]\n", " poststorm_dune_toe_z = profile_features.query(\"profile_type=='poststorm'\").dune_toe_z.tolist()[0]\n", " return prestorm_dune_toe_z - poststorm_dune_toe_z\n", "\n", "toe_ele_change = dune_toe_elevation_change(\"NARRA0018\", df_profile_features)" ] }, { "cell_type": "code", "execution_count": 92, "metadata": { "ExecuteTime": { "end_time": "2018-12-03T03:07:48.328428Z", "start_time": "2018-12-03T03:06:54.182343Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0 of 204\n", "20 of 204\n", "40 of 204\n", "60 of 204\n", "80 of 204\n", "100 of 204\n", "120 of 204\n", "140 of 204\n", "160 of 204\n", "180 of 204\n", "200 of 204\n" ] } ], "source": [ "vols = []\n", "exposed_times = []\n", "toe_ele_changes = []\n", "wave_powers = []\n", "\n", "# Get site ids where we observed collision\n", "observed_site_ids = impacts['observed'].query(\"storm_regime=='collision'\").index.get_level_values('site_id').unique()\n", "\n", "# Get site ids where we forecast swash\n", "forecasted_site_ids = impacts['forecasted']['mean_slope_sto06'].query(\"storm_regime=='swash'\").index.get_level_values('site_id').unique()\n", "\n", "site_ids = set(observed_site_ids).intersection(set(forecasted_site_ids))\n", "\n", "# Calculate for each site\n", "\n", "for n, site_id in enumerate(site_ids):\n", " \n", " if n%20 ==0:\n", " print('{} of {}'.format(n, len(site_ids)))\n", " \n", " try:\n", " query = \"site_id == '{}' and profile_type == '{}'\".format(site_id, 'prestorm')\n", " prestorm_profile = df_profiles.query(query)\n", " profile_features = df_profile_features.query(query)\n", "\n", " vol = collision_protection_vol(x = prestorm_profile.index.get_level_values('x'),\n", " z = prestorm_profile.z,\n", " d_low_x = profile_features.dune_toe_x.tolist()[0],\n", " d_low_z = profile_features.dune_toe_z.tolist()[0],\n", " lower_z = profile_features.dune_toe_z.tolist()[0] - 2, # from mhw\n", " angle = 60, # degrees from the horizontal\n", " )\n", "\n", " df_twl = twls['forecasted']['mean_slope_sto06']\n", " df_twl_site = df_twl.query(\"site_id == '{}'\".format(site_id))\n", "\n", " exposed_time = wl_time(t = df_twl_site.index.get_level_values('datetime'),\n", " wl = df_twl_site.R2.tolist(),\n", " z_lower = profile_features.dune_toe_z.tolist()[0] -2,\n", " z_upper = profile_features.dune_toe_z.tolist()[0],\n", " )\n", "\n", " \n", " power = wave_power(t = df_twl_site.index.get_level_values('datetime'),\n", " wl = df_twl_site.R2.tolist(),\n", " z_lower = profile_features.dune_toe_z.tolist()[0] -2,\n", " z_upper = profile_features.dune_toe_z.tolist()[0],\n", " Hs0=df_twl_site.Hs0.tolist(),\n", " Tp=df_twl_site.Tp.tolist())\n", " \n", " toe_ele_change = dune_toe_elevation_change(site_id, df_profile_features)\n", " except:\n", " continue\n", "\n", " vols.append(vol)\n", " exposed_times.append(exposed_time)\n", " toe_ele_changes.append(toe_ele_change)\n", " wave_powers.append(power)\n", "# if n>100:\n", "# break\n", "\n", " \n", "\n" ] }, { "cell_type": "code", "execution_count": 95, "metadata": { "ExecuteTime": { "end_time": "2018-12-03T03:12:11.598150Z", "start_time": "2018-12-03T03:12:11.590128Z" } }, "outputs": [ { "data": { "text/plain": [ "[0.0,\n", " 62.0,\n", " 35.0,\n", " 53.0,\n", " 4.0,\n", " 41.0,\n", " 31.0,\n", " 98.0,\n", " 17.0,\n", " 103.0,\n", " 81.0,\n", " 1.0,\n", " 31.0,\n", " 0.0,\n", " 24.0,\n", " 35.0,\n", " 67.0,\n", " 53.0,\n", " 77.0,\n", " 130.0,\n", " 66.0,\n", " 89.0,\n", " 31.0,\n", " 13.0,\n", " 41.0,\n", " 118.0,\n", " 0.0,\n", " 29.0,\n", " 82.0,\n", " 52.0,\n", " 84.0,\n", " 31.0,\n", " 2.0,\n", " 43.0,\n", " 84.0,\n", " 421.0,\n", " 29.0,\n", " 22.0,\n", " 54.0,\n", " 150.0,\n", " 115.0,\n", " 46.0,\n", " 29.0,\n", " 86.0,\n", " 103.0,\n", " 82.0,\n", " 77.0,\n", " 50.0,\n", " 104.0,\n", " 81.0,\n", " 67.0,\n", " 81.0,\n", " 141.0,\n", " 64.0,\n", " 17.0,\n", " 95.0,\n", " 46.0,\n", " 80.0,\n", " 77.0,\n", " 4.0,\n", " 159.0,\n", " 115.0,\n", " 33.0,\n", " 13.0,\n", " 46.0,\n", " 110.0,\n", " 162.0,\n", " 16.0,\n", " 77.0,\n", " 77.0,\n", " 86.0,\n", " 46.0,\n", " 6.0,\n", " 42.0,\n", " 77.0,\n", " 89.0,\n", " 29.0,\n", " 172.0,\n", " 68.0,\n", " 204.0,\n", " 8.0,\n", " 13.0,\n", " 78.0,\n", " 0.0,\n", " 81.0,\n", " 156.0,\n", " 57.0,\n", " 0.0,\n", " 62.0,\n", " 46.0,\n", " 37.0,\n", " 52.0,\n", " 310.0,\n", " 287.0,\n", " 4.0,\n", " 89.0,\n", " 2.0,\n", " 22.0,\n", " 31.0,\n", " 20.0,\n", " 45.0,\n", " 54.0,\n", " 46.0,\n", " 43.0,\n", " 0.0,\n", " 89.0,\n", " 122.0,\n", " 5.0,\n", " 46.0,\n", " 24.0,\n", " 0.0,\n", " 77.0,\n", " 51.0,\n", " 0.0,\n", " 43.0,\n", " 53.0,\n", " 151.0,\n", " 52.0,\n", " 29.0,\n", " 103.0,\n", " 35.0,\n", " 68.0,\n", " 17.0,\n", " 29.0,\n", " 34.0,\n", " 211.0,\n", " 55.0,\n", " 85.0,\n", " 21.0,\n", " 14.0,\n", " 103.0,\n", " 227.0,\n", " 208.0,\n", " 78.0,\n", " 43.0,\n", " 17.0,\n", " 104.0,\n", " 50.0,\n", " 37.0,\n", " 54.0,\n", " 78.0,\n", " 349.0,\n", " 80.0,\n", " 49.0,\n", " 29.0,\n", " 17.0,\n", " 82.0,\n", " 91.0,\n", " 1.0,\n", " 75.0,\n", " 46.0,\n", " 210.0,\n", " 205.0,\n", " 16.0,\n", " 35.0,\n", " 82.0,\n", " 49.0,\n", " 0.0,\n", " 29.0,\n", " 58.0,\n", " 57.0,\n", " 103.0,\n", " 29.0,\n", " 0.0,\n", " 46.0,\n", " 48.0,\n", " 1.0,\n", " 17.0,\n", " 48.0,\n", " 29.0,\n", " 17.0,\n", " 165.0,\n", " 45.0,\n", " 17.0,\n", " 426.0,\n", " 30.0]" ] }, "execution_count": 95, "metadata": {}, "output_type": "execute_result" } ], "source": [] }, { "cell_type": "code", "execution_count": 104, "metadata": { "ExecuteTime": { "end_time": "2018-12-03T03:31:19.534072Z", "start_time": "2018-12-03T03:31:19.439822Z" } }, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "01b5e4786fd44414a38db9266df25a12", "version_major": 2, "version_minor": 0 }, "text/html": [ "
Failed to display Jupyter Widget of type FigureWidget
.
\n", " If you're reading this message in the Jupyter Notebook or JupyterLab Notebook, it may mean\n", " that the widgets JavaScript is still loading. If this message persists, it\n", " likely means that the widgets JavaScript library is either not installed or\n", " not enabled. See the Jupyter\n", " Widgets Documentation for setup instructions.\n", "
\n", "\n", " If you're reading this message in another frontend (for example, a static\n", " rendering on GitHub or NBViewer),\n", " it may mean that your frontend doesn't currently support widgets.\n", "
\n" ], "text/plain": [ "FigureWidget({\n", " 'data': [{'marker': {'color': [0.6740000000000004, -0.262, 0.5449999999999999,\n", " 0.35399999999999965, 0.3929999999999998,\n", " 0.35399999999999965, 0.6699999999999999,\n", " -0.22799999999999976, 1.0739999999999998,\n", " -0.125, 0.35199999999999987, 0.5299999999999998,\n", " 0.4109999999999996, -0.1120000000000001,\n", " 0.9180000000000001, 0.621, -0.3579999999999992,\n", " -0.20699999999999985, -0.48400000000000043,\n", " 0.5430000000000001, -0.09600000000000053,\n", " -0.41700000000000026, 0.4340000000000006,\n", " 0.7559999999999998, 0.8720000000000003,\n", " -0.5759999999999996, 0.5119999999999996,\n", " 0.6440000000000001, 0.395, 0.09300000000000042,\n", " 0.923, 0.4339999999999997, 0.5429999999999993,\n", " -0.09699999999999998, -0.20899999999999963,\n", " -0.8130000000000002, 0.706, 0.36399999999999944,\n", " 0.34099999999999975, -0.08499999999999996,\n", " 0.21799999999999953, 0.6240000000000006,\n", " -0.16700000000000026, -0.278,\n", " -0.5110000000000001, 0.726,\n", " -0.35199999999999987, 0.13100000000000023,\n", " -0.3440000000000003, 0.6760000000000002,\n", " -0.5719999999999996, -0.4980000000000002,\n", " -0.2020000000000004, 0.34199999999999964,\n", " 0.41999999999999993, 0.5800000000000001,\n", " 0.7200000000000002, 0.6230000000000002,\n", " -0.49100000000000055, 0.11499999999999977,\n", " -0.242, -0.6240000000000006, 0.9500000000000002,\n", " 1.2479999999999998, 1.2889999999999997,\n", " 0.021999999999999797, -0.3200000000000003,\n", " 0.7829999999999999, -0.32899999999999974,\n", " 0.3389999999999995, -0.01599999999999957,\n", " 0.38399999999999945, -0.04599999999999982,\n", " 0.4849999999999999, 0.02499999999999991,\n", " -0.7749999999999999, 0.617,\n", " -0.06500000000000039, 0.5589999999999997,\n", " -0.1299999999999999, -0.0040000000000000036,\n", " 1.298, 0.05500000000000016, 0.9689999999999999,\n", " -0.28800000000000026, -0.09700000000000042,\n", " -0.03100000000000014, 0.7960000000000003,\n", " 0.33599999999999985, 0.28900000000000015, 0.798,\n", " -0.05900000000000016, -0.0029999999999996696,\n", " -0.15399999999999991, 0.31999999999999984,\n", " -0.613, 0.2650000000000001, 0.258, 1.513,\n", " 0.36199999999999966, 0.7849999999999997,\n", " -0.28900000000000015, 1.0290000000000004,\n", " 0.8260000000000005, -0.03799999999999981,\n", " -0.21300000000000008, -0.492,\n", " 1.6029999999999998, 0.5699999999999998,\n", " 0.9300000000000002, 0.0389999999999997,\n", " -0.3110000000000004, 1.0070000000000001,\n", " 0.11500000000000021, 0.7269999999999999,\n", " 0.30100000000000016, -0.5180000000000002,\n", " 0.40600000000000014, 0.1389999999999998,\n", " -0.08299999999999974, 1.174,\n", " -0.29300000000000015, 0.036000000000000476,\n", " 0.702, 0.49099999999999966, 0.27400000000000047,\n", " 0.6890000000000005, 0.46199999999999974,\n", " 0.7839999999999998, 0.2939999999999996,\n", " 0.013000000000000345, -0.09400000000000075,\n", " 0.022999999999999687, -0.6230000000000002,\n", " 0.9040000000000004, 0.027000000000000135,\n", " -0.1280000000000001, 0.25, 0.488,\n", " -0.1200000000000001, 0.43599999999999994,\n", " 0.04999999999999982, 0.11599999999999966,\n", " 1.5290000000000001, 1.0189999999999997,\n", " 0.5780000000000003, -0.4969999999999999,\n", " -0.39800000000000013, 0.6430000000000002,\n", " -0.04500000000000037, -0.40700000000000003,\n", " 0.07000000000000028, -0.18599999999999994,\n", " 1.4900000000000002, 0.2889999999999997,\n", " -0.3719999999999999, 0.1519999999999997,\n", " 0.10899999999999999, -0.3700000000000001,\n", " 0.7320000000000002, -0.238,\n", " -0.08900000000000041, 0.20699999999999985,\n", " 0.9870000000000001, 1.9829999999999997,\n", " 0.6460000000000004, 0.06100000000000039,\n", " 0.06300000000000017, -0.03500000000000014,\n", " 0.5730000000000004, 0.7640000000000002,\n", " -0.7690000000000006, 0.5750000000000002,\n", " 1.4529999999999998, 0.07500000000000018, 0.649],\n", " 'colorscale': 'Viridis',\n", " 'showscale': True,\n", " 'size': 4},\n", " 'mode': 'markers',\n", " 'text': [0.6740000000000004