Skip to content

Commit

Permalink
fix(mp7particledata): add global_xy to to_prp
Browse files Browse the repository at this point in the history
- default is False which converts vertices to model coordinates
- add simple test
  • Loading branch information
dbrakenhoff committed Dec 23, 2024
1 parent 0001d1f commit a0cfd65
Showing 1 changed file with 49 additions and 16 deletions.
65 changes: 49 additions & 16 deletions flopy/modpath/mp7particledata.py
Original file line number Diff line number Diff line change
Expand Up @@ -362,7 +362,7 @@ def write(self, f=None):
for v in d:
f.write(fmt.format(*v))

def to_coords(self, grid, localz=False) -> Iterator[tuple]:
def to_coords(self, grid, localz=False, global_xy=False) -> Iterator[tuple]:
"""
Compute particle coordinates on the given grid.
Expand Down Expand Up @@ -397,9 +397,12 @@ def cvt_z(p, k, i, j):
span = mx - mn
return mn + span * p

def convert(row) -> tuple[float, float, float]:
def convert(row, global_xy=False) -> tuple[float, float, float]:
verts = grid.get_cell_vertices(row.i, row.j)
xs, ys = list(zip(*verts))
if global_xy:
xs, ys = list(zip(*verts))
else:
xs, ys = grid.get_local_coords(*np.array(verts).T)
return [
cvt_xy(row.localx, xs),
cvt_xy(row.localy, ys),
Expand All @@ -421,19 +424,22 @@ def cvt_z(p, nn):
span = mx - mn
return mn + span * p

def convert(row) -> tuple[float, float, float]:
def convert(row, global_xy=False) -> tuple[float, float, float]:
verts = grid.get_cell_vertices(row.node)
xs, ys = list(zip(*verts))
if global_xy:
xs, ys = list(zip(*verts))
else:
xs, ys = grid.get_local_coords(*np.array(verts).T)
return [
cvt_xy(row.localx, xs),
cvt_xy(row.localy, ys),
row.localz if localz else cvt_z(row.localz, row.node),
]

for t in self.particledata.itertuples():
yield convert(t)
yield convert(t, global_xy=global_xy)

def to_prp(self, grid, localz=False) -> Iterator[tuple]:
def to_prp(self, grid, localz=False, global_xy=False) -> Iterator[tuple]:
"""
Convert particle data to PRT particle release point (PRP)
package data entries for the given grid. A model grid is
Expand All @@ -447,6 +453,8 @@ def to_prp(self, grid, localz=False) -> Iterator[tuple]:
The grid on which to locate particle release points.
localz : bool, optional
Whether to return local z coordinates.
global_xy : bool, optional
Whether to return global x and y coordinates, default is False.
Returns
-------
Expand All @@ -459,7 +467,7 @@ def to_prp(self, grid, localz=False) -> Iterator[tuple]:
for i, (t, c) in enumerate(
zip(
self.particledata.itertuples(index=False),
self.to_coords(grid, localz),
self.to_coords(grid, localz, global_xy=global_xy),
)
):
row = [i] # release point index (irpt)
Expand Down Expand Up @@ -778,10 +786,16 @@ def write(self, f=None):
)


def get_extent(grid, k=None, i=None, j=None, nn=None, localz=False) -> Extent:
def get_extent(
grid, k=None, i=None, j=None, nn=None, localz=False, global_xy=False
) -> Extent:
# get cell coords and span in each dimension
if not (k is None or i is None or j is None):
verts = grid.get_cell_vertices(i, j)
if not global_xy and (
grid.xoffset != 0.0 or grid.yoffset != 0.0 or grid.angrot != 0.0
):
verts = list(zip(*grid.get_local_coords(*np.array(verts).T)))
minz, maxz = (
(0, 1)
if localz
Expand All @@ -792,6 +806,10 @@ def get_extent(grid, k=None, i=None, j=None, nn=None, localz=False) -> Extent:
)
elif nn is not None:
verts = grid.get_cell_vertices(nn)
if not global_xy and (
grid.xoffset != 0.0 or grid.yoffset != 0.0 or grid.angrot != 0.0
):
verts = list(zip(*grid.get_local_coords(*np.array(verts).T)))
if grid.grid_type == "structured":
k, i, j = grid.get_lrc([nn])[0]
minz, maxz = (
Expand Down Expand Up @@ -967,7 +985,14 @@ def get_cell_release_points(subdivisiondata, cellid, extent) -> Iterator[tuple]:


def get_release_points(
subdivisiondata, grid, k=None, i=None, j=None, nn=None, localz=False
subdivisiondata,
grid,
k=None,
i=None,
j=None,
nn=None,
localz=False,
global_xy=False,
) -> Iterator[tuple]:
"""
Get MODPATH 7 release point tuples for the given cell.
Expand All @@ -980,7 +1005,7 @@ def get_release_points(
)

cellid = [k, i, j] if nn is None else [nn]
extent = get_extent(grid, k, i, j, nn, localz)
extent = get_extent(grid, k, i, j, nn, localz, global_xy=global_xy)

if isinstance(subdivisiondata, FaceDataType):
return get_face_release_points(subdivisiondata, cellid, extent)
Expand Down Expand Up @@ -1351,7 +1376,7 @@ def write(self, f=None):
line += "\n"
f.write(line)

def to_coords(self, grid, localz=False) -> Iterator[tuple]:
def to_coords(self, grid, localz=False, global_xy=False) -> Iterator[tuple]:
"""
Compute global particle coordinates on the given grid.
Expand All @@ -1361,6 +1386,8 @@ def to_coords(self, grid, localz=False) -> Iterator[tuple]:
The grid on which to locate particle release points.
localz : bool, optional
Whether to return local z coordinates.
global_xy : bool, optional
Whether to return global x, y coordinates. Default is False.
Returns
-------
Expand All @@ -1369,23 +1396,27 @@ def to_coords(self, grid, localz=False) -> Iterator[tuple]:

for sd in self.subdivisiondata:
for nd in self.nodedata:
for rpt in get_release_points(sd, grid, nn=int(nd[0]), localz=localz):
for rpt in get_release_points(
sd, grid, nn=int(nd[0]), localz=localz, global_xy=global_xy
):
yield (*rpt[1:4],)

def to_prp(self, grid, localz=False) -> Iterator[tuple]:
def to_prp(self, grid, localz=False, global_xy=False) -> Iterator[tuple]:
"""
Convert particle data to PRT particle release point (PRP)
package data entries for the given grid. A model grid is
required because MODPATH supports several ways to specify
particle release locations by cell ID and subdivision info
or local coordinates, but PRT expects global coordinates.
or local coordinates, but PRT expects model coordinates, by default.
Parameters
----------
grid : flopy.discretization.grid.Grid
The grid on which to locate particle release points.
localz : bool, optional
Whether to return local z coordinates.
global_xy : bool, optional
Whether to return global x, y coordinates. Default is False.
Returns
-------
Expand All @@ -1396,7 +1427,9 @@ def to_prp(self, grid, localz=False) -> Iterator[tuple]:
for sd in self.subdivisiondata:
for nd in self.nodedata:
for irpt, rpt in enumerate(
get_release_points(sd, grid, nn=int(nd[0]), localz=localz)
get_release_points(
sd, grid, nn=int(nd[0]), localz=localz, global_xy=global_xy
)
):
row = [irpt]
if grid.grid_type == "structured":
Expand Down

0 comments on commit a0cfd65

Please sign in to comment.