From 22a93e91f1543a403145f00dcf961e785172d52b Mon Sep 17 00:00:00 2001
From: Mike Taves <mwtoews@gmail.com>
Date: Wed, 13 Nov 2024 15:44:01 +1300
Subject: [PATCH] Use ruff format

---
 docs/pestutilslib/docx_to_md.py               | 144 +++---
 examples/clear_output_all.py                  |  15 +-
 ...oring_lowlevel_pypestutils_functions.ipynb | 259 ++++++++---
 examples/exploring_pypestutils_helpers.ipynb  | 207 +++++----
 examples/ppu_helpers.py                       |  23 +-
 .../pypestutils_pstfrom_quadtree_mf6.ipynb    | 268 +++++++----
 .../pypestutils_pstfrom_structured_mf6.ipynb  | 328 ++++++++-----
 examples/run_notebooks.py                     |  29 +-
 examples/swap.py                              | 245 ++++++----
 pypestutils/ctypes_declarations.py            |   1 +
 pypestutils/data.py                           |   4 +-
 pypestutils/enum.py                           |   1 +
 pypestutils/finder.py                         |   1 +
 pypestutils/helpers.py                        | 438 ++++++++++--------
 pypestutils/logger.py                         |   1 +
 pypestutils/pestutilslib.py                   |   1 +
 tests/common.py                               |   1 +
 tests/test_crash.py                           |   1 +
 tests/test_ctypes_declarations.py             |   1 +
 tests/test_data.py                            |   1 +
 tests/test_drivers.py                         |   1 +
 tests/test_enum.py                            |   1 +
 tests/test_helpers.py                         |  43 +-
 tests/test_notebooks.py                       |  21 +-
 tests/test_pestutilslib.py                    |  82 ++--
 25 files changed, 1287 insertions(+), 830 deletions(-)

diff --git a/docs/pestutilslib/docx_to_md.py b/docs/pestutilslib/docx_to_md.py
index 5944137..e8887b8 100644
--- a/docs/pestutilslib/docx_to_md.py
+++ b/docs/pestutilslib/docx_to_md.py
@@ -1,30 +1,30 @@
 import sys
 import os
 
+
 def processFile(inFile, outFile):
-    mdFile = open(inFile, 'r')
+    mdFile = open(inFile, "r")
     toc = []
-    levels = [0,0,0,0,0 ]
-    newFile = open(outFile, 'w')
+    levels = [0, 0, 0, 0, 0]
+    newFile = open(outFile, "w")
     tempFile = []
     tocLoc = 0
     partOfToc = False
-    
+
     for line in mdFile:
-        if partOfToc and line != '\n':
+        if partOfToc and line != "\n":
             continue
         else:
             partOfToc = False
-        if 'Table of Contents' in line:
+        if "Table of Contents" in line:
             tocLoc = len(tempFile) + 1
             partOfToc = True
             line += "\n"
-        elif line[0] == '#':
+        elif line[0] == "#":
             secId = buildToc(line, toc, levels)
-            line = addSectionTag(cleanLine(line), secId) + '\n'
+            line = addSectionTag(cleanLine(line), secId) + "\n"
         tempFile.append(line)
 
-
     for line in toc:
         tempFile.insert(tocLoc, line)
         tocLoc += 1
@@ -36,58 +36,75 @@ def processFile(inFile, outFile):
     mdFile.close()
     newFile.close()
 
+
 def addSectionTag(line, secId):
-    startIndex = line.find(' ')
-    line = line[:startIndex + 1] + '<a id=\'' + secId + '\' />' + line[startIndex + 1:]
+    startIndex = line.find(" ")
+    line = line[: startIndex + 1] + "<a id='" + secId + "' />" + line[startIndex + 1 :]
     return line
 
+
 def buildToc(line, toc, levels):
     line = cleanLine(line)
-    secId = 's'
-    if line[:4] == '####':
-        
-        #raise UserWarning('Header levels greater than 3 not supported')
+    secId = "s"
+    if line[:4] == "####":
+        # raise UserWarning('Header levels greater than 3 not supported')
         levels[4] += 1
-        secId += str(levels[1]) + '-' + str(levels[2]) + '-' + str(levels[3]) + '-' + str(levels[4])
-        toc.append('            - [' + line[5:] + '](#' + secId + ')\n')
-    elif line[:3] == '###':
+        secId += (
+            str(levels[1])
+            + "-"
+            + str(levels[2])
+            + "-"
+            + str(levels[3])
+            + "-"
+            + str(levels[4])
+        )
+        toc.append("            - [" + line[5:] + "](#" + secId + ")\n")
+    elif line[:3] == "###":
         levels[3] += 1
-        secId += str(levels[1]) + '-' + str(levels[2]) + '-' + str(levels[3])
-        toc.append('        - [' + line[4:] + '](#' + secId + ')\n')
-    elif line[:2] == '##':
+        secId += str(levels[1]) + "-" + str(levels[2]) + "-" + str(levels[3])
+        toc.append("        - [" + line[4:] + "](#" + secId + ")\n")
+    elif line[:2] == "##":
         levels[2] += 1
         levels[3] = 0
-        secId += str(levels[1]) + '-' + str(levels[2])
-        toc.append('    - [' + line[3:] + '](#' + secId + ')\n')
-    elif line[:1] == '#':
+        secId += str(levels[1]) + "-" + str(levels[2])
+        toc.append("    - [" + line[3:] + "](#" + secId + ")\n")
+    elif line[:1] == "#":
         levels[1] += 1
         levels[3] = levels[2] = 0
         secId += str(levels[1])
-        toc.append('- [' + line[2:] + '](#' + secId + ')\n')
+        toc.append("- [" + line[2:] + "](#" + secId + ")\n")
     return secId
 
+
 def cleanLine(text):
     text = stripNewline(text)
     text = removeAnchors(text)
     return text
 
+
 def stripNewline(text):
-    return text.replace('\n', '')
+    return text.replace("\n", "")
+
 
 def removeAnchors(text):
-    while ('<' in text and '>' in text):
-        leftTag = text.index('<')
-        rightTag = text.index('>')
-        text = text[0:leftTag] + text[rightTag + 1:]
+    while "<" in text and ">" in text:
+        leftTag = text.index("<")
+        rightTag = text.index(">")
+        text = text[0:leftTag] + text[rightTag + 1 :]
     return text
 
-def clean(docx_file,inFile,outFile,run_pandoc=True):
+
+def clean(docx_file, inFile, outFile, run_pandoc=True):
     if run_pandoc:
-        os.system("pandoc -t gfm --wrap=none --extract-media . -o file.md {0} --mathjax".format(docx_file))
-    num_str = [str(i) for i in range(1,11)]
-    lines = open(inFile,'r').readlines()
+        os.system(
+            "pandoc -t gfm --wrap=none --extract-media . -o file.md {0} --mathjax".format(
+                docx_file
+            )
+        )
+    num_str = [str(i) for i in range(1, 11)]
+    lines = open(inFile, "r").readlines()
 
-    #notoc_lines = []
+    # notoc_lines = []
     i = 0
     # while i < len(lines):
     #     line = lines[i]
@@ -104,31 +121,40 @@ def clean(docx_file,inFile,outFile,run_pandoc=True):
     # notoc_lines = []
     i = 0
     while i < len(lines):
-        #if lines[i].strip().startswith("----"):
+        # if lines[i].strip().startswith("----"):
         #    lines[i-1] = "## " + lines[i-1]
         #    lines[i] = "\n"
         if "<img" in lines[i]:
-            lines[i] = lines[i].replace("#","")
+            lines[i] = lines[i].replace("#", "")
         if "<p>" in lines[i]:
-            lines[i] = lines[i].replace("<p>","").replace("</p>","<br>").replace("#","~")
+            lines[i] = (
+                lines[i].replace("<p>", "").replace("</p>", "<br>").replace("#", "~")
+            )
         if "blockquote" in lines[i]:
-            lines[i] = lines[i].replace("<blockquote>","").replace("</blockquote>","")
+            lines[i] = lines[i].replace("<blockquote>", "").replace("</blockquote>", "")
         if lines[i].strip().startswith("$") and not "bmatrix" in lines[i].lower():
             label = lines[i].split()[-1]
-            eq_str = lines[i].replace("$$","$").split('$')[1]
-            eq_str = r"{0}".format(eq_str).replace("\\\\","\\").replace(" \\ "," ").replace("\\_","_")
-            math_str_pre = "<img src=\"https://latex.codecogs.com/svg.latex?\Large&space;{0}".format(eq_str)
-            math_str_post = "\" title=\"\Large {0}\" />  {1}  <br>".format(eq_str,label)
+            eq_str = lines[i].replace("$$", "$").split("$")[1]
+            eq_str = (
+                r"{0}".format(eq_str)
+                .replace("\\\\", "\\")
+                .replace(" \\ ", " ")
+                .replace("\\_", "_")
+            )
+            math_str_pre = '<img src="https://latex.codecogs.com/svg.latex?\Large&space;{0}'.format(
+                eq_str
+            )
+            math_str_post = '" title="\Large {0}" />  {1}  <br>'.format(eq_str, label)
             lines[i] = math_str_pre + " " + math_str_post
-        if lines[i].strip().startswith("<table>"): # and "pcf" in lines[i].lower():
-            lines[i] = "<div style=\"text-align: left\">" + lines[i] + "</div>"
+        if lines[i].strip().startswith("<table>"):  # and "pcf" in lines[i].lower():
+            lines[i] = '<div style="text-align: left">' + lines[i] + "</div>"
             if "comment" in lines[i].lower():
-                lines[i] = lines[i].replace("~","#")
+                lines[i] = lines[i].replace("~", "#")
         elif lines[i].strip().startswith("[") and "]" in lines[i]:
-            raw = lines[i].strip().split(']')
-            rraw = " ".join(raw[0].split()[:-1])+"]"
+            raw = lines[i].strip().split("]")
+            rraw = " ".join(raw[0].split()[:-1]) + "]"
             new_line = rraw + "".join(raw[1:]) + "\n"
-            print(lines[i],new_line)
+            print(lines[i], new_line)
             lines[i] = new_line
         # elif "bmatrix" in lines[i].lower():
         #     eq_str = lines[i]
@@ -146,20 +172,18 @@ def clean(docx_file,inFile,outFile,run_pandoc=True):
 
         i += 1
 
-    #lines[0] = "# Table of Contents\n"
-    with open(outFile,'w') as f:
-
+    # lines[0] = "# Table of Contents\n"
+    with open(outFile, "w") as f:
         for line in lines:
-            #if line.lower().strip().startswith("[introduction"):
+            # if line.lower().strip().startswith("[introduction"):
             #    f.write("\n# Table of Contents\n\n")
             f.write(line)
 
 
 if __name__ == "__main__":
-
-    clean("new_function_documentation.docx","file.md","fortran_library_documentation.md",True)
-    
-
-
-
-
+    clean(
+        "new_function_documentation.docx",
+        "file.md",
+        "fortran_library_documentation.md",
+        True,
+    )
diff --git a/examples/clear_output_all.py b/examples/clear_output_all.py
index 05843b5..169d7e8 100644
--- a/examples/clear_output_all.py
+++ b/examples/clear_output_all.py
@@ -1,10 +1,17 @@
 import os
+
 dirs = ["."]
 notebook_count = 0
 for d in dirs:
-    nb_files = [os.path.join(d,f) for f in os.listdir(d) if f.lower().endswith(".ipynb")]
+    nb_files = [
+        os.path.join(d, f) for f in os.listdir(d) if f.lower().endswith(".ipynb")
+    ]
     for nb_file in nb_files:
-        print("clearing",nb_file)
-        os.system("jupyter nbconvert --ClearOutputPreprocessor.enabled=True --ClearMetadataPreprocessor.enabled=True --inplace {0}".format(nb_file))
+        print("clearing", nb_file)
+        os.system(
+            "jupyter nbconvert --ClearOutputPreprocessor.enabled=True --ClearMetadataPreprocessor.enabled=True --inplace {0}".format(
+                nb_file
+            )
+        )
         notebook_count += 1
-print(notebook_count," notebooks cleared")
+print(notebook_count, " notebooks cleared")
diff --git a/examples/exploring_lowlevel_pypestutils_functions.ipynb b/examples/exploring_lowlevel_pypestutils_functions.ipynb
index 7fd2e72..419374f 100644
--- a/examples/exploring_lowlevel_pypestutils_functions.ipynb
+++ b/examples/exploring_lowlevel_pypestutils_functions.ipynb
@@ -46,7 +46,7 @@
     "w_d = \"freyberg_explore_lowlevel\"\n",
     "if os.path.exists(w_d):\n",
     "    shutil.rmtree(w_d)\n",
-    "shutil.copytree(org_d,w_d)"
+    "shutil.copytree(org_d, w_d)"
    ]
   },
   {
@@ -88,6 +88,7 @@
    "outputs": [],
    "source": [
     "from pypestutils.pestutilslib import PestUtilsLib\n",
+    "\n",
     "lib = PestUtilsLib()"
    ]
   },
@@ -98,7 +99,9 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "grid_info = lib.install_mf6_grid_from_file(\"structgrid\",os.path.join(w_d,\"freyberg6.dis.grb\"))\n"
+    "grid_info = lib.install_mf6_grid_from_file(\n",
+    "    \"structgrid\", os.path.join(w_d, \"freyberg6.dis.grb\")\n",
+    ")"
    ]
   },
   {
@@ -154,7 +157,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "easting,northing,elev = lib.get_cell_centres_mf6(\"structgrid\",grid_info[\"ncells\"])"
+    "easting, northing, elev = lib.get_cell_centres_mf6(\"structgrid\", grid_info[\"ncells\"])"
    ]
   },
   {
@@ -164,10 +167,10 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "easting = easting[:nrow*ncol]\n",
-    "northing = northing[:nrow*ncol]\n",
-    "Easting = easting.reshape((nrow,ncol))\n",
-    "Northing = northing.reshape((nrow,ncol))"
+    "easting = easting[: nrow * ncol]\n",
+    "northing = northing[: nrow * ncol]\n",
+    "Easting = easting.reshape((nrow, ncol))\n",
+    "Northing = northing.reshape((nrow, ncol))"
    ]
   },
   {
@@ -189,7 +192,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "hds_file = [os.path.join(w_d,f) for f in os.listdir(w_d) if f.endswith(\".hds\")][0]\n",
+    "hds_file = [os.path.join(w_d, f) for f in os.listdir(w_d) if f.endswith(\".hds\")][0]\n",
     "hds_file"
    ]
   },
@@ -208,7 +211,9 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "depvar_file_info = lib.inquire_modflow_binary_file_specs(hds_file,hds_file+\".out\",31,1)"
+    "depvar_file_info = lib.inquire_modflow_binary_file_specs(\n",
+    "    hds_file, hds_file + \".out\", 31, 1\n",
+    ")"
    ]
   },
   {
@@ -236,7 +241,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "df = pd.read_csv(hds_file+\".out\")\n",
+    "df = pd.read_csv(hds_file + \".out\")\n",
     "df.columns = [c.lower() for c in df.columns]\n",
     "df"
    ]
@@ -256,7 +261,9 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "hdsdf = pd.read_csv(os.path.join(\"freyberg_aux_files\",\"gwlevel_obs.csv\"),parse_dates=[\"datetime\"])\n",
+    "hdsdf = pd.read_csv(\n",
+    "    os.path.join(\"freyberg_aux_files\", \"gwlevel_obs.csv\"), parse_dates=[\"datetime\"]\n",
+    ")\n",
     "hdsdf"
    ]
   },
@@ -268,7 +275,7 @@
    "outputs": [],
    "source": [
     "mtop = m.dis.top.array\n",
-    "mtop[mtop<-100] = np.nan"
+    "mtop[mtop < -100] = np.nan"
    ]
   },
   {
@@ -278,15 +285,15 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "fig,axes = plt.subplots(1,2,figsize=(10,10))\n",
-    "for lay,ax in zip([1,3],axes):\n",
+    "fig, axes = plt.subplots(1, 2, figsize=(10, 10))\n",
+    "for lay, ax in zip([1, 3], axes):\n",
     "    ax.set_aspect(\"equal\")\n",
-    "    kdf = hdsdf.loc[hdsdf.layer==lay,:]\n",
+    "    kdf = hdsdf.loc[hdsdf.layer == lay, :]\n",
     "    assert kdf.shape[0] > 0\n",
-    "    ax.pcolormesh(Easting,Northing,mtop)\n",
-    "    ax.scatter(kdf.x,kdf.y,marker=\"^\",c=\"k\",label=\"gw level loc\")\n",
+    "    ax.pcolormesh(Easting, Northing, mtop)\n",
+    "    ax.scatter(kdf.x, kdf.y, marker=\"^\", c=\"k\", label=\"gw level loc\")\n",
     "    ax.legend(loc=\"upper left\")\n",
-    "    ax.set_title(\"gw level locations in layer {0}\".format(lay),loc=\"left\")\n"
+    "    ax.set_title(\"gw level locations in layer {0}\".format(lay), loc=\"left\")"
    ]
   },
   {
@@ -323,8 +330,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "fac_file = os.path.join(w_d,\"obs_interp_fac.bin\")\n",
-    "bln_file = fac_file.replace(\".bin\",\".bln\")"
+    "fac_file = os.path.join(w_d, \"obs_interp_fac.bin\")\n",
+    "bln_file = fac_file.replace(\".bin\", \".bln\")"
    ]
   },
   {
@@ -334,7 +341,15 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "results = lib.calc_mf6_interp_factors(\"structgrid\",usitedf.x.values,usitedf.y.values,usitedf.layer.values,fac_file,\"binary\",bln_file)\n",
+    "results = lib.calc_mf6_interp_factors(\n",
+    "    \"structgrid\",\n",
+    "    usitedf.x.values,\n",
+    "    usitedf.y.values,\n",
+    "    usitedf.layer.values,\n",
+    "    fac_file,\n",
+    "    \"binary\",\n",
+    "    bln_file,\n",
+    ")\n",
     "results"
    ]
   },
@@ -371,7 +386,17 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "head_results = lib.interp_from_mf6_depvar_file(hds_file,fac_file,\"binary\",depvar_file_info[\"ntime\"],\"head\",1e+10,True,-1.0e+30,usitedf.shape[0])"
+    "head_results = lib.interp_from_mf6_depvar_file(\n",
+    "    hds_file,\n",
+    "    fac_file,\n",
+    "    \"binary\",\n",
+    "    depvar_file_info[\"ntime\"],\n",
+    "    \"head\",\n",
+    "    1e10,\n",
+    "    True,\n",
+    "    -1.0e30,\n",
+    "    usitedf.shape[0],\n",
+    ")"
    ]
   },
   {
@@ -401,8 +426,10 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "start_datetime = pd.to_datetime(\"1-1-2018\") \n",
-    "hdsdf.loc[:,\"time\"] = hdsdf.datetime.apply(lambda x: x  - start_datetime).dt.days # we are losing fractional days..oh well...\n",
+    "start_datetime = pd.to_datetime(\"1-1-2018\")\n",
+    "hdsdf.loc[:, \"time\"] = hdsdf.datetime.apply(\n",
+    "    lambda x: x - start_datetime\n",
+    ").dt.days  # we are losing fractional days..oh well...\n",
     "hdsdf.time"
    ]
   },
@@ -423,10 +450,10 @@
    "source": [
     "usite = hdsdf.site.unique()\n",
     "usite.sort()\n",
-    "usite_dict = {s:c for s,c in zip(usite,np.arange(usite.shape[0],dtype=int))}\n",
-    "hdsdf.loc[:,\"isite\"] = hdsdf.site.apply(lambda x: usite_dict[x])\n",
+    "usite_dict = {s: c for s, c in zip(usite, np.arange(usite.shape[0], dtype=int))}\n",
+    "hdsdf.loc[:, \"isite\"] = hdsdf.site.apply(lambda x: usite_dict[x])\n",
     "hdsdf.isite\n",
-    "hdsdf.sort_values(by=[\"isite\",\"time\"],inplace=True)\n",
+    "hdsdf.sort_values(by=[\"isite\", \"time\"], inplace=True)\n",
     "hdsdf"
    ]
   },
@@ -437,7 +464,17 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "ihead_results = lib.interp_to_obstime(head_results[\"nproctime\"],head_results[\"simtime\"],head_results[\"simstate\"],1.e+10,\"L\",35.0,1.0e+30,hdsdf.isite.values,hdsdf.time.values)"
+    "ihead_results = lib.interp_to_obstime(\n",
+    "    head_results[\"nproctime\"],\n",
+    "    head_results[\"simtime\"],\n",
+    "    head_results[\"simstate\"],\n",
+    "    1.0e10,\n",
+    "    \"L\",\n",
+    "    35.0,\n",
+    "    1.0e30,\n",
+    "    hdsdf.isite.values,\n",
+    "    hdsdf.time.values,\n",
+    ")"
    ]
   },
   {
@@ -447,7 +484,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "hdsdf.loc[:,\"simulated\"] = ihead_results"
+    "hdsdf.loc[:, \"simulated\"] = ihead_results"
    ]
   },
   {
@@ -485,7 +522,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "easting,northing,elev = lib.get_cell_centres_mf6(\"structgrid\",grid_info[\"ncells\"])"
+    "easting, northing, elev = lib.get_cell_centres_mf6(\"structgrid\", grid_info[\"ncells\"])"
    ]
   },
   {
@@ -495,7 +532,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "easting.shape,northing.shape"
+    "easting.shape, northing.shape"
    ]
   },
   {
@@ -513,8 +550,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "easting = easting[:nrow*ncol]\n",
-    "northing = northing[:nrow*ncol]"
+    "easting = easting[: nrow * ncol]\n",
+    "northing = northing[: nrow * ncol]"
    ]
   },
   {
@@ -536,18 +573,18 @@
    "source": [
     "# cell area\n",
     "area = np.ones_like(easting)\n",
-    "#active array\n",
-    "active = m.dis.idomain.array[0,:,:].flatten()\n",
+    "# active array\n",
+    "active = m.dis.idomain.array[0, :, :].flatten()\n",
     "# property mean\n",
     "mean = np.ones_like(easting)\n",
     "# property variance\n",
     "var = np.ones_like(easting)\n",
     "# the variogram range\n",
-    "aa = np.ones_like(easting)*1000\n",
+    "aa = np.ones_like(easting) * 1000\n",
     "# anisotropy\n",
-    "anis = np.ones_like(easting)*5\n",
+    "anis = np.ones_like(easting) * 5\n",
     "# bearing\n",
-    "bearing = (np.ones_like(easting) * 55)"
+    "bearing = np.ones_like(easting) * 55"
    ]
   },
   {
@@ -587,9 +624,23 @@
    "source": [
     "transform = \"none\"\n",
     "variogram_type = \"exp\"\n",
-    "power = 1.0 #unused\n",
+    "power = 1.0  # unused\n",
     "num_reals = 10\n",
-    "reals = lib.fieldgen2d_sva(easting,northing,area,active,mean,var,aa,anis,bearing,transform,variogram_type,power,num_reals)\n",
+    "reals = lib.fieldgen2d_sva(\n",
+    "    easting,\n",
+    "    northing,\n",
+    "    area,\n",
+    "    active,\n",
+    "    mean,\n",
+    "    var,\n",
+    "    aa,\n",
+    "    anis,\n",
+    "    bearing,\n",
+    "    transform,\n",
+    "    variogram_type,\n",
+    "    power,\n",
+    "    num_reals,\n",
+    ")\n",
     "reals.shape"
    ]
   },
@@ -600,7 +651,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "m.dis.top = reals[:,0]\n",
+    "m.dis.top = reals[:, 0]\n",
     "m.dis.top.plot()"
    ]
   },
@@ -622,12 +673,26 @@
     "lib2 = PestUtilsLib()\n",
     "transform = \"none\"\n",
     "variogram_type = \"exp\"\n",
-    "power = 1.0 #unused\n",
+    "power = 1.0  # unused\n",
     "num_reals = 10\n",
     "lib2.initialize_randgen(54321)\n",
-    "reals = lib2.fieldgen2d_sva(easting,northing,area,active,mean,var,aa,anis,bearing,transform,variogram_type,power,num_reals)\n",
+    "reals = lib2.fieldgen2d_sva(\n",
+    "    easting,\n",
+    "    northing,\n",
+    "    area,\n",
+    "    active,\n",
+    "    mean,\n",
+    "    var,\n",
+    "    aa,\n",
+    "    anis,\n",
+    "    bearing,\n",
+    "    transform,\n",
+    "    variogram_type,\n",
+    "    power,\n",
+    "    num_reals,\n",
+    ")\n",
     "print(reals)\n",
-    "plt.imshow(reals[:,0].reshape((nrow,ncol)))"
+    "plt.imshow(reals[:, 0].reshape((nrow, ncol)))"
    ]
   },
   {
@@ -645,7 +710,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "bearing = np.add(np.ones((nrow,ncol)),np.atleast_2d(np.arange(ncol)))"
+    "bearing = np.add(np.ones((nrow, ncol)), np.atleast_2d(np.arange(ncol)))"
    ]
   },
   {
@@ -687,8 +752,22 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "reals = lib.fieldgen2d_sva(easting,northing,area,active,mean,var,aa,anis,bearing,transform,variogram_type,power,num_reals)\n",
-    "r = reals[:,0].reshape((nrow,ncol))\n",
+    "reals = lib.fieldgen2d_sva(\n",
+    "    easting,\n",
+    "    northing,\n",
+    "    area,\n",
+    "    active,\n",
+    "    mean,\n",
+    "    var,\n",
+    "    aa,\n",
+    "    anis,\n",
+    "    bearing,\n",
+    "    transform,\n",
+    "    variogram_type,\n",
+    "    power,\n",
+    "    num_reals,\n",
+    ")\n",
+    "r = reals[:, 0].reshape((nrow, ncol))\n",
     "plt.imshow(r)"
    ]
   },
@@ -707,24 +786,24 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "Easting = easting.reshape((nrow,ncol))\n",
-    "Northing = northing.reshape((nrow,ncol))\n",
-    "ppeasting,ppnorthing = [],[]\n",
+    "Easting = easting.reshape((nrow, ncol))\n",
+    "Northing = northing.reshape((nrow, ncol))\n",
+    "ppeasting, ppnorthing = [], []\n",
     "ppval = []\n",
     "pp_space = 20\n",
-    "ib = m.dis.idomain.array[0,:,:]\n",
-    "half_pp_space = int(pp_space/2)\n",
-    "for i in range(half_pp_space,nrow,pp_space):\n",
-    "    for j in range(half_pp_space,ncol,pp_space):\n",
-    "        if ib[i,j] == 0:\n",
+    "ib = m.dis.idomain.array[0, :, :]\n",
+    "half_pp_space = int(pp_space / 2)\n",
+    "for i in range(half_pp_space, nrow, pp_space):\n",
+    "    for j in range(half_pp_space, ncol, pp_space):\n",
+    "        if ib[i, j] == 0:\n",
     "            continue\n",
-    "        ppeasting.append(Easting[i,j])\n",
-    "        ppnorthing.append(Northing[i,j])\n",
-    "        ppval.append(r[i,j])\n",
+    "        ppeasting.append(Easting[i, j])\n",
+    "        ppnorthing.append(Northing[i, j])\n",
+    "        ppval.append(r[i, j])\n",
     "ppeasting = np.asarray(ppeasting)\n",
     "ppnorthing = np.asarray(ppnorthing)\n",
     "ppval = np.asarray(ppval)\n",
-    "ppeasting.shape,ppnorthing.shape,ppval.shape"
+    "ppeasting.shape, ppnorthing.shape, ppval.shape"
    ]
   },
   {
@@ -742,12 +821,12 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "fig,ax = plt.subplots(1,1)\n",
+    "fig, ax = plt.subplots(1, 1)\n",
     "ax.set_aspect(\"equal\")\n",
-    "ax.pcolormesh(Easting,Northing,r)\n",
+    "ax.pcolormesh(Easting, Northing, r)\n",
     "ax.set_title(\"realization\")\n",
     "\n",
-    "ax.scatter(ppeasting,ppnorthing,marker=\".\",s=50,c='k')"
+    "ax.scatter(ppeasting, ppnorthing, marker=\".\", s=50, c=\"k\")"
    ]
   },
   {
@@ -767,15 +846,32 @@
    "source": [
     "max_pts = 50\n",
     "min_pts = 1\n",
-    "search_dist = 1.e+10\n",
-    "aa_pp = aa * pp_space *10 #?\n",
-    "zone_pp = np.ones_like(ppeasting,dtype=int)\n",
-    "fac_file = os.path.join(w_d,\"factors.bin\")\n",
+    "search_dist = 1.0e10\n",
+    "aa_pp = aa * pp_space * 10  # ?\n",
+    "zone_pp = np.ones_like(ppeasting, dtype=int)\n",
+    "fac_file = os.path.join(w_d, \"factors.bin\")\n",
     "from datetime import datetime\n",
+    "\n",
     "s = datetime.now()\n",
-    "ipts = lib.calc_kriging_factors_2d(ppeasting,ppnorthing,zone_pp,easting,northing,ib.flatten(),\n",
-    "                                   \"exp\",\"ordinary\",aa_pp,anis,bearing,search_dist,max_pts,min_pts,fac_file,\"binary\")\n",
-    "\"total points:\",ipts,\" took:\",(datetime.now() - s).total_seconds()"
+    "ipts = lib.calc_kriging_factors_2d(\n",
+    "    ppeasting,\n",
+    "    ppnorthing,\n",
+    "    zone_pp,\n",
+    "    easting,\n",
+    "    northing,\n",
+    "    ib.flatten(),\n",
+    "    \"exp\",\n",
+    "    \"ordinary\",\n",
+    "    aa_pp,\n",
+    "    anis,\n",
+    "    bearing,\n",
+    "    search_dist,\n",
+    "    max_pts,\n",
+    "    min_pts,\n",
+    "    fac_file,\n",
+    "    \"binary\",\n",
+    ")\n",
+    "\"total points:\", ipts, \" took:\", (datetime.now() - s).total_seconds()"
    ]
   },
   {
@@ -793,7 +889,16 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "result = lib.krige_using_file(fac_file,\"binary\",len(easting),\"ordinary\",\"none\",ppval,np.zeros_like(easting),0)"
+    "result = lib.krige_using_file(\n",
+    "    fac_file,\n",
+    "    \"binary\",\n",
+    "    len(easting),\n",
+    "    \"ordinary\",\n",
+    "    \"none\",\n",
+    "    ppval,\n",
+    "    np.zeros_like(easting),\n",
+    "    0,\n",
+    ")"
    ]
   },
   {
@@ -803,19 +908,19 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "rr = result[\"targval\"].reshape(nrow,ncol)\n",
-    "fig,axes = plt.subplots(1,2)\n",
+    "rr = result[\"targval\"].reshape(nrow, ncol)\n",
+    "fig, axes = plt.subplots(1, 2)\n",
     "ax = axes[0]\n",
     "ax.set_aspect(\"equal\")\n",
     "ax.set_title(\"pp interpolated array\")\n",
-    "ax.pcolormesh(Easting,Northing,rr) #the interpolated array\n",
+    "ax.pcolormesh(Easting, Northing, rr)  # the interpolated array\n",
     "ax = axes[1]\n",
     "ax.set_aspect(\"equal\")\n",
     "ax.set_title(\"pp locs with sampled values\")\n",
-    "id_mask = m.dis.idomain.array[0,:,:].copy().astype(float)\n",
-    "id_mask[id_mask!=0] = np.nan\n",
-    "ax.pcolormesh(Easting,Northing,id_mask)\n",
-    "ax.scatter(ppeasting,ppnorthing,marker=\".\",s=50,c=ppval)\n"
+    "id_mask = m.dis.idomain.array[0, :, :].copy().astype(float)\n",
+    "id_mask[id_mask != 0] = np.nan\n",
+    "ax.pcolormesh(Easting, Northing, id_mask)\n",
+    "ax.scatter(ppeasting, ppnorthing, marker=\".\", s=50, c=ppval)"
    ]
   },
   {
diff --git a/examples/exploring_pypestutils_helpers.ipynb b/examples/exploring_pypestutils_helpers.ipynb
index 5b06963..0fa79dc 100644
--- a/examples/exploring_pypestutils_helpers.ipynb
+++ b/examples/exploring_pypestutils_helpers.ipynb
@@ -44,7 +44,7 @@
     "w_d = \"freyberg_highlevel_helpers\"\n",
     "if os.path.exists(w_d):\n",
     "    shutil.rmtree(w_d)\n",
-    "shutil.copytree(org_d,w_d)"
+    "shutil.copytree(org_d, w_d)"
    ]
   },
   {
@@ -62,8 +62,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "delc = np.loadtxt(os.path.join(w_d,\"freyberg6.dis_delc.txt\")).flatten()\n",
-    "delr = np.loadtxt(os.path.join(w_d,\"freyberg6.dis_delr.txt\")).flatten()"
+    "delc = np.loadtxt(os.path.join(w_d, \"freyberg6.dis_delc.txt\")).flatten()\n",
+    "delr = np.loadtxt(os.path.join(w_d, \"freyberg6.dis_delr.txt\")).flatten()"
    ]
   },
   {
@@ -76,8 +76,8 @@
     "nrow = delc.shape[0]\n",
     "ncol = delr.shape[0]\n",
     "nlay = 3\n",
-    "ib = np.loadtxt(os.path.join(w_d,\"freyberg6.dis_idomain_layer1.txt\"),dtype=int)\n",
-    "ib = ib.flatten().reshape(nrow,ncol)\n",
+    "ib = np.loadtxt(os.path.join(w_d, \"freyberg6.dis_idomain_layer1.txt\"), dtype=int)\n",
+    "ib = ib.flatten().reshape(nrow, ncol)\n",
     "plt.imshow(ib)"
    ]
   },
@@ -108,7 +108,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "grb_fname = os.path.join(w_d,\"freyberg6.dis.grb\")\n",
+    "grb_fname = os.path.join(w_d, \"freyberg6.dis.grb\")\n",
     "os.path.exists(grb_fname)"
    ]
   },
@@ -127,7 +127,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "csv_fname = os.path.join(\"freyberg_aux_files\",\"gwlevel_obs.csv\")\n",
+    "csv_fname = os.path.join(\"freyberg_aux_files\", \"gwlevel_obs.csv\")\n",
     "assert os.path.exists(csv_fname)\n",
     "obsdf = pd.read_csv(csv_fname)\n",
     "obsdf"
@@ -148,11 +148,11 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "depvar_fname = os.path.join(w_d,\"freyberg6_freyberg.hds\")\n",
-    "model_type = 31 #structured modflow-6 model type\n",
-    "start_datetime = \"1-1-2018\" # when the simulation starts\n",
-    "depvar_ftype = 1# modflow-6 binary file format: 1 for states, 2 for budgets\n",
-    "depvar_name = \"head\"# the variable name in the output file to process/extract"
+    "depvar_fname = os.path.join(w_d, \"freyberg6_freyberg.hds\")\n",
+    "model_type = 31  # structured modflow-6 model type\n",
+    "start_datetime = \"1-1-2018\"  # when the simulation starts\n",
+    "depvar_ftype = 1  # modflow-6 binary file format: 1 for states, 2 for budgets\n",
+    "depvar_name = \"head\"  # the variable name in the output file to process/extract"
    ]
   },
   {
@@ -162,7 +162,9 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "results = helpers.mod2obs_mf6(grb_fname,depvar_fname,csv_fname,model_type,start_datetime,depvar_ftype)"
+    "results = helpers.mod2obs_mf6(\n",
+    "    grb_fname, depvar_fname, csv_fname, model_type, start_datetime, depvar_ftype\n",
+    ")"
    ]
   },
   {
@@ -185,22 +187,31 @@
     "adf = results[\"all_results\"]\n",
     "idf = results[\"interpolated_results\"]\n",
     "for site in adf.columns:\n",
-    "    aadf = adf.loc[:,site]\n",
-    "    aadf.loc[aadf.values >1e29] = np.nan\n",
-    "    \n",
-    "    iidf = idf.loc[idf.site==site,:]\n",
-    "    iidf.loc[iidf.simulated>1e29,\"simulated\"] = np.nan\n",
-    "    fig,ax = plt.subplots(1,1,figsize=(8,5))\n",
-    "    ax.plot(aadf.index,aadf.values,\"0.5\",lw=0.5,marker=\"o\",mfc='none',label=\"all sim times\")\n",
-    "    ax.scatter(iidf.datetime,iidf.obsval,marker=\"^\",c=\"r\",label=\"observed\")\n",
-    "    ax.scatter(iidf.datetime,iidf.simulated,marker=\"^\",c=\"0.5\",label=\"interp to obs\")\n",
-    "    ax.set_title(site,loc=\"left\")\n",
+    "    aadf = adf.loc[:, site]\n",
+    "    aadf.loc[aadf.values > 1e29] = np.nan\n",
+    "\n",
+    "    iidf = idf.loc[idf.site == site, :]\n",
+    "    iidf.loc[iidf.simulated > 1e29, \"simulated\"] = np.nan\n",
+    "    fig, ax = plt.subplots(1, 1, figsize=(8, 5))\n",
+    "    ax.plot(\n",
+    "        aadf.index,\n",
+    "        aadf.values,\n",
+    "        \"0.5\",\n",
+    "        lw=0.5,\n",
+    "        marker=\"o\",\n",
+    "        mfc=\"none\",\n",
+    "        label=\"all sim times\",\n",
+    "    )\n",
+    "    ax.scatter(iidf.datetime, iidf.obsval, marker=\"^\", c=\"r\", label=\"observed\")\n",
+    "    ax.scatter(\n",
+    "        iidf.datetime, iidf.simulated, marker=\"^\", c=\"0.5\", label=\"interp to obs\"\n",
+    "    )\n",
+    "    ax.set_title(site, loc=\"left\")\n",
     "    ax.legend(loc=\"upper left\")\n",
     "    ax.grid()\n",
     "    plt.tight_layout()\n",
     "    plt.show()\n",
-    "    plt.close(fig)\n",
-    "        "
+    "    plt.close(fig)"
    ]
   },
   {
@@ -218,7 +229,11 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "[f for f in os.listdir(w_d) if f.endswith(\".csv\") and os.path.split(depvar_fname)[1] in f]"
+    "[\n",
+    "    f\n",
+    "    for f in os.listdir(w_d)\n",
+    "    if f.endswith(\".csv\") and os.path.split(depvar_fname)[1] in f\n",
+    "]"
    ]
   },
   {
@@ -252,7 +267,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "sr = helpers.SpatialReference(delc=delc,delr=delr,rotation=-55,xul=0,yul=0)\n",
+    "sr = helpers.SpatialReference(delc=delc, delr=delr, rotation=-55, xul=0, yul=0)\n",
     "sr.rotation"
    ]
   },
@@ -263,7 +278,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "plt.pcolormesh(sr.xcentergrid,sr.ycentergrid,ib)"
+    "plt.pcolormesh(sr.xcentergrid, sr.ycentergrid, ib)"
    ]
   },
   {
@@ -281,7 +296,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "gridspec_fname = os.path.join(w_d,\"grid.spc\")\n",
+    "gridspec_fname = os.path.join(w_d, \"grid.spc\")\n",
     "sr.write_gridspec(gridspec_fname)"
    ]
   },
@@ -300,9 +315,9 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "x,y = sr.xcentergrid.flatten(),sr.ycentergrid.flatten()\n",
-    "grid_df = pd.DataFrame({\"x\":x,\"y\":y,\"layer\":1})\n",
-    "csv_fname = os.path.join(w_d,\"grid.csv\")\n",
+    "x, y = sr.xcentergrid.flatten(), sr.ycentergrid.flatten()\n",
+    "grid_df = pd.DataFrame({\"x\": x, \"y\": y, \"layer\": 1})\n",
+    "csv_fname = os.path.join(w_d, \"grid.csv\")\n",
     "grid_df.to_csv(csv_fname)"
    ]
   },
@@ -322,6 +337,7 @@
    "outputs": [],
    "source": [
     "from pypestutils.pestutilslib import PestUtilsLib\n",
+    "\n",
     "lib = PestUtilsLib()"
    ]
   },
@@ -340,7 +356,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "grb_fname = os.path.join(w_d,\"freyberg6.dis.grb\")\n",
+    "grb_fname = os.path.join(w_d, \"freyberg6.dis.grb\")\n",
     "os.path.exists(grb_fname)"
    ]
   },
@@ -352,7 +368,7 @@
    "outputs": [],
    "source": [
     "grid_info = helpers.get_grid_info_from_mf6_grb(grb_fname)\n",
-    "grid_info['x'].shape"
+    "grid_info[\"x\"].shape"
    ]
   },
   {
@@ -362,8 +378,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "xx = grid_info['x'].reshape((nlay,nrow,ncol))\n",
-    "yy = grid_info['y'].reshape((nlay,nrow,ncol))"
+    "xx = grid_info[\"x\"].reshape((nlay, nrow, ncol))\n",
+    "yy = grid_info[\"y\"].reshape((nlay, nrow, ncol))"
    ]
   },
   {
@@ -373,7 +389,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "xx.shape\n"
+    "xx.shape"
    ]
   },
   {
@@ -383,9 +399,9 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "fig,ax = plt.subplots(1,1)\n",
+    "fig, ax = plt.subplots(1, 1)\n",
     "ax.set_aspect(\"equal\")\n",
-    "ax.pcolormesh(xx[0,:,:],yy[0,:,:],ib)"
+    "ax.pcolormesh(xx[0, :, :], yy[0, :, :], ib)"
    ]
   },
   {
@@ -404,7 +420,7 @@
    "outputs": [],
    "source": [
     "grid_info = helpers.get_2d_grid_info_from_mf6_grb(grb_fname)\n",
-    "grid_info['x'].shape"
+    "grid_info[\"x\"].shape"
    ]
   },
   {
@@ -503,9 +519,11 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "ppdf = helpers.get_2d_pp_info_structured_grid(pp_space=10,gridinfo_fname=gridspec_fname)\n",
-    "plt.pcolormesh(sr.xcentergrid,sr.ycentergrid,ib)\n",
-    "plt.scatter(ppdf.x,ppdf.y)"
+    "ppdf = helpers.get_2d_pp_info_structured_grid(\n",
+    "    pp_space=10, gridinfo_fname=gridspec_fname\n",
+    ")\n",
+    "plt.pcolormesh(sr.xcentergrid, sr.ycentergrid, ib)\n",
+    "plt.scatter(ppdf.x, ppdf.y)"
    ]
   },
   {
@@ -515,9 +533,11 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "ppdf = helpers.get_2d_pp_info_structured_grid(pp_space=10,gridinfo_fname=grb_fname)\n",
-    "plt.pcolormesh(grid_info_2d['x'].reshape((nrow,ncol)),grid_info_2d['y'].reshape((nrow,ncol)),ib)\n",
-    "plt.scatter(ppdf.x,ppdf.y)"
+    "ppdf = helpers.get_2d_pp_info_structured_grid(pp_space=10, gridinfo_fname=grb_fname)\n",
+    "plt.pcolormesh(\n",
+    "    grid_info_2d[\"x\"].reshape((nrow, ncol)), grid_info_2d[\"y\"].reshape((nrow, ncol)), ib\n",
+    ")\n",
+    "plt.scatter(ppdf.x, ppdf.y)"
    ]
   },
   {
@@ -527,9 +547,11 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "ppdf = helpers.get_2d_pp_info_structured_grid(10,gridspec_fname,array_dict={\"zone\":ib})\n",
-    "plt.pcolormesh(sr.xcentergrid,sr.ycentergrid,ib)\n",
-    "plt.scatter(ppdf.x,ppdf.y)"
+    "ppdf = helpers.get_2d_pp_info_structured_grid(\n",
+    "    10, gridspec_fname, array_dict={\"zone\": ib}\n",
+    ")\n",
+    "plt.pcolormesh(sr.xcentergrid, sr.ycentergrid, ib)\n",
+    "plt.scatter(ppdf.x, ppdf.y)"
    ]
   },
   {
@@ -549,8 +571,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "reals = helpers.generate_2d_grid_realizations(gridspec_fname,num_reals=10)\n",
-    "plt.pcolormesh(sr.xcentergrid,sr.ycentergrid,reals[0,:,:])"
+    "reals = helpers.generate_2d_grid_realizations(gridspec_fname, num_reals=10)\n",
+    "plt.pcolormesh(sr.xcentergrid, sr.ycentergrid, reals[0, :, :])"
    ]
   },
   {
@@ -560,8 +582,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "reals = helpers.generate_2d_grid_realizations(grb_fname,num_reals=10)\n",
-    "plt.pcolormesh(sr.xcentergrid,sr.ycentergrid,reals[0,:,:])"
+    "reals = helpers.generate_2d_grid_realizations(grb_fname, num_reals=10)\n",
+    "plt.pcolormesh(sr.xcentergrid, sr.ycentergrid, reals[0, :, :])"
    ]
   },
   {
@@ -579,8 +601,15 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "reals = helpers.generate_2d_grid_realizations(grid_df,num_reals=10)\n",
-    "plt.pcolormesh(sr.xcentergrid,sr.ycentergrid,reals[0,:,].reshape((nrow,ncol)))"
+    "reals = helpers.generate_2d_grid_realizations(grid_df, num_reals=10)\n",
+    "plt.pcolormesh(\n",
+    "    sr.xcentergrid,\n",
+    "    sr.ycentergrid,\n",
+    "    reals[\n",
+    "        0,\n",
+    "        :,\n",
+    "    ].reshape((nrow, ncol)),\n",
+    ")"
    ]
   },
   {
@@ -598,8 +627,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "reals = helpers.generate_2d_grid_realizations(grb_fname,num_reals=10,zone_array=ib)\n",
-    "plt.pcolormesh(sr.xcentergrid,sr.ycentergrid,reals[0,:,:])"
+    "reals = helpers.generate_2d_grid_realizations(grb_fname, num_reals=10, zone_array=ib)\n",
+    "plt.pcolormesh(sr.xcentergrid, sr.ycentergrid, reals[0, :, :])"
    ]
   },
   {
@@ -617,9 +646,9 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "bearing = np.add(np.ones((nrow,ncol)),np.atleast_2d(np.arange(ncol)))\n",
+    "bearing = np.add(np.ones((nrow, ncol)), np.atleast_2d(np.arange(ncol)))\n",
     "plt.imshow(bearing)\n",
-    "#bearing.min(),bearing.max()"
+    "# bearing.min(),bearing.max()"
    ]
   },
   {
@@ -629,8 +658,15 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "reals = helpers.generate_2d_grid_realizations(gridspec_fname,num_reals=10,zone_array=ib,variobearing=bearing,varioaniso=10,variorange=1000)\n",
-    "plt.pcolormesh(sr.xcentergrid,sr.ycentergrid,reals[0,:,:])"
+    "reals = helpers.generate_2d_grid_realizations(\n",
+    "    gridspec_fname,\n",
+    "    num_reals=10,\n",
+    "    zone_array=ib,\n",
+    "    variobearing=bearing,\n",
+    "    varioaniso=10,\n",
+    "    variorange=1000,\n",
+    ")\n",
+    "plt.pcolormesh(sr.xcentergrid, sr.ycentergrid, reals[0, :, :])"
    ]
   },
   {
@@ -648,11 +684,11 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "s = 10**(np.sin(np.linspace(0,np.pi*2,nrow)))\n",
-    "#plt.plot(s) \n",
-    "aniso = np.add(np.ones((nrow,ncol)),np.atleast_2d(s).transpose())\n",
+    "s = 10 ** (np.sin(np.linspace(0, np.pi * 2, nrow)))\n",
+    "# plt.plot(s)\n",
+    "aniso = np.add(np.ones((nrow, ncol)), np.atleast_2d(s).transpose())\n",
     "plt.imshow(aniso)\n",
-    "#aniso.min(),aniso.max()"
+    "# aniso.min(),aniso.max()"
    ]
   },
   {
@@ -662,8 +698,15 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "reals = helpers.generate_2d_grid_realizations(gridspec_fname,num_reals=10,zone_array=ib,variobearing=bearing,varioaniso=aniso,variorange=1000)\n",
-    "plt.pcolormesh(sr.xcentergrid,sr.ycentergrid,reals[0,:,:])"
+    "reals = helpers.generate_2d_grid_realizations(\n",
+    "    gridspec_fname,\n",
+    "    num_reals=10,\n",
+    "    zone_array=ib,\n",
+    "    variobearing=bearing,\n",
+    "    varioaniso=aniso,\n",
+    "    variorange=1000,\n",
+    ")\n",
+    "plt.pcolormesh(sr.xcentergrid, sr.ycentergrid, reals[0, :, :])"
    ]
   },
   {
@@ -681,14 +724,14 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "array_dict={\"zone\":ib,\"value\":reals[0,:,:],\"bearing\":bearing,\"aniso\":aniso}\n",
-    "ppdf = helpers.get_2d_pp_info_structured_grid(10,gridspec_fname,array_dict=array_dict)\n",
-    "fig,axes = plt.subplots(1,2,figsize=(10,10))\n",
-    "axes[0].pcolormesh(sr.xcentergrid,sr.ycentergrid,reals[0])\n",
-    "axes[0].scatter(ppdf.x,ppdf.y,c=\"k\")\n",
-    "axes[1].scatter(ppdf.x,ppdf.y,c=ppdf.value)\n",
+    "array_dict = {\"zone\": ib, \"value\": reals[0, :, :], \"bearing\": bearing, \"aniso\": aniso}\n",
+    "ppdf = helpers.get_2d_pp_info_structured_grid(10, gridspec_fname, array_dict=array_dict)\n",
+    "fig, axes = plt.subplots(1, 2, figsize=(10, 10))\n",
+    "axes[0].pcolormesh(sr.xcentergrid, sr.ycentergrid, reals[0])\n",
+    "axes[0].scatter(ppdf.x, ppdf.y, c=\"k\")\n",
+    "axes[1].scatter(ppdf.x, ppdf.y, c=ppdf.value)\n",
     "for ax in axes:\n",
-    "    ax.set_aspect(\"equal\")\n"
+    "    ax.set_aspect(\"equal\")"
    ]
   },
   {
@@ -726,13 +769,15 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "interp_results = helpers.interpolate_with_sva_pilotpoints_2d(ppdf,grid_df,zone_array=ib)\n",
-    "for tag,arr in interp_results.items():\n",
-    "    fig,ax = plt.subplots(1,1)\n",
-    "    \n",
-    "    cb = ax.pcolormesh(sr.xcentergrid,sr.ycentergrid,arr.reshape((nrow,ncol)))\n",
-    "    plt.colorbar(cb,ax=ax)\n",
-    "    ax.set_title(tag,loc=\"left\")"
+    "interp_results = helpers.interpolate_with_sva_pilotpoints_2d(\n",
+    "    ppdf, grid_df, zone_array=ib\n",
+    ")\n",
+    "for tag, arr in interp_results.items():\n",
+    "    fig, ax = plt.subplots(1, 1)\n",
+    "\n",
+    "    cb = ax.pcolormesh(sr.xcentergrid, sr.ycentergrid, arr.reshape((nrow, ncol)))\n",
+    "    plt.colorbar(cb, ax=ax)\n",
+    "    ax.set_title(tag, loc=\"left\")"
    ]
   },
   {
diff --git a/examples/ppu_helpers.py b/examples/ppu_helpers.py
index a839540..9e4795c 100644
--- a/examples/ppu_helpers.py
+++ b/examples/ppu_helpers.py
@@ -1,33 +1,36 @@
 import os
 import pyemu
 
+
 def setup_pps(d):
     cwd = os.getcwd()
     os.chdir(d)
     apply_pps()
     os.chdir(cwd)
 
+
 def apply_pps():
-    
     from pypestutils import helpers
     import numpy as np
     import pandas as pd
+
     df = pd.read_csv("pp_info.csv")
     gridspec_fname = "org.grb"
-    for model_file,pp_file in zip(df.model_file,df.pp_file):
+    for model_file, pp_file in zip(df.model_file, df.pp_file):
         ppdf = pd.read_csv(pp_file)
-        results = helpers.interpolate_with_sva_pilotpoints_2d(ppdf,gridspec_fname,vartransform="log")
+        results = helpers.interpolate_with_sva_pilotpoints_2d(
+            ppdf, gridspec_fname, vartransform="log"
+        )
         org_arr = np.loadtxt(model_file)
         interp = results["result"]
         interp = interp.reshape(org_arr.shape)
         new_arr = org_arr * interp
         new_arr = new_arr
-        new_arr[new_arr<1.0e-10] = 1.0e-10
-        np.savetxt(model_file,new_arr,fmt="%15.6E")
-        np.savetxt("interp_"+model_file,interp,fmt="%15.6E")
-        np.savetxt("log_"+model_file,np.log10(new_arr),fmt="%15.6E")
-    
+        new_arr[new_arr < 1.0e-10] = 1.0e-10
+        np.savetxt(model_file, new_arr, fmt="%15.6E")
+        np.savetxt("interp_" + model_file, interp, fmt="%15.6E")
+        np.savetxt("log_" + model_file, np.log10(new_arr), fmt="%15.6E")
+
 
 if __name__ == "__main__":
-	setup_pps()
-	
\ No newline at end of file
+    setup_pps()
diff --git a/examples/pypestutils_pstfrom_quadtree_mf6.ipynb b/examples/pypestutils_pstfrom_quadtree_mf6.ipynb
index cd185cb..a3fb36a 100644
--- a/examples/pypestutils_pstfrom_quadtree_mf6.ipynb
+++ b/examples/pypestutils_pstfrom_quadtree_mf6.ipynb
@@ -59,10 +59,10 @@
    "source": [
     "bin_dir = \"quadtree_bin\"\n",
     "if os.path.exists(bin_dir):\n",
-    "   shutil.rmtree(bin_dir)\n",
+    "    shutil.rmtree(bin_dir)\n",
     "os.makedirs(bin_dir)\n",
-    "flopy.utils.get_modflow(bin_dir,downloads_dir=bin_dir)\n",
-    "pyemu.utils.get_pestpp(bin_dir,downloads_dir=bin_dir)"
+    "flopy.utils.get_modflow(bin_dir, downloads_dir=bin_dir)\n",
+    "pyemu.utils.get_pestpp(bin_dir, downloads_dir=bin_dir)"
    ]
   },
   {
@@ -78,7 +78,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "org_model_ws = os.path.join('freyberg_quadtree')\n",
+    "org_model_ws = os.path.join(\"freyberg_quadtree\")\n",
     "os.listdir(org_model_ws)"
    ]
   },
@@ -95,7 +95,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "id_arr = np.loadtxt(os.path.join(org_model_ws,\"freyberg6.disv_idomain_layer3.txt\"))"
+    "id_arr = np.loadtxt(os.path.join(org_model_ws, \"freyberg6.disv_idomain_layer3.txt\"))"
    ]
   },
   {
@@ -114,10 +114,10 @@
     "tmp_model_ws = \"temp_ppu_quadtree\"\n",
     "if os.path.exists(tmp_model_ws):\n",
     "    shutil.rmtree(tmp_model_ws)\n",
-    "shutil.copytree(org_model_ws,tmp_model_ws)\n",
+    "shutil.copytree(org_model_ws, tmp_model_ws)\n",
     "for bin_file in os.listdir(bin_dir):\n",
-    "    shutil.copy2(os.path.join(bin_dir,bin_file),os.path.join(tmp_model_ws,bin_file))\n",
-    "os.listdir(tmp_model_ws)\n"
+    "    shutil.copy2(os.path.join(bin_dir, bin_file), os.path.join(tmp_model_ws, bin_file))\n",
+    "os.listdir(tmp_model_ws)"
    ]
   },
   {
@@ -137,7 +137,9 @@
    "source": [
     "sim = flopy.mf6.MFSimulation.load(sim_ws=tmp_model_ws)\n",
     "m = sim.get_model(\"freyberg6\")\n",
-    "m.modelgrid.set_coord_info(angrot=-55) #for some reason flopy isnt picking up angrot from the disv options...\n"
+    "m.modelgrid.set_coord_info(\n",
+    "    angrot=-55\n",
+    ")  # for some reason flopy isnt picking up angrot from the disv options..."
    ]
   },
   {
@@ -172,10 +174,14 @@
    "outputs": [],
    "source": [
     "template_ws = \"freyberg6_unstruct_template\"\n",
-    "pf = pyemu.utils.PstFrom(original_d=tmp_model_ws, new_d=template_ws,\n",
-    "                 remove_existing=True,\n",
-    "                 longnames=True,\n",
-    "                 zero_based=False,start_datetime=\"1-1-2018\")"
+    "pf = pyemu.utils.PstFrom(\n",
+    "    original_d=tmp_model_ws,\n",
+    "    new_d=template_ws,\n",
+    "    remove_existing=True,\n",
+    "    longnames=True,\n",
+    "    zero_based=False,\n",
+    "    start_datetime=\"1-1-2018\",\n",
+    ")"
    ]
   },
   {
@@ -193,7 +199,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "hdsdf = pd.read_csv(os.path.join(tmp_model_ws,\"heads.csv\"),index_col=0)\n",
+    "hdsdf = pd.read_csv(os.path.join(tmp_model_ws, \"heads.csv\"), index_col=0)\n",
     "hdsdf"
    ]
   },
@@ -210,8 +216,13 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "hds_df = pf.add_observations(\"heads.csv\",insfile=\"heads.csv.ins\",index_cols=\"time\",\n",
-    "                    use_cols=list(hdsdf.columns.values),prefix=\"hds\",)\n",
+    "hds_df = pf.add_observations(\n",
+    "    \"heads.csv\",\n",
+    "    insfile=\"heads.csv.ins\",\n",
+    "    index_cols=\"time\",\n",
+    "    use_cols=list(hdsdf.columns.values),\n",
+    "    prefix=\"hds\",\n",
+    ")\n",
     "hds_df"
    ]
   },
@@ -247,7 +258,12 @@
    "outputs": [],
    "source": [
     "df = pd.read_csv(os.path.join(tmp_model_ws, \"sfr.csv\"), index_col=0)\n",
-    "sfr_df = pf.add_observations(\"sfr.csv\", insfile=\"sfr.csv.ins\", index_cols=\"time\", use_cols=list(df.columns.values))\n",
+    "sfr_df = pf.add_observations(\n",
+    "    \"sfr.csv\",\n",
+    "    insfile=\"sfr.csv.ins\",\n",
+    "    index_cols=\"time\",\n",
+    "    use_cols=list(df.columns.values),\n",
+    ")\n",
     "sfr_df"
    ]
   },
@@ -279,7 +295,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "ppdf = pd.read_csv(os.path.join(\"freyberg_aux_files\",\"pp_info.csv\"),index_col=0)\n",
+    "ppdf = pd.read_csv(os.path.join(\"freyberg_aux_files\", \"pp_info.csv\"), index_col=0)\n",
     "ppdf"
    ]
   },
@@ -289,13 +305,13 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "fig,ax = plt.subplots(1,1)\n",
+    "fig, ax = plt.subplots(1, 1)\n",
     "ax.set_aspect(\"equal\")\n",
-    "#ax.pcolormesh(sr.xcentergrid,sr.ycentergrid,top_arr)\n",
-    "#ax = m.dis.top.plot()\n",
+    "# ax.pcolormesh(sr.xcentergrid,sr.ycentergrid,top_arr)\n",
+    "# ax = m.dis.top.plot()\n",
     "\n",
-    "ax.scatter(m.modelgrid.xcellcenters,m.modelgrid.ycellcenters,alpha=0.5,color=\"0.5\")\n",
-    "ax.scatter(ppdf.x,ppdf.y,marker=\".\",c=\"r\")\n",
+    "ax.scatter(m.modelgrid.xcellcenters, m.modelgrid.ycellcenters, alpha=0.5, color=\"0.5\")\n",
+    "ax.scatter(ppdf.x, ppdf.y, marker=\".\", c=\"r\")\n",
     "plt.tight_layout()"
    ]
   },
@@ -305,7 +321,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "ppdf.loc[:,\"value\"] = 3.0  # HK in layer 1"
+    "ppdf.loc[:, \"value\"] = 3.0  # HK in layer 1"
    ]
   },
   {
@@ -322,14 +338,14 @@
    "outputs": [],
    "source": [
     "id_mask = id_arr.copy()\n",
-    "id_mask[id_mask!=0] = np.nan\n",
-    "fig,axes = plt.subplots(1,3)\n",
-    "for ax,attr in zip(axes,[\"aniso\",\"corrlen\",\"bearing\"]):\n",
-    "    #ax.pcolormesh(sr.xcentergrid,sr.ycentergrid,id_mask)\n",
+    "id_mask[id_mask != 0] = np.nan\n",
+    "fig, axes = plt.subplots(1, 3)\n",
+    "for ax, attr in zip(axes, [\"aniso\", \"corrlen\", \"bearing\"]):\n",
+    "    # ax.pcolormesh(sr.xcentergrid,sr.ycentergrid,id_mask)\n",
     "    m.dis.top.plot(axes=[ax])\n",
-    "    ax.scatter(ppdf.x,ppdf.y,c=ppdf.loc[:,attr])\n",
+    "    ax.scatter(ppdf.x, ppdf.y, c=ppdf.loc[:, attr])\n",
     "    ax.set_aspect(\"equal\")\n",
-    "    ax.set_title(attr,loc=\"left\")\n",
+    "    ax.set_title(attr, loc=\"left\")\n",
     "    ax.set_xticks([])\n",
     "    ax.set_yticks([])\n",
     "plt.tight_layout()"
@@ -348,8 +364,10 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "ppdf.loc[:,\"value\"] = np.random.normal(1.0,0.25,ppdf.shape[0])\n",
-    "results = ppu.interpolate_with_sva_pilotpoints_2d(ppdf,os.path.join(pf.new_d,\"freyberg6.disv.grb\"))"
+    "ppdf.loc[:, \"value\"] = np.random.normal(1.0, 0.25, ppdf.shape[0])\n",
+    "results = ppu.interpolate_with_sva_pilotpoints_2d(\n",
+    "    ppdf, os.path.join(pf.new_d, \"freyberg6.disv.grb\")\n",
+    ")"
    ]
   },
   {
@@ -358,16 +376,16 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "for tag,arr in results.items():\n",
-    "    fig,ax = plt.subplots(1,1)\n",
+    "for tag, arr in results.items():\n",
+    "    fig, ax = plt.subplots(1, 1)\n",
     "    ax.set_aspect(\"equal\")\n",
-    "    ax.set_title(tag,loc=\"left\")\n",
+    "    ax.set_title(tag, loc=\"left\")\n",
     "    m.dis.top = arr\n",
     "    ax = m.dis.top.plot(colorbar=True)\n",
-    "    \n",
-    "    #plt.colorbar(cb)\n",
+    "\n",
+    "    # plt.colorbar(cb)\n",
     "    plt.show()\n",
-    "    plt.close(fig)\n"
+    "    plt.close(fig)"
    ]
   },
   {
@@ -399,7 +417,9 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "hk_arr_files = [f for f in os.listdir(tmp_model_ws) if \"npf_k_\" in f and f.endswith(\".txt\")]\n",
+    "hk_arr_files = [\n",
+    "    f for f in os.listdir(tmp_model_ws) if \"npf_k_\" in f and f.endswith(\".txt\")\n",
+    "]\n",
     "hk_arr_files"
    ]
   },
@@ -420,9 +440,9 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "ppdf.loc[:,\"value\"] = 1.0\n",
+    "ppdf.loc[:, \"value\"] = 1.0\n",
     "pp_v = pyemu.geostats.ExpVario(contribution=1.0, a=2000)\n",
-    "pp_gs = pyemu.geostats.GeoStruct(variograms=pp_v,transform=\"log\")"
+    "pp_gs = pyemu.geostats.GeoStruct(variograms=pp_v, transform=\"log\")"
    ]
   },
   {
@@ -438,29 +458,52 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "pp_files,mod_files = [],[]\n",
+    "pp_files, mod_files = [], []\n",
     "for hk_arr_file in hk_arr_files:\n",
     "    layer = int(hk_arr_file.split(\".\")[1].split(\"layer\")[1])\n",
-    "    base = hk_arr_file.split('.')[1].replace(\"_\",\"\")+\"_attr:\"\n",
-    "    pf.add_parameters(filenames=hk_arr_file,par_type=\"grid\",\n",
-    "                       par_name_base=base+\"gr\",pargp=base+\"gr\",zone_array=np.atleast_2d(ib).transpose(),\n",
-    "                       upper_bound=2.,lower_bound=0.5,ult_ubound=100,ult_lbound=0.01)\n",
+    "    base = hk_arr_file.split(\".\")[1].replace(\"_\", \"\") + \"_attr:\"\n",
+    "    pf.add_parameters(\n",
+    "        filenames=hk_arr_file,\n",
+    "        par_type=\"grid\",\n",
+    "        par_name_base=base + \"gr\",\n",
+    "        pargp=base + \"gr\",\n",
+    "        zone_array=np.atleast_2d(ib).transpose(),\n",
+    "        upper_bound=2.0,\n",
+    "        lower_bound=0.5,\n",
+    "        ult_ubound=100,\n",
+    "        ult_lbound=0.01,\n",
+    "    )\n",
     "    pppdf = ppdf.copy()\n",
-    "    \n",
-    "    pppdf.loc[:,\"name\"] = [n for n in pppdf.ppname.values]\n",
-    "    pppdf.loc[:,\"ppname\"] = pppdf.name.values\n",
-    "    pp_file = os.path.join(pf.new_d,base+\"pp.csv\")\n",
-    "    pppdf.to_csv(pp_file,index=False)\n",
+    "\n",
+    "    pppdf.loc[:, \"name\"] = [n for n in pppdf.ppname.values]\n",
+    "    pppdf.loc[:, \"ppname\"] = pppdf.name.values\n",
+    "    pp_file = os.path.join(pf.new_d, base + \"pp.csv\")\n",
+    "    pppdf.to_csv(pp_file, index=False)\n",
     "    pp_files.append(os.path.split(pp_file)[1])\n",
     "    mod_files.append(hk_arr_file)\n",
-    "    df = pf.add_parameters(os.path.split(pp_file)[1],par_type=\"grid\",index_cols={\"ppname\":\"ppname\",\"x\":\"x\",\"y\":\"y\"},\n",
-    "        use_cols=[\"value\",\"bearing\",\"aniso\",\"corrlen\"],\n",
-    "        par_name_base=[base+\"pp\",base+\"bearing\",base+\"aniso\",base+\"corrlen\"],\n",
-    "        pargp=[base+\"pp\",base+\"bearing\",base+\"aniso\",base+\"corrlen\"],\n",
-    "        upper_bound=[20,pppdf.bearing.max()*1.1,pppdf.aniso.max()*1.1,pppdf.corrlen.max()*1.1],\n",
-    "        lower_bound=[.05,pppdf.bearing.min()*.9,pppdf.aniso.min()*.9,pppdf.corrlen.min()*.9],\n",
-    "        par_style=\"direct\",transform=\"log\",geostruct=pp_gs)\n",
-    "    "
+    "    df = pf.add_parameters(\n",
+    "        os.path.split(pp_file)[1],\n",
+    "        par_type=\"grid\",\n",
+    "        index_cols={\"ppname\": \"ppname\", \"x\": \"x\", \"y\": \"y\"},\n",
+    "        use_cols=[\"value\", \"bearing\", \"aniso\", \"corrlen\"],\n",
+    "        par_name_base=[base + \"pp\", base + \"bearing\", base + \"aniso\", base + \"corrlen\"],\n",
+    "        pargp=[base + \"pp\", base + \"bearing\", base + \"aniso\", base + \"corrlen\"],\n",
+    "        upper_bound=[\n",
+    "            20,\n",
+    "            pppdf.bearing.max() * 1.1,\n",
+    "            pppdf.aniso.max() * 1.1,\n",
+    "            pppdf.corrlen.max() * 1.1,\n",
+    "        ],\n",
+    "        lower_bound=[\n",
+    "            0.05,\n",
+    "            pppdf.bearing.min() * 0.9,\n",
+    "            pppdf.aniso.min() * 0.9,\n",
+    "            pppdf.corrlen.min() * 0.9,\n",
+    "        ],\n",
+    "        par_style=\"direct\",\n",
+    "        transform=\"log\",\n",
+    "        geostruct=pp_gs,\n",
+    "    )"
    ]
   },
   {
@@ -486,10 +529,9 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "with open(os.path.join(template_ws,tpl_files[0]),'r') as f:\n",
+    "with open(os.path.join(template_ws, tpl_files[0]), \"r\") as f:\n",
     "    for _ in range(4):\n",
-    "        print(f.readline().strip())\n",
-    "        "
+    "        print(f.readline().strip())"
    ]
   },
   {
@@ -507,7 +549,9 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "shutil.copy2(os.path.join(pf.new_d,\"freyberg6.disv.grb\"),os.path.join(pf.new_d,\"org.grb\"))"
+    "shutil.copy2(\n",
+    "    os.path.join(pf.new_d, \"freyberg6.disv.grb\"), os.path.join(pf.new_d, \"org.grb\")\n",
+    ")"
    ]
   },
   {
@@ -516,7 +560,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "_ = [print(line.rstrip()) for line in open(\"ppu_helpers.py\",'r').readlines()]"
+    "_ = [print(line.rstrip()) for line in open(\"ppu_helpers.py\", \"r\").readlines()]"
    ]
   },
   {
@@ -533,8 +577,8 @@
    "outputs": [],
    "source": [
     "# save the pp file and model file info\n",
-    "df = pd.DataFrame({\"model_file\":mod_files,\"pp_file\":pp_files})\n",
-    "df.to_csv(os.path.join(pf.new_d,\"pp_info.csv\"))"
+    "df = pd.DataFrame({\"model_file\": mod_files, \"pp_file\": pp_files})\n",
+    "df.to_csv(os.path.join(pf.new_d, \"pp_info.csv\"))"
    ]
   },
   {
@@ -544,6 +588,7 @@
    "outputs": [],
    "source": [
     "import ppu_helpers\n",
+    "\n",
     "ppu_helpers.setup_pps(pf.new_d)"
    ]
   },
@@ -553,7 +598,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "pf.add_py_function(\"ppu_helpers.py\",\"apply_pps()\",is_pre_cmd=True)"
+    "pf.add_py_function(\"ppu_helpers.py\", \"apply_pps()\", is_pre_cmd=True)"
    ]
   },
   {
@@ -569,7 +614,9 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "interp_files = [f for f in os.listdir(pf.new_d) if f.endswith(\".txt\") and f.startswith(\"interp\")]\n",
+    "interp_files = [\n",
+    "    f for f in os.listdir(pf.new_d) if f.endswith(\".txt\") and f.startswith(\"interp\")\n",
+    "]\n",
     "assert len(interp_files) == 3"
    ]
   },
@@ -580,8 +627,8 @@
    "outputs": [],
    "source": [
     "for interp_file in interp_files:\n",
-    "    base = interp_file.replace(\"_\",\"\").replace(\".\",\"\").replace(\"txt\",\"\")\n",
-    "    pf.add_observations(interp_file,prefix=base,obsgp=base)"
+    "    base = interp_file.replace(\"_\", \"\").replace(\".\", \"\").replace(\"txt\", \"\")\n",
+    "    pf.add_observations(interp_file, prefix=base, obsgp=base)"
    ]
   },
   {
@@ -590,9 +637,10 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "\n",
-    "log_files = [f for f in os.listdir(pf.new_d) if f.endswith(\".txt\") and f.startswith(\"log_\")]\n",
-    "assert len(log_files) == 3\n"
+    "log_files = [\n",
+    "    f for f in os.listdir(pf.new_d) if f.endswith(\".txt\") and f.startswith(\"log_\")\n",
+    "]\n",
+    "assert len(log_files) == 3"
    ]
   },
   {
@@ -602,8 +650,8 @@
    "outputs": [],
    "source": [
     "for log_file in log_files:\n",
-    "    base = log_file.replace(\"_\",\"\").replace(\".\",\"\").replace(\"txt\",\"\")\n",
-    "    pf.add_observations(log_file,prefix=base,obsgp=base)"
+    "    base = log_file.replace(\"_\", \"\").replace(\".\", \"\").replace(\"txt\", \"\")\n",
+    "    pf.add_observations(log_file, prefix=base, obsgp=base)"
    ]
   },
   {
@@ -639,7 +687,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "_ = [print(line.rstrip()) for line in open(os.path.join(template_ws,\"forward_run.py\"))]"
+    "_ = [print(line.rstrip()) for line in open(os.path.join(template_ws, \"forward_run.py\"))]"
    ]
   },
   {
@@ -660,8 +708,8 @@
    "outputs": [],
    "source": [
     "pst.control_data.noptmax = 0\n",
-    "pst.write(os.path.join(pf.new_d,\"freyberg.pst\"),version=2)\n",
-    "pyemu.os_utils.run(\"pestpp-ies freyberg.pst\",cwd=pf.new_d)"
+    "pst.write(os.path.join(pf.new_d, \"freyberg.pst\"), version=2)\n",
+    "pyemu.os_utils.run(\"pestpp-ies freyberg.pst\", cwd=pf.new_d)"
    ]
   },
   {
@@ -680,7 +728,7 @@
    "outputs": [],
    "source": [
     "num_reals = 50\n",
-    "pe = pf.draw(num_reals=num_reals,use_specsim=False)"
+    "pe = pf.draw(num_reals=num_reals, use_specsim=False)"
    ]
   },
   {
@@ -705,7 +753,9 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "results = ppu.generate_2d_grid_realizations(os.path.join(pf.new_d,\"org.grb\"),num_reals=num_reals)"
+    "results = ppu.generate_2d_grid_realizations(\n",
+    "    os.path.join(pf.new_d, \"org.grb\"), num_reals=num_reals\n",
+    ")"
    ]
   },
   {
@@ -714,9 +764,9 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "#fig,ax = plt.subplots(1,1)\n",
-    "#ax.set_aspect(\"equal\")\n",
-    "#ax.pcolormesh(m.modelgrid.xcellcenters,m.modelgrid.ycellcenters,results[0])\n",
+    "# fig,ax = plt.subplots(1,1)\n",
+    "# ax.set_aspect(\"equal\")\n",
+    "# ax.pcolormesh(m.modelgrid.xcellcenters,m.modelgrid.ycellcenters,results[0])\n",
     "m.dis.top = results[0]\n",
     "m.dis.top.plot()\n",
     "plt.show()\n",
@@ -736,7 +786,12 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "results = ppu.generate_2d_grid_realizations(os.path.join(pf.new_d,\"org.grb\"),num_reals=num_reals,variobearing=90,varioaniso=3)\n",
+    "results = ppu.generate_2d_grid_realizations(\n",
+    "    os.path.join(pf.new_d, \"org.grb\"),\n",
+    "    num_reals=num_reals,\n",
+    "    variobearing=90,\n",
+    "    varioaniso=3,\n",
+    ")\n",
     "m.dis.top = results[0]\n",
     "m.dis.top.plot(colorbar=True)\n",
     "plt.show()\n",
@@ -766,7 +821,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "grid_groups = [\"npfklayer1_attr:gr\",\"npfklayer2_attr:gr\",\"npfklayer3_attr:gr\"]"
+    "grid_groups = [\"npfklayer1_attr:gr\", \"npfklayer2_attr:gr\", \"npfklayer3_attr:gr\"]"
    ]
   },
   {
@@ -775,22 +830,33 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "for igrp,grp in enumerate(grid_groups):\n",
-    "    grpar = par.loc[par.pargp==grp,:].copy()\n",
+    "for igrp, grp in enumerate(grid_groups):\n",
+    "    grpar = par.loc[par.pargp == grp, :].copy()\n",
     "    assert grpar.shape[0] > 0\n",
     "    grpar[\"i\"] = grpar.i.astype(int)\n",
-    "    \n",
-    "    names,ivals,jvals = grpar.parnme.values,grpar.i.values,grpar.j.values\n",
-    "    results = ppu.generate_2d_grid_realizations(os.path.join(pf.new_d,\"org.grb\"),num_reals=num_reals,\n",
-    "                                            variobearing=1.0,varioaniso=1.0,variorange=1000,variance=0.0125,mean=0.0,random_seed=12345*igrp)\n",
+    "\n",
+    "    names, ivals, jvals = grpar.parnme.values, grpar.i.values, grpar.j.values\n",
+    "    results = ppu.generate_2d_grid_realizations(\n",
+    "        os.path.join(pf.new_d, \"org.grb\"),\n",
+    "        num_reals=num_reals,\n",
+    "        variobearing=1.0,\n",
+    "        varioaniso=1.0,\n",
+    "        variorange=1000,\n",
+    "        variance=0.0125,\n",
+    "        mean=0.0,\n",
+    "        random_seed=12345 * igrp,\n",
+    "    )\n",
     "    print(grpar.shape)\n",
-    "    for ireal,real in enumerate(results):\n",
-    "        real = real[id_arr>0]\n",
-    "        #print(real.shape,ivals.shape,jvals.shape)\n",
-    "        pe.loc[pe.index[ireal],names] = 10**(real)\n",
-    "    print(\"group bound info: \",grpar.parlbnd.unique(),grpar.parubnd.unique())\n",
-    "    print(\"ensemble range info: \",pe.loc[:,names].values.min(),pe.loc[:,names].values.max())\n",
-    "    "
+    "    for ireal, real in enumerate(results):\n",
+    "        real = real[id_arr > 0]\n",
+    "        # print(real.shape,ivals.shape,jvals.shape)\n",
+    "        pe.loc[pe.index[ireal], names] = 10 ** (real)\n",
+    "    print(\"group bound info: \", grpar.parlbnd.unique(), grpar.parubnd.unique())\n",
+    "    print(\n",
+    "        \"ensemble range info: \",\n",
+    "        pe.loc[:, names].values.min(),\n",
+    "        pe.loc[:, names].values.max(),\n",
+    "    )"
    ]
   },
   {
@@ -806,8 +872,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "arr = np.zeros_like(id_arr,dtype=float)\n",
-    "arr[id_arr>0] = pe.loc[pe.index[0],names].values\n",
+    "arr = np.zeros_like(id_arr, dtype=float)\n",
+    "arr[id_arr > 0] = pe.loc[pe.index[0], names].values\n",
     "print(arr)\n",
     "m.dis.top = arr\n",
     "m.dis.top.plot(colorbar=True)"
@@ -829,7 +895,7 @@
    "outputs": [],
    "source": [
     "pe.enforce()\n",
-    "pe.to_csv(os.path.join(template_ws,\"prior.csv\"))"
+    "pe.to_csv(os.path.join(template_ws, \"prior.csv\"))"
    ]
   },
   {
diff --git a/examples/pypestutils_pstfrom_structured_mf6.ipynb b/examples/pypestutils_pstfrom_structured_mf6.ipynb
index 887c129..95e32f8 100644
--- a/examples/pypestutils_pstfrom_structured_mf6.ipynb
+++ b/examples/pypestutils_pstfrom_structured_mf6.ipynb
@@ -52,10 +52,10 @@
    "source": [
     "bin_dir = \"struct_bin\"\n",
     "if os.path.exists(bin_dir):\n",
-    "   shutil.rmtree(bin_dir)\n",
+    "    shutil.rmtree(bin_dir)\n",
     "os.makedirs(bin_dir)\n",
-    "flopy.utils.get_modflow(bin_dir,downloads_dir=bin_dir)\n",
-    "pyemu.utils.get_pestpp(bin_dir,downloads_dir=bin_dir)"
+    "flopy.utils.get_modflow(bin_dir, downloads_dir=bin_dir)\n",
+    "pyemu.utils.get_pestpp(bin_dir, downloads_dir=bin_dir)"
    ]
   },
   {
@@ -71,7 +71,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "org_model_ws = os.path.join('freyberg_monthly')\n",
+    "org_model_ws = os.path.join(\"freyberg_monthly\")\n",
     "os.listdir(org_model_ws)"
    ]
   },
@@ -90,9 +90,9 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "id_arr = np.loadtxt(os.path.join(org_model_ws,\"freyberg6.dis_idomain_layer3.txt\"))\n",
-    "top_arr = np.loadtxt(os.path.join(org_model_ws,\"freyberg6.dis_top.txt\"))\n",
-    "top_arr[id_arr==0] = np.nan\n",
+    "id_arr = np.loadtxt(os.path.join(org_model_ws, \"freyberg6.dis_idomain_layer3.txt\"))\n",
+    "top_arr = np.loadtxt(os.path.join(org_model_ws, \"freyberg6.dis_top.txt\"))\n",
+    "top_arr[id_arr == 0] = np.nan\n",
     "plt.imshow(top_arr)\n",
     "top_arr.shape"
    ]
@@ -113,9 +113,9 @@
     "tmp_model_ws = \"temp_ppu_struct\"\n",
     "if os.path.exists(tmp_model_ws):\n",
     "    shutil.rmtree(tmp_model_ws)\n",
-    "shutil.copytree(org_model_ws,tmp_model_ws)\n",
+    "shutil.copytree(org_model_ws, tmp_model_ws)\n",
     "for bin_file in os.listdir(bin_dir):\n",
-    "    shutil.copy2(os.path.join(bin_dir,bin_file),os.path.join(tmp_model_ws,bin_file))\n",
+    "    shutil.copy2(os.path.join(bin_dir, bin_file), os.path.join(tmp_model_ws, bin_file))\n",
     "os.listdir(tmp_model_ws)"
    ]
   },
@@ -135,7 +135,7 @@
    "outputs": [],
    "source": [
     "sim = flopy.mf6.MFSimulation.load(sim_ws=tmp_model_ws)\n",
-    "m = sim.get_model(\"freyberg6\")\n"
+    "m = sim.get_model(\"freyberg6\")"
    ]
   },
   {
@@ -152,10 +152,15 @@
    "outputs": [],
    "source": [
     "template_ws = \"freyberg6_template\"\n",
-    "pf = pyemu.utils.PstFrom(original_d=tmp_model_ws, new_d=template_ws,\n",
-    "                 remove_existing=True,\n",
-    "                 longnames=True, spatial_reference=m.modelgrid,\n",
-    "                 zero_based=False,start_datetime=\"1-1-2018\")"
+    "pf = pyemu.utils.PstFrom(\n",
+    "    original_d=tmp_model_ws,\n",
+    "    new_d=template_ws,\n",
+    "    remove_existing=True,\n",
+    "    longnames=True,\n",
+    "    spatial_reference=m.modelgrid,\n",
+    "    zero_based=False,\n",
+    "    start_datetime=\"1-1-2018\",\n",
+    ")"
    ]
   },
   {
@@ -173,7 +178,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "hdsdf = pd.read_csv(os.path.join(tmp_model_ws,\"heads.csv\"),index_col=0)\n",
+    "hdsdf = pd.read_csv(os.path.join(tmp_model_ws, \"heads.csv\"), index_col=0)\n",
     "hdsdf"
    ]
   },
@@ -190,8 +195,13 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "hds_df = pf.add_observations(\"heads.csv\",insfile=\"heads.csv.ins\",index_cols=\"time\",\n",
-    "                    use_cols=list(hdsdf.columns.values),prefix=\"hds\",)\n",
+    "hds_df = pf.add_observations(\n",
+    "    \"heads.csv\",\n",
+    "    insfile=\"heads.csv.ins\",\n",
+    "    index_cols=\"time\",\n",
+    "    use_cols=list(hdsdf.columns.values),\n",
+    "    prefix=\"hds\",\n",
+    ")\n",
     "hds_df"
    ]
   },
@@ -227,7 +237,12 @@
    "outputs": [],
    "source": [
     "df = pd.read_csv(os.path.join(tmp_model_ws, \"sfr.csv\"), index_col=0)\n",
-    "sfr_df = pf.add_observations(\"sfr.csv\", insfile=\"sfr.csv.ins\", index_cols=\"time\", use_cols=list(df.columns.values))\n",
+    "sfr_df = pf.add_observations(\n",
+    "    \"sfr.csv\",\n",
+    "    insfile=\"sfr.csv.ins\",\n",
+    "    index_cols=\"time\",\n",
+    "    use_cols=list(df.columns.values),\n",
+    ")\n",
     "sfr_df"
    ]
   },
@@ -266,7 +281,9 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "ppdf = ppu.get_2d_pp_info_structured_grid(10,os.path.join(pf.new_d,\"freyberg6.dis.grb\"))\n",
+    "ppdf = ppu.get_2d_pp_info_structured_grid(\n",
+    "    10, os.path.join(pf.new_d, \"freyberg6.dis.grb\")\n",
+    ")\n",
     "ppdf.head()"
    ]
   },
@@ -276,11 +293,11 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "fig,ax = plt.subplots(1,1)\n",
+    "fig, ax = plt.subplots(1, 1)\n",
     "ax.set_aspect(\"equal\")\n",
-    "#ax.pcolormesh(sr.xcentergrid,sr.ycentergrid,top_arr)\n",
-    "ax.pcolormesh(m.modelgrid.xcellcenters,m.modelgrid.ycellcenters,top_arr)\n",
-    "ax.scatter(ppdf.x,ppdf.y,marker=\".\",c=\"r\")"
+    "# ax.pcolormesh(sr.xcentergrid,sr.ycentergrid,top_arr)\n",
+    "ax.pcolormesh(m.modelgrid.xcellcenters, m.modelgrid.ycellcenters, top_arr)\n",
+    "ax.scatter(ppdf.x, ppdf.y, marker=\".\", c=\"r\")"
    ]
   },
   {
@@ -296,7 +313,11 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "ppdf = ppu.get_2d_pp_info_structured_grid(10,os.path.join(pf.new_d,\"freyberg6.dis.grb\"),array_dict={\"zone\":m.dis.idomain.array[0,:,:]})\n",
+    "ppdf = ppu.get_2d_pp_info_structured_grid(\n",
+    "    10,\n",
+    "    os.path.join(pf.new_d, \"freyberg6.dis.grb\"),\n",
+    "    array_dict={\"zone\": m.dis.idomain.array[0, :, :]},\n",
+    ")\n",
     "ppdf.head()"
    ]
   },
@@ -306,10 +327,10 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "fig,ax = plt.subplots(1,1)\n",
+    "fig, ax = plt.subplots(1, 1)\n",
     "ax.set_aspect(\"equal\")\n",
-    "ax.pcolormesh(m.modelgrid.xcellcenters,m.modelgrid.ycellcenters,top_arr)\n",
-    "ax.scatter(ppdf.x,ppdf.y,marker=\".\",c=\"r\")"
+    "ax.pcolormesh(m.modelgrid.xcellcenters, m.modelgrid.ycellcenters, top_arr)\n",
+    "ax.scatter(ppdf.x, ppdf.y, marker=\".\", c=\"r\")"
    ]
   },
   {
@@ -332,7 +353,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "ppdf.loc[:,\"value\"] = 3.0  # HK in layer 1"
+    "ppdf.loc[:, \"value\"] = 3.0  # HK in layer 1"
    ]
   },
   {
@@ -343,10 +364,10 @@
    "source": [
     "# set aniso to be a function of column value\n",
     "# stronger aniso near the sfr network\n",
-    "jmin,jmax = ppdf.j.min(),float(ppdf.j.max())\n",
-    "ppdf.loc[:,\"aniso\"] = 30 * ppdf.j.values.copy() / jmax\n",
-    "ppdf.loc[ppdf.aniso<1,\"aniso\"] = 1\n",
-    "ppdf.loc[ppdf.aniso>10,\"aniso\"] = 10\n"
+    "jmin, jmax = ppdf.j.min(), float(ppdf.j.max())\n",
+    "ppdf.loc[:, \"aniso\"] = 30 * ppdf.j.values.copy() / jmax\n",
+    "ppdf.loc[ppdf.aniso < 1, \"aniso\"] = 1\n",
+    "ppdf.loc[ppdf.aniso > 10, \"aniso\"] = 10"
    ]
   },
   {
@@ -359,8 +380,8 @@
     "# cl = ppdf.corrlen.min()\n",
     "# ppdf.loc[:,\"corrlen\"] = cl * (ppdf.j.values.copy() / jmax) * sr.delc.max()\n",
     "# ppdf.loc[ppdf.corrlen<cl/20,\"corrlen\"] = cl/20\n",
-    "ppdf.loc[:,\"x\"] = np.round(ppdf.x.values,1)\n",
-    "ppdf.loc[:,\"y\"] = np.round(ppdf.y.values,1)\n",
+    "ppdf.loc[:, \"x\"] = np.round(ppdf.x.values, 1)\n",
+    "ppdf.loc[:, \"y\"] = np.round(ppdf.y.values, 1)\n",
     "ppdf.corrlen *= 10\n",
     "ppdf.corrlen.describe()"
    ]
@@ -371,13 +392,13 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "#set bearing to be a high-freq sin wave in the sfr direction\n",
+    "# set bearing to be a high-freq sin wave in the sfr direction\n",
     "y = np.cumsum(m.dis.delc.array)\n",
-    "phase = np.pi/4\n",
+    "phase = np.pi / 4\n",
     "gain = 40\n",
-    "iy = np.linspace(phase,5*np.pi+phase,y.shape[0])\n",
+    "iy = np.linspace(phase, 5 * np.pi + phase, y.shape[0])\n",
     "ix = 235 + np.sin(iy) * gain\n",
-    "ppdf.loc[:,\"bearing\"] = ix[ppdf.i.values]\n"
+    "ppdf.loc[:, \"bearing\"] = ix[ppdf.i.values]"
    ]
   },
   {
@@ -394,14 +415,14 @@
    "outputs": [],
    "source": [
     "id_mask = id_arr.copy()\n",
-    "id_mask[id_mask!=0] = np.nan\n",
-    "fig,axes = plt.subplots(1,3)\n",
-    "for ax,attr in zip(axes,[\"aniso\",\"corrlen\",\"bearing\"]):\n",
-    "    #ax.pcolormesh(sr.xcentergrid,sr.ycentergrid,id_mask)\n",
-    "    ax.pcolormesh(m.modelgrid.xcellcenters,m.modelgrid.ycellcenters,id_mask)\n",
-    "    ax.scatter(ppdf.x,ppdf.y,c=ppdf.loc[:,attr])\n",
+    "id_mask[id_mask != 0] = np.nan\n",
+    "fig, axes = plt.subplots(1, 3)\n",
+    "for ax, attr in zip(axes, [\"aniso\", \"corrlen\", \"bearing\"]):\n",
+    "    # ax.pcolormesh(sr.xcentergrid,sr.ycentergrid,id_mask)\n",
+    "    ax.pcolormesh(m.modelgrid.xcellcenters, m.modelgrid.ycellcenters, id_mask)\n",
+    "    ax.scatter(ppdf.x, ppdf.y, c=ppdf.loc[:, attr])\n",
     "    ax.set_aspect(\"equal\")\n",
-    "    ax.set_title(attr,loc=\"left\")\n",
+    "    ax.set_title(attr, loc=\"left\")\n",
     "    ax.set_xticks([])\n",
     "    ax.set_yticks([])\n",
     "plt.tight_layout()"
@@ -420,9 +441,11 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "ppdf.to_csv(os.path.join(pf.new_d,\"pp_info.csv\"))\n",
-    "ppdf.loc[:,\"value\"] = np.random.normal(1.0,0.25,ppdf.shape[0])\n",
-    "results = ppu.interpolate_with_sva_pilotpoints_2d(ppdf,os.path.join(pf.new_d,\"freyberg6.dis.grb\"))"
+    "ppdf.to_csv(os.path.join(pf.new_d, \"pp_info.csv\"))\n",
+    "ppdf.loc[:, \"value\"] = np.random.normal(1.0, 0.25, ppdf.shape[0])\n",
+    "results = ppu.interpolate_with_sva_pilotpoints_2d(\n",
+    "    ppdf, os.path.join(pf.new_d, \"freyberg6.dis.grb\")\n",
+    ")"
    ]
   },
   {
@@ -431,21 +454,21 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "fig,ax = plt.subplots(1,1)\n",
-    "ax.pcolormesh(m.modelgrid.xcellcenters,m.modelgrid.ycellcenters,id_mask)\n",
-    "ax.scatter(ppdf.x,ppdf.y,c=ppdf.loc[:,\"value\"])\n",
+    "fig, ax = plt.subplots(1, 1)\n",
+    "ax.pcolormesh(m.modelgrid.xcellcenters, m.modelgrid.ycellcenters, id_mask)\n",
+    "ax.scatter(ppdf.x, ppdf.y, c=ppdf.loc[:, \"value\"])\n",
     "ax.set_aspect(\"equal\")\n",
-    "_= ax.set_title(\"pilot point values\")\n",
+    "_ = ax.set_title(\"pilot point values\")\n",
     "plt.show()\n",
     "plt.close(fig)\n",
-    "for tag,arr in results.items():\n",
-    "    fig,ax = plt.subplots(1,1)\n",
+    "for tag, arr in results.items():\n",
+    "    fig, ax = plt.subplots(1, 1)\n",
     "    ax.set_aspect(\"equal\")\n",
-    "    ax.set_title(tag,loc=\"left\")\n",
-    "    cb = ax.pcolormesh(m.modelgrid.xcellcenters,m.modelgrid.ycellcenters,arr)\n",
+    "    ax.set_title(tag, loc=\"left\")\n",
+    "    cb = ax.pcolormesh(m.modelgrid.xcellcenters, m.modelgrid.ycellcenters, arr)\n",
     "    plt.colorbar(cb)\n",
     "    plt.show()\n",
-    "    plt.close(fig)\n"
+    "    plt.close(fig)"
    ]
   },
   {
@@ -477,7 +500,9 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "hk_arr_files = [f for f in os.listdir(tmp_model_ws) if \"npf_k_\" in f and f.endswith(\".txt\")]\n",
+    "hk_arr_files = [\n",
+    "    f for f in os.listdir(tmp_model_ws) if \"npf_k_\" in f and f.endswith(\".txt\")\n",
+    "]\n",
     "hk_arr_files"
    ]
   },
@@ -498,9 +523,9 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "ppdf.loc[:,\"value\"] = 1.0\n",
+    "ppdf.loc[:, \"value\"] = 1.0\n",
     "pp_v = pyemu.geostats.ExpVario(contribution=1.0, a=2000)\n",
-    "pp_gs = pyemu.geostats.GeoStruct(variograms=pp_v,transform=\"log\")"
+    "pp_gs = pyemu.geostats.GeoStruct(variograms=pp_v, transform=\"log\")"
    ]
   },
   {
@@ -516,29 +541,52 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "pp_files,mod_files = [],[]\n",
+    "pp_files, mod_files = [], []\n",
     "for hk_arr_file in hk_arr_files:\n",
     "    layer = int(hk_arr_file.split(\".\")[1].split(\"layer\")[1])\n",
-    "    base = hk_arr_file.split('.')[1].replace(\"_\",\"\")+\"_attr:\"\n",
-    "    pf.add_parameters(filenames=hk_arr_file,par_type=\"grid\",\n",
-    "                       par_name_base=base+\"gr\",pargp=base+\"gr\",zone_array=ib,\n",
-    "                       upper_bound=2.,lower_bound=0.5,ult_ubound=100,ult_lbound=0.01)\n",
+    "    base = hk_arr_file.split(\".\")[1].replace(\"_\", \"\") + \"_attr:\"\n",
+    "    pf.add_parameters(\n",
+    "        filenames=hk_arr_file,\n",
+    "        par_type=\"grid\",\n",
+    "        par_name_base=base + \"gr\",\n",
+    "        pargp=base + \"gr\",\n",
+    "        zone_array=ib,\n",
+    "        upper_bound=2.0,\n",
+    "        lower_bound=0.5,\n",
+    "        ult_ubound=100,\n",
+    "        ult_lbound=0.01,\n",
+    "    )\n",
     "    pppdf = ppdf.copy()\n",
-    "    \n",
-    "    pppdf.loc[:,\"name\"] = [n for n in pppdf.ppname.values]\n",
-    "    pppdf.loc[:,\"ppname\"] = pppdf.name.values\n",
-    "    pp_file = os.path.join(pf.new_d,base+\"pp.csv\")\n",
-    "    pppdf.to_csv(pp_file,index=False)\n",
+    "\n",
+    "    pppdf.loc[:, \"name\"] = [n for n in pppdf.ppname.values]\n",
+    "    pppdf.loc[:, \"ppname\"] = pppdf.name.values\n",
+    "    pp_file = os.path.join(pf.new_d, base + \"pp.csv\")\n",
+    "    pppdf.to_csv(pp_file, index=False)\n",
     "    pp_files.append(os.path.split(pp_file)[1])\n",
     "    mod_files.append(hk_arr_file)\n",
-    "    df = pf.add_parameters(os.path.split(pp_file)[1],par_type=\"grid\",index_cols={\"ppname\":\"ppname\",\"x\":\"x\",\"y\":\"y\"},\n",
-    "        use_cols=[\"value\",\"bearing\",\"aniso\",\"corrlen\"],\n",
-    "        par_name_base=[base+\"pp\",base+\"bearing\",base+\"aniso\",base+\"corrlen\"],\n",
-    "        pargp=[base+\"pp\",base+\"bearing\",base+\"aniso\",base+\"corrlen\"],\n",
-    "        upper_bound=[20,pppdf.bearing.max()*1.1,pppdf.aniso.max()*1.1,pppdf.corrlen.max()*1.1],\n",
-    "        lower_bound=[.05,pppdf.bearing.min()*.9,pppdf.aniso.min()*.9,pppdf.corrlen.min()*.9],\n",
-    "        par_style=\"direct\",transform=\"log\",geostruct=pp_gs)\n",
-    "    "
+    "    df = pf.add_parameters(\n",
+    "        os.path.split(pp_file)[1],\n",
+    "        par_type=\"grid\",\n",
+    "        index_cols={\"ppname\": \"ppname\", \"x\": \"x\", \"y\": \"y\"},\n",
+    "        use_cols=[\"value\", \"bearing\", \"aniso\", \"corrlen\"],\n",
+    "        par_name_base=[base + \"pp\", base + \"bearing\", base + \"aniso\", base + \"corrlen\"],\n",
+    "        pargp=[base + \"pp\", base + \"bearing\", base + \"aniso\", base + \"corrlen\"],\n",
+    "        upper_bound=[\n",
+    "            20,\n",
+    "            pppdf.bearing.max() * 1.1,\n",
+    "            pppdf.aniso.max() * 1.1,\n",
+    "            pppdf.corrlen.max() * 1.1,\n",
+    "        ],\n",
+    "        lower_bound=[\n",
+    "            0.05,\n",
+    "            pppdf.bearing.min() * 0.9,\n",
+    "            pppdf.aniso.min() * 0.9,\n",
+    "            pppdf.corrlen.min() * 0.9,\n",
+    "        ],\n",
+    "        par_style=\"direct\",\n",
+    "        transform=\"log\",\n",
+    "        geostruct=pp_gs,\n",
+    "    )"
    ]
   },
   {
@@ -564,10 +612,9 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "with open(os.path.join(template_ws,tpl_files[0]),'r') as f:\n",
+    "with open(os.path.join(template_ws, tpl_files[0]), \"r\") as f:\n",
     "    for _ in range(4):\n",
-    "        print(f.readline().strip())\n",
-    "        "
+    "        print(f.readline().strip())"
    ]
   },
   {
@@ -585,7 +632,9 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "shutil.copy2(os.path.join(pf.new_d,\"freyberg6.dis.grb\"),os.path.join(pf.new_d,\"org.grb\"))"
+    "shutil.copy2(\n",
+    "    os.path.join(pf.new_d, \"freyberg6.dis.grb\"), os.path.join(pf.new_d, \"org.grb\")\n",
+    ")"
    ]
   },
   {
@@ -594,7 +643,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "_ = [print(line.rstrip()) for line in open(\"ppu_helpers.py\",'r').readlines()]"
+    "_ = [print(line.rstrip()) for line in open(\"ppu_helpers.py\", \"r\").readlines()]"
    ]
   },
   {
@@ -611,8 +660,8 @@
    "outputs": [],
    "source": [
     "# save the pp file and model file info\n",
-    "df = pd.DataFrame({\"model_file\":mod_files,\"pp_file\":pp_files})\n",
-    "df.to_csv(os.path.join(pf.new_d,\"pp_info.csv\"))"
+    "df = pd.DataFrame({\"model_file\": mod_files, \"pp_file\": pp_files})\n",
+    "df.to_csv(os.path.join(pf.new_d, \"pp_info.csv\"))"
    ]
   },
   {
@@ -622,6 +671,7 @@
    "outputs": [],
    "source": [
     "import ppu_helpers\n",
+    "\n",
     "ppu_helpers.setup_pps(pf.new_d)"
    ]
   },
@@ -631,7 +681,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "pf.add_py_function(\"ppu_helpers.py\",\"apply_pps()\",is_pre_cmd=True)"
+    "pf.add_py_function(\"ppu_helpers.py\", \"apply_pps()\", is_pre_cmd=True)"
    ]
   },
   {
@@ -647,7 +697,9 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "interp_files = [f for f in os.listdir(pf.new_d) if f.endswith(\".txt\") and f.startswith(\"interp\")]\n",
+    "interp_files = [\n",
+    "    f for f in os.listdir(pf.new_d) if f.endswith(\".txt\") and f.startswith(\"interp\")\n",
+    "]\n",
     "assert len(interp_files) == 3"
    ]
   },
@@ -658,8 +710,8 @@
    "outputs": [],
    "source": [
     "for interp_file in interp_files:\n",
-    "    base = interp_file.replace(\"_\",\"\").replace(\".\",\"\").replace(\"txt\",\"\")\n",
-    "    pf.add_observations(interp_file,prefix=base,obsgp=base)"
+    "    base = interp_file.replace(\"_\", \"\").replace(\".\", \"\").replace(\"txt\", \"\")\n",
+    "    pf.add_observations(interp_file, prefix=base, obsgp=base)"
    ]
   },
   {
@@ -668,9 +720,10 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "\n",
-    "log_files = [f for f in os.listdir(pf.new_d) if f.endswith(\".txt\") and f.startswith(\"log_\")]\n",
-    "assert len(log_files) == 3\n"
+    "log_files = [\n",
+    "    f for f in os.listdir(pf.new_d) if f.endswith(\".txt\") and f.startswith(\"log_\")\n",
+    "]\n",
+    "assert len(log_files) == 3"
    ]
   },
   {
@@ -680,8 +733,8 @@
    "outputs": [],
    "source": [
     "for log_file in log_files:\n",
-    "    base = log_file.replace(\"_\",\"\").replace(\".\",\"\").replace(\"txt\",\"\")\n",
-    "    pf.add_observations(log_file,prefix=base,obsgp=base)"
+    "    base = log_file.replace(\"_\", \"\").replace(\".\", \"\").replace(\"txt\", \"\")\n",
+    "    pf.add_observations(log_file, prefix=base, obsgp=base)"
    ]
   },
   {
@@ -717,7 +770,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "_ = [print(line.rstrip()) for line in open(os.path.join(template_ws,\"forward_run.py\"))]"
+    "_ = [print(line.rstrip()) for line in open(os.path.join(template_ws, \"forward_run.py\"))]"
    ]
   },
   {
@@ -738,8 +791,8 @@
    "outputs": [],
    "source": [
     "pst.control_data.noptmax = 0\n",
-    "pst.write(os.path.join(pf.new_d,\"freyberg.pst\"),version=2)\n",
-    "pyemu.os_utils.run(\"pestpp-ies freyberg.pst\",cwd=pf.new_d)"
+    "pst.write(os.path.join(pf.new_d, \"freyberg.pst\"), version=2)\n",
+    "pyemu.os_utils.run(\"pestpp-ies freyberg.pst\", cwd=pf.new_d)"
    ]
   },
   {
@@ -758,7 +811,7 @@
    "outputs": [],
    "source": [
     "num_reals = 30\n",
-    "pe = pf.draw(num_reals=num_reals,use_specsim=True)"
+    "pe = pf.draw(num_reals=num_reals, use_specsim=True)"
    ]
   },
   {
@@ -788,7 +841,9 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "results = ppu.generate_2d_grid_realizations(os.path.join(pf.new_d,\"org.grb\"),num_reals=num_reals)"
+    "results = ppu.generate_2d_grid_realizations(\n",
+    "    os.path.join(pf.new_d, \"org.grb\"), num_reals=num_reals\n",
+    ")"
    ]
   },
   {
@@ -797,9 +852,9 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "fig,ax = plt.subplots(1,1)\n",
+    "fig, ax = plt.subplots(1, 1)\n",
     "ax.set_aspect(\"equal\")\n",
-    "ax.pcolormesh(m.modelgrid.xcellcenters,m.modelgrid.ycellcenters,results[0])\n",
+    "ax.pcolormesh(m.modelgrid.xcellcenters, m.modelgrid.ycellcenters, results[0])\n",
     "plt.show()\n",
     "plt.close(fig)"
    ]
@@ -817,10 +872,15 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "results = ppu.generate_2d_grid_realizations(os.path.join(pf.new_d,\"org.grb\"),num_reals=num_reals,variobearing=90,varioaniso=3)\n",
-    "fig,ax = plt.subplots(1,1)\n",
+    "results = ppu.generate_2d_grid_realizations(\n",
+    "    os.path.join(pf.new_d, \"org.grb\"),\n",
+    "    num_reals=num_reals,\n",
+    "    variobearing=90,\n",
+    "    varioaniso=3,\n",
+    ")\n",
+    "fig, ax = plt.subplots(1, 1)\n",
     "ax.set_aspect(\"equal\")\n",
-    "ax.pcolormesh(m.modelgrid.xcellcenters,m.modelgrid.ycellcenters,results[0])\n",
+    "ax.pcolormesh(m.modelgrid.xcellcenters, m.modelgrid.ycellcenters, results[0])\n",
     "plt.show()\n",
     "plt.close(fig)"
    ]
@@ -839,8 +899,8 @@
    "outputs": [],
    "source": [
     "# some redic bearing that makes a pretty picture...\n",
-    "nrow,ncol = m.dis.nrow.data,m.dis.ncol.data\n",
-    "bearing = np.add(np.ones((nrow,ncol)),np.atleast_2d(np.arange(ncol)))"
+    "nrow, ncol = m.dis.nrow.data, m.dis.ncol.data\n",
+    "bearing = np.add(np.ones((nrow, ncol)), np.atleast_2d(np.arange(ncol)))"
    ]
   },
   {
@@ -849,11 +909,18 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "results = ppu.generate_2d_grid_realizations(os.path.join(pf.new_d,\"org.grb\"),num_reals=num_reals,\n",
-    "                                            variobearing=bearing,varioaniso=7,variorange=2000,variance=0.025,mean=0.0)\n",
-    "fig,ax = plt.subplots(1,1)\n",
+    "results = ppu.generate_2d_grid_realizations(\n",
+    "    os.path.join(pf.new_d, \"org.grb\"),\n",
+    "    num_reals=num_reals,\n",
+    "    variobearing=bearing,\n",
+    "    varioaniso=7,\n",
+    "    variorange=2000,\n",
+    "    variance=0.025,\n",
+    "    mean=0.0,\n",
+    ")\n",
+    "fig, ax = plt.subplots(1, 1)\n",
     "ax.set_aspect(\"equal\")\n",
-    "ax.pcolormesh(m.modelgrid.xcellcenters,m.modelgrid.ycellcenters,results[0])\n",
+    "ax.pcolormesh(m.modelgrid.xcellcenters, m.modelgrid.ycellcenters, results[0])\n",
     "plt.show()\n",
     "plt.close(fig)"
    ]
@@ -881,7 +948,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "grid_groups = [\"npfklayer1_attr:gr\",\"npfklayer2_attr:gr\",\"npfklayer3_attr:gr\"]"
+    "grid_groups = [\"npfklayer1_attr:gr\", \"npfklayer2_attr:gr\", \"npfklayer3_attr:gr\"]"
    ]
   },
   {
@@ -890,19 +957,30 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "for igrp,grp in enumerate(grid_groups):\n",
-    "    grpar = par.loc[par.pargp==grp,:].copy()\n",
+    "for igrp, grp in enumerate(grid_groups):\n",
+    "    grpar = par.loc[par.pargp == grp, :].copy()\n",
     "    assert grpar.shape[0] > 0\n",
     "    grpar[\"i\"] = grpar.i.astype(int)\n",
     "    grpar[\"j\"] = grpar.j.astype(int)\n",
-    "    names,ivals,jvals = grpar.parnme.values,grpar.i.values,grpar.j.values\n",
-    "    results = ppu.generate_2d_grid_realizations(os.path.join(pf.new_d,\"org.grb\"),num_reals=num_reals,\n",
-    "                                            variobearing=1.0,varioaniso=1.0,variorange=1000,variance=0.0125,mean=0.0,random_seed=12345*igrp)\n",
-    "    for ireal,real in enumerate(results):\n",
-    "        pe.loc[pe.index[ireal],names] = 10**(real[ivals,jvals])\n",
-    "    print(\"group bound info: \",grpar.parlbnd.unique(),grpar.parubnd.unique())\n",
-    "    print(\"ensemble range info: \",pe.loc[:,names].values.min(),pe.loc[:,names].values.max())\n",
-    "    "
+    "    names, ivals, jvals = grpar.parnme.values, grpar.i.values, grpar.j.values\n",
+    "    results = ppu.generate_2d_grid_realizations(\n",
+    "        os.path.join(pf.new_d, \"org.grb\"),\n",
+    "        num_reals=num_reals,\n",
+    "        variobearing=1.0,\n",
+    "        varioaniso=1.0,\n",
+    "        variorange=1000,\n",
+    "        variance=0.0125,\n",
+    "        mean=0.0,\n",
+    "        random_seed=12345 * igrp,\n",
+    "    )\n",
+    "    for ireal, real in enumerate(results):\n",
+    "        pe.loc[pe.index[ireal], names] = 10 ** (real[ivals, jvals])\n",
+    "    print(\"group bound info: \", grpar.parlbnd.unique(), grpar.parubnd.unique())\n",
+    "    print(\n",
+    "        \"ensemble range info: \",\n",
+    "        pe.loc[:, names].values.min(),\n",
+    "        pe.loc[:, names].values.max(),\n",
+    "    )"
    ]
   },
   {
@@ -918,9 +996,9 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "arr = np.zeros((nrow,ncol))\n",
-    "arr[ivals,jvals] = pe.loc[pe.index[0],names]\n",
-    "arr[id_arr==0] = np.nan\n",
+    "arr = np.zeros((nrow, ncol))\n",
+    "arr[ivals, jvals] = pe.loc[pe.index[0], names]\n",
+    "arr[id_arr == 0] = np.nan\n",
     "cb = plt.imshow(arr)\n",
     "plt.colorbar(cb)"
    ]
@@ -941,7 +1019,7 @@
    "outputs": [],
    "source": [
     "pe.enforce()\n",
-    "pe.to_csv(os.path.join(template_ws,\"prior.csv\"))"
+    "pe.to_csv(os.path.join(template_ws, \"prior.csv\"))"
    ]
   }
  ],
diff --git a/examples/run_notebooks.py b/examples/run_notebooks.py
index 1a30316..f5bfe74 100644
--- a/examples/run_notebooks.py
+++ b/examples/run_notebooks.py
@@ -7,19 +7,34 @@
 pdf = False
 allow_errors = True
 
-def run_nb(nb_file, nb_dir): 
+
+def run_nb(nb_file, nb_dir):
     os.chdir(nb_dir)
-    worker_dirs = [d for d in os.listdir(".") if os.path.isdir(d) and d.startswith("worker")]
+    worker_dirs = [
+        d for d in os.listdir(".") if os.path.isdir(d) and d.startswith("worker")
+    ]
     for worker_dir in worker_dirs:
         shutil.rmtree(worker_dir)
     if allow_errors:
-        os.system("jupyter nbconvert --execute --ExecutePreprocessor.timeout=180000 --allow-errors --inplace {0}".format(nb_file))
+        os.system(
+            "jupyter nbconvert --execute --ExecutePreprocessor.timeout=180000 --allow-errors --inplace {0}".format(
+                nb_file
+            )
+        )
     else:
-        os.system("jupyter nbconvert --execute --ExecutePreprocessor.timeout=180000 --inplace {0}".format(nb_file))
+        os.system(
+            "jupyter nbconvert --execute --ExecutePreprocessor.timeout=180000 --inplace {0}".format(
+                nb_file
+            )
+        )
     if pdf:
         os.system("jupyter nbconvert --to pdf {0}".format(nb_file))
     if clear:
-        os.system("jupyter nbconvert --ClearOutputPreprocessor.enabled=True --ClearMetadataPreprocessor.enabled=True --allow-errors --inplace {0}".format(nb_file))
+        os.system(
+            "jupyter nbconvert --ClearOutputPreprocessor.enabled=True --ClearMetadataPreprocessor.enabled=True --allow-errors --inplace {0}".format(
+                nb_file
+            )
+        )
     os.chdir(cwd)
     return
 
@@ -30,7 +45,3 @@ def run_nb(nb_file, nb_dir):
 nb_files.sort()
 for nb_file in nb_files:
     run_nb(nb_file, nb_dir)
-
-
-
-
diff --git a/examples/swap.py b/examples/swap.py
index a98ad3d..88374c6 100644
--- a/examples/swap.py
+++ b/examples/swap.py
@@ -21,41 +21,56 @@ def repair_and_prep_quadtree_model():
         if f.startswith("."):
             continue
         try:
-            lines = open(os.path.join(org_d,f),'r').readlines()
+            lines = open(os.path.join(org_d, f), "r").readlines()
         except:
-            print("error for file ",f)
+            print("error for file ", f)
             continue
 
-        with open(os.path.join(new_d,f.replace("project","freyberg6")),'w') as ff:
+        with open(os.path.join(new_d, f.replace("project", "freyberg6")), "w") as ff:
             for line in lines:
-                ff.write(line.lower().replace("project","freyberg6"))
-
+                ff.write(line.lower().replace("project", "freyberg6"))
 
     mn_dir = "freyberg_monthly"
 
-    shutil.copy2(os.path.join(mn_dir,"freyberg6.tdis"),os.path.join(new_d,"freyberg6.tdis"))
+    shutil.copy2(
+        os.path.join(mn_dir, "freyberg6.tdis"), os.path.join(new_d, "freyberg6.tdis")
+    )
 
     wel_files = [f for f in os.listdir(mn_dir) if "wel" in f and f.endswith(".txt")]
-    org_uwell_file = os.path.join(new_d,"freyberg6.wel_stress_period_data_1.txt")
-    org_wel_df = pd.read_csv(org_uwell_file,delim_whitespace=True,header=None,names=["layer","node","flux"])
+    org_uwell_file = os.path.join(new_d, "freyberg6.wel_stress_period_data_1.txt")
+    org_wel_df = pd.read_csv(
+        org_uwell_file,
+        delim_whitespace=True,
+        header=None,
+        names=["layer", "node", "flux"],
+    )
     print(org_wel_df)
     for wel_file in wel_files:
-        df = pd.read_csv(os.path.join(mn_dir,wel_file),delim_whitespace=True,header=None,names=["layer","row","col","flux"])
+        df = pd.read_csv(
+            os.path.join(mn_dir, wel_file),
+            delim_whitespace=True,
+            header=None,
+            names=["layer", "row", "col", "flux"],
+        )
         new_df = org_wel_df.copy()
-        new_df.loc[:,"flux"] = df.flux.values
+        new_df.loc[:, "flux"] = df.flux.values
         print(new_df)
-        new_df.to_csv(os.path.join(new_d,wel_file),sep=" ",index=False,header=False)
-    shutil.copy2(os.path.join(mn_dir,"freyberg6.wel"),os.path.join(new_d,"freyberg6.wel"))
+        new_df.to_csv(os.path.join(new_d, wel_file), sep=" ", index=False, header=False)
+    shutil.copy2(
+        os.path.join(mn_dir, "freyberg6.wel"), os.path.join(new_d, "freyberg6.wel")
+    )
 
     rch_files = [f for f in os.listdir(mn_dir) if "rcha" in f and f.endswith(".txt")]
-    org_urch_file = os.path.join(new_d,"freyberg6.rcha_recharge_1.txt")
+    org_urch_file = os.path.join(new_d, "freyberg6.rcha_recharge_1.txt")
     vals = []
-    with open(org_urch_file,'r') as f:
+    with open(org_urch_file, "r") as f:
         for line in f:
             vals.extend([float(v) for v in line.strip().split()])
     org_urch_arr = np.array(vals)
 
-    shutil.copy(os.path.join(mn_dir,"freyberg6.sto"),os.path.join(new_d,"freyberg6.sto"))
+    shutil.copy(
+        os.path.join(mn_dir, "freyberg6.sto"), os.path.join(new_d, "freyberg6.sto")
+    )
     sto_files = [f for f in os.listdir(mn_dir) if "sto" in f and f.endswith(".txt")]
     print(sto_files)
 
@@ -64,148 +79,192 @@ def repair_and_prep_quadtree_model():
     rch_files.extend(npf_files)
     print(rch_files)
     for rch_file in rch_files:
-        arr = np.loadtxt(os.path.join(mn_dir,rch_file))
+        arr = np.loadtxt(os.path.join(mn_dir, rch_file))
         new_arr = np.zeros_like(org_urch_arr) + arr.mean()
-        np.savetxt(os.path.join(new_d,rch_file),new_arr,fmt="%15.6E")
-    shutil.copy2(os.path.join(mn_dir,"freyberg6.rcha"),os.path.join(new_d,"freyberg6.rcha"))
+        np.savetxt(os.path.join(new_d, rch_file), new_arr, fmt="%15.6E")
+    shutil.copy2(
+        os.path.join(mn_dir, "freyberg6.rcha"), os.path.join(new_d, "freyberg6.rcha")
+    )
 
     # now unwrap all that dumbass wrapped format - what is this 1980?!
-    txt_files = [f for f in os.listdir(os.path.join(new_d)) if f.startswith("freyberg6") and f.endswith(".txt") 
-                and "stress" not in f and "continuous" not in f and "sfr" not in f and "vertices" not in f and "cell" not in f]
+    txt_files = [
+        f
+        for f in os.listdir(os.path.join(new_d))
+        if f.startswith("freyberg6")
+        and f.endswith(".txt")
+        and "stress" not in f
+        and "continuous" not in f
+        and "sfr" not in f
+        and "vertices" not in f
+        and "cell" not in f
+    ]
 
     txt_files.sort()
     for txt_file in txt_files:
-        lines = open(os.path.join(new_d,txt_file),'r').readlines()
+        lines = open(os.path.join(new_d, txt_file), "r").readlines()
         vals = []
         for line in lines:
             vals.extend(line.strip().split())
-        with open(os.path.join(new_d,txt_file),'w') as f:
-            [f.write(v+"\n") for v in vals]
-    #print(txt_files)
-    botm = [32.5,30.0,10.0]
-    for k,b in enumerate(botm):
+        with open(os.path.join(new_d, txt_file), "w") as f:
+            [f.write(v + "\n") for v in vals]
+    # print(txt_files)
+    botm = [32.5, 30.0, 10.0]
+    for k, b in enumerate(botm):
         arr = np.zeros_like(org_urch_arr) + b
-        np.savetxt(os.path.join(new_d,"freyberg6.disv_botm_layer{0}.txt".format(k+1)),arr,fmt="%15.6E")
-    top = np.loadtxt(os.path.join(new_d,"freyberg6.disv_top.txt"))
+        np.savetxt(
+            os.path.join(new_d, "freyberg6.disv_botm_layer{0}.txt".format(k + 1)),
+            arr,
+            fmt="%15.6E",
+        )
+    top = np.loadtxt(os.path.join(new_d, "freyberg6.disv_top.txt"))
     delt = top - botm[0]
 
-    top += np.abs(delt.min()) + 1.1 #?
-    np.savetxt(os.path.join(new_d,"freyberg6.disv_top.txt"),top,fmt="%15.6E")
-    mnsim = flopy.mf6.MFSimulation.load(sim_ws=mn_dir,load_only=["dis"])
+    top += np.abs(delt.min()) + 1.1  # ?
+    np.savetxt(os.path.join(new_d, "freyberg6.disv_top.txt"), top, fmt="%15.6E")
+    mnsim = flopy.mf6.MFSimulation.load(sim_ws=mn_dir, load_only=["dis"])
     mnm = mnsim.get_model()
     mnm.modelgrid.set_coord_info(angrot=0.0)
 
-    usim = flopy.mf6.MFSimulation.load(sim_ws=new_d)#,load_only=["disv"])
+    usim = flopy.mf6.MFSimulation.load(sim_ws=new_d)  # ,load_only=["disv"])
     um = usim.get_model()
     um.modelgrid.set_coord_info(angrot=0.0)
-    print(um.modelgrid,mnm.modelgrid)
-     
+    print(um.modelgrid, mnm.modelgrid)
+
     hobs_csv = "freyberg6.obs_continuous_heads.csv.txt"
-    hds = pd.read_csv(os.path.join(mn_dir,hobs_csv),delim_whitespace=True,header=None,names=["site","otype","layer","row","col"])
-    X,Y = mnm.modelgrid.xcellcenters,mnm.modelgrid.ycellcenters
-    hds.loc[:,"x"] = hds.apply(lambda x: X[x.row-1,x.col-1],axis=1)
-    hds.loc[:,"y"] = hds.apply(lambda x: Y[x.row-1,x.col-1],axis=1)
-    hds.loc[:,"pt"] = hds.apply(lambda x: shapely.geometry.Point(x.x,x.y),axis=1)
+    hds = pd.read_csv(
+        os.path.join(mn_dir, hobs_csv),
+        delim_whitespace=True,
+        header=None,
+        names=["site", "otype", "layer", "row", "col"],
+    )
+    X, Y = mnm.modelgrid.xcellcenters, mnm.modelgrid.ycellcenters
+    hds.loc[:, "x"] = hds.apply(lambda x: X[x.row - 1, x.col - 1], axis=1)
+    hds.loc[:, "y"] = hds.apply(lambda x: Y[x.row - 1, x.col - 1], axis=1)
+    hds.loc[:, "pt"] = hds.apply(lambda x: shapely.geometry.Point(x.x, x.y), axis=1)
 
     i = flopy.utils.GridIntersect(um.modelgrid)
-    #print(i.intersect(shapely.geometry.Point(X[0,0],Y[0,0])))
-    hds.loc[:,"inode"] = hds.pt.apply(lambda x: i.intersect(x)[0][0])
+    # print(i.intersect(shapely.geometry.Point(X[0,0],Y[0,0])))
+    hds.loc[:, "inode"] = hds.pt.apply(lambda x: i.intersect(x)[0][0])
     print(hds.inode)
 
-    hds.loc[:,"site"] = hds.apply(lambda x: "twgw_{0}_{1}".format(x.layer-1,x.inode),axis=1)
-    hds.loc[:,["site","otype","layer","inode"]].to_csv(os.path.join(new_d,hobs_csv),index=False,header=False)
-    shutil.copy2(os.path.join(mn_dir,"freyberg6.obs"),os.path.join(new_d,"freyberg6.obs"))
-    lines = open(os.path.join(new_d,"freyberg6.nam"),'r').readlines()
-    with open(os.path.join(new_d,"freyberg6.nam"),'w') as f:
+    hds.loc[:, "site"] = hds.apply(
+        lambda x: "twgw_{0}_{1}".format(x.layer - 1, x.inode), axis=1
+    )
+    hds.loc[:, ["site", "otype", "layer", "inode"]].to_csv(
+        os.path.join(new_d, hobs_csv), index=False, header=False
+    )
+    shutil.copy2(
+        os.path.join(mn_dir, "freyberg6.obs"), os.path.join(new_d, "freyberg6.obs")
+    )
+    lines = open(os.path.join(new_d, "freyberg6.nam"), "r").readlines()
+    with open(os.path.join(new_d, "freyberg6.nam"), "w") as f:
         for line in lines:
-            
             if "end packages" in line.lower():
                 f.write("obs6  freyberg6.obs  head_obs\n")
                 f.write("sto6  freyberg6.sto\n")
             f.write(line)
 
-    sfr_file = os.path.join(new_d,"freyberg6.sfr_packagedata.txt")
-    df = pd.read_csv(sfr_file,delim_whitespace=True,header=None)
-    df.loc[:,1] = 1 # reset layer to 1
-    #botm = np.loadtxt(os.path.join(new_d,"freyberg6.disv_botm_layer1.txt"))
-    df.loc[:,5] = 5e-5 #gradient
-    df.loc[:,4] = 5 # width
-    cumdist = np.cumsum(df.loc[:,4].values)
+    sfr_file = os.path.join(new_d, "freyberg6.sfr_packagedata.txt")
+    df = pd.read_csv(sfr_file, delim_whitespace=True, header=None)
+    df.loc[:, 1] = 1  # reset layer to 1
+    # botm = np.loadtxt(os.path.join(new_d,"freyberg6.disv_botm_layer1.txt"))
+    df.loc[:, 5] = 5e-5  # gradient
+    df.loc[:, 4] = 5  # width
+    cumdist = np.cumsum(df.loc[:, 4].values)
     up = 34.0
     dn = 33.501
     totdrop = up - dn
     totlen = cumdist[-1]
-    drop_grad = totdrop/totlen
+    drop_grad = totdrop / totlen
     drop_reach = up - (cumdist * drop_grad)
     print(drop_reach)
-    df.loc[:,6] = drop_reach #reach botm elev
-    df.loc[:,7] = 1.0 #thickness of stream bed
-    df.loc[:,8] = 0.1 #stream hk
-    df.loc[:,9] = 0.3 #mannings
+    df.loc[:, 6] = drop_reach  # reach botm elev
+    df.loc[:, 7] = 1.0  # thickness of stream bed
+    df.loc[:, 8] = 0.1  # stream hk
+    df.loc[:, 9] = 0.3  # mannings
     print(df)
-    df.to_csv(os.path.join(new_d,"freyberg6.sfr_packagedata.txt"),index=False,header=False,sep=" ")
+    df.to_csv(
+        os.path.join(new_d, "freyberg6.sfr_packagedata.txt"),
+        index=False,
+        header=False,
+        sep=" ",
+    )
 
-    df = pd.read_csv(os.path.join(new_d,"freyberg6.chd_stress_period_data_1.txt"),header=None,delim_whitespace=True)
-    df.loc[:,2] = 33.5
-    df.to_csv(os.path.join(new_d,"freyberg6.chd_stress_period_data_1.txt"),index=False,header=False)
+    df = pd.read_csv(
+        os.path.join(new_d, "freyberg6.chd_stress_period_data_1.txt"),
+        header=None,
+        delim_whitespace=True,
+    )
+    df.loc[:, 2] = 33.5
+    df.to_csv(
+        os.path.join(new_d, "freyberg6.chd_stress_period_data_1.txt"),
+        index=False,
+        header=False,
+    )
 
     for k in range(3):
-        id_name = os.path.join(new_d,"freyberg6.disv_idomain_layer{0}.txt".format(k+1))
+        id_name = os.path.join(
+            new_d, "freyberg6.disv_idomain_layer{0}.txt".format(k + 1)
+        )
         arr = np.loadtxt(id_name)
-        arr[arr>0] = 1
-        np.savetxt(id_name,arr,fmt="%2.0f")
-
+        arr[arr > 0] = 1
+        np.savetxt(id_name, arr, fmt="%2.0f")
 
-    pyemu.os_utils.run("mf6",cwd=new_d)
+    pyemu.os_utils.run("mf6", cwd=new_d)
 
-    #now write an "obs" twgw file with xy layer info
+    # now write an "obs" twgw file with xy layer info
     mnm.modelgrid.set_coord_info(angrot=-55)
-    X,Y = mnm.modelgrid.xcellcenters,mnm.modelgrid.ycellcenters
-    hds.loc[:,"x"] = hds.apply(lambda x: X[x.row-1,x.col-1],axis=1)
-    hds.loc[:,"y"] = hds.apply(lambda x: Y[x.row-1,x.col-1],axis=1)
+    X, Y = mnm.modelgrid.xcellcenters, mnm.modelgrid.ycellcenters
+    hds.loc[:, "x"] = hds.apply(lambda x: X[x.row - 1, x.col - 1], axis=1)
+    hds.loc[:, "y"] = hds.apply(lambda x: Y[x.row - 1, x.col - 1], axis=1)
 
-    obsdf = pd.read_csv(os.path.join(new_d,"heads.csv"),index_col=0)
+    obsdf = pd.read_csv(os.path.join(new_d, "heads.csv"), index_col=0)
     obsdf.columns = [c.lower() for c in obsdf.columns]
-    obsdf = obsdf.melt(ignore_index=False,var_name="site",value_name="obsval")
-    start_datetime= pd.to_datetime("1-1-2018")
-    obsdf.index = start_datetime + pd.to_timedelta(obsdf.index.values,unit="d")
+    obsdf = obsdf.melt(ignore_index=False, var_name="site", value_name="obsval")
+    start_datetime = pd.to_datetime("1-1-2018")
+    obsdf.index = start_datetime + pd.to_timedelta(obsdf.index.values, unit="d")
     obsdf.index.name = "datetime"
     hds.index = hds.site.values
-    noisex = np.random.uniform(-40,40,hds.x.shape[0])
+    noisex = np.random.uniform(-40, 40, hds.x.shape[0])
     xdict = (hds.x + noisex).to_dict()
-    noisey = np.random.uniform(-40,40,hds.y.shape[0])
+    noisey = np.random.uniform(-40, 40, hds.y.shape[0])
     ydict = (hds.y + noisey).to_dict()
     layerdict = hds.layer.to_dict()
-    obsdf.loc[:,"x"] = obsdf.site.apply(lambda x: xdict[x])
-    obsdf.loc[:,"y"] = obsdf.site.apply(lambda x: ydict[x])
-    obsdf.loc[:,"layer"] = obsdf.site.apply(lambda x: layerdict[x])
-    #now remove some values
-    keep_idx = np.arange(0,obsdf.shape[0],dtype=int)
+    obsdf.loc[:, "x"] = obsdf.site.apply(lambda x: xdict[x])
+    obsdf.loc[:, "y"] = obsdf.site.apply(lambda x: ydict[x])
+    obsdf.loc[:, "layer"] = obsdf.site.apply(lambda x: layerdict[x])
+    # now remove some values
+    keep_idx = np.arange(0, obsdf.shape[0], dtype=int)
     np.random.shuffle(keep_idx)
-    keep_idx = keep_idx[:int(3*keep_idx.shape[0]/4)]
-    kobsdf = obsdf.iloc[keep_idx,:].copy()
+    keep_idx = keep_idx[: int(3 * keep_idx.shape[0] / 4)]
+    kobsdf = obsdf.iloc[keep_idx, :].copy()
     kobsdf.sort_index(inplace=True)
-    print(obsdf.shape,kobsdf.shape)
+    print(obsdf.shape, kobsdf.shape)
     print(kobsdf)
     # now some datetime noise
-    noise = np.random.uniform(-30,30,kobsdf.shape[0])
-    noise = np.round(noise,1)
-    td_noise = pd.to_timedelta(noise,unit='d')
+    noise = np.random.uniform(-30, 30, kobsdf.shape[0])
+    noise = np.round(noise, 1)
+    td_noise = pd.to_timedelta(noise, unit="d")
     kobsdf.index = kobsdf.index + td_noise
-    kobsdf.loc[:,"datetime"] = kobsdf.index.values
-    kobsdf.loc[:,["site","datetime","x","y","layer","obsval"]].to_csv(os.path.join("freyberg_aux_files","gwlevel_obs.csv"),index=False)
+    kobsdf.loc[:, "datetime"] = kobsdf.index.values
+    kobsdf.loc[:, ["site", "datetime", "x", "y", "layer", "obsval"]].to_csv(
+        os.path.join("freyberg_aux_files", "gwlevel_obs.csv"), index=False
+    )
 
     print(kobsdf)
 
+
 def invest_quadtree_results():
     import matplotlib.pyplot as plt
+
     ws = "freyberg_quadtree"
     sim = flopy.mf6.MFSimulation.load(sim_ws=ws)
     m = sim.get_model()
-    hds = flopy.utils.HeadFile(os.path.join(ws,"freyberg6.hds"),model=m)
+    hds = flopy.utils.HeadFile(os.path.join(ws, "freyberg6.hds"), model=m)
     hds.plot(colorbar=True)
     plt.show()
 
+
 if __name__ == "__main__":
     repair_and_prep_quadtree_model()
-    #invest_quadtree_results()
+    # invest_quadtree_results()
diff --git a/pypestutils/ctypes_declarations.py b/pypestutils/ctypes_declarations.py
index 1d7cf46..0f0c98a 100644
--- a/pypestutils/ctypes_declarations.py
+++ b/pypestutils/ctypes_declarations.py
@@ -1,4 +1,5 @@
 """Low-level Fortran-Python ctypes functions."""
+
 from __future__ import annotations
 
 from ctypes import ARRAY, CDLL, POINTER, c_char, c_double, c_int
diff --git a/pypestutils/data.py b/pypestutils/data.py
index 0521e97..374f4bb 100644
--- a/pypestutils/data.py
+++ b/pypestutils/data.py
@@ -1,4 +1,5 @@
 """Data module."""
+
 from __future__ import annotations
 
 from enum import Enum
@@ -157,8 +158,7 @@ def __init__(
             if name in self._names:
                 raise KeyError(f"'{name}' defined more than once")
             self._names.append(name)
-            float_any[name] = ar = np.asarray(
-                float_any[name], np.float64, order="F")
+            float_any[name] = ar = np.asarray(float_any[name], np.float64, order="F")
             if not self.shape and ar.ndim == 1:
                 self.shape = ar.shape
         for name in int_any.keys():
diff --git a/pypestutils/enum.py b/pypestutils/enum.py
index 939cf81..2303ff1 100644
--- a/pypestutils/enum.py
+++ b/pypestutils/enum.py
@@ -1,4 +1,5 @@
 """Enumeration module."""
+
 from __future__ import annotations
 
 from enum import IntEnum
diff --git a/pypestutils/finder.py b/pypestutils/finder.py
index a00dae6..b80353e 100644
--- a/pypestutils/finder.py
+++ b/pypestutils/finder.py
@@ -1,4 +1,5 @@
 """Locate pestutils shared library by any means necessary."""
+
 import ctypes
 import os
 import platform
diff --git a/pypestutils/helpers.py b/pypestutils/helpers.py
index 181b4e3..8335bc4 100644
--- a/pypestutils/helpers.py
+++ b/pypestutils/helpers.py
@@ -7,10 +7,19 @@
 from .pestutilslib import PestUtilsLib
 
 
-def mod2obs_mf6(gridinfo_fname: str,depvar_fname: str,obscsv_fname: str ,model_type: int,start_datetime: str | pd.TimeStamp,depvar_ftype=1,
-                depvar_name="head",interp_thresh=1.0e+30,no_interp_val=1.0e+30,model_timeunit="d",
-                time_extrap=1.0)->dict:
-
+def mod2obs_mf6(
+    gridinfo_fname: str,
+    depvar_fname: str,
+    obscsv_fname: str,
+    model_type: int,
+    start_datetime: str | pd.TimeStamp,
+    depvar_ftype=1,
+    depvar_name="head",
+    interp_thresh=1.0e30,
+    no_interp_val=1.0e30,
+    model_timeunit="d",
+    time_extrap=1.0,
+) -> dict:
     """python implementation of mod2smp and mod2obs using modflow6 binary grid files
     Parameters
     ----------
@@ -45,8 +54,8 @@ def mod2obs_mf6(gridinfo_fname: str,depvar_fname: str,obscsv_fname: str ,model_t
         temporally interpolated simulated results at observation locations (ie mod2obs)
     """
 
-    for fname in [gridinfo_fname,depvar_fname]:
-        assert os.path.exists(fname),"file {0} not found".format(fname)
+    for fname in [gridinfo_fname, depvar_fname]:
+        assert os.path.exists(fname), "file {0} not found".format(fname)
     lib = PestUtilsLib()
     is_mf6 = False
     is_structured = True
@@ -69,68 +78,109 @@ def mod2obs_mf6(gridinfo_fname: str,depvar_fname: str,obscsv_fname: str ,model_t
         raise Exception("unrecognized 'model_type':{0}".format(model_type))
 
     depvar_ftype = int(depvar_ftype)
-    if depvar_ftype not in [1,2]:
+    if depvar_ftype not in [1, 2]:
         raise Exception("unrecognized 'depvar_ftype':{0}".format(depvar_ftype))
 
     if is_mf6:
-        grid_info = lib.install_mf6_grid_from_file("grid",gridinfo_fname)
+        grid_info = lib.install_mf6_grid_from_file("grid", gridinfo_fname)
     else:
         raise NotImplementedError()
-    
-    if isinstance(start_datetime,str):
+
+    if isinstance(start_datetime, str):
         start_datetime = pd.to_datetime(start_datetime)
 
-    depvar_info = lib.inquire_modflow_binary_file_specs(depvar_fname,depvar_fname+".out.csv",model_type,depvar_ftype)    
-    depvar_df = pd.read_csv(depvar_fname+".out.csv")
+    depvar_info = lib.inquire_modflow_binary_file_specs(
+        depvar_fname, depvar_fname + ".out.csv", model_type, depvar_ftype
+    )
+    depvar_df = pd.read_csv(depvar_fname + ".out.csv")
     depvar_df.columns = [c.lower() for c in depvar_df.columns]
-    #print(depvar_df)
+    # print(depvar_df)
 
-    if isinstance(obscsv_fname,str):
+    if isinstance(obscsv_fname, str):
         if not os.path.exists(obscsv_fname):
             raise Exception("obscsv_fname '{0}' not found".format(obscsv_fname))
         # todo: think about supporting a site sample file maybe?
-        obsdf = pd.read_csv(os.path.join(obscsv_fname),parse_dates=["datetime"])
-    elif isinstance(obscsv_fname,pd.DataFrame):
+        obsdf = pd.read_csv(os.path.join(obscsv_fname), parse_dates=["datetime"])
+    elif isinstance(obscsv_fname, pd.DataFrame):
         obsdf = obscsv_fname.copy()
     else:
-        raise Exception("obscsv arg type not recognized (looking for str or pd.DataFrame):'{0}'".format(type(obscsv_fname)))
-    #check obsdf
+        raise Exception(
+            "obscsv arg type not recognized (looking for str or pd.DataFrame):'{0}'".format(
+                type(obscsv_fname)
+            )
+        )
+    # check obsdf
     obsdf.columns = [c.lower() for c in obsdf.columns]
-    for req_col in ["site","x","y","datetime","layer"]:
+    for req_col in ["site", "x", "y", "datetime", "layer"]:
         if req_col not in obsdf.columns:
-            raise Exception("observation dataframe missing column '{0}'".format(req_col))
+            raise Exception(
+                "observation dataframe missing column '{0}'".format(req_col)
+            )
     usitedf = obsdf.groupby("site").first()
     pth = os.path.split(depvar_fname)[0]
-    fac_file = os.path.join(pth,"obs_interp_fac.bin")
-    bln_file = fac_file.replace(".bin",".bln")
-    interp_fac_results = lib.calc_mf6_interp_factors("grid",usitedf.x.values,usitedf.y.values,usitedf.layer.values,fac_file,"binary",bln_file)
+    fac_file = os.path.join(pth, "obs_interp_fac.bin")
+    bln_file = fac_file.replace(".bin", ".bln")
+    interp_fac_results = lib.calc_mf6_interp_factors(
+        "grid",
+        usitedf.x.values,
+        usitedf.y.values,
+        usitedf.layer.values,
+        fac_file,
+        "binary",
+        bln_file,
+    )
     if 0 in interp_fac_results:
-        print("warning: the following site(s) failed to have interpolation factors calculated:")
-        fsites = usitedf.reset_index().site.iloc[interp_fac_results==0].to_list()
+        print(
+            "warning: the following site(s) failed to have interpolation factors calculated:"
+        )
+        fsites = usitedf.reset_index().site.iloc[interp_fac_results == 0].to_list()
         print(fsites)
-    all_results = lib.interp_from_mf6_depvar_file(depvar_fname,fac_file,"binary",depvar_info["ntime"],depvar_name,interp_thresh,True,
-        no_interp_val,usitedf.shape[0])
-    datetimes = start_datetime+pd.to_timedelta(all_results["simtime"],unit=model_timeunit)
-    allresults_df = pd.DataFrame(all_results["simstate"],index=datetimes,columns=usitedf.index)
-    allresults_df.to_csv(depvar_fname+".all.csv")
+    all_results = lib.interp_from_mf6_depvar_file(
+        depvar_fname,
+        fac_file,
+        "binary",
+        depvar_info["ntime"],
+        depvar_name,
+        interp_thresh,
+        True,
+        no_interp_val,
+        usitedf.shape[0],
+    )
+    datetimes = start_datetime + pd.to_timedelta(
+        all_results["simtime"], unit=model_timeunit
+    )
+    allresults_df = pd.DataFrame(
+        all_results["simstate"], index=datetimes, columns=usitedf.index
+    )
+    allresults_df.to_csv(depvar_fname + ".all.csv")
 
     if "totim" in obsdf:
         print("WARNING: replacing existing 'totim' column in observation dataframe")
-    obsdf.loc[:,"totim"] = obsdf.datetime.apply(lambda x: x  - start_datetime).dt.days 
+    obsdf.loc[:, "totim"] = obsdf.datetime.apply(lambda x: x - start_datetime).dt.days
 
     usite = obsdf.site.unique()
     usite.sort()
-    usite_dict = {s:c for s,c in zip(usite,np.arange(usite.shape[0],dtype=int))}
-    obsdf.loc[:,"isite"] = obsdf.site.apply(lambda x: usite_dict[x])
-    obsdf.sort_values(by=["isite","totim"],inplace=True)
-    
-    interp_results = lib.interp_to_obstime(all_results["nproctime"],all_results["simtime"],all_results["simstate"],interp_thresh,"L",
-        time_extrap,no_interp_val,obsdf.isite.values,obsdf.totim.values)
-
-    obsdf.loc[:,"simulated"] = interp_results
-    lib.uninstall_mf6_grid('grid')
+    usite_dict = {s: c for s, c in zip(usite, np.arange(usite.shape[0], dtype=int))}
+    obsdf.loc[:, "isite"] = obsdf.site.apply(lambda x: usite_dict[x])
+    obsdf.sort_values(by=["isite", "totim"], inplace=True)
+
+    interp_results = lib.interp_to_obstime(
+        all_results["nproctime"],
+        all_results["simtime"],
+        all_results["simstate"],
+        interp_thresh,
+        "L",
+        time_extrap,
+        no_interp_val,
+        obsdf.isite.values,
+        obsdf.totim.values,
+    )
+
+    obsdf.loc[:, "simulated"] = interp_results
+    lib.uninstall_mf6_grid("grid")
     lib.free_all_memory()
-    return {"all_results":allresults_df,"interpolated_results":obsdf}
+    return {"all_results": allresults_df, "interpolated_results": obsdf}
+
 
 def get_grid_info_from_gridspec(gridspec_fname: str) -> dict:
     """Read structured grid info from a PEST-style grid specification file
@@ -138,7 +188,7 @@ def get_grid_info_from_gridspec(gridspec_fname: str) -> dict:
     ----------
     gridspec_fname : str
         PEST-style grid specification file
-    
+
     Returns
     -------
     grid_info: dict
@@ -155,7 +205,7 @@ def get_grid_info_from_gridspec(gridspec_fname: str) -> dict:
         "nrow": sr.nrow,
         "ncol": sr.ncol,
         "delr": sr.delr,
-        "delc": sr.delc
+        "delc": sr.delc,
     }
 
 
@@ -165,7 +215,7 @@ def get_grid_info_from_mf6_grb(grb_fname: str) -> dict:
     ----------
     grb_fname: str
         MODFLOW-6 binary grid file
-    
+
     Returns
     -------
     grid_info: dict
@@ -174,8 +224,8 @@ def get_grid_info_from_mf6_grb(grb_fname: str) -> dict:
     if not os.path.exists(grb_fname):
         raise FileNotFoundError(grb_fname)
     lib = PestUtilsLib()
-    data = lib.install_mf6_grid_from_file("grid",grb_fname)
-    data["x"],data["y"],data["z"] = lib.get_cell_centres_mf6("grid",data["ncells"])
+    data = lib.install_mf6_grid_from_file("grid", grb_fname)
+    data["x"], data["y"], data["z"] = lib.get_cell_centres_mf6("grid", data["ncells"])
     lib.uninstall_mf6_grid("grid")
     lib.free_all_memory()
     return data
@@ -191,22 +241,22 @@ def get_2d_grid_info_from_file(fname: str, layer=None) -> dict:
     layer: int (optional)
         the (one-based) layer number to use for 2-D.  If None and
         grid info is 3-D, a value of 1 is used
-    
+
     Returns
     -------
     grid_info: dict
         grid information
-    """ 
+    """
 
     grid_info = None
-    if isinstance(fname,str):
+    if isinstance(fname, str):
         if not os.path.exists(fname):
             raise FileNotFoundError(fname)
         if fname.lower().endswith(".csv"):
             grid_info = pd.read_csv(fname)
             grid_info.columns = [c.lower() for c in grid_info.columns]
-            fname = grid_info # for  checks and processing below
-            
+            fname = grid_info  # for  checks and processing below
+
         else:
             try:
                 grid_info = get_grid_info_from_gridspec(fname)
@@ -214,23 +264,26 @@ def get_2d_grid_info_from_file(fname: str, layer=None) -> dict:
                 try:
                     grid_info = get_2d_grid_info_from_mf6_grb(fname, layer=layer)
                 except Exception as e2:
-                    
-                    raise Exception("error getting grid info from file '{0}'".format(fname))
-        
-    if isinstance(fname,pd.DataFrame):
-        if 'x' not in fname.columns:
+                    raise Exception(
+                        "error getting grid info from file '{0}'".format(fname)
+                    )
+
+    if isinstance(fname, pd.DataFrame):
+        if "x" not in fname.columns:
             raise Exception("required 'x' column not found in grid info dataframe")
-        if 'y' not in fname.columns:
+        if "y" not in fname.columns:
             raise Exception("required 'y' column not found in grid info dataframe")
-        if layer is not None and 'layer' not in fname.columns:
-            print("WARNING: 'layer' arg is not None but 'layer' not found in grid info dataframe...")
+        if layer is not None and "layer" not in fname.columns:
+            print(
+                "WARNING: 'layer' arg is not None but 'layer' not found in grid info dataframe..."
+            )
         # I think these should just be references to column values (not copies)
-        grid_info = {c:fname[c].values for c in fname.columns}
-    
+        grid_info = {c: fname[c].values for c in fname.columns}
+
     return grid_info
 
 
-def get_2d_grid_info_from_mf6_grb(grb_fname: str,layer=None) -> dict:
+def get_2d_grid_info_from_mf6_grb(grb_fname: str, layer=None) -> dict:
     """Read grid info from a MODFLOW-6 binary grid file
     Parameters
     ----------
@@ -239,7 +292,7 @@ def get_2d_grid_info_from_mf6_grb(grb_fname: str,layer=None) -> dict:
     layer: int (optional)
         the layer number to use for 2-D.  If None,
         a value of 1 is used
-    
+
     Returns
     -------
     grid_info: dict
@@ -249,19 +302,23 @@ def get_2d_grid_info_from_mf6_grb(grb_fname: str,layer=None) -> dict:
     nnodes = grid_info["ncells"]
     x = grid_info["x"].copy()
     y = grid_info["y"].copy()
-    nrow,ncol = None,None
+    nrow, ncol = None, None
     if grid_info["idis"] == 1:
         nlay = grid_info["ndim3"]
         if layer is not None:
             assert layer != 0, "Value provided for layer should be 1 based, sorry"
             if layer > nlay:
-                raise Exception("user-supplied 'layer' {0} greater than nlay {1}".format(layer,nlay))
+                raise Exception(
+                    "user-supplied 'layer' {0} greater than nlay {1}".format(
+                        layer, nlay
+                    )
+                )
         else:
             layer = 1
         nrow = grid_info["ndim2"]
         ncol = grid_info["ndim1"]
-        x = x.reshape((nlay,nrow,ncol))[layer-1]
-        y = y.reshape((nlay,nrow,ncol))[layer-1]
+        x = x.reshape((nlay, nrow, ncol))[layer - 1]
+        y = y.reshape((nlay, nrow, ncol))[layer - 1]
         grid_info["nnodes"] = nrow * ncol
         grid_info["x"] = x
         grid_info["y"] = y
@@ -272,12 +329,16 @@ def get_2d_grid_info_from_mf6_grb(grb_fname: str,layer=None) -> dict:
         nlay = grid_info["ndim3"]
         if layer is not None:
             if layer > nlay:
-                raise Exception("user-supplied 'layer' {0} greater than nlay {1}".format(layer,nlay))
+                raise Exception(
+                    "user-supplied 'layer' {0} greater than nlay {1}".format(
+                        layer, nlay
+                    )
+                )
         else:
             layer = 1
         ncpl = grid_info["ndim1"]
-        x = x.reshape((nlay,ncpl))[layer-1]
-        y = y.reshape((nlay,ncpl))[layer-1]
+        x = x.reshape((nlay, ncpl))[layer - 1]
+        y = y.reshape((nlay, ncpl))[layer - 1]
         grid_info["nnodes"] = ncpl
         grid_info["x"] = x
         grid_info["y"] = y
@@ -285,12 +346,9 @@ def get_2d_grid_info_from_mf6_grb(grb_fname: str,layer=None) -> dict:
 
 
 def get_2d_pp_info_structured_grid(
-    pp_space: int,
-    gridinfo_fname: str,
-    array_dict = {},
-    name_prefix="pp"
+    pp_space: int, gridinfo_fname: str, array_dict={}, name_prefix="pp"
 ) -> pd.DataFrame:
-    """Create a grid of pilot point locations for a 
+    """Create a grid of pilot point locations for a
     2-D structured grid
     Parameters
     ----------
@@ -299,14 +357,14 @@ def get_2d_pp_info_structured_grid(
     gridinfo_fname: str
         file contain grid information
     array_dict: dict (optional)
-        a dict of 2-D grid-shape arrays used to populate 
+        a dict of 2-D grid-shape arrays used to populate
         pilot point attributes.  Special values include:
-        "value","zone","bearing","aniso" and "corrlen", 
+        "value","zone","bearing","aniso" and "corrlen",
         although any number of arrays can be passed and will
         sampled at pilot point locations
     name_prefix: str
         pilot point name prefix. Default is "pp"
-    
+
     Returns
     -------
     ppdf: pd.DataaFrame
@@ -317,29 +375,31 @@ def get_2d_pp_info_structured_grid(
     grid_info = get_2d_grid_info_from_file(gridinfo_fname)
     pname, px, py, pval = [], [], [], []
     pi, pj = [], []
-    parr_dict = {k:[] for k in array_dict.keys()}
+    parr_dict = {k: [] for k in array_dict.keys()}
     count = 0
     nrow = grid_info["nrow"]
     ncol = grid_info["ncol"]
-    nlay = grid_info.get("nlay",1)
+    nlay = grid_info.get("nlay", 1)
 
-    zone_array = array_dict.get("zone",None)
+    zone_array = array_dict.get("zone", None)
 
-    x = grid_info['x']
-    y = grid_info['y']
-    x = x.reshape((nlay,nrow,ncol))[0,:,:]
-    y = y.reshape((nlay,nrow,ncol))[0,:,:]
+    x = grid_info["x"]
+    y = grid_info["y"]
+    x = x.reshape((nlay, nrow, ncol))[0, :, :]
+    y = y.reshape((nlay, nrow, ncol))[0, :, :]
     if nrow is None:
-        raise Exception("unstructured grid loaded from gridinfo_fname '{0}'".format(gridspec_fname))
+        raise Exception(
+            "unstructured grid loaded from gridinfo_fname '{0}'".format(gridspec_fname)
+        )
     for i in range(int(pp_space / 2), nrow, pp_space):
         for j in range(int(pp_space / 2), ncol, pp_space):
             if zone_array is not None and zone_array[i, j] <= 0:
                 continue
             px.append(x[i, j])
             py.append(y[i, j])
-            #if zone_array is not None:
+            # if zone_array is not None:
             #    pzone.append(zone_array[i, j])
-            #else:
+            # else:
             #    pzone.append(1)
 
             pname.append(name_prefix + "{0}".format(count))
@@ -356,15 +416,15 @@ def get_2d_pp_info_structured_grid(
         },
         index=pname,
     )
-    df.loc[:,"value"] = 1.0
+    df.loc[:, "value"] = 1.0
     df.loc[:, "bearing"] = 0.0
     df.loc[:, "aniso"] = 1.0
     delx = pp_space * 5 * int((x.max() - x.min()) / float(ncol))
     dely = pp_space * 5 * int((y.max() - y.min()) / float(nrow))
-    df.loc[:, "corrlen"] = max(delx,dely)  # ?
-    df.loc[:,"zone"] = 1
-    for k,arr in array_dict.items():
-        df.loc[:,k] = arr[df.i,df.j]
+    df.loc[:, "corrlen"] = max(delx, dely)  # ?
+    df.loc[:, "zone"] = 1
+    for k, arr in array_dict.items():
+        df.loc[:, k] = arr[df.i, df.j]
     df["zone"] = df.zone.astype(int)
 
     return df
@@ -381,16 +441,16 @@ def interpolate_with_sva_pilotpoints_2d(
     search_dist=1e30,
     zone_array=1,
     verbose=True,
-    layer=None
+    layer=None,
 ) -> dict:
     """Perform 2-D pilot point interpolation using
     spatially varying geostatistical hyper-parameters
     Parameters
     ----------
     pp_info: pandas.DataFrame
-        dataframe with pilot point info.  Required columns 
+        dataframe with pilot point info.  Required columns
         include: "x","y",and "value".  optional columns include:
-        "zone","bearing","aniso",and "corrlen"    
+        "zone","bearing","aniso",and "corrlen"
     gridinfo_fname: str
         file name storing grid information
     vartype: str
@@ -409,7 +469,7 @@ def interpolate_with_sva_pilotpoints_2d(
         search distance to use when looking for nearby pilot points.
         Default is 1.0e+30
     zone_array: int | numpy.ndarray
-        the zone array to match up with "zone" value in `pp_info`.  If 
+        the zone array to match up with "zone" value in `pp_info`.  If
         integer type, a constant zone array of value "zone_array" is used.
         Default is 1
     verbose: bool
@@ -422,7 +482,7 @@ def interpolate_with_sva_pilotpoints_2d(
     Returns
     -------
     results: dict
-        resulting arrays of the various interpolation from pilot 
+        resulting arrays of the various interpolation from pilot
         points to grid-shaped arrays
     """
     # some checks on pp_info
@@ -447,13 +507,13 @@ def interpolate_with_sva_pilotpoints_2d(
     nrow, ncol = None, None
     x, y, area = None, None, None
     grid_info = get_2d_grid_info_from_file(gridinfo_fname, layer)
-    nrow = grid_info.get("nrow",None)
-    ncol = grid_info.get("ncol",None)
-    x = grid_info['x']
-    y = grid_info['y']
-    area = grid_info.get("area",None)
-    idis = grid_info.get("idis",None)
-    nnodes = grid_info.get("nnodes",None)
+    nrow = grid_info.get("nrow", None)
+    ncol = grid_info.get("ncol", None)
+    x = grid_info["x"]
+    y = grid_info["y"]
+    area = grid_info.get("area", None)
+    idis = grid_info.get("idis", None)
+    nnodes = grid_info.get("nnodes", None)
     if area is None:
         area = np.ones_like(x)
     if nnodes is None:
@@ -468,8 +528,7 @@ def interpolate_with_sva_pilotpoints_2d(
         # TODO warn here
         lib.logger.info("casting zone_array from %r to int", zone_array.dtype)
         zone_array = zone_array.astype(int)
-    
-    
+
     hyperfac_ftype = "binary"
     if verbose:
         hyperfac_ftype = "text"
@@ -480,15 +539,26 @@ def interpolate_with_sva_pilotpoints_2d(
     hypertrans = "none"
 
     fac_files = []
-    lib.logger.info("using bearing of %r and aniso of %r for hyperpar interpolation", hyperbearing, hyperaniso)
-    lib.logger.info("using %r variogram with %r transform for hyperpar interpolation",hypervartype,hypertrans)
-       
-    results = {}    
+    lib.logger.info(
+        "using bearing of %r and aniso of %r for hyperpar interpolation",
+        hyperbearing,
+        hyperaniso,
+    )
+    lib.logger.info(
+        "using %r variogram with %r transform for hyperpar interpolation",
+        hypervartype,
+        hypertrans,
+    )
+
+    results = {}
 
     bearing = np.zeros_like(x)
     if "bearing" in pp_info.columns:
         hypernoint = pp_info.bearing.mean()
-        lib.logger.info("using no-interpolation value of %r for 'bearing' hyperpar interpolation", hypernoint)
+        lib.logger.info(
+            "using no-interpolation value of %r for 'bearing' hyperpar interpolation",
+            hypernoint,
+        )
         hyperfac_fname = "tempbearing.fac"
         npts = lib.calc_kriging_factors_auto_2d(
             pp_info.x.values,
@@ -517,17 +587,19 @@ def interpolate_with_sva_pilotpoints_2d(
         fac_files.append(hyperfac_fname)
         if verbose:
             if nrow is not None:
-                np.savetxt("bearing.txt",bearing.reshape(nrow,ncol),fmt="%15.6E")
+                np.savetxt("bearing.txt", bearing.reshape(nrow, ncol), fmt="%15.6E")
             else:
-                np.savetxt("bearing.txt",bearing,fmt="%15.6E")
+                np.savetxt("bearing.txt", bearing, fmt="%15.6E")
     results["bearing"] = bearing
-    
 
     aniso = np.ones_like(x)
     if "aniso" in pp_info.columns:
         hypernoint = pp_info.aniso.mean()
-        lib.logger.info("using no-interpolation value of %r for 'aniso' hyperpar interpolation", hypernoint)
-        
+        lib.logger.info(
+            "using no-interpolation value of %r for 'aniso' hyperpar interpolation",
+            hypernoint,
+        )
+
         hyperfac_fname = "tempaniso.fac"
         npts = lib.calc_kriging_factors_auto_2d(
             pp_info.x.values,
@@ -555,9 +627,9 @@ def interpolate_with_sva_pilotpoints_2d(
         aniso = result["targval"]
         if verbose:
             if nrow is not None:
-                np.savetxt("aniso.txt",aniso.reshape(nrow,ncol),fmt="%15.6E")
+                np.savetxt("aniso.txt", aniso.reshape(nrow, ncol), fmt="%15.6E")
             else:
-                np.savetxt("aniso.txt",aniso,fmt="%15.6E")
+                np.savetxt("aniso.txt", aniso, fmt="%15.6E")
 
     results["aniso"] = aniso
 
@@ -565,7 +637,10 @@ def interpolate_with_sva_pilotpoints_2d(
     corrlen = None  # todo corrlen default where not in ppdf
     if "corrlen" in pp_info.columns:
         hypernoint = pp_info.corrlen.mean()
-        lib.logger.info("using no-interpolation value of %r for 'corrlen' hyperpar interpolation", hypernoint)
+        lib.logger.info(
+            "using no-interpolation value of %r for 'corrlen' hyperpar interpolation",
+            hypernoint,
+        )
         hyperfac_fname = "tempcorrlen.fac"
         npts = lib.calc_kriging_factors_auto_2d(
             pp_info.x.values,
@@ -595,9 +670,9 @@ def interpolate_with_sva_pilotpoints_2d(
         use_auto = False
         if verbose:
             if nrow is not None:
-                np.savetxt("corrlen.txt",corrlen.reshape(nrow,ncol),fmt="%15.6E")
+                np.savetxt("corrlen.txt", corrlen.reshape(nrow, ncol), fmt="%15.6E")
             else:
-                np.savetxt("corrlen.txt",corrlen,fmt="%15.6E")
+                np.savetxt("corrlen.txt", corrlen, fmt="%15.6E")
         results["corrlen"] = corrlen
 
     if not verbose:
@@ -660,8 +735,8 @@ def interpolate_with_sva_pilotpoints_2d(
     results["result"] = result["targval"]
 
     if nrow is not None:
-        for k,v in results.items():
-            results[k] = v.reshape(nrow,ncol)
+        for k, v in results.items():
+            results[k] = v.reshape(nrow, ncol)
 
     return results
 
@@ -679,9 +754,9 @@ def generate_2d_grid_realizations(
     variobearing=0.0,
     random_seed=12345,
     layer=None,
-) ->np.NDArray[float]:
-    """draw 2-D realizations using sequential gaussian 
-    simulations and optionally using spatially varying 
+) -> np.NDArray[float]:
+    """draw 2-D realizations using sequential gaussian
+    simulations and optionally using spatially varying
     geostatistical hyper parameters.
     Parameters
     ----------
@@ -718,24 +793,21 @@ def generate_2d_grid_realizations(
         will be reshaped to NROW X NCOL)
     """
 
-    
-
     nrow, ncol = None, None
     x, y, area = None, None, None
     grid_info = get_2d_grid_info_from_file(gridinfo_fname, layer)
-    nrow = grid_info.get("nrow",None)
-    ncol = grid_info.get("ncol",None)
-    x = grid_info['x']
-    y = grid_info['y']
-    area = grid_info.get("area",None)
-    idis = grid_info.get("idis",None)
-    nnodes = grid_info.get("nnodes",None)
+    nrow = grid_info.get("nrow", None)
+    ncol = grid_info.get("ncol", None)
+    x = grid_info["x"]
+    y = grid_info["y"]
+    area = grid_info.get("area", None)
+    idis = grid_info.get("idis", None)
+    nnodes = grid_info.get("nnodes", None)
     if area is None:
         area = np.ones_like(x)
     if nnodes is None:
         nnodes = x.shape[0]
 
-
     if not isinstance(mean, np.ndarray):
         mean = np.zeros((nnodes)) + mean
     if not isinstance(variance, np.ndarray):
@@ -744,7 +816,7 @@ def generate_2d_grid_realizations(
         delx = x.max() - x.min()
         dely = y.max() - y.min()
 
-        variorange = np.zeros((nnodes)) + max(delx,dely) / 10 # ?
+        variorange = np.zeros((nnodes)) + max(delx, dely) / 10  # ?
     elif not isinstance(variorange, np.ndarray):
         variorange = np.zeros((nnodes)) + variorange
 
@@ -760,8 +832,6 @@ def generate_2d_grid_realizations(
         # TODO warn here
         zone_array = zone_array.astype(int)
 
-    
-
     power = 1.0
 
     lib = PestUtilsLib()
@@ -828,33 +898,29 @@ def __init__(self, delr, delc, xul, yul, rotation=0.0):
         self._ycentergrid = None
 
     @property
-    def xll(self)->float:
-        """lower left x coord
-        """
+    def xll(self) -> float:
+        """lower left x coord"""
         xll = self.xul - (np.sin(self.theta) * self.yedge[0])
         return xll
 
     @property
-    def yll(self)->float:
-        """lower left y coord
-        """
+    def yll(self) -> float:
+        """lower left y coord"""
         yll = self.yul - (np.cos(self.theta) * self.yedge[0])
         return yll
 
     @property
-    def nrow(self)->int:
-        """number of rows
-        """
+    def nrow(self) -> int:
+        """number of rows"""
         return self.delc.shape[0]
 
     @property
-    def ncol(self)->int:
-        """number of cols
-        """
+    def ncol(self) -> int:
+        """number of cols"""
         return self.delr.shape[0]
 
     @classmethod
-    def from_gridspec(cls, gridspec_file)->SpatialReference:
+    def from_gridspec(cls, gridspec_file) -> SpatialReference:
         """instantiate from a pest-style grid specification file
         Parameters
         ----------
@@ -902,71 +968,61 @@ def from_gridspec(cls, gridspec_file)->SpatialReference:
         return cls(np.array(delr), np.array(delc), xul=xul, yul=yul, rotation=rot)
 
     @property
-    def theta(self)->float:
-        """rotation in radians
-        """
+    def theta(self) -> float:
+        """rotation in radians"""
         return -self.rotation * np.pi / 180.0
 
     @property
-    def xedge(self)->np.NDArray[float]:
-        """the xedge array of the grid
-        """
+    def xedge(self) -> np.NDArray[float]:
+        """the xedge array of the grid"""
         return self.get_xedge_array()
 
     @property
-    def yedge(self)->np.NDArray[float]:
-        """the yedge array of the grid
-        """
+    def yedge(self) -> np.NDArray[float]:
+        """the yedge array of the grid"""
         return self.get_yedge_array()
 
     @property
-    def xgrid(self)->np.NDArray[float]:
-        """xgrid array
-        """
+    def xgrid(self) -> np.NDArray[float]:
+        """xgrid array"""
         if self._xgrid is None:
             self._set_xygrid()
         return self._xgrid
 
     @property
-    def ygrid(self)->np.NDArray[float]:
-        """ygrid array
-        """
+    def ygrid(self) -> np.NDArray[float]:
+        """ygrid array"""
         if self._ygrid is None:
             self._set_xygrid()
         return self._ygrid
 
     @property
-    def xcenter(self)->np.NDArray[float]:
-        """grid x center array
-        """
+    def xcenter(self) -> np.NDArray[float]:
+        """grid x center array"""
         return self.get_xcenter_array()
 
     @property
-    def ycenter(self)->np.NDArray[float]:
-        """grid y center array
-        """
+    def ycenter(self) -> np.NDArray[float]:
+        """grid y center array"""
         return self.get_ycenter_array()
 
     @property
-    def ycentergrid(self)->np.NDArray[float]:
-        """grid y center array
-        """
+    def ycentergrid(self) -> np.NDArray[float]:
+        """grid y center array"""
         if self._ycentergrid is None:
             self._set_xycentergrid()
         return self._ycentergrid
 
     @property
-    def xcentergrid(self)->np.NDArray[float]:
-        """grid x center array
-        """
+    def xcentergrid(self) -> np.NDArray[float]:
+        """grid x center array"""
         if self._xcentergrid is None:
             self._set_xycentergrid()
         return self._xcentergrid
 
     @property
-    def areagrid(self)->np.NDArray[float]:
-        """area of grid nodes
-        """
+    def areagrid(self) -> np.NDArray[float]:
+        """area of grid nodes"""
         dr, dc = np.meshgrid(self.delr, self.delc)
         return dr * dc
 
@@ -980,7 +1036,7 @@ def _set_xygrid(self):
         self._xgrid, self._ygrid = np.meshgrid(self.xedge, self.yedge)
         self._xgrid, self._ygrid = self.transform(self._xgrid, self._ygrid)
 
-    def get_xedge_array(self)->np.NDArray[float]:
+    def get_xedge_array(self) -> np.NDArray[float]:
         """
         a numpy one-dimensional float array that has the cell edge x
         coordinates for every column in the grid in model space - not offset
@@ -993,7 +1049,7 @@ def get_xedge_array(self)->np.NDArray[float]:
         xedge = np.concatenate(([0.0], np.add.accumulate(self.delr)))
         return xedge
 
-    def get_yedge_array(self)->np.NDArray[float]:
+    def get_yedge_array(self) -> np.NDArray[float]:
         """
         a numpy one-dimensional float array that has the cell edge y
         coordinates for every row in the grid in model space - not offset or
@@ -1007,7 +1063,7 @@ def get_yedge_array(self)->np.NDArray[float]:
         yedge = np.concatenate(([length_y], length_y - np.add.accumulate(self.delc)))
         return yedge
 
-    def get_xcenter_array(self)->np.NDArray[float]:
+    def get_xcenter_array(self) -> np.NDArray[float]:
         """
         a numpy one-dimensional float array that has the cell center x
         coordinate for every column in the grid in model space - not offset or rotated.
@@ -1019,7 +1075,7 @@ def get_xcenter_array(self)->np.NDArray[float]:
         x = np.add.accumulate(self.delr) - 0.5 * self.delr
         return x
 
-    def get_ycenter_array(self)->np.NDArray[float]:
+    def get_ycenter_array(self) -> np.NDArray[float]:
         """
         a numpy one-dimensional float array that has the cell center x
         coordinate for every row in the grid in model space - not offset of rotated.
@@ -1071,7 +1127,7 @@ def transform(self, x, y, inverse=False):
             y -= self.yll
         return x, y
 
-    def get_extent(self)->tuple[float]:
+    def get_extent(self) -> tuple[float]:
         """
         Get the extent of the rotated and offset grid
 
@@ -1100,7 +1156,7 @@ def get_extent(self)->tuple[float]:
 
         return (xmin, xmax, ymin, ymax)
 
-    def get_vertices(self, i, j)->list[list[float]]:
+    def get_vertices(self, i, j) -> list[list[float]]:
         """Get vertices for a single cell or sequence if i, j locations."""
         pts = []
         xgrid, ygrid = self.xgrid, self.ygrid
@@ -1115,7 +1171,7 @@ def get_vertices(self, i, j)->list[list[float]]:
             vrts = np.array(pts).transpose([2, 0, 1])
             return [v.tolist() for v in vrts]
 
-    def get_ij(self, x, y)->tuple(int):
+    def get_ij(self, x, y) -> tuple(int):
         """Return the row and column of a point or sequence of points
         in real-world coordinates.
 
@@ -1131,7 +1187,6 @@ def get_ij(self, x, y)->tuple(int):
         return r, c
 
     def write_gridspec(self, filename):
-
         """write a PEST-style grid specification file
         Parameters
         ----------
@@ -1157,4 +1212,3 @@ def write_gridspec(self, filename):
             f.write("{0:15.6E} ".format(c))
         f.write("\n")
         return
-
diff --git a/pypestutils/logger.py b/pypestutils/logger.py
index 2c19ac5..d88393d 100644
--- a/pypestutils/logger.py
+++ b/pypestutils/logger.py
@@ -1,4 +1,5 @@
 """Logger module."""
+
 from __future__ import annotations
 
 __all__ = ["get_logger"]
diff --git a/pypestutils/pestutilslib.py b/pypestutils/pestutilslib.py
index e682a09..10403fa 100644
--- a/pypestutils/pestutilslib.py
+++ b/pypestutils/pestutilslib.py
@@ -1,4 +1,5 @@
 """Mid-level pestutilslib module to implement ctypes functions."""
+
 from __future__ import annotations
 
 import logging
diff --git a/tests/common.py b/tests/common.py
index 888529a..3be43e6 100644
--- a/tests/common.py
+++ b/tests/common.py
@@ -1,4 +1,5 @@
 """Common functions for tests. For pytest fixtures, see conftest.py."""
+
 from pathlib import Path
 
 # key is name, value is number of arguments
diff --git a/tests/test_crash.py b/tests/test_crash.py
index a0e5ba6..1907b2d 100644
--- a/tests/test_crash.py
+++ b/tests/test_crash.py
@@ -1,4 +1,5 @@
 """Test previous crash scenarios."""
+
 import numpy as np
 import pytest
 
diff --git a/tests/test_ctypes_declarations.py b/tests/test_ctypes_declarations.py
index 93a494b..f12fcd7 100644
--- a/tests/test_ctypes_declarations.py
+++ b/tests/test_ctypes_declarations.py
@@ -1,4 +1,5 @@
 """Tests for ctypes_declarations module."""
+
 from ctypes import c_int
 
 import pytest
diff --git a/tests/test_data.py b/tests/test_data.py
index 989d4f5..f13fe94 100644
--- a/tests/test_data.py
+++ b/tests/test_data.py
@@ -1,4 +1,5 @@
 """Tests for data module."""
+
 import numpy as np
 import pytest
 
diff --git a/tests/test_drivers.py b/tests/test_drivers.py
index 68123dd..328c032 100644
--- a/tests/test_drivers.py
+++ b/tests/test_drivers.py
@@ -1,4 +1,5 @@
 """John Doherty's driver test programs."""
+
 from pathlib import Path
 
 import numpy as np
diff --git a/tests/test_enum.py b/tests/test_enum.py
index ed66e31..24b4254 100644
--- a/tests/test_enum.py
+++ b/tests/test_enum.py
@@ -1,4 +1,5 @@
 """Tests for enum module."""
+
 import pytest
 
 from pypestutils.enum import KrigType
diff --git a/tests/test_helpers.py b/tests/test_helpers.py
index ce44c7b..9d5b520 100644
--- a/tests/test_helpers.py
+++ b/tests/test_helpers.py
@@ -5,39 +5,48 @@
 import pandas as pd
 
 
-
 def test_mf6_mod2obs():
     import pypestutils.helpers as helpers
+
     start_datetime = pd.to_datetime("1-1-2000")
 
     dts = []
-    for dt in np.linspace(100,300,10):
-        dts.append(start_datetime+pd.to_timedelta(dt,unit="d"))
-    obs_df = pd.DataFrame({"site":["site1" for _ in range(len(dts))],
-        "x":[101 for _ in range(len(dts))],
-        "y":[221 for _ in range(len(dts))],
-        "layer":[1 for _ in range(len(dts))]})
-    obs_df.loc[:,"datetime"] = dts
-
-    ws = os.path.join("tests","data","p09")
+    for dt in np.linspace(100, 300, 10):
+        dts.append(start_datetime + pd.to_timedelta(dt, unit="d"))
+    obs_df = pd.DataFrame(
+        {
+            "site": ["site1" for _ in range(len(dts))],
+            "x": [101 for _ in range(len(dts))],
+            "y": [221 for _ in range(len(dts))],
+            "layer": [1 for _ in range(len(dts))],
+        }
+    )
+    obs_df.loc[:, "datetime"] = dts
+
+    ws = os.path.join("tests", "data", "p09")
     case = "gwf-p09-mf6"
-    obs_csv = os.path.join(ws,"obs.csv")
+    obs_csv = os.path.join(ws, "obs.csv")
     obs_df.to_csv(obs_csv)
 
-    dv_name = os.path.join(ws,case+".hds")
-    grb_name = os.path.join(ws,case+".dis.grb")
+    dv_name = os.path.join(ws, case + ".hds")
+    grb_name = os.path.join(ws, case + ".dis.grb")
 
-    results = helpers.mod2obs_mf6(grb_name,dv_name,obs_csv,31,start_datetime,1,"head")
+    results = helpers.mod2obs_mf6(
+        grb_name, dv_name, obs_csv, 31, start_datetime, 1, "head"
+    )
     iresults = results["interpolated_results"]
     assert iresults.shape[0] == obs_df.shape[0]
 
-    dv_name = os.path.join(ws,"gwt-p09-mf6.ucn")
-    results = helpers.mod2obs_mf6(grb_name,dv_name,obs_csv,31,start_datetime,1,"concentration")
+    dv_name = os.path.join(ws, "gwt-p09-mf6.ucn")
+    results = helpers.mod2obs_mf6(
+        grb_name, dv_name, obs_csv, 31, start_datetime, 1, "concentration"
+    )
     iresults = results["interpolated_results"]
     assert iresults.shape[0] == obs_df.shape[0]
 
 
 if __name__ == "__main__":
     import sys
-    sys.path.insert(0,os.path.join(".."))
+
+    sys.path.insert(0, os.path.join(".."))
     test_mf6_mod2obs()
diff --git a/tests/test_notebooks.py b/tests/test_notebooks.py
index 3b04351..3f37d21 100644
--- a/tests/test_notebooks.py
+++ b/tests/test_notebooks.py
@@ -1,4 +1,5 @@
 """Tests for example notebooks."""
+
 import logging
 import pytest
 import os
@@ -27,17 +28,27 @@
 #     #os.chdir(cwd)
 #     return
 
+
 def test_notebooks():
     pyemu = pytest.importorskip("pyemu")
 
     nb_dir = os.path.join("examples")
     nb_files = [f for f in os.listdir(nb_dir) if f.endswith(".ipynb")]
     for nb_file in nb_files:
-        pyemu.os_utils.run("jupyter nbconvert --execute --ExecutePreprocessor.timeout=180000 --inplace {0}".format(nb_file),cwd=nb_dir)
-        pyemu.os_utils.run("jupyter nbconvert --ClearOutputPreprocessor.enabled=True --ClearMetadataPreprocessor.enabled=True --allow-errors --inplace {0}".format(nb_file),cwd=nb_dir)
+        pyemu.os_utils.run(
+            "jupyter nbconvert --execute --ExecutePreprocessor.timeout=180000 --inplace {0}".format(
+                nb_file
+            ),
+            cwd=nb_dir,
+        )
+        pyemu.os_utils.run(
+            "jupyter nbconvert --ClearOutputPreprocessor.enabled=True --ClearMetadataPreprocessor.enabled=True --allow-errors --inplace {0}".format(
+                nb_file
+            ),
+            cwd=nb_dir,
+        )
     return
 
 
-
-#if __name__ == "__main__":
-    #test_notebook("pypestutils_pstfrom_structured_mf6.ipynb")
+# if __name__ == "__main__":
+# test_notebook("pypestutils_pstfrom_structured_mf6.ipynb")
diff --git a/tests/test_pestutilslib.py b/tests/test_pestutilslib.py
index b0825ee..0f49958 100644
--- a/tests/test_pestutilslib.py
+++ b/tests/test_pestutilslib.py
@@ -1,4 +1,5 @@
 """Tests for pestutilslib module."""
+
 import logging
 from pathlib import PureWindowsPath
 
@@ -34,109 +35,82 @@ def test_create_char_array():
         lib.create_char_array(1, "LENGRIDNAME")
 
 
-def test_inquire_modflow_binary_file_specs():
-    ...
+def test_inquire_modflow_binary_file_specs(): ...
 
 
-def test_retrieve_error_message():
-    ...
+def test_retrieve_error_message(): ...
 
 
-def test_install_structured_grid():
-    ...
+def test_install_structured_grid(): ...
 
 
-def test_get_cell_centres_structured():
-    ...
+def test_get_cell_centres_structured(): ...
 
 
-def test_uninstall_structured_grid():
-    ...
+def test_uninstall_structured_grid(): ...
 
 
-def test_free_all_memory():
-    ...
+def test_free_all_memory(): ...
 
 
-def test_interp_from_structured_grid():
-    ...
+def test_interp_from_structured_grid(): ...
 
 
-def test_interp_to_obstime():
-    ...
+def test_interp_to_obstime(): ...
 
 
-def test_install_mf6_grid_from_file():
-    ...
+def test_install_mf6_grid_from_file(): ...
 
 
-def test_get_cell_centres_mf6():
-    ...
+def test_get_cell_centres_mf6(): ...
 
 
-def test_uninstall_mf6_grid():
-    ...
+def test_uninstall_mf6_grid(): ...
 
 
-def test_calc_mf6_interp_factors():
-    ...
+def test_calc_mf6_interp_factors(): ...
 
 
-def test_interp_from_mf6_depvar_file():
-    ...
+def test_interp_from_mf6_depvar_file(): ...
 
 
-def test_extract_flows_from_cbc_file():
-    ...
+def test_extract_flows_from_cbc_file(): ...
 
 
-def test_calc_kriging_factors_2d():
-    ...
+def test_calc_kriging_factors_2d(): ...
 
 
-def test_calc_kriging_factors_auto_2d():
-    ...
+def test_calc_kriging_factors_auto_2d(): ...
 
 
-def test_calc_kriging_factors_3d():
-    ...
+def test_calc_kriging_factors_3d(): ...
 
 
-def test_krige_using_file():
-    ...
+def test_krige_using_file(): ...
 
 
-def test_build_covar_matrix_2d():
-    ...
+def test_build_covar_matrix_2d(): ...
 
 
-def test_build_covar_matrix_3d():
-    ...
+def test_build_covar_matrix_3d(): ...
 
 
-def test_calc_structural_overlay_factors():
-    ...
+def test_calc_structural_overlay_factors(): ...
 
 
-def test_interpolate_blend_using_file():
-    ...
+def test_interpolate_blend_using_file(): ...
 
 
-def test_ipd_interpolate_2d():
-    ...
+def test_ipd_interpolate_2d(): ...
 
 
-def test_ipd_interpolate_3d():
-    ...
+def test_ipd_interpolate_3d(): ...
 
 
-def test_initialize_randgen():
-    ...
+def test_initialize_randgen(): ...
 
 
-def test_fieldgen2d_sva():
-    ...
+def test_fieldgen2d_sva(): ...
 
 
-def test_fieldgen3d_sva():
-    ...
+def test_fieldgen3d_sva(): ...