chrisl
/
nsw-2016-storm-impact
2
1
Fork 0
Code Issues Pull Requests Releases Activity

Merge branch 'bugfix/fix-mean-slopes' into develop

Browse Source
develop
Chris Leaman 6 years ago
parent 75f372154d 5f59c8f8ee
commit f6c43fda38
10 changed files with 1655 additions and 1474 deletions
Show Stats Download Patch File Download Diff File

2
.gitignore vendored
Unescape Escape View File

@ -1,5 +1,7 @@
# Jupyter NB Checkpoints # Jupyter NB Checkpoints
.ipynb_checkpoints/ .ipynb_checkpoints/
/notebooks/*.png
# exclude data from source control by default # exclude data from source control by default
/data/ /data/

28
Makefile
Unescape Escape View File

@ -11,20 +11,25 @@ CURRENT_DIR = $(shell pwd)
############################### ###############################
# Create python virtual environment # Create python virtual environment
. PHONY: venv_init .PHONY: venv-init
venv-init: ##@environment Setup virtual environment venv-init: ##@environment Setup virtual environment
conda env create -f environment.yml --prefix=.venv python=3.7 conda env create -f environment.yml --prefix=.venv python=3.6
.PHONY: venv-activate
venv-activate: ##@environment Activates the virtual environment venv-activate: ##@environment Activates the virtual environment
activate $(CURRENT_DIR)/.venv activate $(CURRENT_DIR)/.venv
.PHONY: venv-update
venv-update: ##@environment Updates to latest packages venv-update: ##@environment Updates to latest packages
conda update ipykernel && conda update --prefix .venv --all conda update ipykernel && conda update --prefix .venv --all
.PHONY: venv-requirements-install
venv-requirements-install: ##@environment Ensures environment.yml packages are installed venv-requirements-install: ##@environment Ensures environment.yml packages are installed
conda env update conda env update
# The environment.yml file should really be created by hand, but
# this provides a good starting point.
.PHONY: venv-requirements-export
venv-requirements-export: ##@environment Exports current environment to environment.yml venv-requirements-export: ##@environment Exports current environment to environment.yml
conda env export --file environment.yml conda env export --file environment.yml
@ -125,6 +130,17 @@ impacts: ./data/interim/impacts_forecasted_foreshore_slope_sto06.csv ./data/inte
--slope "mean" \ --slope "mean" \
--output-file "./data/interim/twl_mean_slope_sto06.csv" --output-file "./data/interim/twl_mean_slope_sto06.csv"
./data/interim/twl_poststorm_mean_slope_sto06.csv: ./data/interim/waves.csv ./data/interim/tides.csv ./data/interim/profiles.csv ./data/interim/sites.csv ./data/interim/profile_features.csv
activate ./.venv && python ./src/cli.py create-twl-forecast \
--waves-csv "./data/interim/waves.csv" \
--tides-csv "./data/interim/tides.csv" \
--profiles-csv "./data/interim/profiles.csv" \
--profile-features-csv "./data/interim/profile_features.csv" \
--runup-function "sto06" \
--slope "mean" \
--profile-type "poststorm" \
--output-file "./data/interim/twl_poststorm_mean_slope_sto06.csv"
./data/interim/impacts_observed.csv: ./data/interim/profiles.csv ./data/interim/profile_features.csv ./data/interim/impacts_observed.csv: ./data/interim/profiles.csv ./data/interim/profile_features.csv
activate ./.venv && python ./src/cli.py create-observed-impacts \ activate ./.venv && python ./src/cli.py create-observed-impacts \
--profiles-csv "./data/interim/profiles.csv" \ --profiles-csv "./data/interim/profiles.csv" \
@ -143,6 +159,12 @@ impacts: ./data/interim/impacts_forecasted_foreshore_slope_sto06.csv ./data/inte
--forecasted-twl-csv "./data/interim/twl_foreshore_slope_sto06.csv" \ --forecasted-twl-csv "./data/interim/twl_foreshore_slope_sto06.csv" \
--output-file "./data/interim/impacts_forecasted_foreshore_slope_sto06.csv" --output-file "./data/interim/impacts_forecasted_foreshore_slope_sto06.csv"
./data/interim/impacts_forecasted_poststorm_mean_slope_sto06.csv: ./data/interim/profile_features.csv ./data/interim/twl_foreshore_slope_sto06.csv
activate ./.venv && python ./src/cli.py create-forecasted-impacts \
--profile-features-csv "./data/interim/profile_features.csv" \
--forecasted-twl-csv "./data/interim/twl_poststorm_mean_slope_sto06.csv" \
--output-file "./data/interim/impacts_forecasted_poststorm_mean_slope_sto06.csv"
############################### ###############################
# Misc commands # Misc commands

4
README.md
Unescape Escape View File

@ -71,7 +71,7 @@ been corrected for systematic errors, so actual elevations should be taken from
- [ ] Estimate max TWL from elevation where pre storm and post storm profiles are the same. Need to think more about this as runup impacting dune toe will move the dune face back, incorrectly raising the observed twl. Perhaps this estimation of max TWL is only useful for the swash regime. - [ ] Estimate max TWL from elevation where pre storm and post storm profiles are the same. Need to think more about this as runup impacting dune toe will move the dune face back, incorrectly raising the observed twl. Perhaps this estimation of max TWL is only useful for the swash regime.
- [ ] Implement [bayesian change detection algorithm](https://github.com/hildensia/bayesian_changepoint_detection) to help detect dune crests and toes from profiles. Probably low priority at the moment since we are doing manual detection. - [ ] Implement [bayesian change detection algorithm](https://github.com/hildensia/bayesian_changepoint_detection) to help detect dune crests and toes from profiles. Probably low priority at the moment since we are doing manual detection.
- [ ] Implement dune impact calculations as per Palmsten & Holman. Calculation should be done in a new dataframe. - [ ] Implement dune impact calculations as per Palmsten & Holman. Calculation should be done in a new dataframe.
- [ ] Implement data/interim/*.csv file checking using py.test. Check for correct columns, number of nans etc. Testing of code is probably a lower priority than just checking the interim data files at the moment. - [ ] Implement data/interim/*.csv file checking using py.test. Check for correct columns, number of nans etc. Testing of code is probably a lower priority than just checking the interim data files at the moment. Some functions which should be tested are the slope functions in `forecast_twl.py`, as these can be tricky with different profiles.
- [ ] Investigate using [modin](https://github.com/modin-project/modin) to help speed up analysis. - [ ] Investigate using [modin](https://github.com/modin-project/modin) to help speed up analysis.
- [ ] Need to think about how relative imports are handled, see [here](https://chrisyeh96.github.io/2017/08/08/definitive-guide-python-imports.html). Maybe the click CLI interface should be moved to the `./src/` folder and it can import all the other packages? - [X] Need to think about how relative imports are handled, see [here](https://chrisyeh96.github.io/2017/08/08/definitive-guide-python-imports.html). Maybe the click CLI interface should be moved to the `./src/` folder and it can import all the other packages?
- [ ] Simplify runup_models in Stockdon06 - we should really only have one function for each runup model. Need to make it work with individual values or entire dataframe. Use [np.maskedarray](https://docs.scipy.org/doc/numpy-1.15.0/reference/maskedarray.generic.html) - [ ] Simplify runup_models in Stockdon06 - we should really only have one function for each runup model. Need to make it work with individual values or entire dataframe. Use [np.maskedarray](https://docs.scipy.org/doc/numpy-1.15.0/reference/maskedarray.generic.html)

8
environment.yml
Unescape Escape View File

@ -1,5 +1,6 @@
name: C:\Users\z5189959\Desktop\nsw-2016-storm-impact\.venv name: C:\Users\z5189959\Desktop\nsw-2016-storm-impact\.venv
channels: channels:
- plotly
- defaults - defaults
- conda-forge - conda-forge
dependencies: dependencies:
@ -49,6 +50,7 @@ dependencies:
- ipykernel=5.1.0=py36h39e3cac_0 - ipykernel=5.1.0=py36h39e3cac_0
- ipython=7.2.0=py36h39e3cac_0 - ipython=7.2.0=py36h39e3cac_0
- ipython_genutils=0.2.0=py36h3c5d0ee_0 - ipython_genutils=0.2.0=py36h3c5d0ee_0
- ipywidgets=7.4.2=py36_0
- jedi=0.13.1=py36_0 - jedi=0.13.1=py36_0
- jinja2=2.10=py36_0 - jinja2=2.10=py36_0
- jpeg=9b=hb83a4c4_2 - jpeg=9b=hb83a4c4_2
@ -103,6 +105,7 @@ dependencies:
- proj4=4.9.3=hcf24537_7 - proj4=4.9.3=hcf24537_7
- prometheus_client=0.4.2=py36_0 - prometheus_client=0.4.2=py36_0
- prompt_toolkit=2.0.7=py36_0 - prompt_toolkit=2.0.7=py36_0
- psutil=5.4.8=py36he774522_0
- py-boost=1.67.0=py36h8300f20_4 - py-boost=1.67.0=py36h8300f20_4
- pycodestyle=2.4.0=py36_0 - pycodestyle=2.4.0=py36_0
- pycparser=2.19=py36_0 - pycparser=2.19=py36_0
@ -121,6 +124,7 @@ dependencies:
- qt=5.9.6=vc14h1e9a669_2 - qt=5.9.6=vc14h1e9a669_2
- requests=2.20.1=py36_0 - requests=2.20.1=py36_0
- retrying=1.3.3=py36_2 - retrying=1.3.3=py36_2
- scikit-learn=0.20.1=py36hb854c30_0
- scipy=1.1.0=py36h4f6bf74_1 - scipy=1.1.0=py36h4f6bf74_1
- send2trash=1.5.0=py36_0 - send2trash=1.5.0=py36_0
- setuptools=40.6.2=py36_0 - setuptools=40.6.2=py36_0
@ -140,14 +144,18 @@ dependencies:
- wcwidth=0.1.7=py36h3d5aa90_0 - wcwidth=0.1.7=py36h3d5aa90_0
- webencodings=0.5.1=py36_1 - webencodings=0.5.1=py36_1
- wheel=0.32.3=py36_0 - wheel=0.32.3=py36_0
- widgetsnbextension=3.4.2=py36_0
- win_inet_pton=1.0.1=py36_1 - win_inet_pton=1.0.1=py36_1
- wincertstore=0.2=py36h7fe50ca_0 - wincertstore=0.2=py36h7fe50ca_0
- winpty=0.4.3=4 - winpty=0.4.3=4
- xerces-c=3.2.2=ha925a31_0 - xerces-c=3.2.2=ha925a31_0
- xz=5.2.4=h2fa13f4_4 - xz=5.2.4=h2fa13f4_4
- yaml=0.1.7=hc54c509_2 - yaml=0.1.7=hc54c509_2
- yapf=0.25.0=py36_0
- zeromq=4.2.5=he025d50_1 - zeromq=4.2.5=he025d50_1
- zlib=1.2.11=h62dcd97_3 - zlib=1.2.11=h62dcd97_3
- plotly-orca=1.1.1=1
- pip: - pip:
- blackcellmagic==0.0.1
- mat4py==0.4.1 - mat4py==0.4.1
prefix: C:\Users\z5189959\Desktop\nsw-2016-storm-impact\.venv prefix: C:\Users\z5189959\Desktop\nsw-2016-storm-impact\.venv

741
notebooks/01_exploration.ipynb
View File

File diff suppressed because one or more lines are too long

796
notebooks/02_collision_protection_volume.ipynb
View File

File diff suppressed because it is too large Load Diff

762
notebooks/03_dune_to_vs_runup.ipynb
Unescape Escape View File

@ -12,8 +12,8 @@
"execution_count": 1, "execution_count": 1,
"metadata": { "metadata": {
"ExecuteTime": { "ExecuteTime": {
"end_time": "2018-12-03T03:38:44.538853Z", "end_time": "2018-12-03T23:04:57.331037Z",
"start_time": "2018-12-03T03:38:44.189514Z" "start_time": "2018-12-03T23:04:57.006071Z"
} }
}, },
"outputs": [], "outputs": [],
@ -28,8 +28,8 @@
"execution_count": 2, "execution_count": 2,
"metadata": { "metadata": {
"ExecuteTime": { "ExecuteTime": {
"end_time": "2018-12-03T03:38:46.213387Z", "end_time": "2018-12-03T23:04:58.749827Z",
"start_time": "2018-12-03T03:38:44.781382Z" "start_time": "2018-12-03T23:04:57.333943Z"
} }
}, },
"outputs": [], "outputs": [],
@ -61,8 +61,8 @@
"execution_count": 3, "execution_count": 3,
"metadata": { "metadata": {
"ExecuteTime": { "ExecuteTime": {
"end_time": "2018-12-03T03:38:53.297184Z", "end_time": "2018-12-03T23:05:05.800496Z",
"start_time": "2018-12-03T03:38:46.365829Z" "start_time": "2018-12-03T23:04:58.751721Z"
} }
}, },
"outputs": [ "outputs": [
@ -77,7 +77,7 @@
"name": "stderr", "name": "stderr",
"output_type": "stream", "output_type": "stream",
"text": [ "text": [
"C:\\Users\\z5189959\\AppData\\Local\\Continuum\\anaconda3\\lib\\site-packages\\numpy\\lib\\arraysetops.py:472: FutureWarning:\n", "C:\\Users\\z5189959\\Desktop\\nsw-2016-storm-impact\\.venv\\lib\\site-packages\\numpy\\lib\\arraysetops.py:522: FutureWarning:\n",
"\n", "\n",
"elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison\n", "elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison\n",
"\n" "\n"
@ -127,162 +127,158 @@
"cell_type": "markdown", "cell_type": "markdown",
"metadata": {}, "metadata": {},
"source": [ "source": [
"### Compare underpredicted cases" "### 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", "cell_type": "code",
"execution_count": 39, "execution_count": 29,
"metadata": { "metadata": {
"ExecuteTime": { "ExecuteTime": {
"end_time": "2018-12-03T04:05:30.984007Z", "end_time": "2018-12-04T02:20:58.446500Z",
"start_time": "2018-12-03T04:05:30.805508Z" "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": [ "outputs": [
{ {
"data": { "data": {
"text/html": [ "application/vnd.jupyter.widget-view+json": {
"<div>\n", "model_id": "94883b85733444528fe8a73379ce4611",
"<style scoped>\n", "version_major": 2,
" .dataframe tbody tr th:only-of-type {\n", "version_minor": 0
" vertical-align: middle;\n", },
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>dune_toe_z</th>\n",
" <th>R_high</th>\n",
" <th>diff</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>AVOCAn0005</th>\n",
" <td>3.306</td>\n",
" <td>3.260440</td>\n",
" <td>-0.045560</td>\n",
" </tr>\n",
" <tr>\n",
" <th>AVOCAn0008</th>\n",
" <td>3.507</td>\n",
" <td>3.220084</td>\n",
" <td>-0.286916</td>\n",
" </tr>\n",
" <tr>\n",
" <th>BILG0005</th>\n",
" <td>4.807</td>\n",
" <td>3.293445</td>\n",
" <td>-1.513555</td>\n",
" </tr>\n",
" <tr>\n",
" <th>BLUEYS0001</th>\n",
" <td>3.064</td>\n",
" <td>2.800144</td>\n",
" <td>-0.263856</td>\n",
" </tr>\n",
" <tr>\n",
" <th>BLUEYS0002</th>\n",
" <td>2.929</td>\n",
" <td>2.470641</td>\n",
" <td>-0.458359</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [ "text/plain": [
" dune_toe_z R_high diff\n", "FigureWidget({\n",
"AVOCAn0005 3.306 3.260440 -0.045560\n", " 'data': [{'marker': {'color': '#ef8a62'},\n",
"AVOCAn0008 3.507 3.220084 -0.286916\n", " 'name': 'Overpredicted',\n",
"BILG0005 4.807 3.293445 -1.513555\n", " …"
"BLUEYS0001 3.064 2.800144 -0.263856\n",
"BLUEYS0002 2.929 2.470641 -0.458359"
] ]
}, },
"execution_count": 39,
"metadata": {}, "metadata": {},
"output_type": "execute_result" "output_type": "display_data"
} }
], ],
"source": [ "source": [
"# Find site_ids where the forecast has been underpredicted\n", "swash_overpredicted_site_ids = get_site_ids(df_forecasted=impacts['forecasted']['mean_slope_sto06'],\n",
"set1 = set(impacts['forecasted']['mean_slope_sto06'].query(\"storm_regime == 'swash'\").index.get_level_values('site_id'))\n", " df_observed=impacts['observed'],\n",
"set2 = set(impacts['observed'].query(\"storm_regime == 'collision'\").index.get_level_values('site_id'))\n", " forecasted_regime='collision',\n",
"site_ids = list(set1.intersection(set2))\n", " observed_regime='swash')\n",
"\n", "swash_overpredicted_diffs = get_R_high_D_low_diff(site_ids=swash_overpredicted_site_ids,\n",
"# Get dune toes at these sites and predicted max R_high\n", " df_profile_features=df_profile_features,\n",
"df_toes = df_profile_features.loc[site_ids].query('profile_type==\"prestorm\"').dune_toe_z\n", " df_twls=twls['forecasted']['mean_slope_sto06'])\n",
"df_R_highs = twls['forecasted']['mean_slope_sto06'].loc[site_ids].groupby('site_id')['R_high'].max()\n", "\n",
"\n", "swash_correct_site_ids = get_site_ids(df_forecasted=impacts['forecasted']['mean_slope_sto06'],\n",
"# Join into one dataframe\n", " df_observed=impacts['observed'],\n",
"df_twl_toes = pd.concat([df_toes, df_R_highs],axis=1,sort=True)\n", " forecasted_regime='swash',\n",
"df_twl_toes['diff'] = df_twl_toes['R_high'] - df_twl_toes['dune_toe_z']\n", " observed_regime='swash')\n",
"df_twl_toes.head()\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",
"cell_type": "markdown", "\n",
"metadata": {}, "trace1 = go.Histogram(y=swash_correct_diffs.tolist(),\n",
"source": [ " opacity=0.75,\n",
"Now let's plot the comparison between our R_high TWL values and the dune toes to see how far off they were." " 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<br>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", "cell_type": "code",
"execution_count": 41, "execution_count": 54,
"metadata": { "metadata": {
"ExecuteTime": { "ExecuteTime": {
"end_time": "2018-12-03T04:08:15.732169Z", "end_time": "2018-12-04T04:10:47.339268Z",
"start_time": "2018-12-03T04:08:15.656966Z" "start_time": "2018-12-04T04:10:45.796887Z"
} }
}, },
"outputs": [ "outputs": [
{ {
"data": { "data": {
"application/vnd.jupyter.widget-view+json": { "application/vnd.jupyter.widget-view+json": {
"model_id": "35b9331242af473dba2f91761c307022", "model_id": "3933da9295fe446f9413bca8842100c2",
"version_major": 2, "version_major": 2,
"version_minor": 0 "version_minor": 0
}, },
"text/html": [
"<p>Failed to display Jupyter Widget of type <code>FigureWidget</code>.</p>\n",
"<p>\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 <a href=\"https://ipywidgets.readthedocs.io/en/stable/user_install.html\">Jupyter\n",
" Widgets Documentation</a> for setup instructions.\n",
"</p>\n",
"<p>\n",
" If you're reading this message in another frontend (for example, a static\n",
" rendering on GitHub or <a href=\"https://nbviewer.jupyter.org/\">NBViewer</a>),\n",
" it may mean that your frontend doesn't currently support widgets.\n",
"</p>\n"
],
"text/plain": [ "text/plain": [
"FigureWidget({\n", "FigureWidget({\n",
" 'data': [{'type': 'histogram',\n", " 'data': [{'marker': {'color': '#ef8a62'},\n",
" 'uid': '75f0d11f-9242-4fc7-b433-1f04e1e37ba6',\n", " 'name': 'Underpredicted',\n",
" 'y': [-0.045560088746212646, -0.28691603912686325,\n", " …"
" -1.5135547360075963, ..., -0.4613631587476821,\n",
" -0.5212332930925054, -0.3948507473332721]}],\n",
" 'layout': {'bargap': 0.2,\n",
" 'bargroupgap': 0.1,\n",
" 'title': 'D_low - R_high<br>Observed Collision, Forecasted Swash',\n",
" 'xaxis': {'title': 'Count'},\n",
" 'yaxis': {'title': 'z (m AHD)'}}\n",
"})"
] ]
}, },
"metadata": {}, "metadata": {},
@ -290,10 +286,45 @@
} }
], ],
"source": [ "source": [
"trace1 = go.Histogram(y=df_twl_toes['diff'].tolist())\n", "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", "\n",
"layout = go.Layout(\n", "layout = go.Layout(\n",
" title='D_low - R_high<br>Observed Collision, Forecasted Swash',\n", " title='R_high - D_low<br>Collision Regime',\n",
" barmode='overlay',\n",
" yaxis=dict(\n", " yaxis=dict(\n",
" title='z (m AHD)'\n", " title='z (m AHD)'\n",
" ),\n", " ),\n",
@ -301,172 +332,267 @@
" title='Count'\n", " title='Count'\n",
" ),\n", " ),\n",
" bargap=0.2,\n", " bargap=0.2,\n",
" bargroupgap=0.1\n", " bargroupgap=0.1,\n",
" legend=dict(x=.6, y=1)\n",
")\n", ")\n",
"\n", "\n",
"g_plot = go.FigureWidget(data=[trace1], layout=layout)\n", "g_plot_collision = go.FigureWidget(data=[trace2, trace1], layout=layout)\n",
"g_plot" "\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", "cell_type": "markdown",
"metadata": {}, "metadata": {},
"source": [ "source": [
"The above plot shows that the R_high value for most of the incorrectly forecasted collision regimes, was typically underpredicted by less than 0.5 m." "### Does dune toe lower?\n"
] ]
}, },
{ {
"cell_type": "markdown", "cell_type": "code",
"execution_count": null,
"metadata": {}, "metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {
"heading_collapsed": true
},
"source": [ "source": [
"### Compare overpredicted cases" "### 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", "cell_type": "code",
"execution_count": 42, "execution_count": 156,
"metadata": { "metadata": {
"ExecuteTime": { "ExecuteTime": {
"end_time": "2018-12-03T04:08:56.128806Z", "end_time": "2018-12-04T23:11:08.853877Z",
"start_time": "2018-12-03T04:08:55.894182Z" "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": [ "outputs": [
{ {
"data": { "data": {
"text/html": [ "application/vnd.jupyter.widget-view+json": {
"<div>\n", "model_id": "03f2e99d20a347f3922a0e6a36f99ccd",
"<style scoped>\n", "version_major": 2,
" .dataframe tbody tr th:only-of-type {\n", "version_minor": 0
" vertical-align: middle;\n", },
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>dune_toe_z</th>\n",
" <th>R_high</th>\n",
" <th>diff</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>AVOCAn0004</th>\n",
" <td>3.178</td>\n",
" <td>3.416988</td>\n",
" <td>0.238988</td>\n",
" </tr>\n",
" <tr>\n",
" <th>BOOM0004</th>\n",
" <td>3.065</td>\n",
" <td>3.074980</td>\n",
" <td>0.009980</td>\n",
" </tr>\n",
" <tr>\n",
" <th>BOOM0011</th>\n",
" <td>2.771</td>\n",
" <td>6.491824</td>\n",
" <td>3.720824</td>\n",
" </tr>\n",
" <tr>\n",
" <th>BOOM0012</th>\n",
" <td>2.796</td>\n",
" <td>3.148087</td>\n",
" <td>0.352087</td>\n",
" </tr>\n",
" <tr>\n",
" <th>CATHIE0001</th>\n",
" <td>2.780</td>\n",
" <td>3.522792</td>\n",
" <td>0.742792</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [ "text/plain": [
" dune_toe_z R_high diff\n", "FigureWidget({\n",
"AVOCAn0004 3.178 3.416988 0.238988\n", " 'data': [{'line': {'color': 'rgb(205, 0, 0)', 'width': 2},\n",
"BOOM0004 3.065 3.074980 0.009980\n", " 'mode': 'lines',\n",
"BOOM0011 2.771 6.491824 3.720824\n", " …"
"BOOM0012 2.796 3.148087 0.352087\n",
"CATHIE0001 2.780 3.522792 0.742792"
] ]
}, },
"execution_count": 42,
"metadata": {}, "metadata": {},
"output_type": "execute_result" "output_type": "display_data"
} }
], ],
"source": [ "source": [
"# Find site_ids where the forecast has been overpredicted\n", "overpredicted = get_avg_profile(swash_overpredicted_site_ids)\n",
"set1 = set(impacts['forecasted']['mean_slope_sto06'].query(\"storm_regime == 'collision'\").index.get_level_values('site_id'))\n", "correct = get_avg_profile(swash_correct_site_ids)\n",
"set2 = set(impacts['observed'].query(\"storm_regime == 'swash'\").index.get_level_values('site_id'))\n", "\n",
"site_ids = list(set1.intersection(set2))\n", "# Add mean profile\n",
"\n", "trace_overpredicted_mean = go.Scatter(\n",
"# Get dune toes at these sites and predicted max R_high\n", " x=overpredicted['x'],\n",
"df_toes = df_profile_features.loc[site_ids].query('profile_type==\"prestorm\"').dune_toe_z\n", " y=overpredicted['z'],\n",
"df_R_highs = twls['forecasted']['mean_slope_sto06'].loc[site_ids].groupby('site_id')['R_high'].max()\n", " opacity=1,\n",
"\n", " mode='lines',\n",
"# Join into one dataframe\n", " name='Mean overpredicted profile (n={})'.format(overpredicted['n']),\n",
"df_twl_toes = pd.concat([df_toes, df_R_highs],axis=1,sort=True)\n", " line=dict(\n",
"df_twl_toes['diff'] = df_twl_toes['R_high'] - df_twl_toes['dune_toe_z']\n", " color=('rgb(205, 0, 0)'),\n",
"df_twl_toes.head()\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", "cell_type": "code",
"execution_count": 47, "execution_count": 162,
"metadata": { "metadata": {
"ExecuteTime": { "ExecuteTime": {
"end_time": "2018-12-03T04:14:46.601092Z", "end_time": "2018-12-05T02:03:38.394415Z",
"start_time": "2018-12-03T04:14:46.522883Z" "start_time": "2018-12-05T02:00:37.335377Z"
} },
"hidden": true
}, },
"outputs": [ "outputs": [
{ {
"data": { "data": {
"application/vnd.jupyter.widget-view+json": { "application/vnd.jupyter.widget-view+json": {
"model_id": "3ea49a4ac07c4ea19bbb4532326ff94c", "model_id": "1255bccc024e4690b4b8ff4ccc8e9e35",
"version_major": 2, "version_major": 2,
"version_minor": 0 "version_minor": 0
}, },
"text/html": [
"<p>Failed to display Jupyter Widget of type <code>FigureWidget</code>.</p>\n",
"<p>\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 <a href=\"https://ipywidgets.readthedocs.io/en/stable/user_install.html\">Jupyter\n",
" Widgets Documentation</a> for setup instructions.\n",
"</p>\n",
"<p>\n",
" If you're reading this message in another frontend (for example, a static\n",
" rendering on GitHub or <a href=\"https://nbviewer.jupyter.org/\">NBViewer</a>),\n",
" it may mean that your frontend doesn't currently support widgets.\n",
"</p>\n"
],
"text/plain": [ "text/plain": [
"FigureWidget({\n", "FigureWidget({\n",
" 'data': [{'type': 'histogram',\n", " 'data': [{'line': {'color': 'rgb(205, 0, 0)', 'width': 2},\n",
" 'uid': '4a284474-2be1-4fd7-87d5-25364cc78df4',\n", " 'mode': 'lines',\n",
" 'y': [0.23898814460475037, 0.009980312001434566, 3.720823710344608,\n", " …"
" ..., 1.5720238663972683, 0.912998680585452, 1.1419977620500927]}],\n",
" 'layout': {'bargap': 0.2,\n",
" 'bargroupgap': 0.1,\n",
" 'title': 'D_low - R_high<br>Observed Swash, Forecasted Collision',\n",
" 'xaxis': {'title': 'Count'},\n",
" 'yaxis': {'title': 'z (m AHD)'}}\n",
"})"
] ]
}, },
"metadata": {}, "metadata": {},
@ -474,29 +600,114 @@
} }
], ],
"source": [ "source": [
"trace1 = go.Histogram(y=df_twl_toes['diff'].tolist())\n", "underpredicted = get_avg_profile(collision_underpredicted_site_ids)\n",
"correct = get_avg_profile(collision_correct_site_ids)\n",
"\n", "\n",
"layout = go.Layout(\n", "# Add mean profile\n",
" title='D_low - R_high<br>Observed Swash, Forecasted Collision',\n", "trace_underpredicted_mean = go.Scatter(\n",
" yaxis=dict(\n", " x = underpredicted['x'],\n",
" title='z (m AHD)'\n", " y= underpredicted['z'],\n",
" ),\n", " opacity = 1,\n",
" xaxis=dict(\n", " mode='lines',\n",
" title='Count'\n", " name='Mean underpredicted profile (n={})'.format(underpredicted['n']),\n",
" ),\n", " line = dict(\n",
" bargap=0.2,\n", " color = ('rgb(205, 0, 0)'),\n",
" bargroupgap=0.1\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",
"\n", "\n",
"g_plot = go.FigureWidget(data=[trace1], layout=layout)\n", "trace_correct_std_top = go.Scatter(\n",
"g_plot" " 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", "cell_type": "markdown",
"metadata": {}, "metadata": {
"hidden": true
},
"source": [ "source": [
"The errors when we forecast collision but observe swash are much greater than we we forecast swash and observe collision. For this case, errors in excess of 1.0 m common. Why is this?" "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? "
] ]
} }
], ],
@ -506,6 +717,18 @@
"language": "python", "language": "python",
"name": "python3" "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": { "toc": {
"base_numbering": 1, "base_numbering": 1,
"nav_menu": {}, "nav_menu": {},
@ -515,9 +738,14 @@
"title_cell": "Table of Contents", "title_cell": "Table of Contents",
"title_sidebar": "Contents", "title_sidebar": "Contents",
"toc_cell": false, "toc_cell": false,
"toc_position": {}, "toc_position": {
"height": "calc(100% - 180px)",
"left": "10px",
"top": "150px",
"width": "232.391px"
},
"toc_section_display": true, "toc_section_display": true,
"toc_window_display": false "toc_window_display": true
}, },
"varInspector": { "varInspector": {
"cols": { "cols": {

743
notebooks/04_profile_picker.ipynb
View File

File diff suppressed because one or more lines are too long

38
src/analysis/forecast_twl.py
Unescape Escape View File

@ -23,6 +23,7 @@ def forecast_twl(
runup_function, runup_function,
n_processes=MULTIPROCESS_THREADS, n_processes=MULTIPROCESS_THREADS,
slope="foreshore", slope="foreshore",
profile_type='prestorm'
): ):
# Use df_waves as a base # Use df_waves as a base
df_twl = df_waves.copy() df_twl = df_waves.copy()
@ -45,12 +46,14 @@ def forecast_twl(
df_twl["beta"] = pd.concat(results) df_twl["beta"] = pd.concat(results)
elif slope == "mean": elif slope == "mean":
df_temp = df_twl.join(df_profile_features.query("profile_type=='prestorm'").reset_index(level='profile_type') df_temp = df_twl.join(df_profile_features.query("profile_type=='{}'".format(profile_type)).reset_index(
level='profile_type')
, how="inner") , how="inner")
df_temp["mhw"] = 0.5 df_temp["mhw"] = 0.5
with Pool(processes=n_processes) as pool: with Pool(processes=n_processes) as pool:
results = pool.starmap( results = pool.starmap(
mean_slope_for_site_id, [(site_id, df_temp, df_profiles, "dune_toe_z", "mhw") for site_id in site_ids] mean_slope_for_site_id, [(site_id, df_temp, df_profiles, "dune_toe_z", "dune_toe_x", "mhw") for
site_id in site_ids]
) )
df_twl["beta"] = pd.concat(results) df_twl["beta"] = pd.concat(results)
@ -71,7 +74,8 @@ def forecast_twl(
return df_twl return df_twl
def mean_slope_for_site_id(site_id, df_twl, df_profiles, top_elevation_col, btm_elevation_col): def mean_slope_for_site_id(site_id, df_twl, df_profiles, top_elevation_col, top_x_col, btm_elevation_col,
profile_type='prestorm'):
""" """
Calculates the foreshore slope values a given site_id. Returns a series (with same indicies as df_twl) of Calculates the foreshore slope values a given site_id. Returns a series (with same indicies as df_twl) of
foreshore slopes. This function is used to parallelize getting foreshore slopes as it is computationally foreshore slopes. This function is used to parallelize getting foreshore slopes as it is computationally
@ -83,7 +87,7 @@ def mean_slope_for_site_id(site_id, df_twl, df_profiles, top_elevation_col, btm_
""" """
# Get the prestorm beach profile # Get the prestorm beach profile
profile = df_profiles.query("site_id =='{}' and profile_type == 'prestorm'".format(site_id)) profile = df_profiles.query("site_id =='{}' and profile_type == '{}'".format(site_id, profile_type))
profile_x = profile.index.get_level_values("x").tolist() profile_x = profile.index.get_level_values("x").tolist()
profile_z = profile.z.tolist() profile_z = profile.z.tolist()
@ -96,6 +100,7 @@ def mean_slope_for_site_id(site_id, df_twl, df_profiles, top_elevation_col, btm_
top_elevation=row[top_elevation_col], top_elevation=row[top_elevation_col],
btm_elevation=row[btm_elevation_col], btm_elevation=row[btm_elevation_col],
method="end_points", method="end_points",
top_x= row[top_x_col]
), ),
axis=1, axis=1,
) )
@ -191,7 +196,7 @@ def foreshore_slope_from_profile(profile_x, profile_z, tide, runup_function, **k
iteration_count += 1 iteration_count += 1
def slope_from_profile(profile_x, profile_z, top_elevation, btm_elevation, method="end_points"): def slope_from_profile(profile_x, profile_z, top_elevation, btm_elevation, method="end_points", top_x=None, btm_x=None):
""" """
Returns a slope (beta) from a bed profile, given the top and bottom elevations of where the slope should be taken. Returns a slope (beta) from a bed profile, given the top and bottom elevations of where the slope should be taken.
:param x: List of x bed profile coordinates :param x: List of x bed profile coordinates
@ -199,6 +204,9 @@ def slope_from_profile(profile_x, profile_z, top_elevation, btm_elevation, metho
:param top_elevation: Top elevation of where to take the slope :param top_elevation: Top elevation of where to take the slope
:param btm_elevation: Bottom elevation of where to take the slope :param btm_elevation: Bottom elevation of where to take the slope
:param method: Method used to calculate slope (end_points or least_squares) :param method: Method used to calculate slope (end_points or least_squares)
:param top_x: x-coordinate of the top end point. May be needed, as there may be multiple crossings of the
top_elevation.
:param btm_x: x-coordinate of the bottom end point
:return: :return:
""" """
@ -208,7 +216,18 @@ def slope_from_profile(profile_x, profile_z, top_elevation, btm_elevation, metho
end_points = {"top": {"z": top_elevation}, "btm": {"z": btm_elevation}} end_points = {"top": {"z": top_elevation}, "btm": {"z": btm_elevation}}
for end_type in end_points.keys(): for end_type in end_points.keys():
# Add x coordinates if they are specified
if top_x and end_type == 'top':
end_points['top']['x'] = top_x
continue
if btm_x and end_type == 'top':
end_points['btm']['x'] = btm_x
continue
elevation = end_points[end_type]["z"] elevation = end_points[end_type]["z"]
intersection_x = crossings(profile_x, profile_z, elevation) intersection_x = crossings(profile_x, profile_z, elevation)
@ -285,8 +304,10 @@ def crossings(profile_x, profile_z, constant_z):
@click.option("--profile-features-csv", required=True, help="") @click.option("--profile-features-csv", required=True, help="")
@click.option("--runup-function", required=True, help="", type=click.Choice(["sto06"])) @click.option("--runup-function", required=True, help="", type=click.Choice(["sto06"]))
@click.option("--slope", required=True, help="", type=click.Choice(["foreshore", "mean"])) @click.option("--slope", required=True, help="", type=click.Choice(["foreshore", "mean"]))
@click.option("--profile-type", required=True, help="", type=click.Choice(["prestorm", "poststorm"]))
@click.option("--output-file", required=True, help="") @click.option("--output-file", required=True, help="")
def create_twl_forecast(waves_csv, tides_csv, profiles_csv, profile_features_csv, runup_function, slope, output_file): def create_twl_forecast(waves_csv, tides_csv, profiles_csv, profile_features_csv, runup_function, slope,
profile_type,output_file):
logger.info("Creating forecast of total water levels") logger.info("Creating forecast of total water levels")
logger.info("Importing data") logger.info("Importing data")
df_waves = pd.read_csv(waves_csv, index_col=[0, 1]) df_waves = pd.read_csv(waves_csv, index_col=[0, 1])
@ -295,15 +316,16 @@ def create_twl_forecast(waves_csv, tides_csv, profiles_csv, profile_features_csv
df_profile_features = pd.read_csv(profile_features_csv, index_col=[0,1]) df_profile_features = pd.read_csv(profile_features_csv, index_col=[0,1])
logger.info("Forecasting TWL") logger.info("Forecasting TWL")
df_twl_foreshore_slope_sto06 = forecast_twl( df_twl = forecast_twl(
df_tides, df_tides,
df_profiles, df_profiles,
df_waves, df_waves,
df_profile_features, df_profile_features,
runup_function=getattr(runup_models, runup_function), runup_function=getattr(runup_models, runup_function),
slope=slope, slope=slope,
profile_type=profile_type
) )
df_twl_foreshore_slope_sto06.to_csv(output_file) df_twl.to_csv(output_file)
logger.info("Saved to %s", output_file) logger.info("Saved to %s", output_file)
logger.info("Done!") logger.info("Done!")

3
src/cli.py
Unescape Escape View File

@ -11,6 +11,9 @@ import analysis.forecast_twl as forecast_twl
import analysis.forecasted_storm_impacts as forecasted_storm_impacts import analysis.forecasted_storm_impacts as forecasted_storm_impacts
import analysis.observed_storm_impacts as observed_storm_impacts import analysis.observed_storm_impacts as observed_storm_impacts
# Disable numpy warnings
import warnings
warnings.simplefilter(action='ignore', category=FutureWarning)
@click.group() @click.group()
def cli(): def cli():

Write Preview
Loading…
Cancel
Save