TTim plots

docs/index.rst

@@ -78,7 +78,7 @@ Quick Example
             ml.solve()
 
             # Plot head contours at t=2 days
-            ml.contour(win=[-30, 55, -30, 30], ngr=40, t=2, labels=True, decimals=1)
+            ml.plots.contour(win=[-30, 55, -30, 30], ngr=40, t=2, labels=True, decimals=1)
 
 
     .. tab-item:: Result

pyproject.toml

@@ -29,11 +29,9 @@ classifiers = [
     "Operating System :: MacOS",
     "Programming Language :: Python",
     "Programming Language :: Python :: 3",
-    "Programming Language :: Python :: 3.7",
-    "Programming Language :: Python :: 3.8",
-    "Programming Language :: Python :: 3.9",
     "Programming Language :: Python :: 3.10",
     "Programming Language :: Python :: 3.11",
+    "Programming Language :: Python :: 3.12",
     "Programming Language :: Python :: 3 :: Only",
     "Topic :: Scientific/Engineering :: Hydrology",
 ]

ttim/circareasink.py

@@ -123,8 +123,10 @@ def disvecinf(self, x, y, aq=None):
             qy[:] = qr * (y - self.yc) / r
         return qx, qy
 
-    def plot(self):
-        plt.plot(
+    def plot(self, ax=None):
+        if ax is None:
+            _, ax = plt.subplots()
+        ax.plot(
             self.xc + self.R * np.cos(np.linspace(0, 2 * np.pi, 100)),
             self.yc + self.R * np.sin(np.linspace(0, 2 * np.pi, 100)),
             "k",

ttim/element.py

@@ -332,5 +332,5 @@ def run_after_solve(self):
         """
         pass
 
-    def plot(self):
+    def plot(self, ax=None):
         pass

ttim/linedoublet.py

@@ -174,8 +174,10 @@ def disvecinf(self, x, y, aq=None):
         rvy.shape = (self.nparam, aq.naq, self.model.npval)
         return rvx, rvy
 
-    def plot(self):
-        plt.plot([self.x1, self.x2], [self.y1, self.y2], "k")
+    def plot(self, ax=None):
+        if ax is None:
+            _, ax = plt.subplots()
+        ax.plot([self.x1, self.x2], [self.y1, self.y2], "k")
 
 
 class LeakyLineDoublet(LineDoubletHoBase, LeakyWallEquation):
@@ -376,5 +378,7 @@ def disvecinf(self, x, y, aq=None):
             rvy[i * ld.nparam : (i + 1) * ld.nparam, :] = qy
         return rvx, rvy
 
-    def plot(self):
-        plt.plot(self.xldlayout, self.yldlayout, "k")
+    def plot(self, ax=None):
+        if ax is None:
+            _, ax = plt.subplots()
+        ax.plot(self.xldlayout, self.yldlayout, "k")

ttim/linesink.py

@@ -168,8 +168,10 @@ def headinside(self, t):
             :, np.newaxis
         ] * self.discharge(t)
 
-    def plot(self):
-        plt.plot([self.x1, self.x2], [self.y1, self.y2], "k")
+    def plot(self, ax=None):
+        if ax is None:
+            _, ax = plt.subplots()
+        ax.plot([self.x1, self.x2], [self.y1, self.y2], "k")
 
 
 class LineSink(LineSinkBase):
@@ -442,8 +444,8 @@ def headinside(self, t, derivative=0):
             )
         return rv
 
-    def plot(self):
-        plt.plot(self.xlslayout, self.ylslayout, "k")
+    def plot(self, ax):
+        ax.plot(self.xlslayout, self.ylslayout, "k")
 
     def run_after_solve(self):
         for i in range(self.nls):
@@ -998,8 +1000,10 @@ def headinside(self, t):
             :, np.newaxis
         ] * self.discharge(t)
 
-    def plot(self):
-        plt.plot([self.x1, self.x2], [self.y1, self.y2], "k")
+    def plot(self, ax=None):
+        if ax is None:
+            _, ax = plt.subplots()
+        ax.plot([self.x1, self.x2], [self.y1, self.y2], "k")
 
 
 class HeadLineSinkHo(LineSinkHoBase, HeadEquationNores):

ttim/model.py

@@ -1,5 +1,6 @@
 # from .invlap import *
 import inspect  # Used for storing the input
+from warnings import warn
 
 import numpy as np
 
@@ -8,10 +9,10 @@
 
 # from .bessel import *
 from .invlapnumba import compute_laplace_parameters_numba, invlap, invlapcomp
-from .util import PlotTtim
+from .plots import PlotTtim
 
 
-class TimModel(PlotTtim):
+class TimModel:
     def __init__(
         self,
         kaq=[1, 1],
@@ -68,6 +69,25 @@ def __init__(
         if self.timmlmodel is not None:
             self.timmlmodel.solve()
 
+        self.plots = PlotTtim(self)
+        self.plot = self.plots.topview
+
+        # NOTE: reinstate later, after deprecation below is removed?
+        # self.xsection = self.plots.xsection
+
+    def xsection(*args, **kwargs):
+        raise DeprecationWarning(
+            "This method is deprecated. Use `ml.plots.head_along_line()` instead."
+        )
+
+    def contour(self, *args, **kwargs):
+        warn(
+            category=DeprecationWarning,
+            message="This method is deprecated. Use `ml.plots.contour()` instead.",
+            stacklevel=1,
+        )
+        self.plots.contour(*args, **kwargs)
+
     def __repr__(self):
         return "Model"
 
@@ -349,7 +369,7 @@ def velocomp(self, x, y, z, t, aq=None, layer_ltype=None):
                 )
                 qz = (h[1, 0] - h[0, 0]) / aq.c[
                     layer
-                ]  # TO DO include storage in leaky layer
+                ]  # TODO: include storage in leaky layer
             vz = qz / aq.porll[layer]
         else:  # in aquifer layer
             h = self.head(x, y, t, layers=layer, aq=aq, neglect_steady=True)
@@ -387,7 +407,7 @@ def velocomp(self, x, y, z, t, aq=None, layer_ltype=None):
                     )[:, 0]
             # this works because c[0] = 1e100 for impermeable top
             qztop = (h[1] - h[0]) / self.aq.c[layer]
-            # TO DO modify for infiltration in top aquifer
+            # TODO: modify for infiltration in top aquifer
             # if layer == 0:
             #    qztop += self.qztop(x, y)
             if layer < aq.naq - 1: