Merge pull request #99 from langevin-usgs/patch-gwt

update gwt-demo notebook
DOI-USGS · Jan 31, 2024 · 6f4d84a · 6f4d84a
2 parents 676ac02 + 9f9761f
commit 6f4d84a
Showing 1 changed file with 50 additions and 60 deletions.
diff --git a/notebooks/part1_flopy/09-gwt-voronoi-demo.ipynb b/notebooks/part1_flopy/09-gwt-voronoi-demo.ipynb
@@ -314,6 +314,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
+    "head = gwf.output.head().get_data()\n",
     "spdis = gwf.output.budget().get_data(text=\"DATA-SPDIS\")[0]\n",
     "qx, qy, qz = flopy.utils.postprocessing.get_specific_discharge(\n",
     "    spdis, gwf, head=gwf.output.head().get_data(),\n",
@@ -333,73 +334,61 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "263dd83d-0579-4340-9107-fdb9c7bc600a",
+   "id": "36492a7c",
    "metadata": {},
    "outputs": [],
    "source": [
-    "fig, ax = plt.subplots(\n",
-    "    ncols=1,\n",
-    "    nrows=1,\n",
-    "    figsize=(4, 8),\n",
-    "    constrained_layout=True,\n",
-    "    subplot_kw=dict(aspect=\"equal\"),\n",
-    ")\n",
+    "fig, ax = plt.subplots(figsize=(4, 6), constrained_layout=True)\n",
+    "ax.set_aspect(1)\n",
+    "ax.set_xlabel(r'x')\n",
+    "ax.set_ylabel(r'y')\n",
+    "title = ax.set_title(f\"Time = {times[0]} days\")\n",
     "\n",
-    "radius = 100.\n",
-    "vmin, vmax = 1e-1, 100.\n",
+    "# plot persistent items\n",
+    "vmin, vmax = 1e-3, 100.\n",
     "norm = colors.LogNorm(vmin=vmin, vmax=vmax)\n",
     "\n",
-    "mm = flopy.plot.PlotMapView(model=gwf, ax=ax, layer=0)\n",
-    "cb = mm.plot_array(np.zeros((nlay, ncpl)),\n",
-    "                   norm=norm,\n",
-    "                   masked_values=[0], vmin=vmin, vmax=vmax)\n",
-    "mm.plot_vector(qx, qy, qz, normalize=True)\n",
-    "mm.plot_shapefile(river_shp, edgecolor=\"cyan\", facecolor=\"red\", lw=2)\n",
-    "mm.plot_shapefile(wells_shp, edgecolor=\"red\", facecolor=\"red\", radius=radius)\n",
-    "\n",
-    "plt.colorbar(cb, ax=ax, shrink=0.5)\n",
+    "pmv = flopy.plot.PlotMapView(gwt, ax=ax)\n",
+    "pmv.plot_grid(lw=0.5, color=\"0.5\")\n",
+    "pmv.contour_array(\n",
+    "    head, \n",
+    "    levels=np.linspace(0, 30, 30), \n",
+    "    tri_mask=True,\n",
+    "    linestyles=\"-\",\n",
+    "    colors=\"blue\",\n",
+    "    linewidths=0.5,\n",
+    ")\n",
+    "ca_dict = {\n",
+    "    \"vmin\": vmin,\n",
+    "    \"vmax\": vmax,\n",
+    "    \"norm\": norm,\n",
+    "    \"masked_values\": [0],\n",
+    "}\n",
+    "conc_alldata = gwt.output.concentration().get_alldata()\n",
+    "c = conc_alldata[0]\n",
+    "c[c < vmin] = 0.\n",
+    "cont = pmv.plot_array(c, **ca_dict)\n",
+    "clb = fig.colorbar(\n",
+    "    cont, \n",
+    "    shrink=0.5, \n",
+    ")\n",
     "\n",
-    "time = 0.\n",
-    "ax.set_title(f\"{time=} days\")\n",
+    "def animate(i):\n",
+    "    c = conc_alldata[i].flatten()\n",
+    "    c[c < vmin] = 0.\n",
+    "    cont.set_array(c)\n",
+    "    title = ax.set_title(f\"Time = {times[i]} days\")\n",
+    "    return cont\n",
     "\n",
+    "import matplotlib.animation\n",
+    "ani = matplotlib.animation.FuncAnimation(fig, animate, frames=conc_alldata.shape[0])\n",
+    "plt.close()\n",
     "\n",
-    "for time in times:\n",
-    "    ax.cla()\n",
-    "    \n",
-    "    conc = gwt.output.concentration().get_data(totim=time)\n",
-    "    conc[conc < 1e-3] = 0.\n",
-    "    cb = mm.plot_array(conc,\n",
-    "                       norm=norm,\n",
-    "                       masked_values=[0], vmin=vmin, vmax=vmax)\n",
-    "    mm.plot_vector(qx, qy, qz, normalize=True)\n",
-    "    mm.plot_grid(lw=0.5, color=\"0.5\")\n",
-    "    mm.plot_shapefile(river_shp, edgecolor=\"cyan\", facecolor=\"red\", lw=2)\n",
-    "    mm.plot_shapefile(wells_shp, edgecolor=\"red\", facecolor=\"red\", radius=radius)\n",
-    "    cs = mm.contour_array(\n",
-    "        conc,\n",
-    "        colors=\"red\",\n",
-    "        levels=(0.01,),\n",
-    "        linestyles=\":\",\n",
-    "        linewidths=1.0,\n",
-    "        tri_mask=True,\n",
-    ")\n",
+    "from IPython.display import HTML\n",
+    "HTML(ani.to_jshtml())\n",
     "\n",
-    "    \n",
-    "    ax.set_title(f\"{time=} days\")\n",
-    "    \n",
-    "    display(fig)\n",
-    "    clear_output(wait=True)\n",
-    "    plt.pause(0.1)    \n"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "7f2f2f31-0cc6-4917-88dc-4be55da84ff3",
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "line = [(0, 8100), (2000, 8300), (5000, 8100)]"
+    "# can use this command to write animation to file\n",
+    "#ani.save(\"voronoi-conc-animation.avi\")"
    ]
   },
   {
@@ -409,8 +398,9 @@
    "metadata": {},
    "outputs": [],
    "source": [
+    "line = [(0, 8100), (2000, 8300), (5000, 8100)]\n",
     "xs = flopy.plot.PlotCrossSection(model=gwt, line={\"line\": line}) \n",
-    "cb = xs.plot_array(conc,\n",
+    "cb = xs.plot_array(conc_alldata[-1],\n",
     "                   head=gwf.output.head().get_data(),\n",
     "                   norm=norm,\n",
     "                   masked_values=[0], vmin=vmin, vmax=vmax)\n",
@@ -421,7 +411,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "50d87d5f-488d-4422-85f8-594553fea796",
+   "id": "eb2c863a",
    "metadata": {},
    "outputs": [],
    "source": []
@@ -443,7 +433,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.9.13"
+   "version": "3.10.11"
   }
  },
  "nbformat": 4,