{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "## Investigate how dune toe compares to R_high" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "ExecuteTime": { "end_time": "2018-12-03T23:04:57.331037Z", "start_time": "2018-12-03T23:04:57.006071Z" } }, "outputs": [], "source": [ "%matplotlib inline\n", "%reload_ext autoreload\n", "%autoreload" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "ExecuteTime": { "end_time": "2018-12-03T23:04:58.749827Z", "start_time": "2018-12-03T23:04:57.333943Z" } }, "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": 3, "metadata": { "ExecuteTime": { "end_time": "2018-12-03T23:05:05.800496Z", "start_time": "2018-12-03T23:04:58.751721Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Importing profiles.csv\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "C:\\Users\\z5189959\\Desktop\\nsw-2016-storm-impact\\.venv\\lib\\site-packages\\numpy\\lib\\arraysetops.py:522: 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": [ "### Compare predicted R_high with D_low\n", "Let's see what the distribution of R_high is compared with D_low. How far off are the predicted water levels compared with the dune toes?" ] }, { "cell_type": "code", "execution_count": 29, "metadata": { "ExecuteTime": { "end_time": "2018-12-04T02:20:58.446500Z", "start_time": "2018-12-04T02:20:58.439480Z" } }, "outputs": [], "source": [ "def get_site_ids(df_forecasted, df_observed, forecasted_regime, observed_regime):\n", " \"\"\"\n", " Returns list of site_ids which match the given forecasted and observed regime\n", " \"\"\"\n", " set1 = set(df_forecasted.query(\"storm_regime == '{}'\".format(\n", " forecasted_regime)).index.get_level_values('site_id'))\n", " set2 = set(df_observed.query(\"storm_regime == '{}'\".format(\n", " observed_regime)).index.get_level_values('site_id'))\n", " return sorted(list(set1.intersection(set2)))\n", "\n", "\n", "def get_R_high_D_low_diff(site_ids, df_profile_features, df_twls):\n", " \"\"\"\n", " Returns a dataframe of the difference between the R_high and D_low differences. \n", " Positive values indicate R_high is larger than D_low.\n", " \"\"\"\n", " # Get dune toes at these sites and predicted max R_high\n", " df_toes = df_profile_features.loc[site_ids].query(\n", " 'profile_type==\"prestorm\"').dune_toe_z\n", " df_R_highs = df_twls.loc[site_ids].groupby('site_id')['R_high'].max()\n", "\n", " # Join into one dataframe\n", " df_twl_toes = pd.concat([df_toes, df_R_highs], axis=1, sort=True)\n", " df_twl_toes['diff'] = df_twl_toes['R_high'] - df_twl_toes['dune_toe_z']\n", " return df_twl_toes['diff']\n" ] }, { "cell_type": "code", "execution_count": 53, "metadata": { "ExecuteTime": { "end_time": "2018-12-04T03:55:51.858020Z", "start_time": "2018-12-04T03:55:50.879155Z" } }, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "94883b85733444528fe8a73379ce4611", "version_major": 2, "version_minor": 0 }, "text/plain": [ "FigureWidget({\n", " 'data': [{'marker': {'color': '#ef8a62'},\n", " 'name': 'Overpredicted',\n", " …" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "swash_overpredicted_site_ids = get_site_ids(df_forecasted=impacts['forecasted']['mean_slope_sto06'],\n", " df_observed=impacts['observed'],\n", " forecasted_regime='collision',\n", " observed_regime='swash')\n", "swash_overpredicted_diffs = get_R_high_D_low_diff(site_ids=swash_overpredicted_site_ids,\n", " df_profile_features=df_profile_features,\n", " df_twls=twls['forecasted']['mean_slope_sto06'])\n", "\n", "swash_correct_site_ids = get_site_ids(df_forecasted=impacts['forecasted']['mean_slope_sto06'],\n", " df_observed=impacts['observed'],\n", " forecasted_regime='swash',\n", " observed_regime='swash')\n", "swash_correct_diffs = get_R_high_D_low_diff(site_ids=swash_correct_site_ids,\n", " df_profile_features=df_profile_features,\n", " df_twls=twls['forecasted']['mean_slope_sto06'])\n", "\n", "\n", "trace1 = go.Histogram(y=swash_correct_diffs.tolist(),\n", " opacity=0.75,\n", " name='Correctly predicted',\n", " marker=dict(\n", " color='#67a9cf',\n", " ),\n", " ybins=dict(\n", " size=0.1\n", "),)\n", "trace2 = go.Histogram(y=swash_overpredicted_diffs.tolist(),\n", " opacity=0.75,\n", " name='Overpredicted',\n", " marker=dict(\n", " color='#ef8a62',\n", "),\n", " ybins=dict(\n", " size=0.1\n", "),)\n", "\n", "layout = go.Layout(\n", " title='R_high - D_low
Swash Regime',\n", " barmode='overlay',\n", " yaxis=dict(\n", " title='z (m AHD)'\n", " ),\n", " xaxis=dict(\n", " title='Count'\n", " ),\n", " bargap=0.2,\n", " bargroupgap=0.1,\n", " legend=dict(x=.6, y=1)\n", ")\n", "\n", "g_plot_swash = go.FigureWidget(data=[trace2, trace1], layout=layout)\n", "\n", "# To output to file\n", "img_bytes = pio.write_image(g_plot_swash, 'g_plot_swash.png',format='png', width=600, height=400, scale=5)\n", "\n", "g_plot_swash\n", "\n" ] }, { "cell_type": "code", "execution_count": 54, "metadata": { "ExecuteTime": { "end_time": "2018-12-04T04:10:47.339268Z", "start_time": "2018-12-04T04:10:45.796887Z" } }, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "3933da9295fe446f9413bca8842100c2", "version_major": 2, "version_minor": 0 }, "text/plain": [ "FigureWidget({\n", " 'data': [{'marker': {'color': '#ef8a62'},\n", " 'name': 'Underpredicted',\n", " …" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "collision_underpredicted_site_ids = get_site_ids(df_forecasted=impacts['forecasted']['mean_slope_sto06'],\n", " df_observed=impacts['observed'],\n", " forecasted_regime='swash',\n", " observed_regime='collision')\n", "collision_underpredicted_diffs = get_R_high_D_low_diff(site_ids=collision_underpredicted_site_ids,\n", " df_profile_features=df_profile_features,\n", " df_twls=twls['forecasted']['mean_slope_sto06'])\n", "\n", "collision_correct_site_ids = get_site_ids(df_forecasted=impacts['forecasted']['mean_slope_sto06'],\n", " df_observed=impacts['observed'],\n", " forecasted_regime='collision',\n", " observed_regime='collision')\n", "collision_correct_diffs = get_R_high_D_low_diff(site_ids=collision_correct_site_ids,\n", " df_profile_features=df_profile_features,\n", " df_twls=twls['forecasted']['mean_slope_sto06'])\n", "\n", "\n", "trace1 = go.Histogram(y=collision_correct_diffs.tolist(),\n", " opacity=0.75,\n", " name='Correctly predicted',\n", " marker=dict(\n", " color='#67a9cf',\n", " ),\n", " ybins=dict(\n", " size=0.1\n", "),)\n", "trace2 = go.Histogram(y=collision_underpredicted_diffs.tolist(),\n", " opacity=0.75,\n", " name='Underpredicted',\n", " marker=dict(\n", " color='#ef8a62',\n", " ),\n", " ybins=dict(\n", " size=0.1\n", "),)\n", "\n", "layout = go.Layout(\n", " title='R_high - D_low
Collision Regime',\n", " barmode='overlay',\n", " yaxis=dict(\n", " title='z (m AHD)'\n", " ),\n", " xaxis=dict(\n", " title='Count'\n", " ),\n", " bargap=0.2,\n", " bargroupgap=0.1,\n", " legend=dict(x=.6, y=1)\n", ")\n", "\n", "g_plot_collision = go.FigureWidget(data=[trace2, trace1], layout=layout)\n", "\n", "# To output to file\n", "img_bytes = pio.write_image(g_plot_collision, 'g_plot_collision.png',format='png', width=600, height=400, scale=5)\n", "\n", "g_plot_collision" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Does dune toe lower?\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": { "heading_collapsed": true }, "source": [ "### What do over predicted and underpredicted profiles look like?" ] }, { "cell_type": "markdown", "metadata": { "hidden": true }, "source": [ "Define a function for getting the average beach profile for a number of given site_ids:" ] }, { "cell_type": "code", "execution_count": 156, "metadata": { "ExecuteTime": { "end_time": "2018-12-04T23:11:08.853877Z", "start_time": "2018-12-04T23:11:08.846876Z" }, "hidden": true }, "outputs": [], "source": [ "def get_avg_profile(site_ids, debug=False):\n", " rows = []\n", " for n,site_id in enumerate(site_ids):\n", " profile = df_profiles.query(\"site_id == '{}' and profile_type == 'prestorm'\".format(site_id))\n", " profile_z = np.array(profile.z.tolist())\n", " profile_x = np.array(profile.index.get_level_values('x').tolist())\n", " \n", " # Let's center the profile based on the z=0 location\n", " idx_last_z_val = max(np.argwhere(~np.isnan(profile_z)==True))[0]\n", " x_last_val = profile_x[idx_last_z_val]\n", " profile_x = [x - x_last_val for x in profile_x]\n", " \n", " # Put values into a dictionary\n", " for x,z in zip(profile_x, profile_z):\n", " rows.append({'x':x, 'z': z})\n", "\n", " # Return early for debugging\n", " if debug and n>3:\n", " break\n", " \n", " # Create dataframe from rows\n", " df = pd.DataFrame(rows)\n", " avg_profile = df.groupby('x').agg({'z': [np.nanmean, np.nanstd]}).reset_index()\n", "\n", " return {\n", " 'x': avg_profile.x.tolist(),\n", " 'z': avg_profile.z.nanmean.tolist(),\n", " 'std': avg_profile.z.nanstd.tolist(),\n", " 'n': n+1 # number of profiles\n", " }" ] }, { "cell_type": "markdown", "metadata": { "hidden": true }, "source": [ "Now, let's look at whether there is a difference between the average beach profile of correctly forecasted site_ids and incorrectly forecasted site_ids. First, looking at sites where we observed swash regime." ] }, { "cell_type": "code", "execution_count": 161, "metadata": { "ExecuteTime": { "end_time": "2018-12-05T02:00:36.853374Z", "start_time": "2018-12-05T01:58:21.839366Z" }, "code_folding": [], "hidden": true }, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "03f2e99d20a347f3922a0e6a36f99ccd", "version_major": 2, "version_minor": 0 }, "text/plain": [ "FigureWidget({\n", " 'data': [{'line': {'color': 'rgb(205, 0, 0)', 'width': 2},\n", " 'mode': 'lines',\n", " …" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "overpredicted = get_avg_profile(swash_overpredicted_site_ids)\n", "correct = get_avg_profile(swash_correct_site_ids)\n", "\n", "# Add mean profile\n", "trace_overpredicted_mean = go.Scatter(\n", " x=overpredicted['x'],\n", " y=overpredicted['z'],\n", " opacity=1,\n", " mode='lines',\n", " name='Mean overpredicted profile (n={})'.format(overpredicted['n']),\n", " line=dict(\n", " color=('rgb(205, 0, 0)'),\n", " width=2)\n", ")\n", "\n", "trace_overpredited_std_top = go.Scatter(\n", " x=overpredicted['x'],\n", " y=np.add(overpredicted['z'], overpredicted['std']),\n", " opacity=1,\n", " hoverinfo='none',\n", " showlegend=False,\n", " mode='lines',\n", " line=dict(\n", " color=('rgb(205, 0, 0)'),\n", " width=0.5,\n", " dash='dash')\n", ")\n", "\n", "trace_overpredited_std_btm = go.Scatter(\n", " x=overpredicted['x'],\n", " y=np.subtract(overpredicted['z'], overpredicted['std']),\n", " opacity=1,\n", " hoverinfo='none',\n", " mode='lines',\n", " showlegend=False,\n", " line=dict(\n", " color=('rgb(205, 0, 0)'),\n", " width=0.5,\n", " dash='dash')\n", ")\n", "\n", "trace_correct_mean = go.Scatter(\n", " x=avg_correct_x,\n", " y=avg_correct_z,\n", " opacity=1,\n", " mode='lines',\n", " name='Mean correct profile (n={})'.format(correct['n']),\n", " line=dict(\n", " color=('rgb(0, 205, 0)'),\n", " width=2)\n", ")\n", "\n", "trace_correct_std_top = go.Scatter(\n", " x=avg_correct_x,\n", " y=np.add(avg_correct_z, avg_correct_std),\n", " opacity=1,\n", " hoverinfo='none',\n", " showlegend=False,\n", " mode='lines',\n", " line=dict(\n", " color=('rgb(0, 205, 0)'),\n", " width=0.5,\n", " dash='dash')\n", ")\n", "\n", "trace_correct_std_btm = go.Scatter(\n", " x=avg_correct_x,\n", " y=np.subtract(avg_correct_z, avg_correct_std),\n", " opacity=1,\n", " hoverinfo='none',\n", " mode='lines',\n", " showlegend=False,\n", " line=dict(\n", " color=('rgb(0, 205, 0)'),\n", " width=0.5,\n", " dash='dash')\n", ")\n", "\n", "layout = dict(showlegend=True,\n", " title='Observed Swash Impact Regime',\n", " legend=dict(x=.6, y=1),\n", " xaxis=dict(\n", " range=[-150, 0]),\n", " yaxis=dict(\n", " range=[0, 10]))\n", "\n", "fig = go.FigureWidget(data=[trace_overpredicted_mean,\n", " trace_overpredited_std_top,\n", " trace_overpredited_std_btm,\n", " trace_correct_mean,\n", " trace_correct_std_top,\n", " trace_correct_std_btm],\n", " layout=layout)\n", "\n", "# To output to file\n", "img_bytes = pio.write_image(\n", " fig, 'mean_profiles_swash.png', format='png', width=600, height=600, scale=5)\n", "\n", "fig" ] }, { "cell_type": "markdown", "metadata": { "hidden": true }, "source": [ "We can see that the difference is pretty minimal. For cases where we predicted collision, but observed swash (overprediction), we see that overpredicted profiles are slightly more concave than correctly predicted sites." ] }, { "cell_type": "code", "execution_count": 162, "metadata": { "ExecuteTime": { "end_time": "2018-12-05T02:03:38.394415Z", "start_time": "2018-12-05T02:00:37.335377Z" }, "hidden": true }, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "1255bccc024e4690b4b8ff4ccc8e9e35", "version_major": 2, "version_minor": 0 }, "text/plain": [ "FigureWidget({\n", " 'data': [{'line': {'color': 'rgb(205, 0, 0)', 'width': 2},\n", " 'mode': 'lines',\n", " …" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "underpredicted = get_avg_profile(collision_underpredicted_site_ids)\n", "correct = get_avg_profile(collision_correct_site_ids)\n", "\n", "# Add mean profile\n", "trace_underpredicted_mean = go.Scatter(\n", " x = underpredicted['x'],\n", " y= underpredicted['z'],\n", " opacity = 1,\n", " mode='lines',\n", " name='Mean underpredicted profile (n={})'.format(underpredicted['n']),\n", " line = dict(\n", " color = ('rgb(205, 0, 0)'),\n", " width = 2)\n", ")\n", "\n", "trace_underpredicted_std_top = go.Scatter(\n", " x = underpredicted['x'],\n", " y= np.add(underpredicted['z'],underpredicted['std']),\n", " opacity = 1,\n", " hoverinfo='none',\n", " showlegend=False,\n", " mode='lines',\n", " line = dict(\n", " color = ('rgb(205, 0, 0)'),\n", " width = 0.5,\n", " dash = 'dash')\n", ") \n", "\n", "trace_underpredicted_std_btm = go.Scatter(\n", " x = underpredicted['x'],\n", " y= np.subtract(underpredicted['z'],underpredicted['std']),\n", " opacity = 1,\n", " hoverinfo='none',\n", " mode='lines',\n", " showlegend=False,\n", " line = dict(\n", " color = ('rgb(205, 0, 0)'),\n", " width = 0.5,\n", " dash = 'dash')\n", ") \n", "\n", "trace_correct_mean = go.Scatter(\n", " x = avg_correct_x,\n", " y= avg_correct_z,\n", " opacity = 1,\n", " mode='lines',\n", " name='Mean correct profile (n={})'.format(correct['n']),\n", " line = dict(\n", " color = ('rgb(0, 205, 0)'),\n", " width = 2)\n", ")\n", "\n", "trace_correct_std_top = go.Scatter(\n", " x = avg_correct_x,\n", " y= np.add(avg_correct_z, avg_correct_std),\n", " opacity = 1,\n", " hoverinfo='none',\n", " showlegend=False,\n", " mode='lines',\n", " line = dict(\n", " color = ('rgb(0, 205, 0)'),\n", " width = 0.5,\n", " dash = 'dash')\n", ") \n", "\n", "trace_correct_std_btm = go.Scatter(\n", " x = avg_correct_x,\n", " y= np.subtract(avg_correct_z, avg_correct_std),\n", " opacity = 1,\n", " hoverinfo='none',\n", " mode='lines',\n", " showlegend=False,\n", " line = dict(\n", " color = ('rgb(0, 205, 0)'),\n", " width = 0.5,\n", " dash = 'dash')\n", ") \n", " \n", "layout = dict(showlegend=True,\n", " title='Observed Collision Impact Regime',\n", " legend=dict(x=.6, y=1),\n", " xaxis=dict(\n", " range=[-150,0]),\n", " yaxis=dict(\n", " range=[0,10]))\n", " \n", "fig=go.FigureWidget(data=[trace_underpredicted_mean, \n", " trace_underpredicted_std_top,\n", " trace_underpredicted_std_btm, \n", " trace_correct_mean, \n", " trace_correct_std_top, \n", " trace_correct_std_btm], \n", " layout=layout)\n", "\n", "# To output to file\n", "img_bytes = pio.write_image(fig, 'mean_profiles_collision.png',format='png', width=600, height=600, scale=5)\n", "\n", "fig\n", "\n" ] }, { "cell_type": "markdown", "metadata": { "hidden": true }, "source": [ "This plot is a bit more interesting. It shows that we are correctly forecasting collision when the profile is more accreted/convex, but when the profile is more eroded/concave, the water level is underpredicted. Why is this? " ] } ], "metadata": { "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.7" }, "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": { "height": "calc(100% - 180px)", "left": "10px", "top": "150px", "width": "232.391px" }, "toc_section_display": true, "toc_window_display": true }, "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 }