Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

fix(mp7particledata): add global_xy to to_prp #2405

Merged
merged 3 commits into from
Dec 23, 2024
Merged

fix(mp7particledata): add global_xy to to_prp #2405

Show file tree
Hide file tree
Changes from 2 commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
0001d1f
fix(mp7particledata): add global_xy to `to_prp`
dbrakenhoff Dec 23, 2024
a0cfd65
fix(mp7particledata): add global_xy to `to_prp`
dbrakenhoff Dec 23, 2024
b9ce6b0
use _has_ref_coordinates
dbrakenhoff Dec 23, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 3 additions & 1 deletion autotest/test_grid_cases.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,7 +10,7 @@

class GridCases:
@staticmethod
def structured_small():
def structured_small(xoff=0.0, yoff=0.0):
nlay, nrow, ncol = 3, 2, 3
delc = 1.0 * np.ones(nrow, dtype=float)
delr = 1.0 * np.ones(ncol, dtype=float)
Expand All @@ -29,6 +29,8 @@ def structured_small():
top=top,
botm=botm,
idomain=idomain,
xoff=xoff,
yoff=yoff,
)

@staticmethod
Expand Down
27 changes: 27 additions & 0 deletions autotest/test_particledata.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -245,6 +245,33 @@ def test_particledata_to_prp_dis_1():
assert np.isclose(rpt[6], minz + (exp[ci][5] * (maxz - minz)))


def test_particledata_to_prp_dis_1_global_xy():
# model grid
xoff = 100.0
yoff = 300.0
grid = GridCases().structured_small(xoff=xoff, yoff=yoff)

# particle data
cells = [(0, 1, 1), (0, 1, 2)]
part_data = ParticleData(partlocs=cells, structured=True)

# convert to global coordinates
rpts_prt_model_coords = flatten(list(part_data.to_prp(grid)))
rpts_prt_global_coords = flatten(list(part_data.to_prp(grid, global_xy=True)))

# check global and model coords
# x
assert np.all(
np.array(rpts_prt_global_coords)[:, -3] - xoff
== np.array(rpts_prt_model_coords)[:, -3]
)
# y
assert np.all(
np.array(rpts_prt_global_coords)[:, -2] - yoff
== np.array(rpts_prt_model_coords)[:, -2]
)


def test_particledata_to_prp_dis_9():
# minimal structured grid
grid = GridCases().structured_small()
Expand Down
65 changes: 49 additions & 16 deletions flopy/modpath/mp7particledata.py
View file
Open in desktop
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
dbrakenhoff marked this conversation as resolved.
Show resolved Hide resolved
):
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,
wpbonelli marked this conversation as resolved.
Show resolved Hide resolved
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
Loading