bug: ParticleData.to_prp() returns global coordinates

[dbrakenhoff]

Describe the bug
Not sure if this can yet be qualified as bug, or is simply still under active development and I'm a bit early in posting this, but I noticed that ParticleData.to_prp() returns global coordinates. PRT6 currently expects model coordinates for the input. This isn't caught in the tests as as far as I could see, as no PRT6 models are run using a grid with offsets (in both flopy and modflow6 repos).

I saw other options for specifying model coordinates are under review in Modflow6, which might also solve this issue by allowing computation with global coordinates: MODFLOW-USGS/modflow6#1901. Nevertheless, it might be nice to have both global and model coordinate options in to_prp().

EDIT: just noticed issue #2388. If vertices are specified in model coordinates, that would also solve this issue.

To Reproduce

import flopy as fp
import numpy as np

gr = fp.discretization.StructuredGrid(
    10 * np.ones(10),
    10 * np.ones(10),
    top=np.ones((10, 10)),
    botm=np.zeros((1, 10, 10)),
    nlay=1,
    xoff=100,
    yoff=300,
)
pd = fp.modpath.ParticleData([(0, 0, 0), (0, 5, 5), (0, 9, 9)], structured=True)
list(pd.to_prp(gr))

Expected behavior
Perhaps the option of returning global or model coordinates could be added to to_prp()? Something like the code below could do the trick. I used the term model_coords since in modpath local signifies a location within a cell.

If this is something you would want to support, Id be happy to submit a PR with this small change, let me know!

# class ParticleData:
# ...
    def to_coords(self, grid, localz=False, model_coords=True) -> Iterator[tuple]:  # <-- add model_coords=True
        # ...
        def cvt_xy(p, vs):
            mn, mx = min(vs), max(vs)
            span = mx - mn
            return mn + span * p

        if grid.grid_type == "structured":
            if not hasattr(self.particledata, "k"):
                raise ValueError(
                    "Particle representation is not structured but grid is"
                )

            def cvt_z(p, k, i, j):
                mn, mx = (
                    grid.botm[k, i, j],
                    grid.top[i, j] if k == 0 else grid.botm[k - 1, i, j],
                )
                span = mx - mn
                return mn + span * p

            def convert(row, model_coords=True) -> tuple[float, float, float]:  # <-- add model_coords
                verts = grid.get_cell_vertices(row.i, row.j)
                if model_coords:
                    xs, ys = grid.get_local_coords(*np.array(verts).T)  # <--convert global coords to model coords
                else:
                    xs, ys = list(zip(*verts))
                return [
                    cvt_xy(row.localx, xs),
                    cvt_xy(row.localy, ys),
                    row.localz if localz else cvt_z(row.localz, row.k, row.i, row.j),
                ]

        else:
            if hasattr(self.particledata, "k"):
                raise ValueError(
                    "Particle representation is structured but grid is not"
                )

            def cvt_z(p, nn):
                k, j = grid.get_lni([nn])[0]
                mn, mx = (
                    grid.botm[k, j],
                    grid.top[j] if k == 0 else grid.botm[k - 1, j],
                )
                span = mx - mn
                return mn + span * p

            def convert(row, model_coords=True) -> tuple[float, float, float]:  # <-- add model_coords
                verts = grid.get_cell_vertices(row.node)
                if model_coords:
                    xs, ys = grid.get_local_coords(*np.array(verts).T)  # <--convert global coords to model coords
                else:
                    xs, ys = list(zip(*verts))
                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, model_coords=model_coords)

    def to_prp(self, grid, localz=False, model_coords=True) -> Iterator[tuple]:  # <-- add model_coords
        # ...
        for i, (t, c) in enumerate(
            zip(
                self.particledata.itertuples(index=False),
                self.to_coords(grid, localz, model_coords=model_coords),  # <-- pass model_coords
            )
        ):
            row = [i]  # release point index (irpt)
            if "node" in self.particledata:
                k, j = grid.get_lni([t.node])[0]
                row.extend([(k, j)])
            else:
                row.extend([t.k, t.i, t.j])
            row.extend(c)
            yield tuple(row)

[wpbonelli]

@dbrakenhoff thanks a lot for catching this — this should probably convert to model coords by default, since that's PRT's expectation... so maybe we can invert the flag? maybe global_xy to match MODFLOW-USGS/modflow6#1901? since the z coordinate is considered separately. please do feel free to submit!

[wpbonelli]

this is all early stage and actively developed but feedback and/or contributions always extremely welcome.