Spatially-aware package/result access (query by geometry, not index)

[mjreno]

From user feedback: a recurring friction point in GIS-driven MF6 workflows is that relating model inputs/outputs back to real-world features requires knowing internal indices — reach numbers, UZF info numbers, grid cell numbers — rather than being able to query by geometry or feature ID. This makes it harder to change grids or package specs without rebuilding that
indexing by hand. Worth exploring what a "spatially aware" FloPy 4 API could look like — e.g. package data addressable by feature/geometry, consistent CRS propagation throughout. Opened to track discussion; feedback and related issues welcome.

[wpbonelli]

We've put a fair bit of thought into this, but it's not all cross-referenced. This seems as good a place as any.

Essentially, we want to be able to go smoothly between

• spatial coordinates, model or global
• grid cell (or other relevant grid feature) indices
• boundary condition features, by index (ifno etc) or name

Given any of the above, the others corresponding to it should be easily discoverable.

There was some (early) discussion in #63. #120 got the ball rolling, modulo xattree specifics which will go away.

#278 sketches a design for an xarray extension, the main entry point for which would be a .grid accessor on xarray data arrays/sets. Basically what xugrid does but scoped specifically for MF6 rather than UGRID. Needs some more attention, I let this drop after making some initial progress. #207 and #208 describe the core pieces, scoped to structured grids for now, though we'll want to make them generic. We already have a stupid simple custom xarray Index for structured grids in flopy4/mf6/indexes.py . An accessor could be built on something like this. #324 lays some groundwork. #220 just makes explicit the idea that this MF6-scoped extension can live side-by-side with xugrid, so you can choose .grid or .ugrid depending which representation you want.

#160 will add spatial awareness for package data in 2 dimensions. User can choose between geopandas/xarray as needed, they each have their pros and cons.

#199 is a typical user story: you have a grid cell index, you want to know geometry/hydraulic/BC data in that cell.

#345 and #346 are loosely related to all this insofar as features (e.g. connections, cross section points) must be addressable in some way in order to correlate them with anything else.

[lmioduszewski]

Thanks for pulling this together, @mjreno and @wpbonelli . The coordinate <--> cell/feature index <--> BC-feature (ifno/name) triangle is the right framing, and it's what lets you swap a grid without hand-rebuilding index bookkeeping.

The piece I'd push hardest on is unstructured grids, and the BC-feature layer specifically. Redirecting .grid to .ugrid for DISV covers array topology, but SFR reaches / LAK lakes / UZF cells / connections need their own stable feature-name/id (cells map that UGRID doesn't model) ideally keyed off source geometry so it survives a grid rebuild.

Roughly the shape I'd love to be able to write:

# features are defined from GIS geometry, named from a source field (not grid
# indices), so 'main_stem' / 'north_pond' survive a grid rebuild. Each call
# returns a package spec; the model is assembled from a flat list.
# package specs are standalone and independent of a specific model context or grid

sfr = Sfr(
    streams="streams.gpkg", name_field="name",
    width=18.0, roughness=0.030,                          # scalar --> all reaches
    streambed_k=0.05, streambed_thickness=1.5,
    gradient={"main_stem": 0.0012, "north_trib": 0.004},  # per-reach, by name
    connections=[("north_trib", "main_stem"),             # network wired by name
                 ("south_trib", "main_stem")],
    inflow={0: {"north_trib": 15_000.0}},                 # per period, by name
)

lak = Lak(
    lakes="lakes.gpkg", name_field="name",
    bed_leakance=0.11,                                    # scalar --> all lakes
    starting_stage={"north_pond": 101.0},                 # per-lake, by name
    lake_bottom={"north_pond": 96.0},
)

# movers reference features by name:
mvr = Mvr(moves=[("sfr", "main_stem", "lak", "north_pond")])

gwf = Gwf("valley", grid=grid, packages=[sfr, lak, mvr])

# after building, the same features are addressable by name or geometry:
node  = gwf.grid.locate(x=1250.0, y=980.0)   # coordinate --> (layer, node)
reach = gwf.sfr.nearest((1250.0, 980.0))     # geometry   --> reach (+ ifno, cell)
lake  = gwf.lak["north_pond"]                # name       --> cells, ifno, geometry

[wpbonelli]

Thanks for this @lmioduszewski. Super helpful for us!