Add rainfall intensity calculation function

physics/rainfall_intensity.py

@@ -0,0 +1,143 @@
+"""
+Rainfall Intensity
+==================
+This module contains functions to calculate the intensity of
+a rainfall event for a given duration and return period.
+
+This function uses the Sherman intensity-duration-frequency curve.
+
+References
+----------
+- Aparicio, F. (1997): Fundamentos de Hidrología de Superficie.
+    Balderas, México, Limusa. 303 p.
+- https://en.wikipedia.org/wiki/Intensity-duration-frequency_curve
+"""
+
+
+def rainfall_intensity(
+    coefficient_k: float,
+    coefficient_a: float,
+    coefficient_b: float,
+    coefficient_c: float,
+    return_period: float,
+    duration: float,
+) -> float:
+    """
+    Calculate the intensity of a rainfall event for a given duration and return period.
+    It's based on the Sherman intensity-duration-frequency curve:
+
+    I = K * T^a / (D + b)^c
+
+    where:
+        I = Intensity of the rainfall event [mm/h]
+        k, a, b, c = Coefficients obtained through statistical distribution adjust
+        T = Return period in years
+        D = Rainfall event duration in minutes
+
+    Parameters
+    ----------
+    coefficient_k : float
+        Coefficient obtained through statistical distribution adjust.
+    coefficient_a : float
+        Coefficient obtained through statistical distribution adjust.
+    coefficient_b : float
+        Coefficient obtained through statistical distribution adjust.
+    coefficient_c : float
+        Coefficient obtained through statistical distribution adjust.
+    return_period : float
+        Return period in years.
+    duration : float
+        Rainfall event duration in minutes.
+
+    Returns
+    -------
+    intensity : float
+        Intensity of the rainfall event in mm/h.
+
+    Raises
+    ------
+    ValueError
+        If any of the parameters are not positive.
+
+    Examples
+    --------
+
+    >>> rainfall_intensity(1000, 0.2, 11.6, 0.81, 10, 60)
+    49.83339231138578
+
+    >>> rainfall_intensity(1000, 0.2, 11.6, 0.81, 10, 30)
+    77.36319588106228
+
+    >>> rainfall_intensity(1000, 0.2, 11.6, 0.81, 5, 60)
+    43.382487747633625
+
+    >>> rainfall_intensity(0, 0.2, 11.6, 0.81, 10, 60)
+    Traceback (most recent call last):
+    ...
+    ValueError: All parameters must be positive.
+
+    >>> rainfall_intensity(1000, -0.2, 11.6, 0.81, 10, 60)
+    Traceback (most recent call last):
+    ...
+    ValueError: All parameters must be positive.
+
+    >>> rainfall_intensity(1000, 0.2, -11.6, 0.81, 10, 60)
+    Traceback (most recent call last):
+    ...
+    ValueError: All parameters must be positive.
+
+    >>> rainfall_intensity(1000, 0.2, 11.6, -0.81, 10, 60)
+    Traceback (most recent call last):
+    ...
+    ValueError: All parameters must be positive.
+
+    >>> rainfall_intensity(1000, 0, 11.6, 0.81, 10, 60)
+    Traceback (most recent call last):
+    ...
+    ValueError: All parameters must be positive.
+
+    >>> rainfall_intensity(1000, 0.2, 0, 0.81, 10, 60)
+    Traceback (most recent call last):
+    ...
+    ValueError: All parameters must be positive.
+
+    >>> rainfall_intensity(1000, 0.2, 11.6, 0, 10, 60)
+    Traceback (most recent call last):
+    ...
+    ValueError: All parameters must be positive.
+
+    >>> rainfall_intensity(0, 0.2, 11.6, 0.81, 10, 60)
+    Traceback (most recent call last):
+    ...
+    ValueError: All parameters must be positive.
+
+    >>> rainfall_intensity(1000, 0.2, 11.6, 0.81, 0, 60)
+    Traceback (most recent call last):
+    ...
+    ValueError: All parameters must be positive.
+
+    >>> rainfall_intensity(1000, 0.2, 11.6, 0.81, 10, 0)
+    Traceback (most recent call last):
+    ...
+    ValueError: All parameters must be positive.
+
+    """
+    if (
+        coefficient_k <= 0
+        or coefficient_a <= 0
+        or coefficient_b <= 0
+        or coefficient_c <= 0
+        or return_period <= 0
+        or duration <= 0
+    ):
+        raise ValueError("All parameters must be positive.")
+    intensity = (coefficient_k * (return_period**coefficient_a)) / (
+        (duration + coefficient_b) ** coefficient_c
+    )
+    return intensity
+
+
+if __name__ == "__main__":
+    import doctest
+
+    doctest.testmod()