Merge aeb5534 into fcde9ff

examples/07_diagnostic_plot.py

@@ -0,0 +1,15 @@
+"""
+Diagnostic plot
+-------------------
+
+A diagnostic plot is a simultaneous plot of the drawdown and the
+logarithmic derivative of the drawdown in a log-log plot.
+Often, this plot is used to identify the right approach for the aquifer estimations.
+
+"""
+
+
+import welltestpy as wtp
+
+campaign = wtp.load_campaign("Cmp_UFZ-campaign.cmp")
+campaign.diagnostic_plot("well_0", "well_1")

welltestpy/data/campaignlib.py

@@ -108,7 +108,7 @@ class Campaign:
     """Class for a well based campaign.
 
     This is a class for a well based test campaign on a field site.
-    It has a name, a descrition and a timeframe.
+    It has a name, a description and a timeframe.
 
     Parameters
     ----------
@@ -118,7 +118,7 @@ class Campaign:
         The field site.
         Default: ``"Fieldsite"``
     wells : :class:`dict`, optional
-        The wells within the fild site. Keys are the well names and values
+        The wells within the field site. Keys are the well names and values
         are an instance of :class:`Well`.
         Default: ``None``
     wells : :class:`dict`, optional
@@ -437,6 +437,40 @@ def plot_wells(self, **kwargs):
         """
         return plotter.campaign_well_plot(self, **kwargs)
 
+    def diagnostic_plot(self, pumping_test, observation_well, **kwargs):
+        """Generate a diagnostic plot.
+
+        Parameters
+        ----------
+        pumping_test : :class:`str`
+            The pumping well that is saved in the campaign.
+
+        observation_well : :class:`str`
+            Observation point to make the diagnostic plot.
+
+        **kwargs
+            Keyword-arguments forwarded to :any:`campaign_well_plot`.
+        """
+        # check if this is a pumping test
+        if pumping_test in self.tests:
+            if not isinstance(self.tests[pumping_test], testslib.PumpingTest):
+                raise ValueError(
+                    f"diagnostic_plot: test '{pumping_test}' is not of instance PumpingTest!"
+                )
+            # check if the well is present
+            if observation_well in self.wells:
+                return self.tests[pumping_test].diagnostic_plot(
+                    observation_well=observation_well, **kwargs
+                )
+            else:
+                raise ValueError(
+                    f"diagnostic_plot: well '{observation_well}' could not be found!"
+                )
+        else:
+            raise ValueError(
+                f"diagnostic_plot: test '{pumping_test}' could not be found!"
+            )
+
     def save(self, path="", name=None):
         """Save the campaign to file.
 

welltestpy/data/testslib.py

@@ -14,7 +14,7 @@
 
 import numpy as np
 
-from ..tools import plotter
+from ..tools import plotter, diagnostic_plots
 from . import varlib, data_io
 from ..process import processlib
 
@@ -490,7 +490,7 @@ def plot(self, wells, exclude=None, fig=None, ax=None, **kwargs):
         ax : :class:`Axes`
             Axes where the plot should be done.
         wells : :class:`dict`
-            Dictonary containing the well classes sorted by name.
+            Dictionary containing the well classes sorted by name.
         exclude: :class:`list`, optional
             List of wells that should be excluded from the plot.
             Default: ``None``
@@ -508,6 +508,30 @@ def plot(self, wells, exclude=None, fig=None, ax=None, **kwargs):
             **kwargs,
         )
 
+    def diagnostic_plot(self, observation_well, **kwargs):
+        """Make a diagnostic plot.
+
+
+         Parameters
+         ----------
+        observation_well : :class:`str`
+             The observation well for the data to make the diagnostic plot.
+
+         Notes
+         -----
+         This will be used by the Campaign class.
+        """
+        if observation_well in self.observations:
+            observation = self.observations[observation_well]
+            rate = self.pumpingrate()
+            return diagnostic_plots.diagnostic_plot_pump_test(
+                observation=observation, rate=rate, **kwargs
+            )
+        else:
+            raise ValueError(
+                f"diagnostic_plot: well '{observation_well}' could not be found!"
+            )
+
     def save(self, path="", name=None):
         """Save a pumping test to file.
 

welltestpy/process/__init__.py

@@ -14,17 +14,20 @@
     combinepumptest
     filterdrawdown
     cooper_jacob_correction
+    smoothing_derivative
 """
 from .processlib import (
     normpumptest,
     combinepumptest,
     filterdrawdown,
     cooper_jacob_correction,
+    smoothing_derivative,
 )
 
 __all__ = [
     "normpumptest",
     "combinepumptest",
     "filterdrawdown",
     "cooper_jacob_correction",
+    "smoothing_derivative",
 ]

welltestpy/process/processlib.py

@@ -11,6 +11,7 @@
    combinepumptest
    filterdrawdown
    cooper_jacob_correction
+   smoothing_derivative
 """
 from copy import deepcopy as dcopy
 import numpy as np
@@ -23,6 +24,7 @@
     "combinepumptest",
     "filterdrawdown",
     "cooper_jacob_correction",
+    "smoothing_derivative",
 ]
 
 
@@ -273,3 +275,50 @@ def cooper_jacob_correction(observation, sat_thickness):
     observation(observation=head)
 
     return observation
+
+
+def smoothing_derivative(head, time, method="bourdet"):
+    """Calculate the derivative of the drawdown curve.
+
+    Parameters
+    ----------
+    head : :class: 'array'
+        An array with the observed head values.
+    time: :class: 'array'
+        An array with the time values for the observed head values.
+    method : :class:`str`, optional
+        Method to calculate the time derivative.
+        Default: "bourdet"
+
+    Returns
+    ---------
+    the derivative of the observed heads.
+
+    """
+    # create arrays for the input of head and time.
+    derhead = np.zeros(len(head))
+    t = np.arange(len(time))
+    if method == "bourdet":
+        for i in t[1:-1]:
+            # derivative approximation by Bourdet (1989)
+            dh = (
+                (
+                    (head[i] - head[i - 1])
+                    / (np.log(time[i]) - np.log(time[i - 1]))
+                    * (np.log(time[i + 1]) - np.log(time[i]))
+                )
+                + (
+                    (head[i + 1] - head[i])
+                    / (np.log(time[i + 1]) - np.log(time[i]))
+                    * (np.log(time[i]) - np.log(time[i - 1]))
+                )
+            ) / (
+                (np.log(time[i]) - np.log(time[i - 1]))
+                + (np.log(time[i + 1]) - np.log(time[i]))
+            )
+            derhead[i] = dh
+        return derhead
+    else:
+        raise ValueError(
+            f"smoothing_derivative: method '{method}' is unknown!"
+        )

welltestpy/tools/__init__.py

@@ -25,8 +25,9 @@
    plotparainteract
    plotparatrace
    plotsensitivity
+   diagnostic_plot_pump_test
 """
-from . import plotter, trilib
+from . import plotter, trilib, diagnostic_plots
 
 from .trilib import triangulate, sym
 
@@ -42,6 +43,8 @@
     plotsensitivity,
 )
 
+from .diagnostic_plots import diagnostic_plot_pump_test
+
 __all__ = [
     "triangulate",
     "sym",
@@ -54,5 +57,6 @@
     "plotparainteract",
     "plotparatrace",
     "plotsensitivity",
+    "diagnostic_plot_pump_test",
 ]
-__all__ += ["plotter", "trilib"]
+__all__ += ["plotter", "trilib", "diagnostic_plots"]

welltestpy/tools/diagnostic_plots.py

@@ -0,0 +1,123 @@
+"""
+welltestpy subpackage to make diagnostic plots.
+
+.. currentmodule:: welltestpy.tools.diagnostic_plots
+
+The following classes and functions are provided
+
+.. autosummary::
+   diagnostic_plot_pump_test
+
+"""
+# pylint: disable=C0103
+
+
+import numpy as np
+
+from ..process import processlib
+from . import plotter
+
+import matplotlib.pyplot as plt
+
+
+def diagnostic_plot_pump_test(
+    observation,
+    rate,
+    method="bourdet",
+    linthresh_time=1.0,
+    linthresh_head=1e-5,
+    fig=None,
+    ax=None,
+    plotname=None,
+    style="WTP",
+):
+    """plot the derivative with the original data.
+
+    Parameters
+    ----------
+    observation : :class:`welltestpy.data.Observation`
+         The observation to calculate the derivative.
+    rate : :class:`float`
+        Pumping rate.
+    method : :class:`str`, optional
+        Method to calculate the time derivative.
+        Default: "bourdet"
+    linthresh_time : :class: 'float'
+        Range of time around 0 that behaves linear.
+        Default: 1
+    linthresh_head : :class: 'float'
+        Range of head values around 0 that behaves linear.
+        Default: 1e-5
+    fig : Figure, optional
+        Matplotlib figure to plot on.
+        Default: None.
+    ax : :class:`Axes`
+        Matplotlib axes to plot on.
+        Default: None.
+    plotname : str, optional
+        Plot name if the result should be saved.
+        Default: None.
+    style : str, optional
+        Plot style.
+        Default: "WTP".
+
+     Returns
+     ---------
+     Diagnostic plot
+    """
+    head, time = observation()
+    head = np.array(head, dtype=float).reshape(-1)
+    time = np.array(time, dtype=float).reshape(-1)
+    if rate < 0:
+        head = head * -1
+    derivative = processlib.smoothing_derivative(
+        head=head, time=time, method=method
+    )
+    # setting variables
+    x = time
+    y = head
+    sx = time
+    sy = head
+    dx = time[1:-1]
+    dy = derivative[1:-1]
+
+    # plotting
+    if style == "WTP":
+        style = "ggplot"
+        font_size = plt.rcParams.get("font.size", 10.0)
+        keep_fs = True
+    with plt.style.context(style):
+        if keep_fs:
+            plt.rcParams.update({"font.size": font_size})
+        fig, ax = plotter._get_fig_ax(fig, ax)
+        ax.scatter(x, y, color="red", label="drawdown")
+        ax.plot(sx, sy, c="red")
+        ax.plot(dx, dy, c="black", linestyle="dashed")
+        ax.scatter(dx, dy, c="black", marker="+", label="time derivative")
+        ax.set_xscale("symlog", linthresh=linthresh_time)
+        ax.set_yscale("symlog", linthresh=linthresh_head)
+        ax.set_xlabel("$t$ in [s]", fontsize=16)
+        ax.set_ylabel("$h$ and $dh/dx$ in [m]", fontsize=16)
+        lgd = ax.legend(loc="upper left", facecolor="w")
+        min_v = min(np.min(head), np.min(dy))
+        max_v = max(np.max(head), np.max(dy))
+        max_e = int(np.ceil(np.log10(max_v)))
+        if min_v < linthresh_head:
+            min_e = -np.inf
+        else:
+            min_e = int(np.floor(np.log10(min_v)))
+        ax.set_ylim(10.0 ** min_e, 10.0 ** max_e)
+        yticks = [0 if min_v < linthresh_head else 10.0 ** min_e]
+        thresh_e = int(np.floor(np.log10(linthresh_head)))
+        first_e = thresh_e if min_v < linthresh_head else (min_e + 1)
+        yticks += list(10.0 ** np.arange(first_e, max_e + 1))
+        ax.set_yticks(yticks)
+        fig.tight_layout()
+        if plotname is not None:
+            fig.savefig(
+                plotname,
+                format="pdf",
+                bbox_extra_artists=(lgd,),
+                bbox_inches="tight",
+            )
+    return ax

welltestpy/tools/plotter.py

@@ -141,7 +141,7 @@ def campaign_plot(campaign, select_test=None, fig=None, style="WTP", **kwargs):
     fig : Figure, optional
         Matplotlib figure to plot on. The default is None.
     style : str, optional
-        Plot stlye. The default is "WTP".
+        Plot style. The default is "WTP".
     **kwargs : TYPE
         Keyword arguments forwarded to the tests plotting routines.
 
@@ -895,3 +895,4 @@ def plotsensitivity(
                 bbox_inches="tight",
             )
     return ax
+