12312733157

Committed 13 Dec 2024 09:08AM UTC coverage: 77.784% (+2.4%) from 75.342%

Build # 12312733157

Build Type

Pull #81

github

Committed by

web-flow

Commit Message

Merge fcd788ee4 into f9283bd02

Pull Request Pull Request #81: Add cross-section models

Run Details

503 of 632 new or added lines in 12 files covered. (79.59%)

42 existing lines in 3 files now uncovered.

3011 of 3871 relevant lines covered (77.78%)

2.33 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

76.54

/ttim/fit.py

import warnings
from typing import Iterable

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from scipy.optimize import least_squares

# import lmfit


class Calibrate:
    def __init__(self, model):
        """Initialize Calibration class.

        Parameters
        ----------
        model : ttim.Model
            model to calibrate
        """
        self.model = model
        self.parameters = pd.DataFrame(
            columns=[
                "layers",
                "optimal",
                "std",
                "perc_std",
                "pmin",
                "pmax",
                "initial",
                "inhoms",
                "parray",
            ]
        )
        self.seriesdict = {}
        self.seriesinwelldict = {}

    def set_parameter(
        self, name=None, layers=None, initial=0, pmin=-np.inf, pmax=np.inf, inhoms=None
    ):
        """Set parameter to be optimized.

        Parameters
        ----------
        name : str
            parameter name, can include layer information.
            name can be 'kaq', 'Saq' or 'c'. A number after the parameter
            name denotes the layer number, i.e. 'kaq0' refers to the hydraulic
            conductivity of layer 0.
            name also supports layer ranges, entered by adding a '_' and a
            layer number, i.e. 'kaq0_3' denotes conductivity for layers 0 up to
            and including 3.
        initial : float, optional
            initial value for the parameter (the default is 0)
        pmin : float, optional
            lower bound for parameter value (the default is -np.inf)
        pmax : float, optional
            upper bound for paramater value (the default is np.inf)
        """
        assert isinstance(name, str), "Error: name must be string"
        # find numbers in name str for support layer ranges
        # layers_from_name = re.findall(r"\d+", name)

        if isinstance(layers, Iterable):
            from_lay = min(layers)
            to_lay = max(layers)
            if from_lay == to_lay:
                to_lay += 1
            if (np.diff(layers) > 1).any():
                warnings.warn(
                    "Non-consecutive layers are not supported. "
                    f"Setting parameter '{name}' for layers {from_lay} - {to_lay}.",
                    stacklevel=1,
                )
        elif isinstance(layers, int):
            from_lay = layers
            to_lay = layers + 1
        else:
            raise DeprecationWarning(
                "Setting layers in the parameter name is no longer supported. "
                f"Set layers keyword argument for parameter '{name}'."
            )

        # get aquifer information and create list if necessary
        if inhoms is None:
            aq = [self.model.aq]
        elif not isinstance(inhoms, (list, tuple)):
            aq = [inhoms]
        else:
            aq = inhoms

        # convert aquifer names to aquifer objects
        for i, iaq in enumerate(aq):
            if isinstance(iaq, str):
                aq[i] = self.model.aq.inhomdict[iaq]

        p = None
        # from_lay, to_lay = [int(i) for i in layers_from_name]
        if name[:3] == "kaq":
            p = aq[0].kaq[from_lay : to_lay + 1]
            param = name[:3]
        elif name[:3] == "Saq":
            p = aq[0].Saq[from_lay : to_lay + 1]
            param = name[:3]
        elif name[0] == "c":
            p = aq[0].c[from_lay : to_lay + 1]
            param = name[0]
        elif name[:3] == "Sll":
            p = aq[0].Sll[from_lay : to_lay + 1]
            param = name[:3]
        elif name[0:8] == "kzoverkh":
            p = aq[0].kzoverkh[from_lay : to_lay + 1]
            param = name[:8]

        if p is None:  # no parameter set
            print("parameter name not recognized or no parameter ref supplied")
            return

        # set all other aquifer parameter references to same array
        if len(aq) > 1:
            for iaq in aq:
                setattr(iaq, param, p)

        if inhoms is None:
            pname = name
        else:
            pname = f"{name}_{'_'.join([iaq.name for iaq in aq])}"
        self.parameters.loc[pname] = {
            "layers": layers,
            "optimal": initial,
            "std": None,
            "perc_std": None,
            "pmin": pmin,
            "pmax": pmax,
            "initial": initial,
            "inhoms": aq if inhoms is not None else None,
            "parray": p[:],
        }

    def set_parameter_by_reference(
        self, name=None, parameter=None, initial=0, pmin=-np.inf, pmax=np.inf
    ):
        """Set parameter to be optimized.

        Parameters
        ----------
        name : str
            parameter name
        parameter : np.array
            array reference containing the parameter to be optimized. must be
            specified as reference, i.e. w.rc[0:]
        initial : float, optional
            initial value for the parameter (the default is 0)
        pmin : float, optional
            lower bound for parameter value (the default is -np.inf)
        pmax : float, optional
            upper bound for paramater value (the default is np.inf)
        """
        assert isinstance(name, str), "Error: name must be string"
        if parameter is not None:
            assert isinstance(
                parameter, np.ndarray
            ), "Error: parameter needs to be numpy array"
            p = parameter
        self.parameters.loc[name] = {
            "optimal": initial,
            "std": None,
            "perc_std": None,
            "pmin": pmin,
            "pmax": pmax,
            "initial": initial,
            "parray": p[:],
        }

    def series(self, name, x, y, layer, t, h, weights=None):
        """Method to add observations to Calibration object.

        Parameters
        ----------
        name : str
            name of series
        x : float
            x-coordinate
        y : float
            y-coordinate
        layer : int
            layer number, 0-indexed
        t : np.array
            array containing timestamps of timeseries
        h : np.array
            array containing timeseries values, i.e. head observations
        """
        s = Series(x, y, layer, t, h, weights=weights)
        self.seriesdict[name] = s

    def seriesinwell(self, name, element, t, h):
        """Method to add observations to Calibration object.

        Parameters
        ----------
        name : str
            name of series
        element: element object with headinside function
        t : np.array
            array containing timestamps of timeseries
        h : np.array
            array containing timeseries values, i.e. head observations
        """
        e = SeriesInWell(element, t, h)
        self.seriesinwelldict[name] = e

    def residuals(self, p, printdot=False, weighted=True, layers=None, series=None):
        """Method to calculate residuals given certain parameters.

        Parameters
        ----------
        p : np.array
            array containing parameter values
        printdot : bool, optional
            print dot for each function call

        Returns
        -------
        np.array
            array containing all residuals
        """
        if printdot:
            print(".", end="")
        # set the values of the variables

        if printdot == 7:
            print(p)

        if layers is None:
            layers = range(self.model.aq.naq)

        for i, k in enumerate(self.parameters.index):
            # [:] needed to do set value in array
            self.parameters.loc[k, "parray"][:] = p[i]

        self.model.solve(silent=True)

        rv = np.empty(0)
        cal_series = self.seriesdict.keys() if series is None else series
        for key in cal_series:
            s = self.seriesdict[key]
            if s.layer not in layers:
                continue
            h = self.model.head(s.x, s.y, s.t, layers=s.layer)
            w = s.weights if ((s.weights is not None) and weighted) else np.ones_like(h)
            rv = np.append(rv, (s.h - h) * w)
        for key in self.seriesinwelldict:
            s = self.seriesinwelldict[key]
            h = s.element.headinside(s.t)[0]
            rv = np.append(rv, s.h - h)
        return rv

    def residuals_lmfit(self, lmfitparams, printdot=False):
        vals = lmfitparams.valuesdict()
        p = np.array([vals[k] for k in self.parameters.index])
        # p = np.array([vals[k] for k in vals])
        return self.residuals(p, printdot)

    def fit_least_squares(self, report=True, diff_step=1e-4, xtol=1e-8, method="lm"):
        self.fitresult = least_squares(
            self.residuals,
            self.parameters.initial.values,
            args=(True,),
            bounds=(self.parameters.pmin.values, self.parameters.pmax.values),
            method=method,
            diff_step=diff_step,
            xtol=xtol,
            x_scale="jac",
        )
        print("", flush=True)
        # Call residuals to specify optimal values for model
        res = self.residuals(self.fitresult.x)
        for ipar in self.parameters.index:
            self.parameters.loc[ipar, "optimal"] = self.parameters.loc[ipar, "parray"][
                0
            ]
        nparam = len(self.fitresult.x)
        H = self.fitresult.jac.T @ self.fitresult.jac
        sigsq = np.var(res, ddof=nparam)
        self.covmat = np.linalg.inv(H) * sigsq
        self.sig = np.sqrt(np.diag(self.covmat))
        D = np.diag(1 / self.sig)
        self.cormat = D @ self.covmat @ D
        self.parameters["std"] = self.sig
        self.parameters["perc_std"] = self.sig / self.parameters["optimal"] * 100
        if report:
            print(self.parameters)
            print(self.sig)
            print(self.covmat)
            print(self.cormat)

    def fit_lmfit(self, report=False, printdot=True):
        import lmfit

        self.lmfitparams = lmfit.Parameters()
        for name in self.parameters.index:
            p = self.parameters.loc[name]
            self.lmfitparams.add(name, value=p["initial"], min=p["pmin"], max=p["pmax"])
        fit_kws = {"epsfcn": 1e-4}
        self.fitresult = lmfit.minimize(
            self.residuals_lmfit,
            self.lmfitparams,
            method="leastsq",
            kws={"printdot": printdot},
            **fit_kws,
        )
        print("", flush=True)
        print(self.fitresult.message)
        if self.fitresult.success:
            for name in self.parameters.index:
                self.parameters.loc[name, "optimal"] = (
                    self.fitresult.params.valuesdict()[name]
                )
            if hasattr(self.fitresult, "covar"):
                self.parameters["std"] = np.sqrt(np.diag(self.fitresult.covar))
                self.parameters["perc_std"] = (
                    100 * self.parameters["std"] / np.abs(self.parameters["optimal"])
                )
            else:
                self.parameters["std"] = np.nan
                self.parameters["perc_std"] = np.nan
        if report:
            print(lmfit.fit_report(self.fitresult))

    def fit(self, report=False, printdot=True):
        # current default fitting routine is lmfit
        # return self.fit_least_squares(report) # does not support bounds by default
        return self.fit_lmfit(report, printdot)

    def rmse(self, weighted=True, layers=None):
        """Calculate root-mean-squared-error.

        Returns
        -------
        float
            return rmse value
        """
        r = self.residuals(
            self.parameters["optimal"].values, weighted=weighted, layers=layers
        )
        return np.sqrt(np.mean(r**2))

    def topview(self, ax=None, layers=None, labels=True):
        """Plot topview of model with calibration points.

        Parameters
        ----------
        ax : matplotlib.axes.Axes, optional
            axes to plot on (the default is None, which creates a new figure)
        """
        if ax is None:
            _, ax = plt.subplots()
            # self.model.plots.topview(ax=ax)
        for key, s in self.seriesdict.items():
            if layers is None or s.layer in layers:
                ax.plot(s.x, s.y, "ko")
                if labels:
                    ax.text(s.x, s.y, key, ha="left", va="bottom")
        return ax

    def xsection(self, ax=None, labels=True):
        """Plot cross-section of model with calibration points.

        Parameters
        ----------
        ax : matplotlib.axes.Axes, optional
            axes to plot on (the default is None, which creates a new figure)
        """
        if ax is None:
            _, ax = plt.subplots()
        for key, s in self.seriesdict.items():
            aq = self.model.aq.find_aquifer_data(s.x, s.y)
            ztop = aq.z[0]
            zpb_top = aq.zaqtop[s.layer]
            zpb_bot = aq.zaqbot[s.layer]
            ax.plot([s.x, s.x], [zpb_top, zpb_bot], c="k", ls="dotted")
            ax.plot([s.x, s.x], [ztop + 1, zpb_top], c="k", ls="solid", lw=1.0)
            if labels:
                ax.text(s.x, ztop + 1, key, ha="left", va="bottom", rotation=45)
        return ax


class Series:
    def __init__(self, x, y, layer, t, h, weights=None):
        self.x = x
        self.y = y
        self.layer = layer
        self.t = t
        self.h = h
        self.weights = weights


class SeriesInWell:
    def __init__(self, element, t, h):
        self.element = element
        self.t = t
        self.h = h

1	import warnings	3✔
2	from typing import Iterable	3✔
3
4	import matplotlib.pyplot as plt	3✔
5	import numpy as np	3✔
6	import pandas as pd	3✔
7	from scipy.optimize import least_squares	3✔
8
9	# import lmfit
10
11
12	class Calibrate:	3✔
13	def __init__(self, model):	3✔
14	"""Initialize Calibration class.
15
16	Parameters
17	----------
18	model : ttim.Model
19	model to calibrate
20	"""
21	self.model = model	3✔
22	self.parameters = pd.DataFrame(	3✔
23	columns=[
24	"layers",
25	"optimal",
26	"std",
27	"perc_std",
28	"pmin",
29	"pmax",
30	"initial",
31	"inhoms",
32	"parray",
33	]
34	)
35	self.seriesdict = {}	3✔
36	self.seriesinwelldict = {}	3✔
37
38	def set_parameter(	3✔
39	self, name=None, layers=None, initial=0, pmin=-np.inf, pmax=np.inf, inhoms=None
40	):
41	"""Set parameter to be optimized.
42
43	Parameters
44	----------
45	name : str
46	parameter name, can include layer information.
47	name can be 'kaq', 'Saq' or 'c'. A number after the parameter
48	name denotes the layer number, i.e. 'kaq0' refers to the hydraulic
49	conductivity of layer 0.
50	name also supports layer ranges, entered by adding a '_' and a
51	layer number, i.e. 'kaq0_3' denotes conductivity for layers 0 up to
52	and including 3.
53	initial : float, optional
54	initial value for the parameter (the default is 0)
55	pmin : float, optional
56	lower bound for parameter value (the default is -np.inf)
57	pmax : float, optional
58	upper bound for paramater value (the default is np.inf)
59	"""
60	assert isinstance(name, str), "Error: name must be string"	3✔
61	# find numbers in name str for support layer ranges
62	# layers_from_name = re.findall(r"\d+", name)
63
64	if isinstance(layers, Iterable):	3✔
65	from_lay = min(layers)	3✔
66	to_lay = max(layers)	3✔
67	if from_lay == to_lay:	3✔
NEW 68	to_lay += 1	×
69	if (np.diff(layers) > 1).any():	3✔
NEW 70	warnings.warn(	×
71	"Non-consecutive layers are not supported. "
72	f"Setting parameter '{name}' for layers {from_lay} - {to_lay}.",
73	stacklevel=1,
74	)
75	elif isinstance(layers, int):	3✔
76	from_lay = layers	3✔
77	to_lay = layers + 1	3✔
78	else:
NEW 79	raise DeprecationWarning(	×
80	"Setting layers in the parameter name is no longer supported. "
81	f"Set layers keyword argument for parameter '{name}'."
82	)
83
84	# get aquifer information and create list if necessary
85	if inhoms is None:	3✔
86	aq = [self.model.aq]	3✔
NEW 87	elif not isinstance(inhoms, (list, tuple)):	×
NEW 88	aq = [inhoms]	×
89	else:
NEW 90	aq = inhoms	×
91
92	# convert aquifer names to aquifer objects
93	for i, iaq in enumerate(aq):	3✔
94	if isinstance(iaq, str):	3✔
NEW 95	aq[i] = self.model.aq.inhomdict[iaq]	×
96
97	p = None	3✔
98	# from_lay, to_lay = [int(i) for i in layers_from_name]
99	if name[:3] == "kaq":	3✔
100	p = aq[0].kaq[from_lay : to_lay + 1]	3✔
101	param = name[:3]	3✔
102	elif name[:3] == "Saq":	3✔
103	p = aq[0].Saq[from_lay : to_lay + 1]	3✔
104	param = name[:3]	3✔
105	elif name[0] == "c":	3✔
106	p = aq[0].c[from_lay : to_lay + 1]	3✔
107	param = name[0]	3✔
NEW 108	elif name[:3] == "Sll":	×
NEW 109	p = aq[0].Sll[from_lay : to_lay + 1]	×
NEW 110	param = name[:3]	×
NEW 111	elif name[0:8] == "kzoverkh":	×
NEW 112	p = aq[0].kzoverkh[from_lay : to_lay + 1]	×
NEW 113	param = name[:8]	×
114
115	if p is None: # no parameter set	3✔
116	print("parameter name not recognized or no parameter ref supplied")	×
117	return	×
118
119	# set all other aquifer parameter references to same array
120	if len(aq) > 1:	3✔
NEW 121	for iaq in aq:	×
NEW 122	setattr(iaq, param, p)	×
123
124	if inhoms is None:	3✔
125	pname = name	3✔
126	else:
NEW 127	pname = f"{name}_{'_'.join([iaq.name for iaq in aq])}"	×
128	self.parameters.loc[pname] = {	3✔
129	"layers": layers,
130	"optimal": initial,
131	"std": None,
132	"perc_std": None,
133	"pmin": pmin,
134	"pmax": pmax,
135	"initial": initial,
136	"inhoms": aq if inhoms is not None else None,
137	"parray": p[:],
138	}
139
140	def set_parameter_by_reference(	3✔
141	self, name=None, parameter=None, initial=0, pmin=-np.inf, pmax=np.inf
142	):
143	"""Set parameter to be optimized.
144
145	Parameters
146	----------
147	name : str
148	parameter name
149	parameter : np.array
150	array reference containing the parameter to be optimized. must be
151	specified as reference, i.e. w.rc[0:]
152	initial : float, optional
153	initial value for the parameter (the default is 0)
154	pmin : float, optional
155	lower bound for parameter value (the default is -np.inf)
156	pmax : float, optional
157	upper bound for paramater value (the default is np.inf)
158	"""
159	assert isinstance(name, str), "Error: name must be string"	3✔
160	if parameter is not None:	3✔
161	assert isinstance(	3✔
162	parameter, np.ndarray
163	), "Error: parameter needs to be numpy array"
164	p = parameter	3✔
165	self.parameters.loc[name] = {	3✔
166	"optimal": initial,
167	"std": None,
168	"perc_std": None,
169	"pmin": pmin,
170	"pmax": pmax,
171	"initial": initial,
172	"parray": p[:],
173	}
174
175	def series(self, name, x, y, layer, t, h, weights=None):	3✔
176	"""Method to add observations to Calibration object.
177
178	Parameters
179	----------
180	name : str
181	name of series
182	x : float
183	x-coordinate
184	y : float
185	y-coordinate
186	layer : int
187	layer number, 0-indexed
188	t : np.array
189	array containing timestamps of timeseries
190	h : np.array
191	array containing timeseries values, i.e. head observations
192	"""
193	s = Series(x, y, layer, t, h, weights=weights)	3✔
194	self.seriesdict[name] = s	3✔
195
196	def seriesinwell(self, name, element, t, h):	3✔
197	"""Method to add observations to Calibration object.
198
199	Parameters
200	----------
201	name : str
202	name of series
203	element: element object with headinside function
204	t : np.array
205	array containing timestamps of timeseries
206	h : np.array
207	array containing timeseries values, i.e. head observations
208	"""
209	e = SeriesInWell(element, t, h)	3✔
210	self.seriesinwelldict[name] = e	3✔
211
212	def residuals(self, p, printdot=False, weighted=True, layers=None, series=None):	3✔
213	"""Method to calculate residuals given certain parameters.
214
215	Parameters
216	----------
217	p : np.array
218	array containing parameter values
219	printdot : bool, optional
220	print dot for each function call
221
222	Returns
223	-------
224	np.array
225	array containing all residuals
226	"""
227	if printdot:	3✔
228	print(".", end="")	3✔
229	# set the values of the variables
230
231	if printdot == 7:	3✔
232	print(p)	×
233
234	if layers is None:	3✔
235	layers = range(self.model.aq.naq)	3✔
236
237	for i, k in enumerate(self.parameters.index):	3✔
238	# [:] needed to do set value in array
239	self.parameters.loc[k, "parray"][:] = p[i]	3✔
240
241	self.model.solve(silent=True)	3✔
242
243	rv = np.empty(0)	3✔
244	cal_series = self.seriesdict.keys() if series is None else series	3✔
245	for key in cal_series:	3✔
246	s = self.seriesdict[key]	3✔
247	if s.layer not in layers:	3✔
248	continue	×
249	h = self.model.head(s.x, s.y, s.t, layers=s.layer)	3✔
250	w = s.weights if ((s.weights is not None) and weighted) else np.ones_like(h)	3✔
251	rv = np.append(rv, (s.h - h) * w)	3✔
252	for key in self.seriesinwelldict:	3✔
253	s = self.seriesinwelldict[key]	3✔
254	h = s.element.headinside(s.t)[0]	3✔
255	rv = np.append(rv, s.h - h)	3✔
256	return rv	3✔
257
258	def residuals_lmfit(self, lmfitparams, printdot=False):	3✔
259	vals = lmfitparams.valuesdict()	3✔
260	p = np.array([vals[k] for k in self.parameters.index])	3✔
261	# p = np.array([vals[k] for k in vals])
262	return self.residuals(p, printdot)	3✔
263
264	def fit_least_squares(self, report=True, diff_step=1e-4, xtol=1e-8, method="lm"):	3✔
265	self.fitresult = least_squares(	3✔
266	self.residuals,
267	self.parameters.initial.values,
268	args=(True,),
269	bounds=(self.parameters.pmin.values, self.parameters.pmax.values),
270	method=method,
271	diff_step=diff_step,
272	xtol=xtol,
273	x_scale="jac",
274	)
275	print("", flush=True)	3✔
276	# Call residuals to specify optimal values for model
277	res = self.residuals(self.fitresult.x)	3✔
278	for ipar in self.parameters.index:	3✔
279	self.parameters.loc[ipar, "optimal"] = self.parameters.loc[ipar, "parray"][	3✔
280	0
281	]
282	nparam = len(self.fitresult.x)	3✔
283	H = self.fitresult.jac.T @ self.fitresult.jac	3✔
284	sigsq = np.var(res, ddof=nparam)	3✔
285	self.covmat = np.linalg.inv(H) * sigsq	3✔
286	self.sig = np.sqrt(np.diag(self.covmat))	3✔
287	D = np.diag(1 / self.sig)	3✔
288	self.cormat = D @ self.covmat @ D	3✔
289	self.parameters["std"] = self.sig	3✔
290	self.parameters["perc_std"] = self.sig / self.parameters["optimal"] * 100	3✔
291	if report:	3✔
292	print(self.parameters)	3✔
293	print(self.sig)	3✔
294	print(self.covmat)	3✔
295	print(self.cormat)	3✔
296
297	def fit_lmfit(self, report=False, printdot=True):	3✔
298	import lmfit	3✔
299
300	self.lmfitparams = lmfit.Parameters()	3✔
301	for name in self.parameters.index:	3✔
302	p = self.parameters.loc[name]	3✔
303	self.lmfitparams.add(name, value=p["initial"], min=p["pmin"], max=p["pmax"])	3✔
304	fit_kws = {"epsfcn": 1e-4}	3✔
305	self.fitresult = lmfit.minimize(	3✔
306	self.residuals_lmfit,
307	self.lmfitparams,
308	method="leastsq",
309	kws={"printdot": printdot},
310	**fit_kws,
311	)
312	print("", flush=True)	3✔
313	print(self.fitresult.message)	3✔
314	if self.fitresult.success:	3✔
315	for name in self.parameters.index:	3✔
316	self.parameters.loc[name, "optimal"] = (	3✔
317	self.fitresult.params.valuesdict()[name]
318	)
319	if hasattr(self.fitresult, "covar"):	3✔
320	self.parameters["std"] = np.sqrt(np.diag(self.fitresult.covar))	3✔
321	self.parameters["perc_std"] = (	3✔
322	100 * self.parameters["std"] / np.abs(self.parameters["optimal"])
323	)
324	else:
325	self.parameters["std"] = np.nan	×
326	self.parameters["perc_std"] = np.nan	×
327	if report:	3✔
328	print(lmfit.fit_report(self.fitresult))	3✔
329
330	def fit(self, report=False, printdot=True):	3✔
331	# current default fitting routine is lmfit
332	# return self.fit_least_squares(report) # does not support bounds by default
333	return self.fit_lmfit(report, printdot)	3✔
334
335	def rmse(self, weighted=True, layers=None):	3✔
336	"""Calculate root-mean-squared-error.
337
338	Returns
339	-------
340	float
341	return rmse value
342	"""
343	r = self.residuals(	3✔
344	self.parameters["optimal"].values, weighted=weighted, layers=layers
345	)
346	return np.sqrt(np.mean(r**2))	3✔
347
348	def topview(self, ax=None, layers=None, labels=True):	3✔
349	"""Plot topview of model with calibration points.
350
351	Parameters
352	----------
353	ax : matplotlib.axes.Axes, optional
354	axes to plot on (the default is None, which creates a new figure)
355	"""
NEW 356	if ax is None:	×
NEW 357	_, ax = plt.subplots()	×
358	# self.model.plots.topview(ax=ax)
NEW 359	for key, s in self.seriesdict.items():	×
NEW 360	if layers is None or s.layer in layers:	×
NEW 361	ax.plot(s.x, s.y, "ko")	×
NEW 362	if labels:	×
NEW 363	ax.text(s.x, s.y, key, ha="left", va="bottom")	×
NEW 364	return ax	×
365
366	def xsection(self, ax=None, labels=True):	3✔
367	"""Plot cross-section of model with calibration points.
368
369	Parameters
370	----------
371	ax : matplotlib.axes.Axes, optional
372	axes to plot on (the default is None, which creates a new figure)
373	"""
NEW 374	if ax is None:	×
NEW 375	_, ax = plt.subplots()	×
NEW 376	for key, s in self.seriesdict.items():	×
NEW 377	aq = self.model.aq.find_aquifer_data(s.x, s.y)	×
NEW 378	ztop = aq.z[0]	×
NEW 379	zpb_top = aq.zaqtop[s.layer]	×
NEW 380	zpb_bot = aq.zaqbot[s.layer]	×
NEW 381	ax.plot([s.x, s.x], [zpb_top, zpb_bot], c="k", ls="dotted")	×
NEW 382	ax.plot([s.x, s.x], [ztop + 1, zpb_top], c="k", ls="solid", lw=1.0)	×
NEW 383	if labels:	×
NEW 384	ax.text(s.x, ztop + 1, key, ha="left", va="bottom", rotation=45)	×
NEW 385	return ax	×
386
387
388	class Series:	3✔
389	def __init__(self, x, y, layer, t, h, weights=None):	3✔
390	self.x = x	3✔
391	self.y = y	3✔
392	self.layer = layer	3✔
393	self.t = t	3✔
394	self.h = h	3✔
395	self.weights = weights	3✔
396
397
398	class SeriesInWell:	3✔
399	def __init__(self, element, t, h):	3✔
400	self.element = element	3✔
401	self.t = t	3✔
402	self.h = h	3✔

mbakker7 / ttim / 12312733157

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous