Merge pull request #2 from martinvonk/dev

Update main to v0.0.3

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)

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

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,23 +74,23 @@
             "theta_s": 0.4,
             "alpha": 0.02,
             "beta": 2.3,
-            "brook": 0.5,
+            "brook": 10.0,
         },
         "p_min": {
             "k_s": 0.001,
             "theta_r": 1e-5,
             "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,

src/pedon/_version.py

@@ -1 +1 @@
-__version__ = "0.0.2"
+__version__ = "0.0.3"

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()

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)