{ "cells": [ { "cell_type": "markdown", "id": "81c052ce-0ca3-4f8c-a6e3-f626433018b3", "metadata": {}, "source": [ "# 05: Unstructured grid and general mesh generation with FloPy\n", "\n", "In this notebook we'll cover how to build Structured and Unstructured Meshes (Grids) with FloPy for use in MODFLOW models. \n", "\n", "This notebook will start with generating a simple Structured Grid and we will use this grid object to begin developing Unstructured Meshes. The grids generated in theis notebook include the following:\n", "\n", " - Regular structured grid (review)\n", " - Irregular structured grid with variable row and column spacings\n", " - Local Grid Refinement (LGR): A locally refined mesh that can be subsetted into a second model within a MODFLOW simulation\n", " - Quadtree Mesh: A mesh with a four fold (quad) refinement for each grid cell in a given refinement area\n", " - Triangular Mesh: A triangular mesh generated using Delaunay triangulation\n", " - Voronoi Mesh: A mesh generated from the connected centers of the Delaunay triangulation circumcircles\n" ] }, { "cell_type": "code", "execution_count": null, "id": "296be634-fa9e-4373-805c-15952a73baa4", "metadata": {}, "outputs": [], "source": [ "import warnings\n", "from pathlib import Path\n", "\n", "import flopy\n", "import geopandas as gpd\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "import pandas as pd\n", "import contextily as ctx\n", "from dataretrieval import nldi, nwis\n", "from flopy.discretization import StructuredGrid, VertexGrid\n", "from flopy.plot import PlotMapView\n", "from flopy.utils.gridgen import Gridgen\n", "from flopy.utils.gridintersect import GridIntersect\n", "from flopy.utils.lgrutil import Lgr\n", "from flopy.utils.rasters import Raster\n", "from flopy.utils.triangle import Triangle\n", "from flopy.utils.voronoi import VoronoiGrid\n", "from shapely.errors import ShapelyDeprecationWarning\n", "\n", "warnings.simplefilter(\"ignore\", DeprecationWarning)\n", "warnings.simplefilter(\"ignore\", ShapelyDeprecationWarning)" ] }, { "cell_type": "markdown", "id": "96fbc012-0c27-4335-9cc6-2314b5e445e5", "metadata": {}, "source": [ "## Creating a regular structured (rectilinear) grid\n", "\n", "This activity will start by creating a simple sturctured (rectilinear) grid from a 1/3 arc second Digital Elevation Model of the Sagehen Creek watershed as an illustration of how we can create an initial mesh that can be used to develop a model from.\n", "\n", "This is a quick review for more information, see the previous activity `04_Modelgrid_and_intersection.ipynb`" ] }, { "cell_type": "markdown", "id": "a140b317-5a20-4af5-9e1e-07392acdeb2a", "metadata": {}, "source": [ "Load the Digital Elevation Model and visualize it" ] }, { "cell_type": "code", "execution_count": null, "id": "32163ef9-ede0-4cdd-a238-c64b8bef43da", "metadata": {}, "outputs": [], "source": [ "data_ws = Path(\"./data\")\n", "geospatial_ws = data_ws / \"unstructured_grid_data\"\n", "dem = geospatial_ws / \"dem.img\"\n", "rstr = Raster.load(dem)\n", "rstr.plot();" ] }, { "cell_type": "markdown", "id": "3e09da28-5aac-43c4-b248-98488e783bb2", "metadata": {}, "source": [ "Get the DEM boundaries" ] }, { "cell_type": "code", "execution_count": null, "id": "11714a91-b05b-4586-b1c2-278faa625bdb", "metadata": {}, "outputs": [], "source": [ "xmin, xmax, ymin, ymax = rstr.bounds\n", "epsg = rstr.crs.to_epsg()\n", "print(f\"{xmin=}, {xmax=}\")\n", "print(f\"{epsg=}\")" ] }, { "cell_type": "markdown", "id": "5b234f85-67ed-4ef0-b6b7-6f4d6cddb9de", "metadata": {}, "source": [ "And begin developing an initial mesh for the model" ] }, { "cell_type": "code", "execution_count": null, "id": "29b63ada-d57b-4254-beb0-e03d7ec9322e", "metadata": {}, "outputs": [], "source": [ "# live code modelgrid creation\n", "cellsize = 200\n" ] }, { "cell_type": "markdown", "id": "998ebaa1-7304-4cbb-b506-7536a91a0a34", "metadata": {}, "source": [ "Finally create a fake top, bottom, and idomain for the model for now. These will be refined in a subsequent step." ] }, { "cell_type": "code", "execution_count": null, "id": "562186f7-e3a8-47ef-9fc8-7202f98be7ae", "metadata": {}, "outputs": [], "source": [ "# live code modelgrid creation\n" ] }, { "cell_type": "markdown", "id": "7d9a5b38-de02-45ed-85ea-57d596c7e0eb", "metadata": {}, "source": [ "And now we can generate a `StructuredGrid` object (called `sgrid`) that can be used for developing model inputs" ] }, { "cell_type": "code", "execution_count": null, "id": "7b52fe9c-991a-4c1e-a088-cd080c13af5c", "metadata": {}, "outputs": [], "source": [ "# live code modelgrid creation\n" ] }, { "cell_type": "markdown", "id": "9232a243-d892-4aac-8af1-191346429ecd", "metadata": {}, "source": [ "We can plot our `SturcturedGrid` object over the existing raster data to check that it is correctly orientied in space" ] }, { "cell_type": "code", "execution_count": null, "id": "88b9c624-fae6-4ebc-8686-11d4bd159da2", "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots(figsize=(6, 6))\n", "rstr.plot(ax=ax)\n", "pmv = PlotMapView(modelgrid=sgrid, ax=ax)\n", "pmv.plot_grid(lw=0.5);" ] }, { "cell_type": "markdown", "id": "b5d6c8c1-bdc9-488d-936d-11cace8c2024", "metadata": {}, "source": [ "### Resampling data to finish creating our initial grid.\n", "\n", "Once we have an initial grid, raster and vector data can be resampled or intersected with the grid to begin building a model domain. Although the focus of this notebook is to generate different types of meshes for modeling, the resampling and intersection process in this notebook can be applied many different data types and be used to generate boundary conditions within a MODFLOW model." ] }, { "cell_type": "markdown", "id": "c3e392d0-6c07-4d97-b3be-f29902094081", "metadata": {}, "source": [ "#### Raster resampling\n", "We'll use FloPy's built in raster resampling method to generate values for land surface elevation (`top`) using the existing `StructuredGrid` object. For more information on raster resampling, please see this [example](https://flopy.readthedocs.io/en/latest/Notebooks/raster_intersection_example.html)" ] }, { "cell_type": "code", "execution_count": null, "id": "369ee997-04f2-4618-94a5-e6fd51458627", "metadata": {}, "outputs": [], "source": [ "top = rstr.resample_to_grid(\n", " sgrid, rstr.bands[0], method=\"min\", extrapolate_edges=True\n", ")" ] }, { "cell_type": "code", "execution_count": null, "id": "7d639180-2dda-48f0-abff-f9b90ae718cf", "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots(figsize=(7, 7))\n", "pmv = PlotMapView(modelgrid=sgrid)\n", "pc = pmv.plot_array(top)\n", "pmv.plot_grid(lw=0.5)\n", "plt.colorbar(pc, shrink=0.7);" ] }, { "cell_type": "markdown", "id": "e714343a-2d99-44e8-b245-5a3a6dea62d4", "metadata": {}, "source": [ "### Vector (shapefile) intersection review\n", "\n", "We can also intersect lines, points, and polygons with an existing modelgrid object. In this example, we'll load up the watershed boundary and get the active and inactive areas of the grid. For more information on intersecting shapefiles with a modelgrid/model see this [example](https://flopy.readthedocs.io/en/latest/Notebooks/grid_intersection_example.html)" ] }, { "cell_type": "markdown", "id": "0aaa3503-a6f0-47c5-9092-454baeedf233", "metadata": {}, "source": [ "Our first step in this process is to get the watershed boundary. We can use the `dataretrieval` package to begin downloading information about the watershed from NWIS.\n", "\n", "**Note**: this information has been stored in the repository as shapefiles if internet access is not available." ] }, { "cell_type": "code", "execution_count": null, "id": "a53d258e", "metadata": {}, "outputs": [], "source": [ "# the get_info() call accepts a decimal lat/lon bounding box\n", "griddf = sgrid.geo_dataframe.set_crs(epsg=26911)\n", "grid_wgs = griddf.to_crs(epsg=4326)\n", "\n", "wxmin, wxmax, wymin, wymax = grid_wgs.total_bounds\n", "wgs_bounds = [wxmin, wymin, wxmax, wymax]\n", "\n", "try:\n", " info, metadata = nwis.get_info(bBox=[f\"{i :.2f}\" for i in wgs_bounds])\n", "except (ValueError, ConnectionError, NameError):\n", " info = gpd.read_file(geospatial_ws / \"gage_info.shp\")\n", "\n", "info.head()" ] }, { "cell_type": "markdown", "id": "d83035a3", "metadata": {}, "source": [ "Reproject the gage data into UTM Zone 11 N" ] }, { "cell_type": "code", "execution_count": null, "id": "3abe3530", "metadata": {}, "outputs": [], "source": [ "info = info.to_crs(epsg=epsg)" ] }, { "cell_type": "markdown", "id": "80e102e3", "metadata": {}, "source": [ "And plot the gage data on the raster" ] }, { "cell_type": "code", "execution_count": null, "id": "0c451f3e", "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots()\n", "\n", "rstr.plot(ax=ax)\n", "info.plot(\n", " column=\"site_no\",\n", " ax=ax,\n", " cmap=\"tab20\",\n", " legend=True,\n", " legend_kwds={\"fontsize\": 6},\n", ");" ] }, { "cell_type": "markdown", "id": "1d1ba528", "metadata": {}, "source": [ "Grab the site record for gage number 10343500" ] }, { "cell_type": "code", "execution_count": null, "id": "36ed5138", "metadata": {}, "outputs": [], "source": [ "sitedf = info.loc[info.site_no == \"10343500\"]\n", "sitedf" ] }, { "cell_type": "markdown", "id": "28eb81e2", "metadata": {}, "source": [ "Get the upstream basin from nldi using the NWIS stream gage station code" ] }, { "cell_type": "code", "execution_count": null, "id": "27c167f8", "metadata": {}, "outputs": [], "source": [ "try:\n", " basindf = nldi.get_basin(\n", " feature_source=\"nwissite\",\n", " feature_id=f\"USGS-{sitedf.site_no.values[0]}\",\n", " )\n", " basindf = basindf.to_crs(epsg=epsg)\n", "except (ValueError, ConnectionError, NameError):\n", " basindf = gpd.read_file(geospatial_ws / \"basin_boundary.shp\")" ] }, { "cell_type": "markdown", "id": "e5a29e7b", "metadata": {}, "source": [ "Plot up the basin boundary with the existing raster" ] }, { "cell_type": "code", "execution_count": null, "id": "099710ee", "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots()\n", "rstr.plot(ax=ax)\n", "basindf.plot(ax=ax, facecolor=\"None\", edgecolor=\"k\", lw=1.5)\n", "info.plot(\n", " column=\"site_no\",\n", " ax=ax,\n", " cmap=\"tab20\",\n", " legend=True,\n", " legend_kwds={\"fontsize\": 6},\n", ");" ] }, { "cell_type": "markdown", "id": "cce0613c-1187-4739-b316-61ae03bf0237", "metadata": {}, "source": [ "Doing the intersection" ] }, { "cell_type": "code", "execution_count": null, "id": "8720cab6-8f03-401d-8150-216a4f6b8a71", "metadata": {}, "outputs": [], "source": [ "ix = GridIntersect(sgrid, method=\"vertex\")\n", "result = ix.intersect(basindf.geometry.values[0], contains_centroid=True)\n", "rowcol = result[\"cellids\"]" ] }, { "cell_type": "code", "execution_count": null, "id": "ce6c3715-0a4a-456c-a833-b2830a6cf764", "metadata": {}, "outputs": [], "source": [ "row, col = list(zip(*rowcol))" ] }, { "cell_type": "markdown", "id": "0a9e6f47-527f-4356-8873-a28c59faf095", "metadata": {}, "source": [ "Building idomain from the intersection" ] }, { "cell_type": "code", "execution_count": null, "id": "12287849-a3ef-4efd-9950-7c544b5e0c79", "metadata": {}, "outputs": [], "source": [ "idomain = np.zeros((nlay, nrow, ncol), dtype=int)\n", "idomain[:, row, col] = 1" ] }, { "cell_type": "markdown", "id": "de193428-71a0-4e32-b6de-e5ba9c04e221", "metadata": {}, "source": [ "#### Putting it all together into a completed `StructuredGrid`" ] }, { "cell_type": "code", "execution_count": null, "id": "2b911cd2-cabe-4ff3-8a48-5ae23c9dbc67", "metadata": {}, "outputs": [], "source": [ "sgrid = StructuredGrid(\n", " delc=sgrid.delc,\n", " delr=sgrid.delr,\n", " top=top,\n", " botm=sgrid.botm,\n", " idomain=idomain,\n", " nlay=sgrid.nlay,\n", " nrow=sgrid.nrow,\n", " ncol=sgrid.ncol,\n", " xoff=sgrid.xoffset,\n", " yoff=sgrid.yoffset,\n", " crs=sgrid.crs,\n", ")\n", "sgrid" ] }, { "cell_type": "markdown", "id": "a329bace", "metadata": {}, "source": [ "Plotting the completed grid with gage locations and the watershed boundary" ] }, { "cell_type": "code", "execution_count": null, "id": "db12ad0b-d15e-4cd3-87c3-6c13cef9fdf9", "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots(figsize=(5, 5))\n", "pmv = PlotMapView(modelgrid=sgrid)\n", "pc = pmv.plot_array(sgrid.top)\n", "ib = pmv.plot_inactive()\n", "pmv.plot_grid(lw=0.5)\n", "plt.colorbar(pc, shrink=0.7)\n", "\n", "basindf.plot(ax=ax, facecolor=\"None\", edgecolor=\"r\", lw=1.5)\n", "sitedf.plot(\n", " column=\"site_no\",\n", " ax=ax,\n", " cmap=\"tab20\",\n", " legend=True,\n", " legend_kwds={\"fontsize\": 6},\n", " zorder=6\n", ");" ] }, { "cell_type": "markdown", "id": "7c512728-fc6e-4e83-8fd8-b8e13b0e844a", "metadata": {}, "source": [ "Because this is a watershed, we can also overlay a shapefile of the streams and get the stream cells using `GridIntersect`" ] }, { "cell_type": "markdown", "id": "1b4155d8", "metadata": {}, "source": [ "NHD-Plus stream lines can be downloaded using streamgage information with the `dataretrieval` package" ] }, { "cell_type": "code", "execution_count": null, "id": "9557e678-89de-475d-b2ac-f138373a0648", "metadata": {}, "outputs": [], "source": [ "try:\n", " nhddf = nldi.get_flowlines(\n", " navigation_mode=\"UT\",\n", " distance=999,\n", " feature_source=\"nwissite\",\n", " feature_id=f\"USGS-{sitedf.site_no.values[0]}\",\n", " )\n", " nhddf = nhddf.to_crs(epsg=epsg)\n", "except (ValueError, ConnectionError, NameError):\n", " nhddf = gpd.read_file(geospatial_ws / \"nhd_streams.shp\")\n", "\n", "nhddf.head()" ] }, { "cell_type": "markdown", "id": "8da12304", "metadata": {}, "source": [ "Intersect stream lines with the existing modelgrid" ] }, { "cell_type": "code", "execution_count": null, "id": "f3285925-ce26-48ec-b9ce-ffeb2c1f6249", "metadata": {}, "outputs": [], "source": [ "rowcols = []\n", "ix = GridIntersect(sgrid, method=\"vertex\")\n", "for geom in nhddf.geometry.values:\n", " rcs = ix.intersects(geom)[\"cellids\"]\n", " rowcols.extend(list(rcs))" ] }, { "cell_type": "code", "execution_count": null, "id": "c3a12d60-e5d7-4dad-bb7e-edead8697191", "metadata": {}, "outputs": [], "source": [ "row, col = list(zip(*rowcols))\n", "strms = np.zeros((nrow, ncol))\n", "strms[row, col] = 1\n", "strms" ] }, { "cell_type": "markdown", "id": "12c3ae47-a03b-4ee2-b879-47a39bf9c1b8", "metadata": {}, "source": [ "Plot the intersection results" ] }, { "cell_type": "code", "execution_count": null, "id": "1d8b1e04-1590-4acd-a89b-ee3283220e43", "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots(figsize=(5, 5))\n", "pmv = PlotMapView(modelgrid=sgrid, ax=ax)\n", "pc = pmv.plot_array(sgrid.top, cmap=\"viridis\")\n", "pmv.plot_array(strms, masked_values=[0,], alpha=0.3, cmap=\"Reds_r\")\n", "pmv.plot_inactive()\n", "pmv.plot_grid(lw=0.5)\n", "plt.colorbar(pc, shrink=0.7)\n", "nhddf.plot(ax=ax);" ] }, { "cell_type": "markdown", "id": "3548e80e-aacc-4df5-8850-f8643431bf6f", "metadata": {}, "source": [ "Information from the grid we've created here can now be used for generating more complex meshes." ] }, { "cell_type": "markdown", "id": "2afa63f8-4086-437f-80e6-d49c976ff7aa", "metadata": {}, "source": [ "## Irregular structured Grid (DIS)\n", "\n", "In this example, an irregular structured grid is created to demonstrate how to add local refinement to an area of the model using a structured (rectilinear) grid. " ] }, { "cell_type": "code", "execution_count": null, "id": "13244acf-3649-4d98-9e9c-5cc9621509cb", "metadata": {}, "outputs": [], "source": [ "lx = xmax - xmin\n", "ly = ymax - ymin\n", "cellsize = 200\n", "rcellsize = 100\n", "\n", "smooth = np.linspace(100, 200, 10).tolist()\n", "delr = 9 * [cellsize] + smooth[::-1] + 8 * [rcellsize] + smooth + 9 * [cellsize] \n", "delc = 8 * [cellsize] + smooth[::-1] + 7 * [rcellsize] + smooth + 8 * [cellsize]\n", "nlay = 1\n", "nrow = len(delc)\n", "ncol = len(delr)\n", "\n", "# create some fake data for the moment\n", "top = np.ones((nrow, ncol))\n", "botm = np.zeros((nlay, nrow, ncol))\n", "idomain = np.ones(botm.shape, dtype=int)" ] }, { "cell_type": "code", "execution_count": null, "id": "e10e99f1-cd5b-40f4-aab9-9b7c9a8e95cf", "metadata": {}, "outputs": [], "source": [ "rfgrid = StructuredGrid(\n", " delc=np.array(delc),\n", " delr=np.array(delr),\n", " nlay=nlay,\n", " top=top,\n", " botm=botm,\n", " idomain=idomain,\n", " xoff=xmin,\n", " yoff=ymin\n", ")\n", "print(nlay, nrow, ncol)\n", "print(rfgrid.nrow, rfgrid.ncol)" ] }, { "cell_type": "markdown", "id": "737b140f-a786-473c-9632-f76f3b2ae440", "metadata": {}, "source": [ "Intersect the basin boundary with the grid to set active and inactive cells (idomain)" ] }, { "cell_type": "code", "execution_count": null, "id": "4820e68f-dd08-4530-88fb-0f841ae58ab9", "metadata": {}, "outputs": [], "source": [ "ix = GridIntersect(rfgrid, method=\"vertex\")\n", "result = ix.intersect(basindf.geometry.values[0], contains_centroid=True)\n", "rowcol = result[\"cellids\"]\n", "\n", "row, col = list(zip(*rowcol))\n", "idomain = np.zeros((nlay, nrow, ncol), dtype=int)\n", "idomain[:, row, col] = 1\n", "\n", "# tricky way so we don't have to recreate the grid\n", "rfgrid._idomain = idomain" ] }, { "cell_type": "markdown", "id": "f2d24d71-6b47-4caf-b9a0-dbc651cc7c94", "metadata": {}, "source": [ "Resample the top elevation" ] }, { "cell_type": "code", "execution_count": null, "id": "89b48168-809f-422f-9943-6bdc796927dd", "metadata": {}, "outputs": [], "source": [ "top = rstr.resample_to_grid(\n", " rfgrid, rstr.bands[0], method=\"min\", extrapolate_edges=True\n", ")\n", "\n", "# tricky way so we don't have to recreate the grid\n", "rfgrid._top = top" ] }, { "cell_type": "markdown", "id": "0814ada3-16ac-481b-b129-49d241b3ea26", "metadata": {}, "source": [ "Now we can plot the grid" ] }, { "cell_type": "code", "execution_count": null, "id": "4d450173-15f0-43b7-abbd-c9946cd885bb", "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots(figsize=(5, 5))\n", "pmv = PlotMapView(modelgrid=rfgrid)\n", "pc = pmv.plot_array(rfgrid.top, cmap=\"viridis\")\n", "ib = pmv.plot_inactive()\n", "pmv.plot_grid(lw=0.5)\n", "plt.colorbar(pc, shrink=0.7)\n", "\n", "basindf.plot(ax=ax, facecolor=\"None\", edgecolor=\"r\", lw=1.5)\n", "sitedf.plot(\n", " column=\"site_no\",\n", " ax=ax,\n", " cmap=\"tab20\",\n", " legend=True,\n", " legend_kwds={\"fontsize\": 6},\n", " zorder=6\n", ");" ] }, { "cell_type": "markdown", "id": "e31adf58-9a3e-4b03-ab42-1d1a247bc8ce", "metadata": {}, "source": [ "## Local Grid Refinement (LGR) mesh\n", "\n", "A Local Grid Refinement (LGR) mesh can be created in an area of interest to get finer discretization than the existing/surrounding mesh. FloPy's `Lgr` utility, creates the refined mesh as a seperate modelgrid, and calculates the exchanges between the multi-model simulation mesh that results from refinement. \n", "\n", "The `Lgr` utility has a number of input parameters:\n", " - `nlayp`: the number of parent model layers\n", " - `nrowp`: the number of parent model rows\n", " - `ncolp`: the number of parent model columns\n", " - `delrp`: the parent model delr array\n", " - `delcp`: the parent model delc array\n", " - `topp`: the parent model top array\n", " - `botmp`: the parent model botm array\n", " - `idomainp`: an idomain array that is used to create the child grid. Values of 1 indicate parent model cells, values of 0 indicate child model cells. The child model must be a rectangular region that is continuous.\n", " - `ncpp`: number of child cells along the face of a parent cell\n", " - `ncppl`: number of child cells per parent layer\n", " - `xllp`: (optional) x location of parent grid lower left corner\n", " - `yllp`: (optional) y-location of parent grid lower left corner" ] }, { "cell_type": "markdown", "id": "6729b854-0fd2-4595-8fbc-5411f80adf86", "metadata": {}, "source": [ "**Load up the refinement area polygon**" ] }, { "cell_type": "code", "execution_count": null, "id": "e2764262-3604-4314-96fd-0f6e7b16b209", "metadata": {}, "outputs": [], "source": [ "lgr_gdf = gpd.read_file(geospatial_ws / \"lgr_bound.shp\")" ] }, { "cell_type": "markdown", "id": "5c81271d-5554-448f-bc4d-464ce5617b8a", "metadata": {}, "source": [ "Plot it to see if it lines up properly with our model" ] }, { "cell_type": "code", "execution_count": null, "id": "5b143c4c-a4b5-4fe6-825b-6e540568ac6f", "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots(figsize=(5, 5))\n", "pmv = PlotMapView(modelgrid=sgrid, ax=ax)\n", "pc = pmv.plot_array(sgrid.top, cmap=\"viridis\")\n", "pmv.plot_array(strms, masked_values=[0,], alpha=0.3, cmap=\"Reds_r\")\n", "pmv.plot_inactive()\n", "pmv.plot_grid(lw=0.5)\n", "nhddf.plot(ax=ax)\n", "lgr_gdf.plot(ax=ax, alpha=0.5)\n", "plt.colorbar(pc, shrink=0.7);" ] }, { "cell_type": "markdown", "id": "058053a4-d54c-4dea-bdae-52cbe110c6af", "metadata": {}, "source": [ "Creating the parent/child idomainp array for `Lgr` using `GridIntersect`" ] }, { "cell_type": "code", "execution_count": null, "id": "55b19ebb-3baa-4718-83a9-a965b3f79ead", "metadata": {}, "outputs": [], "source": [ "lgr_gdf" ] }, { "cell_type": "markdown", "id": "fe0e95b8", "metadata": {}, "source": [ "## Activity: Create the `idomainp` (parent, child) array from the polygon provided by the lgr geodataframe\n", "\n", "`idomainp` should be the same shape as the existing modelgrid's `idomain` array. The child model extent should be represented with the value 0 and the parent 1.\n", "\n", "*Hint*: `GridIntersect()` can be used to find the cellids where the lgr polygon intersects with the existing modelgrid" ] }, { "cell_type": "code", "execution_count": null, "id": "ccd7f65c-911e-409c-8382-4fe1655b5caf", "metadata": {}, "outputs": [], "source": [ "# live code this/give it as an activity\n", "idomainp = np.ones(sgrid.idomain.shape, dtype=int)\n", "\n", "gix = GridIntersect(sgrid, method=\"vertex\")\n", "result = gix.intersect(lgr_gdf.geometry.values[0])\n", "cid = result[\"cellids\"]\n", "\n", "crow, ccol = list(zip(*cid))\n", "idomainp[:, crow, ccol] = 0" ] }, { "cell_type": "markdown", "id": "21d37cdb-96f0-4deb-95cf-755416422091", "metadata": {}, "source": [ "Creating the `Lgr` object/child grid" ] }, { "cell_type": "code", "execution_count": null, "id": "7552cda5-806e-4e86-aaa5-0e9351fc875b", "metadata": {}, "outputs": [], "source": [ "lgr = Lgr(\n", " nlayp=sgrid.nlay,\n", " nrowp=sgrid.nrow,\n", " ncolp=sgrid.ncol,\n", " delrp=sgrid.delr,\n", " delcp=sgrid.delc,\n", " topp=sgrid.top,\n", " botmp=sgrid.botm,\n", " idomainp=idomainp,\n", " ncpp=4,\n", " xllp=sgrid.xoffset,\n", " yllp=sgrid.yoffset,\n", ")" ] }, { "cell_type": "markdown", "id": "36bab4ba-d52d-4e15-8f68-4f389775fd34", "metadata": {}, "source": [ "Get a `StructuredGrid` of the child model" ] }, { "cell_type": "code", "execution_count": null, "id": "6860ee97-09eb-41d7-a614-d413217cc031", "metadata": {}, "outputs": [], "source": [ "childgrid = lgr.child.modelgrid\n", "childgrid" ] }, { "cell_type": "markdown", "id": "8742de05-8964-4bc6-aa9a-4c8df0cc1825", "metadata": {}, "source": [ "Resample top elevations to the child grid" ] }, { "cell_type": "code", "execution_count": null, "id": "fb559103-4392-44a2-8aba-bf68ae03e8c9", "metadata": {}, "outputs": [], "source": [ "childtop = rstr.resample_to_grid(childgrid, band=rstr.bands[0], method=\"min\")\n", "print(childtop.shape)\n", "print(childgrid.shape)" ] }, { "cell_type": "markdown", "id": "b820dced-8c11-428a-adec-95a76887d8f8", "metadata": {}, "source": [ "Reset the top elevations in the `childgrid` object" ] }, { "cell_type": "code", "execution_count": null, "id": "a8621ac6-067d-47b4-bffd-9f5ff71a42d8", "metadata": {}, "outputs": [], "source": [ "childgrid._top = childtop" ] }, { "cell_type": "markdown", "id": "8b21c6fa-cf67-4ba0-bc9d-225a403a5d5a", "metadata": {}, "source": [ "Intersect streams with the Child mesh to define stream cells" ] }, { "cell_type": "code", "execution_count": null, "id": "f2bef21c-e3cd-4377-b4e4-ac3545ea626e", "metadata": {}, "outputs": [], "source": [ "ix = GridIntersect(childgrid, method=\"vertex\")\n", "cstrmcells = []\n", "for geom in nhddf.geometry.values:\n", " rcs = ix.intersects(geom)[\"cellids\"]\n", " if len(rowcol) > 0:\n", " cstrmcells.extend(list(rcs))\n", "\n", "cstrms = np.zeros((childgrid.nrow, childgrid.ncol), dtype=int)\n", "row, col = list(zip(*cstrmcells))\n", "cstrms[row, col] = 1" ] }, { "cell_type": "markdown", "id": "a6deba6b-ff4b-4fe0-ad63-e90bed6623a4", "metadata": {}, "source": [ "And now we can visualize the new LGR mesh" ] }, { "cell_type": "code", "execution_count": null, "id": "99f30580-98eb-44c4-bb24-837373be5bdb", "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots(figsize=(10, 10))\n", "vmin, vmax = np.min(sgrid.top), np.max(sgrid.top)\n", "\n", "# plot the parent\n", "pmvp = PlotMapView(modelgrid=sgrid, ax=ax)\n", "ptop = sgrid.top.copy()\n", "ptop[crow, ccol] = np.nan\n", "strms[crow, ccol] = np.nan\n", "pc = pmvp.plot_array(ptop, cmap=\"viridis\", vmin=vmin, vmax=vmax)\n", "pmvp.plot_array(strms, masked_values=[0,], alpha=0.3, cmap=\"Reds_r\")\n", "pmvp.plot_inactive()\n", "pmvp.plot_grid(lw=0.3)\n", "nhddf.plot(ax=ax)\n", "\n", "# plot the child\n", "pmvc = PlotMapView(modelgrid=childgrid, ax=ax, extent=sgrid.extent)\n", "pmvc.plot_array(childtop, cmap=\"viridis\", vmin=vmin, vmax=vmax)\n", "pmvc.plot_array(cstrms, masked_values=[0,], alpha=0.3, cmap=\"Reds_r\")\n", "pmvc.plot_grid(lw=0.3)\n", "plt.colorbar(pc, shrink=0.7);" ] }, { "attachments": { "Quadmesh.png": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXMAAAFjCAIAAABADnxeAAAhQ0lEQVR42u3dd1wTZwPA8SSEKQoIDlBQVKxb68BVrXXVXSuOWvesq9W6qtY966zWPftat9bV2jrr3oJbQEVRwQEie2Xc+QYUDFMSEDD8fp/+Yckd4TK+uXtyeSJ5TUSU1Um4CYgIWYgIWYgIWYiIkIWIkIWIkIWICFmIyCBk8fb2Pk15rFOnTrHh2bx6zt7mt2/fFkUx+2QRBMHe3t7IyMiY8kwymUwqlebNDZdIJHK5XL/VpfHpfe2adTV/QI5suGaTXVxcIiMjs08WtVpta2s7e/bs+5Rn6tGjR9u2bfPghk+ZMsXZ2fn69ev6rV6nTh03Nzf91r13756rq+uwYcNyZMMPHTpUsmTJ8PDw7JZl/fr1HE/mnTSP72+//TYPbviyZcs++eQT/Z5gmpo0adK3b1/91tUciXz++eca2nJkwz09PZGFkAVZkIWQBVmQBVkIWZAFWQhZskkWlSLa6/b1O3cfCmm/Nfs+WcSn97wCQ6KQBVmQBVnin/Zqxeyxw2pW+kRuYb3j7B39ZPHxOFnaJt/Gw7eRBVmQBVninwWKKN+HfmpFRN1SlutPXtVDFkX4y06tP5dIJH8cvoMsyIIsyPK22LCgucO/MbOwPOr+QA9ZNsybsHHbxiyU5eTeDWPGjBmbbqNHjz58zRNZCFly8ThLdOS5v7d/6lykxeCJusry1ONE007dF/8yUSPLqMXbMy+LqI7uUb+85H0ZmVufvfsCWQhZcqksishgj5uegij+uWRG30kLdZXlyfUzo4b/0MWtjebZ7jZiVuZlCfG94VjA+L2yVP1quDJjHwVCFkKWHJDl8cU9+SwtB09ZOHhAX0//QD2OhjTduXRQ82zfeNQz87K4H9ksyUC/Hr7POAshS+6VRVApbl/zuHjx8tOXr9JZPX1ZYmMivb29w6IVmZdl7cxRb+ywsLKpVqPWZw20+6yQhZHmIjOHSkEKNbIQsuTqcZaMlF1nygk/dvxcbmI3ZeHaoNBwlVqtOUxLTBXmV8rSTCNL4+5jMj4nArIQsuR1WYTY8Ca1q28/cSPVIRT3A7/LpXG7M3M3H8r4tSMLIUtel0WtiLpxwzON/RHx535tNawYW9l7vwhDFkIWZMmCM+WE2NDqJYtoZCnXsItKl/nhkIWQBVnS7Mn14/Gjt5Ihs1bptCKyELIgS5qtmfq9hhWpicWx20+RhZAlF8liZ2c3Y/r0eXPn6vFfqVKlatasqd+6mpydnSdNmqT3Hy8K6q/rxJ2YW6hahxjhNbIQsuSWli5damZmVqtWrdp6VahQIQ1MtfUtf/7848aN0/uPj/K7amscN0P4iNX/6roushCyZMfRkKhXvXv37tmzp6hvDRs2zMTRkLjse7e4E+SKlHsSHosshCyGM87SJ16WHBlnUUc+q2hXQCNLg07DBN2/NQhZCFmQJZVuHN1pEnckJJmybq8eqyMLIQuypLL2tIGd476DzbLwjSfByELIgixZIIuoiKhbtphGlpK12yj0+QJVZCFkySFZokODjh47dfDgmbNXbqc1kJG+LIJKcdP9/OHjZ8SsluW510VLo7gPC/X6ebF+G44shCw5I8vRjQvj5y2QNu+b5meI05FFGR3SvXn9jsMnvwiJyPJ9lp1Lpsf/abLtJ+8gCyFLild1tTo2VqFWC7lOFlEY+123MZOnb9y688GjF7ofDYkLxvdp3P2H8Jg053jLjCzftf1MA4uJQ7VnMWpkIWRJfNqqTx3a823Hr5xLOBUpUsS+mFPTtp13HzojiGIukUUREVizZGGZTGZhVXDxzkO6yqIMuudc0KpZW7cyzo5T12zJWlmUL70Km8o1svRcsEPvDUcWMjRZFFGhY3u1T2VOV5m84/il0SohN8gS/8wXXjy89U3zmhXqd47VcZwl8uFFmwKWB919hzWtbOvaMmvHWZ5fOTdjxozpM2b6vIpAFkKWt0+oWcN7pD33vMnPK7Zl545Lmt9kJqgfPX6iFsQg3+st2vVQ6iiLMuCOg5Xl5kPXlo1yc6jXRsjqo6HMhyxkULJ4nz5gYiQ1MjYuXbFa4y8aFS9incwWExtnv4jYHJclNiygVgnHBq26zJk3/4iHlx7jLKtnjWzfo2/bxq7Lth1+jSzIgiwfdIdlQq929hVcD168IcS/kRsZ5Nfis2rJcNl83Cs3HA0pYqJfvQpWpTu6nO67zmJEeFhUrCrN2wJZCFmyJFWoX/3Klc54P9H+4Uu/u5VLJtlzWb3jbG6QJSPl4Nn9yELI8jZBpXj85EmKIQvx/K75Rlqy/HH0NrIgCyFLZnt67VjBd7CYXg2IyLarRhZkIYOVxffSAcsEVxyqNY3NxhPnkAVZyEBlEcW9c75POKXFaO7GA9l55ciCLGSYsgjK6PafOr2BpfwX3V/FqJAFWQhZMtu9KwdN4lkp6Fjhsm9gNl+7RhZra+t+/foN0iuNSmXLlh2kb/b29jVq1BiUE3Xp0qV48eIhISHIQgYpiziqR/P4j+3KF+w5mf1Xr5HFzMysZs2a+k2RrVmxRvXqmZlhWyNLZ71q1qyZVCpt2LChfqu3aNHC0dExNDQUWcgAZXl07ZStWdwuy9e9x8XmxCeeM3k0lKnDQEGoX7/+2rVr9Vv92rVrxsbGR48e5WiIkCVJseEBrepWkhrJO4+ZERqtyJENRxZkIYOSRVQrJvRuptlbaTro52Tnzr187v8iOJue6siCLGRQsuxdETcrWoXGX/uHRCZ9/Cm6tGy84/IjZEEWQhZd9lZE9cktv1nnM7VwqLhg+ZrNiW3atH7tmn6d2lkXq+gXI+RxWUTNzaRWpzOhRDqyiFqraa5FTO23IAsZmizel47Y5ov/nJCJSaqTtJRo3DHbZmjJnbIIasWq+bOGDeq9cN0utaiDLBqObp4/Nn/59jf/+8DjxODvBs5fuUWVYopwZCGDkkUR8aJxpTKSdGs/YPLrvC3LtUPbPnFt+ezRDSdbq4NXfTIsi3jqz80NqpRp1TfuW+jD/L3L2Fuv+vOIazmX30/cRRYyWFlU0aH9O30peV9T1u7O1bKIYkxk2D3PW7t2btM8TTZt2nLp6vXA4FBdJ/FNR5ZVc8e5uLaKUsQ2KFd4wvxtOh0NTRzw9RtZTu1dLZFI9533aVm+VIXuI5GFDFUW8cCi4UaS97ft1P1cK4soin9vXV2jYkmpVGJsnr9SlSp2+fNp/uYCRR3Hz10Zq8skvunIcnbPGmk+h92Hzn/qZDVp8U4dZenwRpb9y8cnylLs6yHIQgYrS2TIS98MFKVU51ZZ1L/PGW1qEjdvfuES1U+7e0ZFRz+87V69gnP8l//IO4+envHz/dKRRRUbtWX9il8mjjI2Mz/gfl8/WU7sXB4ny4U4Wcp3H4EsZMjjLLmtjMsiCqrts0eaxH8zoabF+95NfHflwCbTtz+W9J+/Scy0LJpdo6c+NzrUrTLop3m6jeAK6h+/adq85xjNoVnIo5slChdYufVQvXKlVx65jSyELLlRltBn90pZm76brPfKw8SLFIF3SxayePvtA9aVg5RiJmW5637yh9Hj/znrrhLS/FWpynJ017ruXTp2+rb7rhOX4vg490/XXj2nLVqj5r0hQpbcKYv3mb3ag0G/bD72bk8hOrCMY6GEL02Se4YpMymLIKjfOx6cwTPl1CqVwPkshCy5Vpbr/6zXluWTzzspE3YExKgXpR3tEmQp6BulzqQsGYlzcAlZDEGWB1cOab+xJZWZzFq7740tQT7uRSzffunjJ80GKDI20IIshCzI8loZHtiykmPSb4u1+Kr30Psvgrf+OunNTxyr1b39NDiDV40shCzIEpfX8c22ZvJkZ984lK5cumjcVw/YOZS+ct8/41eNLIQsyPJmQEW9ecFU09TO7ivoVOrCncc6XTWyELIgS4ItonBw02L7ghYpcXFt0zn5jBAZkOW3334L1aszZ87I5fJ9+/bpt/rly5eRhZAlt8iiscXzyomGlQqn+rmEJgPGKtWiTrIUKFDAvmhRPf7TPEnj9pVsbPRb3c7OzsHBITg4GFkIWXL+aGjHqlm2+eN2WOp82eWnQT1MZUlkkRqb/LT+gE6yfPPNN4v1asyYMTKZbMiQIfqtPn78eGbYJmTJeVlUsREzv2tr/OYsfiOr03efaY6MPI5sLGNvmfQNo0I3AyIzLgvjLIQseVoWz/92GyfwYVPry5DYt6fDndy8LOnbRdL1+y8gC7IQsmRIluXThifi0bj7NFXCKfOCMqbrl7W0aVm4429kQRZClgzJMq1vi0Q7vhqxWnuuhEt//aF9eu6vO/9FFmQhZMmQLPNGdk+0o82IFdqyBNw+UyhxSl+p0V53P2RBluSJoqhUKnWdeRBZDF6WYztWJszBInFq0DFc8e5jh088jlsl7LSUqNsxNBOzKIiCOuip75JpE896PdQ8FJ/73ps+eoR7akPCqcqiiI489tf2havipqETVTF71vza0c1t+qJ1SmZRyMkEtfeNy8P7dC1bqvSNl1HIgixJnhcx4WO6NZe91cVk9d/nEiwQ5o/q+1accp/dfBGW0YdbarIEB/ht/32hhUSy7YyHKipoxfzJconkgO+rjMkiXDiwp0XtSm/mlLt58q/xM5ZMHt5DJjE++CAQWXKmR14eA7p+ld/C7M1DxCMgAlmQJVnKqFf9OjR9+/ZQkbJ7jl8MDXn1x9IZ5vHzV1o7uVzw9Mv4Vad5NKQOKCmVamSJO87yvmKtgyxxJc5WqYrb9X798v5VayPLc/7ByJKtKWOj3c8dH9G/h7Vlku++QRZkSTVFdMT6pfO+av1liWL2tnaF7B3sCxWyb9S81aQ5C739A3U6hP7Qsrw9iNu1qnXfcRwNZeuIiq+nx9fNP0v1NG1kQZb3HTqrIyLCw8LDFSqVngffH14WMTa8c5smDwJDX/j7IUs2sfLfhtmudeoOGPJD/56dbfLJkQVZsrk0ZYnxKyyVbDtzVfNP/zsXLOJkCc6wLOK4Xq1a9ZsYPzAUNrnbl2Vrf/Zd3+6jl+xBlmxKrVKq1PHD+6K4Y8EUZEGW3CBLZHDgnxtW9O/ff+5vG0Mjw7atWar595xlG18ERWZElhvnj44fOXzE2Annb92/duJg/4RO+AYhSw50/8hmZEGW3LLPkrE4n+Uj6NGpbciCLMiCLMiCLMiCLMiCLMiCLMiCLMiCLMiCLMiCLMiCLMiCLMiCLHlHFlNT08J5L0tLS5lMZmNjUzAnMjEx0fwB+q1rZWWl+csLFCig9+rMsI0s2SFLu3btHuW9pk+frnnpPn/+/PWcyMPDQ7ProffqV65c0XvdvXv3ss+CLBwNGeDRUM7G0RCyIAuyIAuyIAuyIAuyIAuyIAuyIAuyIAuyfMyynHocgizIgizIkql8jmxJJst+r0BkQRZkQZZM9e//FieTZevZe8iCLMiCLJlIUPRp3TCZLJ1HzUYWZEEWZNG5yGc+SxfOGzt8WMOaVY2l0mSyyEzN6jZrPWHy1C37DyELsiALsmQ0ZeSrC2dPv7dbdx8gC7IgC7IQsiALshCyIAuyELIgC7IgCyELsiALIQuyIAshC7IgC7IQsiALshCyIAuyELIgC7IgCyELsiALIQuyIAshC7IgC7IQsiALshCyIAuyELIgC7IgC7J8wNQqxQOvm3d9/d6/qCgqYqIe+tw7e/bMmbPnPe/eD4uKEcWPVRZRFP197125dVePTRBF4Yb7pdOnT7vfvosshCxJUkaF/71ry5f1a5nIJAOmL0h/4eiQwBkTfixbuoRM9m6iv6JOZYaOmxkarfi4ZFHGRJ079k//rm5WFiYVO/ZR6U7LpQObLIyNNLdAtc6DkYWQJeElV1Af2rmmikuxRCYGpitLbFhA+8aukjRyGzkly/dcPpAsmv2UR7cvNKlX1ShhwyvpLosqKsi1tMOb1T9FFkIWTcEBTzevXdqkbo1kMwynI0uQ//0WrmUlaSeTWyzce1bM3bJER4T8vfOP9q0aWxgZaf/xusqijAge2K5u4urIQnldFlGtOrhpZdkSRVPVIR1ZJvT4SrOAVZHiPfoPnjzxp4Y1qqRc3b5qo8BYde6URVQr/9qwpEYF51Q3XFdZti0cK9NaHVkor8uijgmbOmH01t37r111nzWyv1HGZHnqdd7GOn+fn2b6+AeI8c9ARXTkkskjzWVJd3qsnG8+Cc2dsmg2fOK40Rs277h85fKG36YXtTDVW5aAe1cdbQtIkIWQJfWXcWVY5UIFMiCLuGBs32kbdqUYplEtn9RH+6VbYu180y+XypKs8d3b6ieLKupVm5qVJElFRRZCliSPphblC2ZknyUiNEwpCCl/HvngvF1+k8TVi5Rv9CJa9VHIsmxMN31kEZQzBnwtkcgbNaqLLIQsmZUlrRR+HoWszN6sKzU2n7Tlv1w+gpspWUTxwMZF5nJp7a96bpo7ElkIWT6MLKLw95IJxtK3rvwwc32Wny6Xq2R5fOdicStzu0r1HrwI371oNLIQsnwQWUL871Yoahm3mpGk/dgZCrXw2pBlEfq1qGVm7XDS+7nmf5CFkOUDyCKKQf6+nRrXS1zRoVTF/ScvGeo+i6COXTtpqFwun7B855ulkIWQJWtlEfauW9a1Q2u7ApbJTgkxymf9w2+bolSCockiiruXzTQxkjR0GxyTsF+GLIQsWSzLuB4tTdM6D9dIvnL3KQOTJerlo5K25hIbpxt+rxJ/iCyELFl9NCSK4S+frl06r7SDXUpbKjT/ViWIBiOLqFaM7t7GxKbwpiPXtC9HFkKWDzDOEt8rP58fv20ilyWlxazok/BYw5BFUMUsH9tLLpfPWLcn2WXJZKnWBVkIWbJIljfPvcHtaifbbbnyPMwwZPE8s99IJpWb2AwdMTJZ7RrV1F698CdVR44aNfbnSUGxArIQsmRWFk3P73oUtzAxSFmuH9kq0SUza9tHkSpkIWTJAllEUdyzaKzW77B5EBqDLMhCyJIpWTSFPPRwsnj7Eb0C5evEZN0bzzktyxZdZfFFFkKWlLIMmDZfj98iRARUKlXkzW9oP2hiFs4sl32ydOid8i2t0IDHf6bR0B7ttFcvUb+F5od79v8VoWKchZAlhSxtRs5MTQ713k3Lvhvy/ZaDF1SpqaEK9a/gZBO3vkn+A5ezcq7pbJPFvEEnhVoHEXnXmZAlnTGSmHoO1trPkJYjpqdYRti6aKRZ/AxPUiOTToPHhkYlm0lb9Dq+rYBcIpEZd520PFYtfhSyzB7QUXvDjeu56fS5J2QhZEmz0GfeBfMlOZ+2crNu6qQHM6Ja2atNkndYqzdocfjibWXC89Dn6n/ViltLjEyGjf9NkdWfSfxQsqgV7RtVSzobnsuDV1HIgizIktliIoJ/6tYhxen5Jsv+PplkwEEU9y0ak2wpuYlF03ZdV2/YMG/SaPuCBTSr9Zu4QClk/XcOfQhZFFFh68Z/b5psnk2JpP3wSVFKNbIgC7Lol7ByxsiWzZuVcXJyLJ5Kjo6O1V3rtm3b4Xbo27c5BGX0yjmTShQtbJT02Whqnq+Yc5luA74/dvmmWvggX2WWpbIISyYNa9u6VaWyLo5pbHnlqtXbtm2356T7e3/Xv+vnaK/aeugEZCH2WfQpOiLM89aN06dO/Kvp4MFTZ876+D5RqNQf9Er59lVkIQOXJUdCFmQhZEEWZCFkQRZkQRZkQRZkQRZCFmRBFkIWZEEWZEEWZEEWZCFkQRZkIWRBFmQhZEEWZEEWQhZkQRZCFmRBFkIWZEEWZCFkQRZkIWRBFmQhZEEWZEEWQhZkQRZCFmRBFkIWZEEWZCFkQRZkIWRBFmQhZEEWZEEWQhZkQRZCFmRBFkIWZEEWZEEWZEEWZCFkQRZkIWRBlrwui0ql0shibGxsRnkmo/jM816ax7lUKs2DG6650+3t7YODg7N7n6VXr17rKM/UuHHjqlWrzsl7zZ49e/r06Xlww0eOHOnk5BQWFsbREHE0RIyzELIQsiALsiALsiALIQshCyELIQuyIAshC7IQshCyELIQshCyELIgCyELIQshCyELIQshC7IQshCyELIQshCyELIgCyELIQshCyELIQshC7IQshCyELIQshCyELIgCyELsiALIQshCyELIQuyIAuyIAuyELIQshCyELIgC7IQsiALIQshCyELIQshCyELshCyELIQshCyELIQsiALIQshCyELIQshCyELshCyELIQshCyELIQsiALIQuyIAshCyELIQshC7IQsiALshCyELIQshCyIAuyELIgCyELIQshCyELsiALIQuyELIQshCyELIQshCyIAshCyELIQshCyELIQuyELIQshCyELIQshCyIAshCyELIQshCyELIQuyELIgC7IQshCyELIQsiALIQuyIAshCyELIQshS56RRbx54cSaTbti1GL6ywmCoHpfIrJ8DL14dHf10qWPgyP1e8DERITedL+8f8+uVSuW7zp8AlmQJemmKWJOHdrXt3P7/OamjnW+CFYI6S+/YcHPJdOvfI0H4Qpkyb2vIaLw8M7VqeNGFLe1kkgkJ3ye6rS6IiJ4z6Z1Hdu0KFncwdzMrFSFT4eO+um/K7eQBVnePrx8bl+dN/HHcqUcJAk5vU8WISa4QTknSbq5TVwtiiKy5MK7PCo85PDurc3qV5XJpIn3V8ZlCXh8b2jfb/NbmGnWsrAu0uW7EaeveaoFQTSUWwhZMlt40NNfJo3Kb2GaDIX3yuL93/b8cmk6rMit7N2fBDHOktsS1Mq/tqytWq5UyrssQ7KIgtf5f8qXKPJmFdviZY6732OcBVm0HyGqX0b0KePsWKlqVUtzE51kEdWKid0apb/D4tqut8JQXsMMRRbx3ql9rRq6upSv6OJcTA9ZNLu3mxbPKpTwOlTIseZJDx+DHHhClkztEt+4dTNKodL80+/mGWcby4zL8sL7op2J5oFVtlva7Tl9zWAeZwazzxIa+OTOgydxh6iC6rdx3+kmiyj8sXBq4sIlqje8Hxjy2kBDlqxSRjG0foWMy7Jl6RTNMoOmrnidNzLIEdwH/+2R6SLL/fO7CpvL3y4qk288eMmA73FkybKj72luVTIsi9izaU2JPN/+S16aW8MwxmjzoCyBHgeMpBmVJcTfu1apIolLOtdsHa0WkQVZslKWV14nrORGmmWKOBSvUu3TBo2bDhw24vctu7zuP1GpBWQxRFnE8b2/0l5y6tr9hv1agizZL4u47McOqQ7ZmlrYtOk24MGzV8hiYLKoI54VMTXSWjD/pcdBgkoVEvQyICAwVqlCFmTJrCzRL32rFbNO5y2hgsVL/bxoXUhULLIYjCwnd6zQWkpqUbzSwN7dHe3t8lnEVdDOvnmbLjv/PaMynNNZkCXbZbn01/9kkvf35beDYwzoyCiPyzK2b/uUd7GRLMkDQWZi2m/+7zEqAVmQRXdZRHHKgE4ZgEUilcmmrNgpIosByKKO/KJCcW1DWvQccuLcJa9b1+aM/M5c6xReidx08qqdyIIsOssiCup9OzatWrVq+bJlc2fPGDSwT+PP6hayLSiXpXIyrpmN/eGbfsjyscsS/eSKg432Kdr5dp25m/hKM3NEzyRjbXal/SOUyIIsOh8NpdyLiQp9uW/TqjKO9qmchtv1J4HPDX3ksvie3G2tfdxjZP33pUeJl/qc25dsb3XtPxeRBVkyLUsCMMEv/NYsmulc1C7Jo8zK5XmMGlk+alm8Dm1Jcna2eZFTnoGJl0YFPXbKL9e+fPDMtciCLFkly9tePfSqWTLJB1IuPItEFgOWRVRFN6zsrH15+4ETkAVZslgWTRcObbcwfrf3fORRKLJ81LI8PrfXWnunxNj2n2taw2eiqlu9ytq/pFmfIciCLFkvi6iKGeFWN/GUqvsGMZ6Xl2WJfX6nmG0+7VHa9f9d177Dh7Surv1Lug6diSzIkvWyaHp2fLNp/GPWpdUANSO4H7ksr0Vlp88qaS82YPl+rTtV/P7rutqX/rhoM7IgS5qyOLo2Co7VU5Yo7+NmxkYSicmyIzcM49bJI7Icv5f6mXLLJw7WXqxR13HKdy8Y4uA21bVPZPrz9F1kQZZ3mzW1fRJZ5OXq+Kc2ha3PtVMdWzetXqNmm07dtvxzXCWkos+NfcuMZdKKX7iFxKqQJffK4p5cls2XH6S65Is75yy0TlmyKFknIDLxsSEMavnpuw8UFS3/MlaNLMiSuFVRnaq4JHnDuESNRyExKXaM1X2a1NZerLlb3yveT7RPWgnyvVWrTJH8tsXPefobzM1jkLJcO7At2QmOv5/yTH1RUZw2+BvtJTceezvUIigi65ZLPENXPnXVbsO4cZAla3r5+JZDQfOkHwux/Nc9xSuYqOpWu2qy0+HMCxSasmSdIv5TQmGB/i0/rSCR5Fu1+6Ah3T4GKcvy2SOT3ZWjFvyR1sIRL31bVSmcuGTrgePePnJ8rlonTIf8xTeDwg1lLxVZMpXfg7se7ld2b93Utm7NlGfQlqxR739bd7l7eDx+nnj+gnhiz1r7pPNaxu3gSKX2JUvXdq1Z2K5go1ZfH3H3FgxrOiiDkUUVE3Hzqsexf/6aNXFMsfz5kt2PFoWLzV6y8uzFy3fupTwsEiMCfAZ2bik3ijulwMjYYtKC5aeP/9v6C9e4d6JN8/UYMTkkWmkw9ziyZCbx2J7/TclARy8mGYgNefpoza9zurRv/WmVKi4uLmXKuFSoWKlFO7efZ849635bbYhTzBmMLDHBfotmvf8eX7Ntb+qPGEF97tiBUUMG1KleTXPXu5T9xLVew5/nLLpx18/AphZElpyESa1SRmmKjIpVKAx7ykq+fTXFs0AVd9dHR6sFw5xFEFkIWQhZCFkIWZCFkAVZkIWQhZCFkIWQBVmQhdsBWZCFkIWQhZCFkAVZkIWQBVkIWQhZCFkIWQhZCFmQhZCFkIWQhZCFkIWQBVkIWQhZCFkIWQhZCFmQhZCFkIWQhZCFkIWQBVkIWZAFWQhZCFkIWQhZkIWQBVmQhZCFkIWQhZAFWZCFkAVZCFkIWQhZCFmQBVkIWZCFkIWQhZCFkIWQhZAFWQhZCFkIWQhZCFkIWZCFkIWQhZCFkIWQhZAFWQhZCFkIWQhZCFkIWZCFkAVZkIWQhZCFkIWQBVkIWZAFWQhZCFkIWQhZkAVZCFmQhZCFkIWQhZAFWZCFkAVZCFkIWQhZCFkIWQhZkIWQhZCFkIWQhZCFkAVZCFkIWQhZCFkIWQhZkIWQhZCFkIWQhZCFkAVZCFmQBVkIWQhZCFkIWZCFkAVZskmWJk2afE95pipVqpQtW5bbIe/UvXv37JZFFMXx48d3pbxUp06dOnTowO2Qp5owYYJSqcw+WYiIkIWIkIWIkIWICFmICFmICFmIiJCFiJCFiJCFiJCFiCjL+j/g69xcEyj9XQAAAABJRU5ErkJggg==" } }, "cell_type": "markdown", "id": "6af4c1c2-3fea-4aac-824d-422deebf8899", "metadata": {}, "source": [ "## Quadtree Mesh\n", "\n", "A quadtree mesh is an unstructured refinement to a rectilinear modelgrid. In a quadtree mesh each node is split into 4 child nodes per level of refinement. The figure below is from [Lien and others, 2019](https://pubs.usgs.gov/of/2014/1109/pdf/ofr2014-1109.pdf) and shows multiple levels of recursive refinement.\n", "\n", "![Quadmesh.png](attachment:Quadmesh.png)\n", "\n", "For this example, we will be using GRIDGEN and the FloPy utility `Gridgen` to create quadtree refinement around the streams in our example." ] }, { "cell_type": "code", "execution_count": null, "id": "f7ab3319-2555-40cb-bc40-35760cb059e8", "metadata": {}, "outputs": [], "source": [ "gridgen_dir = data_ws / \"sagehen_gridgen\"\n", "gridgen_dir.mkdir(exist_ok=True)" ] }, { "cell_type": "markdown", "id": "b7f86d9b-9e30-44ab-8a3b-ae8160ed76cb", "metadata": {}, "source": [ "The `Gridgen` class a number of input parameters and further documentation/examples of the class can be found [here](https://flopy.readthedocs.io/en/latest/Notebooks/gridgen_example.html). For our example, we will be refining the area around the stream channels." ] }, { "cell_type": "code", "execution_count": null, "id": "950682df-493e-4d17-b7eb-5ad49341da0a", "metadata": {}, "outputs": [], "source": [ "g = Gridgen(sgrid, model_ws=gridgen_dir)\n", "for geom in nhddf.geometry.values:\n", " xy = list(zip(*geom.coords.xy))\n", " g.add_refinement_features([xy,], \"line\", 2, [0,])" ] }, { "cell_type": "code", "execution_count": null, "id": "45024b09-8e86-4ea3-a77d-23fcd5870bd4", "metadata": {}, "outputs": [], "source": [ "g.build()" ] }, { "cell_type": "markdown", "id": "afe3071c-6fa8-4411-8461-de554810f186", "metadata": {}, "source": [ "After gridgen builds the mesh, it can be used to create a `VertexGrid` (or DISV) that represents the basin" ] }, { "cell_type": "code", "execution_count": null, "id": "3d06162e-b420-4435-b616-d4b0f08759a3", "metadata": {}, "outputs": [], "source": [ "gridprops = g.get_gridprops_disv()\n", "gridprops.pop(\"nvert\")\n", "\n", "quadgrid = VertexGrid(**gridprops)\n", "quadgrid.extent" ] }, { "cell_type": "markdown", "id": "845ea4b5-2e51-4d45-a311-45e727f1c9dc", "metadata": {}, "source": [ "Resample top elevations to the quadtree mesh and intersect the basin boundary and stream vectors" ] }, { "cell_type": "code", "execution_count": null, "id": "cb60bb16-8381-4c72-a97b-61e0c362170f", "metadata": {}, "outputs": [], "source": [ "# create the top array\n", "top = rstr.resample_to_grid(quadgrid, band=rstr.bands[0], method=\"min\")\n", "quadgrid._top = top" ] }, { "cell_type": "code", "execution_count": null, "id": "8f977f97-38dd-40aa-8cf1-b9f89824f0ca", "metadata": {}, "outputs": [], "source": [ "# create the idomain array\n", "idomain = np.zeros(quadgrid.shape, dtype=int)\n", "ix = GridIntersect(quadgrid, method=\"vertex\")\n", "nodes = ix.intersect(basindf.geometry.values[0], contains_centroid=True)[\"cellids\"]\n", "\n", "idomain[:, list(nodes)] = 1\n", "quadgrid._idomain = idomain" ] }, { "cell_type": "code", "execution_count": null, "id": "9dfaa463-e79f-451d-855a-d033e1957552", "metadata": {}, "outputs": [], "source": [ "# create an array of stream cells\n", "qstr = []\n", "ix = GridIntersect(quadgrid, method=\"vertex\")\n", "for geom in nhddf.geometry.values:\n", " nodes = ix.intersects(geom)[\"cellids\"]\n", " qstr.extend(list(nodes))" ] }, { "cell_type": "code", "execution_count": null, "id": "0d004d92-1e98-43d6-8034-7e3f4e73451a", "metadata": {}, "outputs": [], "source": [ "qstrms = np.zeros((quadgrid.ncpl,), dtype=int)\n", "qstrms[qstr] = 1" ] }, { "cell_type": "markdown", "id": "18574f1f-46f1-4343-82b7-c61debc16655", "metadata": {}, "source": [ "Plot the quadtree mesh with top elevations and streams" ] }, { "cell_type": "code", "execution_count": null, "id": "73679f36-5d8d-427c-9706-75b57bf5d4f8", "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots(figsize=(5, 5))\n", "\n", "pmv = PlotMapView(modelgrid=quadgrid)\n", "pc = pmv.plot_array(quadgrid.top)\n", "pmv.plot_array(qstrms, masked_values=[0,], alpha=0.3, cmap=\"Reds_r\")\n", "pmv.plot_inactive()\n", "pmv.plot_grid(lw=0.3)\n", "nhddf.plot(ax=ax)\n", "plt.colorbar(pc, shrink=0.7);" ] }, { "cell_type": "markdown", "id": "aa833f5b-ce9a-4d68-9e4b-832aea1c09f9", "metadata": {}, "source": [ "## Triangular Mesh\n", "\n", "An unstructured/vertex based triangular mesh can be generated using the code \"triangle\"([Shewchuk, 2002](https://www.sciencedirect.com/science/article/pii/S0925772101000475)). And the FloPy utility, `Triangle` that creates inputs for and translates the outputs from \"triangle\". \n", "\n", "For this example, we will be using triangle to generate a mesh and add refinement around the streams." ] }, { "cell_type": "code", "execution_count": null, "id": "dc3ab2b7-378e-4eb8-8a9a-fa1ea7ae9bdd", "metadata": {}, "outputs": [], "source": [ "tri_dir = data_ws / \"sagehen_tri\"\n", "tri_dir.mkdir(exist_ok=True)" ] }, { "cell_type": "markdown", "id": "1e8987ed-1970-480a-a796-7bd71c86af09", "metadata": {}, "source": [ "For this example, we will use the watershed boundary we downloaded from NLDI and create a dissovled polygon of stream segments to generate our triangulation regions." ] }, { "cell_type": "code", "execution_count": null, "id": "6ffe2440-b82c-448e-8735-d89651b562ea", "metadata": {}, "outputs": [], "source": [ "basindf.head()" ] }, { "cell_type": "markdown", "id": "1fa96714-8444-427e-97df-a435b9d4a7f9", "metadata": {}, "source": [ "We will also dissolve and buffer the stream network into a single polygon that defines our area of refinement" ] }, { "cell_type": "code", "execution_count": null, "id": "e75ec798-9cc4-4e3b-8712-030d54e168dd", "metadata": {}, "outputs": [], "source": [ "bufgdf = nhddf.dissolve()\n", "bufgdf[\"geometry\"] = bufgdf.geometry.buffer(100, cap_style=2, join_style=3)\n", "bufgdf.plot();" ] }, { "cell_type": "markdown", "id": "e8571d37-71fd-4324-9905-13565b2c8835", "metadata": {}, "source": [ "These two polygons can now be used to generate a triangular mesh with the `Triangle` utility. The `Triangle` class has a number of input parameters. For more detail on the inputs and additional examples please see this [notebook](https://flopy.readthedocs.io/en/latest/Notebooks/dis_voronoi_example.html)" ] }, { "cell_type": "code", "execution_count": null, "id": "1668863c-03fe-46d7-9f1e-ede287eacb36", "metadata": {}, "outputs": [], "source": [ "# define point locations within the watershed and the stream for adding regions\n", "wsloc = (220000, 4368000)\n", "stloc = (219250, 4370000)" ] }, { "cell_type": "code", "execution_count": null, "id": "3da8bd04-c9cf-4b59-af7a-b54073aa9be9", "metadata": {}, "outputs": [], "source": [ "tri = Triangle(angle=30, model_ws=tri_dir)\n", "\n", "# define the model/mesh boundary\n", "tri.add_polygon(basindf.geometry.values[0])\n", "tri.add_region(wsloc, 0, maximum_area=100 * 300)\n", "\n", "# define the stream refinement area\n", "tri.add_polygon(bufgdf.geometry.values[0])\n", "tri.add_region(stloc, 1, maximum_area=40 * 40)" ] }, { "cell_type": "code", "execution_count": null, "id": "a8a020c5-0878-4a19-aab0-3aa2ff0bac68", "metadata": {}, "outputs": [], "source": [ "tri.build()" ] }, { "cell_type": "markdown", "id": "a5bbe12a-d44c-479f-b76a-6d17f1f88128", "metadata": {}, "source": [ "Visualize the triangular mesh" ] }, { "cell_type": "code", "execution_count": null, "id": "0cd510b8-901c-4648-af86-1986434017fa", "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots(figsize=(6, 6))\n", "tri.plot(ax=ax)\n", "nhddf.plot(ax=ax);" ] }, { "cell_type": "markdown", "id": "22c159b8-185b-41b4-85a9-5dd57b65d7aa", "metadata": {}, "source": [ "Create a `VertexGrid` from the triangular mesh and resample the land surface elevations" ] }, { "cell_type": "code", "execution_count": null, "id": "647735f0-4544-4e3a-b207-be3970661a71", "metadata": {}, "outputs": [], "source": [ "cell2d = tri.get_cell2d()\n", "vertices = tri.get_vertices()\n", "ncpl = len(cell2d)\n", "nlay = 1\n", "idomain = np.ones((nlay, ncpl), dtype=int)\n", "\n", "# set fake values for top and botm for now)\n", "top = np.ones((ncpl,))\n", "botm = np.zeros((nlay, ncpl))" ] }, { "cell_type": "code", "execution_count": null, "id": "0934cc74-7c65-4546-9707-51643f0b1d85", "metadata": {}, "outputs": [], "source": [ "trigrid = VertexGrid(\n", " vertices=vertices,\n", " cell2d=cell2d,\n", " ncpl=ncpl,\n", " nlay=nlay,\n", " idomain=idomain,\n", " top=top,\n", " botm=botm,\n", " crs=\"EPSG:26911\",\n", ")\n", "trigrid.extent" ] }, { "cell_type": "code", "execution_count": null, "id": "3fe7d3e9-4bfa-454d-9ff3-2e54aaa18a39", "metadata": {}, "outputs": [], "source": [ "top = rstr.resample_to_grid(\n", " trigrid, band=rstr.bands[0], method=\"min\", extrapolate_edges=True\n", ")\n", "trigrid._top = top" ] }, { "cell_type": "markdown", "id": "8522ed38-a10a-42ff-82cf-2dba5d91e56f", "metadata": {}, "source": [ "And perform intersection to identify stream cells" ] }, { "cell_type": "code", "execution_count": null, "id": "d870c5cf-8ce1-4692-93cc-49e18b90b05f", "metadata": {}, "outputs": [], "source": [ "tristr = []\n", "ix = GridIntersect(trigrid, method=\"vertex\")\n", "for geom in nhddf.geometry.values:\n", " nodes = ix.intersects(geom)[\"cellids\"]\n", " tristr.extend(list(nodes))\n", "\n", "tristrms = np.zeros((trigrid.ncpl,), dtype=int)\n", "tristrms[tristr] = 1" ] }, { "cell_type": "markdown", "id": "54e52bc2-ee4a-42c5-ad20-68656c78896e", "metadata": {}, "source": [ "Now to plot up the triangular grid" ] }, { "cell_type": "code", "execution_count": null, "id": "cb064b04-fd24-4830-9eb3-69a67a2f61eb", "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots(figsize=(7, 7))\n", "\n", "pmv = PlotMapView(modelgrid=trigrid)\n", "pc = pmv.plot_array(trigrid.top)\n", "pmv.plot_array(tristrms, masked_values=[0,], alpha=0.3, cmap=\"Reds_r\")\n", "pmv.plot_inactive()\n", "pmv.plot_grid(lw=0.3)\n", "nhddf.plot(ax=ax)\n", "plt.colorbar(pc, shrink=0.7);" ] }, { "cell_type": "markdown", "id": "34736781-43ea-483c-acbf-962bb0b507ef", "metadata": {}, "source": [ "## Voronoi Mesh\n", "\n", "A voronoi mesh can be generated from the circumcirle centroids of each triangular node in a Delaunay triangulation. The centroid of each circumcirle produces a vertex for the voronoi mesh and these vertices are connected via adjacency relationships of the triangulation. \n", "\n", "FloPy has a built in utility `VoronoiGrid` that produces voronoi meshes from Delaunay triangulations. For more information and additional examples on `VoronoiGrid`'s usage, see this [notebook](https://flopy.readthedocs.io/en/latest/Notebooks/dis_voronoi_example.html)" ] }, { "cell_type": "markdown", "id": "6c0852aa-8d11-4c7c-b449-30921ca11216", "metadata": {}, "source": [ "To generate a voronoi mesh, we just need to pass our triangluation object (`tri`) to the `VoronoiGrid` utility" ] }, { "cell_type": "code", "execution_count": null, "id": "ed25dcc4-8b4a-4162-b603-35f605fe45e2", "metadata": {}, "outputs": [], "source": [ "vor = VoronoiGrid(tri)" ] }, { "cell_type": "markdown", "id": "345a997b-9584-4a4f-8f5d-4b5c657efbf2", "metadata": {}, "source": [ "Build a `VertexGrid`" ] }, { "cell_type": "code", "execution_count": null, "id": "9bba5fd7-8a3e-4ff0-8002-1ee521c040f3", "metadata": {}, "outputs": [], "source": [ "gridprops = vor.get_gridprops_vertexgrid()\n", "nlay = 1\n", "idomain = np.ones(gridprops[\"ncpl\"], dtype=int)\n", "top = idomain.copy().astype(float)\n", "botm = np.zeros((nlay, gridprops[\"ncpl\"]))" ] }, { "cell_type": "code", "execution_count": null, "id": "c87424f9-d9e1-4798-b5a1-99ec2b35f87d", "metadata": {}, "outputs": [], "source": [ "vorgrid = VertexGrid(\n", " nlay=nlay, idomain=idomain, top=top, botm=botm, **gridprops\n", ")" ] }, { "cell_type": "markdown", "id": "d2cd1136-11ba-46d2-ad45-80db9408368e", "metadata": {}, "source": [ "Resample the land surface elevation from the DEM to the voronoi grid" ] }, { "cell_type": "code", "execution_count": null, "id": "c42a08d0-71b2-4f38-95d4-be4b9e3356bf", "metadata": {}, "outputs": [], "source": [ "top = rstr.resample_to_grid(\n", " vorgrid, band=rstr.bands[0], method=\"min\", extrapolate_edges=True\n", ")\n", "vorgrid._top = top" ] }, { "cell_type": "markdown", "id": "997be18a-d2b7-42a1-ad8f-49285d36de49", "metadata": {}, "source": [ "Identify stream cells with `GridIntersect`" ] }, { "cell_type": "code", "execution_count": null, "id": "93dd9221-4f32-4b90-9adb-13426ef4dcaa", "metadata": {}, "outputs": [], "source": [ "vorst = []\n", "ix = GridIntersect(vorgrid, method=\"vertex\")\n", "for geom in nhddf.geometry.values:\n", " nodes = ix.intersects(geom)[\"cellids\"]\n", " vorst.extend(list(nodes))\n", "\n", "vorstrms = np.zeros((vorgrid.ncpl,), dtype=int)\n", "vorstrms[vorst] = 1" ] }, { "cell_type": "markdown", "id": "0f749c17-6db8-49d3-9e35-7309f183f8db", "metadata": {}, "source": [ "Finally plot the voronoi mesh with the DEM data and stream locations" ] }, { "cell_type": "code", "execution_count": null, "id": "a6ba684f-498a-4e7e-b10c-ad61ba8a9f1f", "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots(figsize=(7, 7))\n", "\n", "pmv = PlotMapView(modelgrid=vorgrid)\n", "pc = pmv.plot_array(vorgrid.top)\n", "pmv.plot_array(vorstrms, masked_values=[0], alpha=0.3, cmap=\"Reds_r\")\n", "pmv.plot_inactive()\n", "pmv.plot_grid(lw=0.3)\n", "nhddf.plot(ax=ax)\n", "plt.colorbar(pc, shrink=0.7);" ] }, { "cell_type": "markdown", "id": "82758f91-a8dc-43f4-a729-83023ef5c4a8", "metadata": {}, "source": [ "### plot the meshes we produced side by side for comparison!" ] }, { "cell_type": "code", "execution_count": null, "id": "c106ccf9-8950-4164-a398-d15db0daecac", "metadata": {}, "outputs": [], "source": [ "from flopy.plot import styles\n", "from matplotlib.gridspec import GridSpec\n", "\n", "with styles.USGSMap():\n", " fig, axm = plt.subplots(nrows=3, ncols=2, figsize=(10,12))\n", " axs = axm.ravel()\n", " \n", " # structured grid\n", " pmv = PlotMapView(modelgrid=sgrid, ax=axs[0])\n", " pmv.plot_array(sgrid.top, cmap=\"viridis\")\n", " pmv.plot_array(strms, masked_values=[0,], alpha=0.3, cmap=\"Reds_r\")\n", " pmv.plot_inactive()\n", " pmv.plot_grid(lw=0.3)\n", " nhddf.plot(ax=axs[0])\n", " styles.heading(axs[0], \"A. \", heading=\"Structured Grid\", x=0.08)\n", "\n", " # non-uniform DIS\n", " pmv = PlotMapView(modelgrid=rfgrid, ax=axs[1])\n", " pmv.plot_array(rfgrid.top, cmap=\"viridis\")\n", " pmv.plot_inactive()\n", " pmv.plot_grid(lw=0.3)\n", " nhddf.plot(ax=axs[1])\n", " styles.heading(axs[1], \"B.\", heading=\"Structured Grid: non-uniform\", x=0.08)\n", "\n", " # LGR grid\n", " pmvp = PlotMapView(modelgrid=sgrid, ax=axs[2])\n", " ptop = sgrid.top.copy()\n", " ptop[crow, ccol] = np.nan\n", " strms[crow, ccol] = np.nan\n", " pc = pmvp.plot_array(ptop, cmap=\"viridis\", vmin=vmin, vmax=vmax)\n", " pmvp.plot_array(strms, masked_values=[0], alpha=0.3, cmap=\"Reds_r\")\n", " pmvp.plot_inactive()\n", " pmvp.plot_grid(lw=0.3)\n", " nhddf.plot(ax=axs[2])\n", "\n", " pmvc = PlotMapView(modelgrid=childgrid, ax=axs[2])\n", " pmvc.plot_array(childtop, cmap=\"viridis\", vmin=vmin, vmax=vmax)\n", " pmvc.plot_array(cstrms, masked_values=[0], alpha=0.3, cmap=\"Reds_r\")\n", " pmvc.plot_grid(lw=0.3)\n", " styles.heading(axs[2], \"C. \", heading=\"LGR Grid\", x=0.08)\n", "\n", " # Quadtree grid\n", " pmv = PlotMapView(modelgrid=quadgrid, ax=axs[3])\n", " pc = pmv.plot_array(quadgrid.top)\n", " pmv.plot_array(qstrms, masked_values=[0,], alpha=0.3, cmap=\"Reds_r\")\n", " pmv.plot_inactive()\n", " pmv.plot_grid(lw=0.2)\n", " nhddf.plot(ax=axs[3])\n", " styles.heading(axs[3], \"D. \", heading=\"Quadtree Grid\", x=0.08)\n", "\n", " # Triangular grid\n", " pmv = PlotMapView(modelgrid=trigrid, ax=axs[4])\n", " pc = pmv.plot_array(trigrid.top)\n", " pmv.plot_array(tristrms, masked_values=[0,], alpha=0.3, cmap=\"Reds_r\")\n", " pmv.plot_inactive()\n", " pmv.plot_grid(lw=0.3)\n", " nhddf.plot(ax=axs[4])\n", " styles.heading(axs[4], \"E. \", heading=\"Triangular Grid\", x=0.08)\n", "\n", " # Voronoi grid\n", " pmv = PlotMapView(modelgrid=vorgrid, ax=axs[5])\n", " pc = pmv.plot_array(vorgrid.top)\n", " pmv.plot_array(vorstrms, masked_values=[0,], alpha=0.3, cmap=\"Reds_r\")\n", " pmv.plot_inactive()\n", " pmv.plot_grid(lw=0.3)\n", " nhddf.plot(ax=axs[5])\n", " styles.heading(axs[5], \"F. \", heading=\"Voronoi Grid\", x=0.08)\n", "\n", " # colorbar\n", "\n", " \n", " # ax5.spines[[\"left\", \"right\", \"top\", \"bottom\"]].set_visible(False)\n", " fig.colorbar(\n", " pc,\n", " ax=axm,\n", " location=\"top\",\n", " label=\"Land surface elevation, in meters\",\n", " shrink=0.5,\n", " fraction=.05\n", " )\n" ] }, { "cell_type": "markdown", "id": "be66dfa8-8e48-4614-b1f3-86baafc73b14", "metadata": {}, "source": [ "### More information \n", "For more inforation and and examples on how to use the utilities presented in this notebook please see:\n", "\n", "#### Examples\n", " - [FloPy Examples Gallery](https://flopy.readthedocs.io/en/latest/examples.html)\n", "\n", "#### API reference\n", " - [FloPy utils API reference](https://flopy.readthedocs.io/en/latest/code.html#flopy-utilities)\n", " - [FloPy lgrutil API reference](https://flopy.readthedocs.io/en/latest/source/flopy.utils.lgrutil.html)\n", " - [FloPy gridgen API reference](https://flopy.readthedocs.io/en/latest/source/flopy.utils.gridgen.html)\n", " - [FloPy triangle API reference](https://flopy.readthedocs.io/en/latest/source/flopy.utils.triangle.html)\n", " - [FloPy voronoi API reference](https://flopy.readthedocs.io/en/latest/source/flopy.utils.voronoi.html)" ] }, { "cell_type": "markdown", "id": "67c62f41", "metadata": {}, "source": [ "## Activity (If time permits); Build a Voronoi mesh/Vertex Grid:\n", "\n", "Define a voronoi grid for the basin upstream of the Black Earth Creek streamgauge near Cross Plains, WI. The upstream basin information is supplied by NLDI and has been saved to the class repo in case of internet connectivity or NLDI service issues.\n", "\n", "Some of the previous code can be re-used to create a voronoi mesh and subsequent `VertexGrid` representation of the basin. See the Triangular and Voronoi grid sections of the notebook.\n", "\n", "**Optional**, use the Upstream Main segment of Black Earth Creek to refine the voronoi mesh around the stream location. Refinement control points have been pre-defined below." ] }, { "cell_type": "code", "execution_count": null, "id": "d58d2db7", "metadata": {}, "outputs": [], "source": [ "epsg = 5070 # NAD83 CONUS Albers\n", "# try:\n", " # todo: update this for black earth creek near cross plains WI site-no 05406479\n", "be_cr_nr_cp = \"05406479\"\n", "basindf = nldi.get_basin(\n", " feature_source=\"nwissite\", feature_id=f\"USGS-{be_cr_nr_cp}\"\n", ")\n", "basindf = basindf.to_crs(epsg=epsg)\n", "basindf.to_file(geospatial_ws / \"black_earth_basin.shp\")\n", "# except (ValueError, ConnectionError, NameError):\n", "# basindf = gpd.read_file(geospatial_ws / \"black_earth_basin.shp\")" ] }, { "cell_type": "code", "execution_count": null, "id": "47c81ac1-e436-4685-a943-dffe9e023f1a", "metadata": {}, "outputs": [], "source": [ "try:\n", " nhddf = nldi.get_flowlines(\n", " navigation_mode=\"UM\",\n", " distance=9999,\n", " feature_source=\"nwissite\",\n", " feature_id=f\"USGS-{be_cr_nr_cp}\",\n", " )\n", " nhddf = nhddf.to_crs(epsg=epsg)\n", " nhddf.to_file(geospatial_ws / \"black_earth_creek_main.shp\")\n", "except (ValueError, ConnectionError, NameError):\n", " nhddf = gpd.read_file(geospatial_ws / \"black_earth_creek_main.shp\")" ] }, { "cell_type": "markdown", "id": "dcb74970-e19d-40b0-8443-f68f11e5a615", "metadata": {}, "source": [ "Visualize the example \"study area\"" ] }, { "cell_type": "code", "execution_count": null, "id": "32c1f501-88e3-41e4-90fb-1f5d4299e984", "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots()\n", "\n", "basindf.plot(ax=ax, facecolor=\"None\", edgecolor=\"k\")\n", "nhddf.plot(ax=ax)\n", "ctx.add_basemap(\n", " ax=ax,\n", " crs=basindf.crs,\n", " source=ctx.providers.USGS.USTopo\n", ")" ] }, { "cell_type": "markdown", "id": "cc6b79b0-198a-410c-9428-13fa814d3052", "metadata": {}, "source": [ "Grid refinement control points are defined here for the optional part of this exercise" ] }, { "cell_type": "code", "execution_count": null, "id": "4bacc5c8-6b0e-436f-8c67-42ae6e393bf5", "metadata": {}, "outputs": [], "source": [ "wsloc = basindf.centroid.values[0]\n", "stloc = list(zip(*nhddf.geometry.values[10].xy))[0]" ] }, { "cell_type": "code", "execution_count": null, "id": "6b6d8a83-e226-4020-8889-0e6c35bd3af3", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "id": "f5bb2150-fe6e-4fe1-918b-7e09239d3694", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "id": "ccecfea2-121b-442f-bfca-0e345ccc5173", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" } }, "nbformat": 4, "nbformat_minor": 5 }