Current working stack for setting up PEST interface



In [ ]:

    
%matplotlib inline
import os
import shutil
import platform
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import flopy
import pyemu

load the existing model just to see what is going on



In [ ]:

    
nam_file = "freyberg.nam"
org_model_ws = "freyberg_sfr_update"
m = flopy.modflow.Modflow.load(nam_file,model_ws=org_model_ws,check=False)



In [ ]:

    
m.dis.nper #stress periods



In [ ]:

    
m.dis.botm[2].plot()
m.export("shape.shp")

change the working dir and write a new copy of model files to keep the others safe



In [ ]:

    
m.change_model_ws("temp",reset_external=True)
m.external_path = '.'
m.exe_name="mfnwt"
m.write_input()
EXE_DIR = os.path.join("..","bin")
if "window" in platform.platform().lower():
    EXE_DIR = os.path.join(EXE_DIR,"win")
elif "darwin" in platform.platform().lower():
    EXE_DIR = os.path.join(EXE_DIR,"mac")
else:
    EXE_DIR = os.path.join(EXE_DIR,"linux")

[shutil.copy2(os.path.join(EXE_DIR,f),os.path.join('temp',f)) for f in os.listdir(EXE_DIR)]

m.run_model()

build up args for which properties and outputs we want to include in the interface



In [ ]:

    
props,hds = [],[]
for k in range(m.nlay):
    props.append(["upw.hk",k])
    props.append(["upw.vka",k])
    for kper in range(m.nper):
        hds.append([kper,k])
    
props.append(["rch.rech",0])
props.append(["rch.rech",1])
props

Here's where the cool stuff happens: this call will build a pest interace entirely using mulitplier parameters - a mixture of uniform (constant) and grid-scale parameters for all props listed above, plus multipliers for all wells in all stress periods and SFR components.

For observations, we will get the MODFLOW list file budget values, sfr outputs and headsave file array values (all of them!). All observations will be given simulated outputs - this is very useful for error checking...

The interface package will be written to "template" and includes a python forward run script - "template" is the whole package....



In [ ]:

    
ph = pyemu.helpers.PstFromFlopyModel(nam_file,org_model_ws="temp",
                                     new_model_ws="template",grid_props=props,pp_space=3,
                                    const_props=props,spatial_list_props=[["wel.flux",2]],
                                    sfr_pars=True,hds_kperk=hds,
                                    remove_existing=True, sfr_obs=True,
                                    model_exe_name="mfnwt",build_prior=False)
EXE_DIR = os.path.join("..","bin")
if "window" in platform.platform().lower():
    EXE_DIR = os.path.join(EXE_DIR,"win")
elif "darwin" in platform.platform().lower():
    EXE_DIR = os.path.join(EXE_DIR,"mac")
else:
    EXE_DIR = os.path.join(EXE_DIR,"linux")

[shutil.copy2(os.path.join(EXE_DIR,f),os.path.join('temp',f)) for f in os.listdir(EXE_DIR)]

let's inspect what just happened...



In [ ]:

    
pst = ph.pst



In [ ]:

    
pst.npar,pst.nobs

WAT!



In [ ]:

    
pst.parameter_data.head()



In [ ]:

    
pst.observation_data.tail() # the last few observations

write parameter and observtion summary LaTeX tables



In [ ]:

    
pst.write_par_summary_table()
pst.write_obs_summary_table()

plot the prior distribution of the parameter implied by the parameter bounds (assuming bounds represent 95% credible limits)



In [ ]:

    
figs = pst.plot(kind="prior",unique_only=True,echo=False)

but we can do better! We can use geostatistics to build a prior parameter covariance matrix for spatially distributed parameters...



In [ ]:

    
cov = ph.build_prior(fmt="none")



In [ ]:

    
plt.imshow(np.ma.masked_where(cov.x==0,cov.x))
plt.show()

Let's run pestpp once to get residuals...which should be nearly zero since all observations are set to current simulated outputs and all parameters are multipliers - this is a good check!



In [ ]:

    
pyemu.helpers.run("pestpp-glm freyberg.pst",cwd="template")

Reload the pst instance to update the residuals



In [ ]:

    
pst = pyemu.Pst(os.path.join("template","freyberg.pst"))



In [ ]:

    
pst.phi_components

all residuals are nearly zero - this is good!



In [ ]:

    
pst.plot(kind='phi_pie')

Let's do some Monte Carlo!

Generate a parameter ensemble and run it in parallel



In [ ]:

    
pe = pyemu.ParameterEnsemble.from_gaussian_draw(pst=ph.pst,num_reals=100,
                                               cov=cov)



In [ ]:

    
# save to a csv file
pe.to_csv(os.path.join("template","sweep_in.csv"))



In [ ]:

    
# run with sweep using 20 workers
pyemu.helpers.start_workers("template","pestpp-swp","freyberg.pst",num_workers=20,
                          master_dir="sweep_master")

more eye candy using the plot helpers...

Every 50th parameter



In [ ]:

    
pyemu.plot_utils.ensemble_helper(pe.iloc[:,::50],facecolor='b',
                                 deter_vals=pst.parameter_data.parval1.to_dict(),
                                 filename=None) # you can also pass pdf filename

load the output csv file....



In [ ]:

    
df = pd.read_csv(os.path.join("sweep_master","sweep_out.csv"))
df.columns = df.columns.map(str.lower)
df = df.loc[:,[o for o in pst.obs_names if "hds_00" in o ]]

and plot...

every 20th observation



In [ ]:

    
pyemu.plot_utils.ensemble_helper(df.iloc[:,:20],filename=None)



In [ ]: