diff --git a/README.md b/README.md index 6a938b1..d71fef8 100644 --- a/README.md +++ b/README.md @@ -2,6 +2,19 @@ *from Greek: πέδον, pedon -> soil* +[![PyPI](https://img.shields.io/pypi/v/pedon?style=flat-square)](https://pypi.org/project/pedon/) +[![PyPi Supported Python Versions](https://img.shields.io/pypi/pyversions/pedon?style=flat-square)](https://pypi.org/project/pedon/) +[![Code Size](https://img.shields.io/github/languages/code-size/martinvonk/pedon?style=flat-square)](https://pypi.org/project/pedon/) +[![PyPi Downloads](https://img.shields.io/pypi/dm/pedon?style=flat-square) ![License](https://img.shields.io/pypi/l/pedon?style=flat-square)](https://pypi.org/project/pedon/) + +[![Tests](https://img.shields.io/github/actions/workflow/status/martinvonk/pedon/tests.yaml?style=flat-square)](https://github.com/martinvonk/pedon/actions/workflows/tests.yaml) +[![MyPy](https://img.shields.io/badge/type_checker-mypy-2A6DB2?style=flat-square)](https://mypy-lang.org/) +[![Format: isort](https://img.shields.io/badge/imports-isort-ef8336?style=flat-square)](https://pycqa.github.io/isort/index.html) +[![Format: Black](https://img.shields.io/badge/code_style-black-black?style=flat-square)](https://github.com/psf/black) +[![Linter: flake8](https://img.shields.io/badge/linter-flake8-yellowgreen?style=flat-square)](https://flake8.pycqa.org/) +[![Linter: ruff](https://img.shields.io/badge/linter-ruff-red?style=flat-square)](https://github.com/charliermarsh/ruff) + + Python package for (unsaturated) soil properties including pedotransfer functions. This package takes an object-oriented approach to soils, soil samples and soil models. Soil models that are available in this package are: - [Mualem-Van Genuchten](http://www.soilphysics.okstate.edu/teaching/soil-6583/references-folder/van%20Genuchten%201980.pdf) - [Brooks-Corey](https://www.wipp.energy.gov/library/cra/2009_cra/references/Others/Brooks_Corey_1964_Hydraulic_Properties_ERMS241117.pdf) diff --git a/src/pedon/__init__.py b/src/pedon/__init__.py index 3677d6f..5bc4911 100644 --- a/src/pedon/__init__.py +++ b/src/pedon/__init__.py @@ -2,5 +2,9 @@ # flake8: noqa from . import _params, soil, soilmodel from ._version import __version__ -from .soil import Soil, SoilSample -from .soilmodel import Brooks, Fredlund, Gardner, Genuchten, Panday, plot_hcf, plot_swrc +from .soil import Soil as Soil +from .soil import SoilSample as SoilSample +from .soilmodel import Brooks as Brooks +from .soilmodel import Genuchten as Genuchten +from .soilmodel import Panday as Panday +from .soilmodel import plot_hcf, plot_swrc diff --git a/src/pedon/_params.py b/src/pedon/_params.py index e61abaa..cd1eefb 100644 --- a/src/pedon/_params.py +++ b/src/pedon/_params.py @@ -22,7 +22,7 @@ "k_s": 100000.0, "theta_r": 0.2, "theta_s": 0.9, - "alpha": 0.15, + "alpha": 0.20, "n": 12, "l": 8, }, @@ -74,7 +74,7 @@ "theta_s": 0.4, "alpha": 0.02, "beta": 2.3, - "brook": 0.5, + "brook": 10.0, }, "p_min": { "k_s": 0.001, @@ -82,15 +82,15 @@ "theta_s": 0.2, "alpha": 0.001, "beta": 1.0, - "brook": -7, + "brook": 1.0, }, "p_max": { "k_s": 100000.0, "theta_r": 0.2, "theta_s": 0.5, - "alpha": 0.15, + "alpha": 0.30, "beta": 12, - "brook": 8, + "brook": 50.0, }, "swrc": { "k_s": False, diff --git a/src/pedon/_version.py b/src/pedon/_version.py index 3b93d0b..27fdca4 100644 --- a/src/pedon/_version.py +++ b/src/pedon/_version.py @@ -1 +1 @@ -__version__ = "0.0.2" +__version__ = "0.0.3" diff --git a/src/pedon/soil.py b/src/pedon/soil.py index 8e98506..b2828ed 100644 --- a/src/pedon/soil.py +++ b/src/pedon/soil.py @@ -9,7 +9,7 @@ from ._params import get_params from ._typing import FloatArray -from .soilmodel import Brooks, Genuchten, SoilModel +from .soilmodel import Brooks, Genuchten, SoilModel, get_soilmodel @dataclass @@ -116,14 +116,22 @@ def fit( W1: float | None = None, W2: float | None = None, return_res: bool = False, + k_s: float | None = None, ) -> SoilModel: """Same method as RETC""" + + theta = self.theta + k = self.k + if pbounds is None: pbounds = get_params(sm.__name__) - pbounds.loc["k_s", "p_ini"] = max(self.k) - pbounds.loc["k_s", "p_max"] = max(self.k) * 10 - pbounds.loc["theta_s", "p_ini"] = max(self.theta) - pbounds.loc["theta_s", "p_max"] = max(self.theta) + 0.02 + if k_s is not None: + pbounds = pbounds.drop("k_s") + else: + pbounds.loc["k_s", "p_ini"] = max(k) + pbounds.loc["k_s", "p_max"] = max(k) * 10 + pbounds.loc["theta_s", "p_ini"] = max(theta) + pbounds.loc["theta_s", "p_max"] = max(theta) + 0.02 if weights is None: weights = ones(self.h.shape) @@ -132,14 +140,17 @@ def fit( W1 = 0.1 if W2 is None: - M = len(self.k) + len(self.theta) - N = len(self.theta) - W2 = (M - N) * sum(weights * self.theta) / (N * sum(weights * self.k)) + M = len(k) + len(theta) + N = len(theta) + W2 = (M - N) * sum(weights * theta) / (N * sum(weights * k)) def fit_staring(p: FloatArray) -> FloatArray: - sml = sm(**dict(zip(pbounds.index, p))) - theta_diff = sml.theta(h=self.h) - self.theta - k_diff = log(sml.k(h=self.h)) - log(self.k) + est_pars = dict(zip(pbounds.index, p)) + if k_s is not None: + est_pars["k_s"] = k_s + sml = sm(**est_pars) + theta_diff = sml.theta(h=self.h) - theta + k_diff = log(sml.k(h=self.h)) - log(k) diff = append(weights * theta_diff, weights * W1 * W2 * k_diff) return diff @@ -152,6 +163,8 @@ def fit_staring(p: FloatArray) -> FloatArray: ), ) opt_pars = dict(zip(pbounds.index, res.x)) + if k_s is not None: + opt_pars["k_s"] = k_s opt_sm = sm(**opt_pars) if return_res: return opt_sm, {"res": res} @@ -406,13 +419,25 @@ class Soil: source: str | None = None description: str | None = None - def from_name(self, sm: Type[SoilModel]) -> "Soil": + def from_name(self, sm: Type[SoilModel] | SoilModel | str) -> "Soil": + if isinstance(sm, SoilModel): + if hasattr(sm, "__name__"): + smn = sm.__name__ + else: + smn = sm.__class__.__name__ + sm = type(sm) + elif isinstance(sm, str): + smn = sm + sm = get_soilmodel(smn) + + else: + raise ValueError( + f"Argument must either be Type[SoilModel] | SoilModel | str," + f"not {type(sm)}" + ) + path = Path(__file__).parent / "datasets/Soil_Parameters.xlsx" - ser = ( - read_excel(path, sheet_name=sm.__name__, index_col=0) - .loc[self.name] - .to_dict() - ) + ser = read_excel(path, sheet_name=smn, index_col=0).loc[self.name].to_dict() # path = Path(__file__).parent / f"datasets/{sm.__name__}.csv" # ser = read_csv(path, index_col=["name"]).loc[self.name].to_dict() self.__setattr__("type", ser.pop("soil type")) @@ -421,7 +446,31 @@ def from_name(self, sm: Type[SoilModel]) -> "Soil": self.__setattr__("model", sm(**ser)) return self + @staticmethod + def list_names(sm: Type[SoilModel] | SoilModel | str) -> list[str]: + if isinstance(sm, SoilModel): + if hasattr(sm, "__name__"): + smn = sm.__name__ + else: + smn = sm.__class__.__name__ + sm = type(sm) + elif isinstance(sm, str): + smn = sm + sm = get_soilmodel(smn) + + else: + raise ValueError( + f"Argument must either be Type[SoilModel] | SoilModel | str," + f"not {type(sm)}" + ) + + path = Path(__file__).parent / "datasets/Soil_Parameters.xlsx" + names = read_excel(path, sheet_name=smn).loc[:, "name"].to_list() + return names + def from_staring(self, year: str = "2018") -> "Soil": + if year not in ("2001", 2001, "2018", 2018): + raise ValueError(f"Year must either be '2001' or '2018', not {year}") path = Path(__file__).parent / f"datasets/Staring_{year}.xlsx" parameters = read_excel(path, sheet_name="parameters", index_col=0) ser = parameters.loc[self.name].to_dict() diff --git a/src/pedon/soilmodel.py b/src/pedon/soilmodel.py index 2ff1716..dc4ec18 100644 --- a/src/pedon/soilmodel.py +++ b/src/pedon/soilmodel.py @@ -1,5 +1,5 @@ -from dataclasses import dataclass -from typing import Protocol +from dataclasses import dataclass, field +from typing import Protocol, Type, runtime_checkable import matplotlib.pyplot as plt from numpy import abs as npabs @@ -8,6 +8,7 @@ from ._typing import FloatArray +@runtime_checkable class SoilModel(Protocol): def theta(self, h: FloatArray) -> FloatArray: """Method to calculate the soil moisture content from the pressure head h""" @@ -17,11 +18,15 @@ def s(self, h: FloatArray) -> FloatArray: """Method to calculate the effective saturation from the pressure head h""" ... + def k_r(self, h: FloatArray, s: FloatArray | None = None) -> FloatArray: + """Method to calcualte the relative permeability from the pressure head h""" + ... + def k(self, h: FloatArray, s: FloatArray | None = None) -> FloatArray: """Method to calcualte the permeability from the pressure head h""" ... - def plot(self): + def plot(self, ax: plt.Axes | None = None) -> plt.Axes: """Method to plot the soil water retention curve""" ... @@ -55,13 +60,16 @@ def theta(self, h: FloatArray) -> FloatArray: def s(self, h: FloatArray) -> FloatArray: return (self.theta(h) - self.theta_r) / (self.theta_s - self.theta_r) - def k(self, h: FloatArray, s: FloatArray | None = None) -> FloatArray: + def k_r(self, h: FloatArray, s: FloatArray | None = None) -> FloatArray: if s is None: s = self.s(h) - return self.k_s * s**self.l * (1 - (1 - s ** (1 / self.m)) ** self.m) ** 2 + return s**self.l * (1 - (1 - s ** (1 / self.m)) ** self.m) ** 2 - def plot(self): - return plot_swrc(self) + def k(self, h: FloatArray, s: FloatArray | None = None) -> FloatArray: + return self.k_s * self.k_r(h=h, s=s) + + def plot(self, ax: plt.Axes | None = None) -> plt.Axes: + return plot_swrc(self, ax=ax) @dataclass @@ -103,13 +111,16 @@ def s(self, h: FloatArray) -> FloatArray: s[h >= self.h_b] = (h[h >= self.h_b] / self.h_b) ** -self.l return s - def k(self, h: FloatArray, s: FloatArray | None = None) -> FloatArray: + def k_r(self, h: FloatArray, s: FloatArray | None = None) -> FloatArray: if s is None: s = self.s(h) - return self.k_s * s ** (3 + 2 / self.l) + return s ** (3 + 2 / self.l) - def plot(self): - return plot_swrc(self) + def k(self, h: FloatArray, s: FloatArray | None = None) -> FloatArray: + return self.k_s * self.k_r(h=h, s=s) + + def plot(self, ax: plt.Axes | None = None) -> plt.Axes: + return plot_swrc(self, ax=ax) @dataclass @@ -132,14 +143,17 @@ def theta(self, h: FloatArray) -> FloatArray: def s(self, h: FloatArray) -> FloatArray: return (self.theta(h) - self.theta_r) / (self.theta_s - self.theta_r) - def k(self, h: FloatArray, s: FloatArray | None = None) -> FloatArray: + def k_r(self, h: FloatArray, s: FloatArray | None = None) -> FloatArray: if s is not None: theta = s * (self.theta_s - self.theta_r) + self.theta_r - return self.k_s * self.a * theta**self.m - return self.k_s * (self.a / (self.b + npabs(h) ** self.m)) + return self.a * theta**self.m + return self.a / (self.b + npabs(h) ** self.m) + + def k(self, h: FloatArray, s: FloatArray | None = None) -> FloatArray: + return self.k_s * self.k_r(h=h, s=s) - def plot(self): - return plot_swrc(self) + def plot(self, ax: plt.Axes | None = None) -> plt.Axes: + return plot_swrc(self, ax=ax) @dataclass @@ -150,29 +164,37 @@ class Panday: """ k_s: float - theta_r: float - theta_s: float + theta_r: float = field(repr=False) + theta_s: float = field(repr=False) alpha: float # alpha beta: float # n brook: float # brooks-corey l + h_b: float = field(default=0.0, repr=False) + sr: float = field(init=False, repr=True) + gamma: float = field(init=False, repr=False) # 1 - 1 / beta + sy: float = field(init=False, repr=False) def __post_init__(self): - self.gamma = 1 - 1 / self.beta # m 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) def theta(self, h: FloatArray) -> FloatArray: return (self.sr + self.s(h) * (1 - self.sr)) * self.theta_s def s(self, h: FloatArray) -> FloatArray: - return (1 + npabs(self.alpha * h) ** self.beta) ** -self.gamma + return (1 + npabs(self.alpha * (h - self.h_b)) ** self.beta) ** -self.gamma - def k(self, h: FloatArray, s: FloatArray | None = None) -> FloatArray: + def k_r(self, h: FloatArray, s: FloatArray | None = None) -> FloatArray: if s is None: s = self.s(h) - return self.k_s * s**self.brook + return s**self.brook - def plot(self): - return plot_swrc(self) + def k(self, h: FloatArray, s: FloatArray | None = None) -> FloatArray: + return self.k_s * self.k_r(h=h, s=s) + + def plot(self, ax: plt.Axes | None = None) -> plt.Axes: + return plot_swrc(self, ax=ax) @dataclass @@ -195,7 +217,7 @@ def theta(self, h: FloatArray) -> FloatArray: def s(self, h: FloatArray) -> FloatArray: return self.theta(h) / self.theta_s - def k(self, h: FloatArray, s: FloatArray | None = None): + def k_r(self, h: FloatArray, s: FloatArray | None = None): if s is not None: raise NotImplementedError( "Can only calculate the hydraulic conductivity" @@ -232,11 +254,24 @@ def theta_d( / exp(y) * theta_d(exp(y), self.a, self.n, self.m, self.theta_s) ) + return teller / noemer + + def k(self, h: FloatArray, s: FloatArray | None = None) -> FloatArray: + return self.k_s * self.k_r(h=h, s=s) - return self.k_s * (teller / noemer) + def plot(self, ax: plt.Axes | None = None) -> plt.Axes: + return plot_swrc(self, ax=ax) - def plot(self): - return plot_swrc(self) + +def get_soilmodel(soilmodel_name: str) -> Type[SoilModel]: + sms = { + "Genuchten": Genuchten, + "Brooks": Brooks, + "Gardner": Gardner, + "Panday": Panday, + "Fredlund": Fredlund, + } + return sms[soilmodel_name] def plot_swrc( @@ -278,7 +313,6 @@ def plot_hcf( ax.set_xscale("log") h = logspace(-6, 10, num=1000) - k = sm.k(h=h) ax.plot(k, h, label=sm.__class__.__name__, **kwargs)