Merge pull request #69 from ebgoldstein/master

Added Passarella 2018

.isort.cfg

@@ -2,4 +2,4 @@
 line_length = 88
 multi_line_output = 3
 include_trailing_comma = True
-known_third_party = matplotlib,numpy,pandas,pytest,sklearn
+known_third_party = joblib,matplotlib,numpy,pandas,pkg_resources,pytest,sklearn

README.rst

@@ -60,6 +60,8 @@ The following `wave runup models`_ are available for use:
 - ``models.Ruggiero2001``: Based on runup measurements from dissipative Orgeon beaches.
 - ``models.Vousdoukas2012``: Based on runup from European Atlantic coast
 - ``models.Senechal2011``: Based on extreme storm condition at Truc Vert, France
+- ``models.Beuzen2019``: Gaussian Process (GP) runup model
+- ``models.Passarella2018``: Genetic Programming (infragravity and total) swash model 
 
 To get calculate runup, setup and swash, define your offshore conditions in your
 selected runup model then you can access each parameter:

docs/models.rst

@@ -14,4 +14,6 @@ Wave runup models
  Nielsen2009
  Ruggiero2001
  Vousdoukas2012
- Senechal2011
+ Senechal2011
+ Beuzen2019
+ Passarella2018

py_wave_runup/datasets.py

@@ -116,7 +116,6 @@ def load_random_sto06(
  2000-01-01 04:00:00 2.906808 6.937742 0.099069 ... 2.001251 1.476527 1.687905
  <BLANKLINE>
  [5 rows x 8 columns]
-
  """
 
  # Create the time index for our dataframe

py_wave_runup/ensembles.py

@@ -61,10 +61,10 @@ def estimate(self, param):
  >>> ensemble = EnsembleRaw(Hs=4, Tp=11, beta=0.1)
  >>> ensemble_r2 = ensemble.estimate('R2')
  >>> ensemble_r2
- Stockdon2006_R2 Power2018_R2 ... Atkinson2017_R2 Senechal2011_R2
- 0  2.542036 NaN ...  3.1775  1.972371
+ Stockdon2006_R2 Power2018_R2 ... Senechal2011_R2 Beuzen2019_R2
+ 0 2.542036 NaN ... 1.972371 2.181613
  <BLANKLINE>
- [1 rows x 8 columns]
+ [1 rows x 9 columns]
  """
  return self._model_estimates(param)
 
@@ -125,8 +125,8 @@ def estimate(self, param):
  >>> ensemble = EnsembleMean(Hs=[3,4], Tp=[10,11], beta=[0.09,0.1])
  >>> ensemble_r2 = ensemble.estimate('R2')
  >>> ensemble_r2
- 0 2.01381
- 1 2.64924
+ 0 1.969634
+ 1 2.590786
  Name: mean_R2, dtype: float64
  """
  df = self._model_estimates(param)

py_wave_runup/models.py

@@ -4,9 +4,12 @@
 typically based on Hs, Tp, and beta.
 """
 
+import warnings
 from abc import ABCMeta, abstractmethod
 
+import joblib
 import numpy as np
+from pkg_resources import resource_filename
 
 
 class RunupModel(metaclass=ABCMeta):
@@ -73,27 +76,6 @@ def _return_one_or_array(self, val):
  else:
  return val
 
- @property
- @abstractmethod
- def R2(self):
- raise NotImplementedError
-
- @property
- def setup(self):
- raise NotImplementedError
-
- @property
- def sinc(self):
- raise NotImplementedError
-
- @property
- def sig(self):
- raise NotImplementedError
-
- @property
- def swash(self):
- raise NotImplementedError
-
 
 class Stockdon2006(RunupModel):
  """
@@ -557,10 +539,111 @@ def R2(self):
  def sig(self):
  """
  Returns:
- Infragravity componennt of swash:
+ Infragravity component of swash:
 
  .. math: S_{ig} = 0.05 * (H_{s} L_{p})^{0.5}
  """
  result = 0.05 * np.sqrt(self.Hs * self.Lp)
  result = self._return_one_or_array(result)
  return result
+
+
+class Beuzen2019(RunupModel):
+ """
+ Implements the GP runup model from Beuzen et al (2019).
+
+ Beuzen, T., Goldstein, E. B., & Splinter, K. D., 2019. Ensemble models from
+ machine learning: an example of wave runup and coastal dune erosion.
+ https://doi.org/10.5194/nhess-19-2295-2019
+
+ Examples:
+ Calculate 2% exceedence runup level given Hs=4m, Tp=11s, beta=0.1
+
+ >>> from py_wave_runup.models import Beuzen2019
+ >>> beu19 = Beuzen2019(Hs=4, Tp=11, beta=0.1)
+ >>> beu19.R2
+ 2.181613070940485
+ """
+
+ @property
+ def R2(self):
+ """
+ Returns:
+ The 2% exceedence runup level from a pre-trained Gaussian process model
+ """
+ model_path = resource_filename(
+ "py_wave_runup", "datasets/gp_runup_model.joblib"
+ )
+
+ # Ignore the warning when unpickling GaussianProcessRegressor from version
+ # 0.22.1.
+ with warnings.catch_warnings():
+ warnings.simplefilter("ignore")
+ with open(model_path, "rb") as f:
+ model = joblib.load(f)
+
+ result = np.squeeze(
+ model.predict(np.column_stack((self.Hs, self.Tp, self.beta)))
+ )
+ result = self._return_one_or_array(result)
+ return result
+
+
+class Passarella2018(RunupModel):
+ """
+ Implements the Infragravity Swash model from Passarella et al (2018)
+
+ Passarella, M., E. B. Goldstein, S. De Muro, G. Coco, 2018.
+ The use of genetic programming to develop a predictor of swash excursion on sandy beaches.
+ Nat. Hazards Earth Syst. Sci., 18, 599-611,
+ https://doi.org/10.5194/nhess-18-599-2018
+
+ Examples:
+ Calculate IG swash given Hs=4m, Tp=11s, beta=0.1
+
+ >>> from py_wave_runup.models import Passarella2018
+ >>> pas18 = Passarella2018(Hs=4, Tp=11, beta=0.1)
+ >>> pas18.sig
+ 1.5687930560916425
+ """
+
+ @property
+ def sig(self):
+ """
+ Returns:
+ Infragravity component of swash using Eqn (14):
+
+ .. math::
+
+ S_{ig} = \\frac{\\beta}{0.028+\\beta} +
+ \\frac{-1}{2412.255 \\beta - 5.521 \\beta L_{p}} +
+ \\frac{H_{s} -0.711}{0.465 + 173.470 (\\frac{H_{s}}{L_{p}})}
+
+ """
+ result = (
+ (self.beta / (0.028 + self.beta))
+ + (-1 / ((2412.255 * self.beta) - (5.521 * self.beta * self.Lp)))
+ + ((self.Hs - 0.711) / (0.465 + (173.470 * (self.Hs / self.Lp))))
+ )
+ result = self._return_one_or_array(result)
+ return result
+
+ @property
+ def swash(self):
+ """
+ Returns:
+ Total amount of swash using Eqn (12):
+
+ .. math::
+
+ S = 146.737\\beta^{2} + \\frac{T_{p}H_{s}^{3}}{5.800+10.595H_{
+ s}^{3}} - 4397.838\\beta^4
+
+ """
+ result = (
+ (146.737 * (self.beta ** 2))
+ + ((self.Tp * (self.Hs ** 3)) / (5.800 + (10.595 * (self.Hs ** 3))))
+ - 4397.838 * (self.beta ** 4)
+ )
+ result = self._return_one_or_array(result)
+ return result

requirements.txt

@@ -1,3 +1,4 @@
+joblib
 black
 codecov
 coverage

tests/test_models.py

@@ -129,3 +129,25 @@ def test_dissipative(self):
  model = models.Senechal2011(Hs=4, Tp=11, beta=0.001)
  assert model.R2 == approx(1.97, abs=0.01)
  assert model.sig == approx(1.37, abs=0.01)
+
+
+class TestBeuzenl2019(object):
+ def test_reflective(self):
+ model = models.Beuzen2019(Hs=4, Tp=11, beta=0.1)
+ assert model.R2 == approx(2.18, abs=0.01)
+
+ def test_dissipative(self):
+ model = models.Beuzen2019(Hs=4, Tp=11, beta=0.001)
+ assert model.R2 == approx(2.12, abs=0.01)
+
+
+class TestPassarella2018(object):
+ def test_reflective(self):
+ model = models.Passarella2018(Hs=4, Tp=11, beta=0.1)
+ assert model.sig == approx(1.57, abs=0.01)
+ assert model.swash == approx(2.05, abs=0.01)
+
+ def test_dissipative(self):
+ model = models.Passarella2018(Hs=4, Tp=11, beta=0.001)
+ assert model.sig == approx(0.10, abs=0.01)
+ assert model.swash == approx(1.02, abs=0.01)