diff --git a/src/pedon/_version.py b/src/pedon/_version.py index 27fdca4..81f0fde 100644 --- a/src/pedon/_version.py +++ b/src/pedon/_version.py @@ -1 +1 @@ -__version__ = "0.0.3" +__version__ = "0.0.4" diff --git a/src/pedon/soilmodel.py b/src/pedon/soilmodel.py index dc4ec18..a0f6286 100644 --- a/src/pedon/soilmodel.py +++ b/src/pedon/soilmodel.py @@ -3,7 +3,7 @@ import matplotlib.pyplot as plt from numpy import abs as npabs -from numpy import exp, full, linspace, log, logspace +from numpy import exp, full, linspace, log, logspace, log10 from ._typing import FloatArray @@ -26,6 +26,10 @@ def k(self, h: FloatArray, s: FloatArray | None = None) -> FloatArray: """Method to calcualte the permeability from the pressure head h""" ... + def h(self, theta: FloatArray) -> FloatArray: + """Method to calcualte the pressure head h from the water content""" + ... + def plot(self, ax: plt.Axes | None = None) -> plt.Axes: """Method to plot the soil water retention curve""" ... @@ -33,7 +37,7 @@ def plot(self, ax: plt.Axes | None = None) -> plt.Axes: @dataclass class Genuchten: - """Mualem- van Genuchten Soil Model + """Mualem-van Genuchten Soil Model van Genuchten, M. Th. (1970) - A Closed-form Equation for Predicting the Hydraulic Conductivity of Unsaturated Soil @@ -68,6 +72,11 @@ def k_r(self, h: FloatArray, s: FloatArray | None = None) -> FloatArray: def k(self, h: FloatArray, s: FloatArray | None = None) -> FloatArray: return self.k_s * self.k_r(h=h, s=s) + def h(self, theta: FloatArray) -> FloatArray: + se = (theta - self.theta_r) / (self.theta_s - self.theta_r) + h = 1 / self.alpha * ((1 / se) ** (1 / self.m) - 1) ** (1 / self.n) + return h + def plot(self, ax: plt.Axes | None = None) -> plt.Axes: return plot_swrc(self, ax=ax) @@ -119,6 +128,22 @@ def k_r(self, h: FloatArray, s: FloatArray | None = None) -> FloatArray: def k(self, h: FloatArray, s: FloatArray | None = None) -> FloatArray: return self.k_s * self.k_r(h=h, s=s) + def h(self, theta: FloatArray) -> FloatArray: + if isinstance(theta, float): + if theta >= self.theta_r: + return self.h_b * ((theta - self.theta_r) / (self.s(theta))) ** ( + -1 / self.l + ) + else: + return self.h_b + else: + h = full(theta.shape, self.h_b) + mask = theta >= self.theta_r + h[mask] = self.h_b * ( + (theta[mask] - self.theta_r) / (self.s(theta[mask])) + ) ** (-1 / self.l) + return h + def plot(self, ax: plt.Axes | None = None) -> plt.Axes: return plot_swrc(self, ax=ax) @@ -152,6 +177,9 @@ def k_r(self, h: FloatArray, s: FloatArray | None = None) -> FloatArray: def k(self, h: FloatArray, s: FloatArray | None = None) -> FloatArray: return self.k_s * self.k_r(h=h, s=s) + def h(self, theta: FloatArray) -> FloatArray: + return (npabs(theta) / self.a) ** (-1 / self.b) + def plot(self, ax: plt.Axes | None = None) -> plt.Axes: return plot_swrc(self, ax=ax) @@ -177,7 +205,11 @@ class Panday: def __post_init__(self): self.sr = self.theta_r / self.theta_s # theta_r / theta_s self.gamma = 1 - 1 / self.beta # m - self.sy = self.theta_s - self.theta_r - self.theta(10**2) + theta_fc = ( + self.beta ** -(0.60 * (2 + log10(self.k_s))) * (self.theta_s - self.theta_r) + + self.theta_r + ) + self.sy = self.theta_s - theta_fc def theta(self, h: FloatArray) -> FloatArray: return (self.sr + self.s(h) * (1 - self.sr)) * self.theta_s @@ -193,6 +225,11 @@ def k_r(self, h: FloatArray, s: FloatArray | None = None) -> FloatArray: def k(self, h: FloatArray, s: FloatArray | None = None) -> FloatArray: return self.k_s * self.k_r(h=h, s=s) + def h(self, theta: FloatArray) -> FloatArray: + se = (theta - self.theta_r) / (self.theta_s - self.theta_r) + h = 1 / self.alpha * ((1 / se) ** (1 / self.gamma) - 1) ** (1 / self.beta) + return h + def plot(self, ax: plt.Axes | None = None) -> plt.Axes: return plot_swrc(self, ax=ax) @@ -259,6 +296,11 @@ def theta_d( def k(self, h: FloatArray, s: FloatArray | None = None) -> FloatArray: return self.k_s * self.k_r(h=h, s=s) + def h(self, theta: FloatArray) -> FloatArray: + return self.a * (exp(self.theta_s / theta) ** (1 / self.m) - exp(1)) ** ( + 1 / self.n + ) + def plot(self, ax: plt.Axes | None = None) -> plt.Axes: return plot_swrc(self, ax=ax) @@ -275,10 +317,7 @@ def get_soilmodel(soilmodel_name: str) -> Type[SoilModel]: def plot_swrc( - sm: SoilModel, - saturation: bool = False, - ax: plt.Axes | None = None, - **kwargs: dict, + sm: SoilModel, saturation: bool = False, ax: plt.Axes | None = None, **kwargs ) -> plt.Axes: """Plot soil water retention curve""" @@ -294,7 +333,12 @@ def plot_swrc( else: sw = sm.theta(h=h) - ax.plot(sw, -h, label=sm.__class__.__name__, **kwargs) + if "label" in kwargs: + label = kwargs.pop("label") + else: + label = getattr(getattr(sm, "__class__"), "__name__") + + ax.plot(sw, -h, label=label, **kwargs) ax.set_ylim(1e-3, 1e6) ax.grid(True) return ax @@ -303,7 +347,7 @@ def plot_swrc( def plot_hcf( sm: SoilModel, ax: plt.Axes | None = None, - **kwargs: dict, + **kwargs, ) -> plt.Axes: """Plot the hydraulic conductivity function""" @@ -315,7 +359,12 @@ def plot_hcf( h = logspace(-6, 10, num=1000) k = sm.k(h=h) - ax.plot(k, h, label=sm.__class__.__name__, **kwargs) + if "label" in kwargs: + label = kwargs.pop("label") + else: + label = getattr(getattr(sm, "__class__"), "__name__") + + ax.plot(k, h, label=label, **kwargs) ax.set_ylim(1e-3, 1e6) ax.set_xlim() ax.grid(True) diff --git a/tests/test_datasets.py b/tests/test_datasets.py index 44fab47..b8c64d9 100644 --- a/tests/test_datasets.py +++ b/tests/test_datasets.py @@ -1,17 +1,21 @@ import pedon as pe -def test_sample_staring_2018(): +def test_sample_staring_2018() -> None: pe.soil.SoilSample().from_staring("B01", year="2018") -def test_sample_staring_2001(): +def test_sample_staring_2001() -> None: pe.soil.SoilSample().from_staring("B02", year="2001") -def test_soil_from_name(): +def test_soil_from_name() -> None: pe.soil.Soil("VS2D_Del Monte Sand").from_name(pe.Brooks) -def test_soil_from_staring(): +def test_soil_from_staring() -> None: pe.soil.Soil("O01").from_staring() + + +def test_soil_list_genuchten() -> None: + pe.soil.Soil.list_names(pe.Genuchten)