Skip to content

Commit

Permalink
feat(disl grids): support for 1d vertex grids. fix for writing gridli…
Browse files Browse the repository at this point in the history 
…nes to shapefile (#799)

* feat(disl grids): support for 1d vertex grids. fix for writing gridlines to shapefile

* fix(parameter ordering): the adding of a new parameter in the VertexGrid constructor was breaking some notebook examples. new parameter moved to end of parameter list to support backward compatibility.
  • Loading branch information
@spaulins-usgs
spaulins-usgs committed Feb 4, 2020
1 parent ab8120b commit 469727b
Show file tree
Hide file tree
Showing 8 changed files with 138 additions and 51 deletions.
96 changes: 66 additions & 30 deletions flopy/discretization/vertexgrid.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -32,13 +32,14 @@ class for a vertex model grid
returns vertices for a single cell at cellid.
"""

def __init__(self, vertices=None, cell2d=None, top=None, botm=None, idomain=None,
lenuni=None, epsg=None, proj4=None, prj=None, xoff=0.0,
yoff=0.0, angrot=0.0, grid_type='vertex',
nlay=None, ncpl=None):
def __init__(self, vertices=None, cell2d=None, top=None,
botm=None, idomain=None, lenuni=None, epsg=None, proj4=None,
prj=None, xoff=0.0, yoff=0.0, angrot=0.0, grid_type='vertex',
nlay=None, ncpl=None, cell1d=None):
super(VertexGrid, self).__init__(grid_type, top, botm, idomain, lenuni,
epsg, proj4, prj, xoff, yoff, angrot)
self._vertices = vertices
self._cell1d = cell1d
self._cell2d = cell2d
self._top = top
self._botm = botm
Expand All @@ -52,26 +53,32 @@ def __init__(self, vertices=None, cell2d=None, top=None, botm=None, idomain=None

@property
def is_valid(self):
if self._vertices is not None and self._cell2d is not None:
if self._vertices is not None and (self._cell2d is not None or
self._cell1d is not None):
return True
return False

@property
def is_complete(self):
if self._vertices is not None and self._cell2d is not None and \
if self._vertices is not None and (self._cell2d is not None or
self._cell1d is not None) and \
super(VertexGrid, self).is_complete:
return True
return False

@property
def nlay(self):
if self._botm is not None:
if self._cell1d is not None:
return 1
elif self._botm is not None:
return len(self._botm)
else:
return self._nlay

@property
def ncpl(self):
if self._cell1d is not None:
return len(self._cell1d)
if self._botm is not None:
return len(self._botm[0])
else:
Expand Down Expand Up @@ -233,34 +240,62 @@ def _build_grid_geometry_info(self):
cache_index_cc = 'cellcenters'
cache_index_vert = 'xyzgrid'

vertexdict = {v[0]: [v[1], v[2]]
for v in self._vertices}
xcenters = []
ycenters = []
xvertices = []
yvertices = []

# build xy vertex and cell center info
for cell2d in self._cell2d:
cell2d = tuple(cell2d)
xcenters.append(cell2d[1])
ycenters.append(cell2d[2])

vert_number = []
for i in cell2d[4:]:
if i is not None:
vert_number.append(int(i))
if self._cell1d is not None:
zcenters = []
zvertices = []
vertexdict = {v[0]: [v[1], v[2], v[3]]
for v in self._vertices}
for cell1d in self._cell1d:
cell1d = tuple(cell1d)
xcenters.append(cell1d[1])
ycenters.append(cell1d[2])
zcenters.append(cell1d[3])

vert_number = []
for i in cell1d[3:]:
if i is not None:
vert_number.append(int(i))

xcellvert = []
ycellvert = []
zcellvert = []
for ix in vert_number:
xcellvert.append(vertexdict[ix][0])
ycellvert.append(vertexdict[ix][1])
zcellvert.append(vertexdict[ix][2])
xvertices.append(xcellvert)
yvertices.append(ycellvert)
zvertices.append(zcellvert)

xcellvert = []
ycellvert = []
for ix in vert_number:
xcellvert.append(vertexdict[ix][0])
ycellvert.append(vertexdict[ix][1])
xvertices.append(xcellvert)
yvertices.append(ycellvert)

# build z cell centers
zvertices, zcenters = self._zcoords()
else:
vertexdict = {v[0]: [v[1], v[2]]
for v in self._vertices}
# build xy vertex and cell center info
for cell2d in self._cell2d:
cell2d = tuple(cell2d)
xcenters.append(cell2d[1])
ycenters.append(cell2d[2])

vert_number = []
for i in cell2d[4:]:
if i is not None:
vert_number.append(int(i))

xcellvert = []
ycellvert = []
for ix in vert_number:
xcellvert.append(vertexdict[ix][0])
ycellvert.append(vertexdict[ix][1])
xvertices.append(xcellvert)
yvertices.append(ycellvert)

# build z cell centers
zvertices, zcenters = self._zcoords()

if self._has_ref_coordinates:
# transform x and y
Expand All @@ -270,7 +305,8 @@ def _build_grid_geometry_info(self):
# vertices are a list within a list
for xcellvertices, ycellvertices in zip(xvertices, yvertices):
xcellvertices, \
ycellvertices = self.get_coords(xcellvertices, ycellvertices)
ycellvertices = self.get_coords(xcellvertices,
ycellvertices)
xvertxform.append(xcellvertices)
yvertxform.append(ycellvertices)
xvertices = xvertxform
Expand Down
4 changes: 2 additions & 2 deletions flopy/export/shapefile_utils.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -51,9 +51,9 @@ def write_gridlines_shapefile(filename, mg):
" instead.",
category=DeprecationWarning)
else:
grid_lines = mg.grid_lines()
grid_lines = mg.grid_lines
for i, line in enumerate(grid_lines):
wr.poly([line])
wr.line([line])
wr.record(i)

wr.close()
Expand Down
4 changes: 4 additions & 0 deletions flopy/mf6/coordinates/modeldimensions.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -336,6 +336,10 @@ def _create_model_grid(self, grid_type):
elif grid_type == DiscretizationType.DISU:
self._model_grid = UnstructuredModelGrid(self.model_name,
self.simulation_data)
elif grid_type == DiscretizationType.DISL:
self._model_grid = ModelGrid(self.model_name,
self.simulation_data,
DiscretizationType.DISL)
else:
self._model_grid = ModelGrid(self.model_name,
self.simulation_data,
Expand Down
44 changes: 28 additions & 16 deletions flopy/mf6/coordinates/modelgrid.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -401,6 +401,10 @@ def get_grid_type(simulation_data, model_name):
'disu{}'.format(structure.get_version_string()),
0) is not None:
return DiscretizationType.DISU
elif package_recarray.search_data(
'disl{}'.format(structure.get_version_string()),
0) is not None:
return DiscretizationType.DISL

return DiscretizationType.UNDEFINED

Expand All @@ -411,6 +415,9 @@ def get_idomain(self):
elif self._grid_type == DiscretizationType.DISV:
return self._simulation_data.mfdata[
(self._model_name, 'disv', 'griddata', 'idomain')].get_data()
elif self._grid_type == DiscretizationType.DISL:
return self._simulation_data.mfdata[
(self._model_name, 'disl', 'griddata', 'idomain')].get_data()
elif self._grid_type == DiscretizationType.DISU:
except_str = 'ERROR: Can not return idomain for model {}. This ' \
'model uses a DISU grid that does not ' \
Expand Down Expand Up @@ -447,10 +454,11 @@ def get_horizontal_cross_section_dim_arrays(self):
np.arange(1, self.num_columns() + 1, 1, np.int32)]
elif self.grid_type() == DiscretizationType.DISV:
return [np.arange(1, self.num_cells_per_layer() + 1, 1, np.int32)]
elif self.grid_type() == DiscretizationType.DISU:
elif self.grid_type() == DiscretizationType.DISU or \
self.grid_type() == DiscretizationType.DISL:
except_str = 'ERROR: Can not get horizontal plane arrays for ' \
'model "{}" DISU grid. DISU grids do not support ' \
'individual layers.'.format(self._model_name)
'model "{}" grid. DISU and DISL grids do not ' \
'support individual layers.'.format(self._model_name)
print(except_str)
raise MFGridException(except_str)

Expand All @@ -459,7 +467,8 @@ def get_model_dim(self):
return [self.num_layers(), self.num_rows(), self.num_columns()]
elif self.grid_type() == DiscretizationType.DISV:
return [self.num_layers(), self.num_cells_per_layer()]
elif self.grid_type() == DiscretizationType.DISU:
elif self.grid_type() == DiscretizationType.DISU or \
self.grid_type() == DiscretizationType.DISL:
return [self.num_cells()]

def get_model_dim_arrays(self):
Expand All @@ -470,7 +479,8 @@ def get_model_dim_arrays(self):
elif self.grid_type() == DiscretizationType.DISV:
return [np.arange(1, self.num_layers() + 1, 1, np.int32),
np.arange(1, self.num_cells_per_layer() + 1, 1, np.int32)]
elif self.grid_type() == DiscretizationType.DISU:
elif self.grid_type() == DiscretizationType.DISU or \
self.grid_type() == DiscretizationType.DISL:
return [np.arange(1, self.num_cells() + 1, 1, np.int32)]

def get_row_array(self):
Expand All @@ -487,7 +497,8 @@ def get_horizontal_cross_section_dim_names(self):
return ['row', 'column']
elif self.grid_type() == DiscretizationType.DISV:
return ['layer_cell_num']
elif self.grid_type() == DiscretizationType.DISU:
elif self.grid_type() == DiscretizationType.DISU or \
self.grid_type() == DiscretizationType.DISL:
except_str = 'ERROR: Can not get layer dimension name for model ' \
'"{}" DISU grid. DISU grids do not support ' \
'layers.'.format(self._model_name)
Expand All @@ -499,23 +510,19 @@ def get_model_dim_names(self):
return ['layer', 'row', 'column']
elif self.grid_type() == DiscretizationType.DISV:
return ['layer', 'layer_cell_num']
elif self.grid_type() == DiscretizationType.DISU:
elif self.grid_type() == DiscretizationType.DISU or \
self.grid_type() == DiscretizationType.DISL:
return ['node']

def get_num_spatial_coordinates(self):
if self.grid_type() == DiscretizationType.DIS:
return 3
elif self.grid_type() == DiscretizationType.DISV:
return 2
elif self.grid_type() == DiscretizationType.DISU:
elif self.grid_type() == DiscretizationType.DISU or \
self.grid_type() == DiscretizationType.DISL:
return 1

def change_grid_spacing(self, spacing_factor):
self.test = 1

def change_discretization_type(self, new_dis_type):
self.test = 1

def num_rows(self):
if self.grid_type() != DiscretizationType.DIS:
except_str = 'ERROR: Model "{}" does not have rows. Can not ' \
Expand Down Expand Up @@ -567,7 +574,8 @@ def num_layers(self):
elif self.grid_type() == DiscretizationType.DISV:
return self._simulation_data.mfdata[
(self._model_name, 'disv', 'dimensions', 'nlay')].get_data()
elif self.grid_type() == DiscretizationType.DISU:
elif self.grid_type() == DiscretizationType.DISU or \
self.grid_type() == DiscretizationType.DISL:
return None

def num_cells(self):
Expand All @@ -578,6 +586,9 @@ def num_cells(self):
elif self.grid_type() == DiscretizationType.DISU:
return self._simulation_data.mfdata[
(self._model_name, 'disu', 'dimensions', 'nodes')].get_data()
elif self.grid_type() == DiscretizationType.DISL:
return self._simulation_data.mfdata[
(self._model_name, 'disl', 'dimensions', 'nodes')].get_data()

def get_all_model_cells(self):
model_cells = []
Expand All @@ -592,7 +603,8 @@ def get_all_model_cells(self):
for layer_cellid in range(0, self.num_rows()):
model_cells.append((layer + 1, layer_cellid + 1))
return model_cells
elif self.grid_type() == DiscretizationType.DISU:
elif self.grid_type() == DiscretizationType.DISU or \
self.grid_type() == DiscretizationType.DISL:
for node in range(0, self.num_cells()):
model_cells.append(node + 1)
return model_cells
Expand Down
3 changes: 2 additions & 1 deletion flopy/mf6/data/mfdatastorage.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1535,7 +1535,8 @@ def _verify_list(self, data):
model_grid = self.data_dimensions.get_model_grid()
cellid_size = model_grid.\
get_num_spatial_coordinates()
if len(data_line[index]) != cellid_size:
if cellid_size != 1 and \
len(data_line[index]) != cellid_size:
message = 'Cellid "{}" contains {} integer(s). ' \
'Expected a cellid containing {} ' \
'integer(s) for grid type' \
Expand Down
6 changes: 4 additions & 2 deletions flopy/mf6/data/mfstructure.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -65,8 +65,9 @@ def __init__(self):
# FIX: Transport - multi packages are hard coded
self.multi_package = {'gwfmvr': 0, 'exggwfgwf': 0, 'gwfchd': 0,
'gwfrch': 0, 'gwfwel': 0,
'gwfdrn': 0, 'gwfriv': 0, 'utlobs': 0,
'utlts': 0, 'utltas': 0}
'gwfdrn': 0, 'gwfriv': 0, 'lnfcgeo': 0,
'lnfrgeo': 0, 'lnfngeo':0,
'utlobs': 0, 'utlts': 0, 'utltas': 0}

def get_file_list(self):
file_order = ['sim-nam', # dfn completed tex updated
Expand All @@ -77,6 +78,7 @@ def get_file_list(self):
'gwf-dis', # dfn completed tex updated
'gwf-disv', # dfn completed tex updated
'gwf-disu', # dfn completed tex updated
'lnf-disl', # dfn completed tex updated
'gwf-ic', # dfn completed tex updated
'gwf-npf', # dfn completed tex updated
'gwf-sto', # dfn completed tex updated
Expand Down
31 changes: 31 additions & 0 deletions flopy/mf6/mfmodel.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -368,6 +368,33 @@ def modelgrid(self):
epsg=self._modelgrid.epsg, xoff=self._modelgrid.xoffset,
yoff=self._modelgrid.yoffset, angrot=self._modelgrid.angrot,
nodes=dis.nodes.get_data())
elif self.get_grid_type() == DiscretizationType.DISL:
dis = self.get_package('disl')
if not hasattr(dis, '_init_complete'):
if not hasattr(dis, 'cell1d'):
# disv package has not yet been initialized
return self._modelgrid
else:
# disv package has been partially initialized
self._modelgrid = VertexGrid(vertices=dis.vertices.array,
cell1d=dis.cell1d.array,
top=None,
botm=None,
idomain=None,
lenuni=None,
proj4=self._modelgrid.proj4,
epsg=self._modelgrid.epsg,
xoff=self._modelgrid.xoffset,
yoff=self._modelgrid.yoffset,
angrot=self._modelgrid.angrot)
else:
self._modelgrid = VertexGrid(
vertices=dis.vertices.array, cell1d=dis.cell1d.array,
top=dis.top.array, botm=dis.botm.array,
idomain=dis.idomain.array, lenuni=dis.length_units.array,
proj4=self._modelgrid.proj4, epsg=self._modelgrid.epsg,
xoff=self._modelgrid.xoffset, yoff=self._modelgrid.yoffset,
angrot=self._modelgrid.angrot)
else:
return self._modelgrid

Expand Down Expand Up @@ -671,6 +698,10 @@ def get_grid_type(self):
'disu{}'.format(structure.get_version_string()),
0) is not None:
return DiscretizationType.DISU
elif package_recarray.search_data(
'disl{}'.format(structure.get_version_string()),
0) is not None:
return DiscretizationType.DISL

return DiscretizationType.UNDEFINED

Expand Down
1 change: 1 addition & 0 deletions flopy/mf6/utils/mfenums.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,3 +9,4 @@ class DiscretizationType(Enum):
DIS = 1
DISV = 2
DISU = 3
DISL = 4

0 comments on commit 469727b

Please sign in to comment.