<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 <IMG SRC="https
Step1:  2. Importing the time series
In this codeblock the time series are imported and collected into a Pastas Project. This Project is a Pastas Class that contains methods that aim to ease the work when dealing with multiple time series. The following time series are imported
Step2:  3/4/5. Creating and optimizing the Time Series Model
For each time series of groundwater head observations a TFN model is constructed with the following model components
Step3:  Make plots for publication
In the next codeblocks the Figures used in the Pastas paper are created. The first codeblock sets the matplotlib parameters to obtain publication-quality figures. The following figures are created
Step4:  Figure of the drawdown estimated for each observations well
Step5:  Example figure of a TFN model
Step6:  Figure of the pumping rate of the well field
<ASSISTANT_TASK:>
Python Code:
# Import the packages
import pandas as pd
import pastas as ps
import numpy as np
import os
import matplotlib.pyplot as plt
%matplotlib inline
# This notebook has been developed using Pastas version 0.9.9 and Python 3.7
print("Pastas version: {}".format(ps.__version__))
print("Pandas version: {}".format(pd.__version__))
print("Numpy version: {}".format(np.__version__))
print("Python version: {}".format(os.sys.version))
ps.set_log_level("ERROR")
Explanation: <IMG SRC="https://raw.githubusercontent.com/pastas/pastas/master/doc/_static/logo.png" WIDTH=250 ALIGN="right">
Example 2: Analysis of groundwater monitoring networks using Pastas
This notebook is supplementary material to the following paper submitted to Groundwater:
R.A. Collenteur, M. Bakker, R. Calje, S. Klop, F. Schaars, (2019) Pastas: open source software for the analysis of groundwater time series, Manuscript under review.
In this second example, it is demonstrated how scripts can be used to analyze a large number of time series. Consider a pumping well field surrounded by a number of observations wells. The pumping wells are screened in the middle aquifer of a three-aquifer system. The objective is to estimate the drawdown caused by the groundwater pumping in each observation well.
1. Import the packages
End of explanation
# Start an Pastas Project
pr = ps.Project(name="Example2")
# Load a metadata-file with xy-coordinates from the groundwater heads
metadata = pd.read_csv("data/metadata_heads.csv", index_col=0)
# Add the groundwater head observations to the database
for fname in os.listdir("./data/heads/"):
    fname = os.path.join("./data/heads/", fname)
    obs = pd.read_csv(fname, parse_dates=True, index_col=0, squeeze=True)
    meta = metadata.loc[obs.name].to_dict()
    pr.add_series(obs, kind="oseries", metadata=meta)
# Load a metadata-file with xy-coordinates from the explanatory variables
metadata = pd.read_csv("data/metadata_stresses.csv", index_col=0)
# Import the precipitation time series
rain = pd.read_csv("data/rain.csv", parse_dates=True, index_col=0, squeeze=True)
pr.add_series(rain, kind="prec", metadata=metadata.loc[rain.name].to_dict())
# Import the evaporation time series
evap = pd.read_csv("data/evap.csv", parse_dates=True, index_col=0, squeeze=True)
pr.add_series(evap, kind="evap", metadata=metadata.loc["Bilt"].to_dict())
# Import the well abstraction time series 
well = pd.read_csv("data/well.csv", parse_dates=True, index_col=0, squeeze=True)
pr.add_series(well, kind="well", metadata=metadata.loc[well.name].to_dict())
# Use the plotting method of the Pastas Project to plot the stresses
pr.plots.stresses(["prec", "evap", "well"], cols=1, figsize=(10,5), sharex=True);
plt.xlim("1960", "2018");
Explanation: 2. Importing the time series
In this codeblock the time series are imported and collected into a Pastas Project. This Project is a Pastas Class that contains methods that aim to ease the work when dealing with multiple time series. The following time series are imported:
44 time series with head observations [m] from the monitoring network;
precipitation [m/d] from KNMI station Oudenbosch;
potential evaporation [m/d] from KNMI station de Bilt;
Total pumping rate [m3/d] from well field Seppe.
End of explanation
# Create folder to save the model figures
mlpath = "models"
if not os.path.exists(mlpath):
    os.mkdir(mlpath)
# Choose the calibration period
tmin = "1970"
tmax = "2017-09"
num = 0
for oseries in pr.oseries.index:
    # Create a Model for each time series and add a StressModel2 for the recharge
    ml = pr.add_model(oseries)
    # Add the RechargeModel to simulate the effect of rainfall and evaporation
    rm = ps.RechargeModel(rain, evap, rfunc=ps.Gamma, name="recharge")
    ml.add_stressmodel(rm)
    
    # Add a StressModel to simulate the effect of the groundwater extractions
    sm = ps.StressModel(well, rfunc=ps.Hantush, name="well", settings="well", up=False)
    ml.add_stressmodel(sm)
    # Since we are dealing with different measurement frequencies for the head, 
    # change the initial parameter for the noise model a little.
#     alpha_init = ml.oseries.series.index.to_series().diff() / pd.Timedelta(1, 'd')
#     ml.set_initial("noise_alpha", alpha_init.mean())
    ml.set_initial("noise_alpha", 3)
    
    # Estimate the model parameters
    ml.solve(tmin=tmin, tmax=tmax, report=False, solver=ps.LmfitSolve)
    # Check if the estimated effect of the groundwater extraction is significant.
    # If not, delete the stressmodel and calibrate the model again.
    gain, stderr = ml.parameters.loc["well_A", ["optimal", "stderr"]]
    if 1.96 * stderr > -gain:
        num += 1
        ml.del_stressmodel("well")
        ml.solve(tmin=tmin, tmax=tmax, report=False)
        
    # Plot the results and store the plot
    ml.plots.results()
    path = os.path.join(mlpath, ml.name + ".png")
    plt.savefig(path, bbox_inches="tight")
    plt.close()
    
print("The number of models where the well is dropped from the model is:", str(num))
Explanation: 3/4/5. Creating and optimizing the Time Series Model
For each time series of groundwater head observations a TFN model is constructed with the following model components:
- A Constant
- A NoiseModel
- A RechargeModel object to simulate the effect of recharge
- A StressModel object to simulate the effect of groundwater extraction
Calibrating all models can take a couple of minutes!!
End of explanation
params = {
   'axes.labelsize': 18,
   'axes.axisbelow': True,
   'font.size': 16,
   'font.family': 'serif',
   'legend.fontsize': 16,
   'xtick.labelsize': 16,
   'ytick.labelsize': 16,
   'text.usetex': False,
   'figure.figsize': [8.2, 5],
   'lines.linewidth' : 2,
   }
plt.rcParams.update(params)
# Save figures or not
savefig = True
figpath = "figures"
if not os.path.exists(figpath):
    os.mkdir(figpath)
Explanation: Make plots for publication
In the next codeblocks the Figures used in the Pastas paper are created. The first codeblock sets the matplotlib parameters to obtain publication-quality figures. The following figures are created:
Figure of the drawdown estimated for each observations well; 
Figure of the decomposition of the different contributions;
Figure of the pumping rate of the well field.
End of explanation
try:
    from timml import ModelMaq, Well
    plot_timml = True
    # Values from REGIS II v2.2 (Site id B49F0240)
    z = [9, -25, -83, -115, -190] # Reference to NAP
    kv = np.array([1e-3, 5e-3]) # Min-Max of Vertical hydraulic conductivity for both leaky layer
    D1 = z[0]-z[1] # Estimated thickness of leaky layer
    c1 = D1/kv # Estimated resistance
    D2 = z[2] - z[3]
    c2 = D2 / kv
    kh1 = np.array([1e0, 2.5e0]) # Min-Max of Horizontal hydraulic conductivity for aquifer 1
    kh2 = np.array([1e1, 2.5e1]) # Min-Max of Horizontal hydraulic conductivity for aquifer 2
    mlm = ModelMaq(kaq=[kh1.mean(), 35], z=z, c=[c1.max(), c2.mean()], \
                  topboundary='semi', hstar=0)
    w = Well(mlm, 0, 0, 34791, layers=1)
    mlm.solve()
    x = np.linspace(100, 5000, 100)
    h = mlm.headalongline(x, 0)
except:
    plot_timml = False
# Get the parameters and distances to plot
param_A = pr.get_parameters(["well_A"])*well.loc["2007":].mean()
#param_A.fillna(0.0, inplace=True)
distances = pr.get_distances(oseries=param_A.index, kind="well").squeeze()
# Select model per aquifer
shallow = pr.oseries.z.loc[(pr.oseries.z<96)].index
aquifer = pr.oseries.z.loc[(pr.oseries.z<186) & (pr.oseries.z>96)].index
deep = pr.oseries.z.loc[pr.oseries.z>186].index
# Make the plot
fig = plt.figure(figsize=(8,5))
plt.grid(zorder=-10)
display_error_bars = True
if display_error_bars:
    std = pr.get_parameters(["well_A"], param_value="stderr")*well.loc["2007":].mean()
    plt.errorbar(distances[shallow], param_A[shallow], yerr=1.96*std[shallow], linestyle="", 
                 elinewidth=2, marker="", markersize=10, capsize=4)
    plt.errorbar(distances[aquifer], param_A[aquifer], yerr=1.96*std[aquifer], linestyle="", 
                 elinewidth=2, marker="", capsize=4)   
#     plt.errorbar(distances[deep], param_A[deep], yerr=1.96*std[deep], linestyle="", 
#                  elinewidth=2, marker="", capsize=4)
    
plt.scatter(distances[shallow], param_A[shallow], marker="^", s=80)
plt.scatter(distances[aquifer], param_A[aquifer], marker="s", s=80)
#plt.scatter(distances[deep], param_A[deep], marker="v", s=100)
# Plot two-layer TimML model for comparison
if plot_timml:
    plt.plot(x, h[0], color="C0", linestyle="--" )
    plt.plot(x, h[1], color="C1", linestyle="--" )
    legend = ["TimML L1", "TimML L2", "aquifer 1", "aquifer 2"]
else:
    legend = ["aquifer 1", "aquifer 2"]
plt.ylabel("steady drawdown (m)")
plt.xlabel("radial distance from the center of the well field (m)")
plt.xlim(0, 4501)
plt.ylim(-10,0)
plt.legend(legend, loc=4)
if savefig:
    path = os.path.join(figpath, "drawdown.eps")
    plt.savefig(path, bbox_inches="tight", dpi=300)
Explanation: Figure of the drawdown estimated for each observations well
End of explanation
# Select a model to plot
ml = pr.models["B49F0232_5"]
# Create the figure
[ax1, ax2, ax3] = ml.plots.decomposition(split=False, figsize=(7,6), ytick_base=1, tmin="1985")
plt.xticks(rotation=0)
ax1.set_yticks([2, 0, -2])
ax1.set_ylabel("head (m)")
ax1.legend().set_visible(False)
ax3.set_yticks([-4, -6])
ax2.set_ylabel("contributions (m)                ") # Little trick to get the label right
ax3.set_xlabel("year")
ax3.set_title("pumping well")
if savefig:
    path = os.path.join(figpath, ml.name + ".eps")
    plt.savefig(path, bbox_inches="tight", dpi=300)
Explanation: Example figure of a TFN model
End of explanation
fig, ax = plt.subplots(1,1, figsize=(8,2.5), sharex=True)
ax.plot(well, color="k")
ax.set_ylabel("pumping rate\n[m$^3$/day]")
ax.set_xlabel("year")
ax.set_xlim("1951", "2018")
if savefig:
    path = os.path.join(figpath, "extraction.eps")
    plt.savefig(path, bbox_inches="tight", dpi=300)
Explanation: Figure of the pumping rate of the well field
End of explanation
<END_TASK>