FloPy Creating Layered Quadtree Grids with GRIDGEN

FloPy has a module that can be used to drive the GRIDGEN program. This notebook shows how it works.

The Flopy GRIDGEN module requires that the gridgen executable can be called using subprocess.



In [1]:

    
%matplotlib inline
import os
import numpy as np
import matplotlib.pyplot as plt
import flopy
from flopy.utils.gridgen import Gridgen

Setup Base MODFLOW Grid

GRIDGEN works off of a base MODFLOW grid. The following information defines the basegrid.



In [2]:

    
Lx = 100.
Ly = 100.
nlay = 2
nrow = 51
ncol = 51
delr = Lx / ncol
delc = Ly / nrow
h0 = 10
h1 = 5
top = h0
botm = np.zeros((nlay, nrow, ncol), dtype=np.float32)
botm[1, :, :] = -10.



In [3]:

    
ms = flopy.modflow.Modflow(rotation=20.)
dis = flopy.modflow.ModflowDis(ms, nlay=nlay, nrow=nrow, ncol=ncol, delr=delr,
                               delc=delc, top=top, botm=botm)

Create the Gridgen Object



In [4]:

    
model_ws = os.path.join('.', 'data')
g = Gridgen(dis, model_ws=model_ws)

Add an Optional Active Domain

Cells outside of the active domain will be clipped and not numbered as part of the final grid. If this step is not performed, then all cells will be included in the final grid.



In [5]:

    
# setup the active domain
adshp = os.path.join(model_ws, 'ad0')
adpoly = [[[(0, 0), (0, 60), (40, 80), (60, 0), (0, 0)]]]
# g.add_active_domain(adpoly, range(nlay))

Refine the Grid



In [6]:

    
x = Lx * np.random.random(10)
y = Ly * np.random.random(10)
wells = zip(x, y)
g.add_refinement_features(wells, 'point', 3, range(nlay))



In [7]:

    
river = [[[(-20, 10), (60, 60)]]]
g.add_refinement_features(river, 'line', 3, range(nlay))



In [8]:

    
g.add_refinement_features(adpoly, 'polygon', 1, range(nlay))

Build the Grid



In [9]:

    
g.build()

Plot the Grid



In [10]:

    
fig = plt.figure(figsize=(15, 15))
ax = fig.add_subplot(1, 1, 1, aspect='equal')
g.plot(ax, linewidth=0.5)
flopy.plot.plot_shapefile(adshp, ax=ax, facecolor='yellow', alpha=0.1)









    Out[10]:





<matplotlib.collections.PatchCollection at 0x10d458e90>

Create a Flopy ModflowDisu Object



In [11]:

    
mu = flopy.modflow.Modflow(model_ws=model_ws, modelname='mfusg')
disu = g.get_disu(mu)
# disu.write_file()
# print(disu)

Intersect Features with the Grid



In [12]:

    
adpoly_intersect = g.intersect(adpoly, 'polygon', 0)
print(adpoly_intersect.dtype.names)
print(adpoly_intersect)
print(adpoly_intersect.nodenumber)









    



('nodenumber', 'polyid', 'totalarea', 'SHAPEID')
[(713, 0, 0.961169, 0) (712, 0, 0.961169, 0) (716, 0, 0.961169, 0) ...,
 (6885, 0, 0.792964, 0) (6888, 0, 0.0768934, 0) (6887, 0, 0.956363, 0)]
[ 713  712  716 ..., 6885 6888 6887]



In [13]:

    
well_intersect = g.intersect(wells, 'point', 0)
print(well_intersect.dtype.names)
print(well_intersect)
print(well_intersect.nodenumber)









    



('nodenumber', 'pointid', 'SHAPEID')
[(290, 0, 0) (995, 1, 1) (237, 2, 2) (769, 3, 3) (602, 4, 4) (45, 5, 5)
 (5744, 6, 6) (955, 7, 7) (1952, 8, 8) (383, 9, 9)]
[ 290  995  237  769  602   45 5744  955 1952  383]



In [14]:

    
river_intersect = g.intersect(river, 'line', 0)
print(river_intersect.dtype.names)
# print(river_intersect)
# print(river_intersect.nodenumber)









    



('nodenumber', 'arcid', 'length', 'starting_distance', 'ending_distance', 'SHAPEID')

Plot Intersected Features



In [15]:

    
a = np.zeros((g.nodes), dtype=np.int)
a[adpoly_intersect.nodenumber] = 1
a[well_intersect.nodenumber] = 2
a[river_intersect.nodenumber] = 3
fig = plt.figure(figsize=(15, 15))
ax = fig.add_subplot(1, 1, 1, aspect='equal')
g.plot(ax, a=a, masked_values=[0], edgecolor='none', cmap='jet')









    Out[15]:





<matplotlib.collections.PatchCollection at 0x10d799c10>



In [ ]: