From 81c1a1f8f8a3196d5a930dd8f431f08acf17a472 Mon Sep 17 00:00:00 2001 From: Onno Ebbens Date: Mon, 9 Feb 2026 13:35:47 +0100 Subject: [PATCH 1/3] update to flopy 3.10.0 --- nlmod/dims/grid.py | 61 +++++++++++++++++--------------------- nlmod/dims/layers.py | 2 +- nlmod/gwf/surface_water.py | 4 +-- nlmod/mfoutput/mfoutput.py | 4 +-- nlmod/plot/dcs.py | 8 ++--- pyproject.toml | 2 +- 6 files changed, 38 insertions(+), 43 deletions(-) diff --git a/nlmod/dims/grid.py b/nlmod/dims/grid.py index 2fa5c4e8..f13fa543 100644 --- a/nlmod/dims/grid.py +++ b/nlmod/dims/grid.py @@ -157,7 +157,7 @@ def get_icell2d_from_xy(x, y, ds, gi=None, rotated=True): assert ds.gridtype == "vertex", msg if gi is None: gi = flopy.utils.GridIntersect( - modelgrid_from_ds(ds, rotated=rotated), method="vertex" + modelgrid_from_ds(ds, rotated=rotated) ) cellids = gi.intersects(Point(x, y))["cellids"] if len(cellids) < 1: @@ -622,7 +622,7 @@ def get_cellids_from_xy(x, y, ds=None, modelgrid=None, gi=None): # build geometries and cellids for grid if gi is None: - gi = flopy.utils.GridIntersect(modelgrid, method="vertex") + gi = flopy.utils.GridIntersect(modelgrid) # spatial join points with grid and add resulting cellid to obs_ds spatial_join = gpd.GeoDataFrame(geometry=gpd.points_from_xy(x, y)).sjoin( @@ -1345,18 +1345,16 @@ def polygon_to_area(modelgrid, polygon, da, gridtype="structured"): f'input geometry should by of type "Polygon" not {polygon.geom_type}' ) - ix = GridIntersect(modelgrid, method="vertex") - opp_cells = ix.intersect(polygon) + ix = GridIntersect(modelgrid) + df = ix.intersect(polygon, geo_dataframe=True) if gridtype == "structured": area_array = util.get_da_from_da_ds(da, dims=("y", "x"), data=0) - for cellid, area in zip(opp_cells["cellids"], opp_cells["areas"]): - area_array[cellid[0], cellid[1]] = area + for row, col, area in zip(df["row"], df["col"], df["areas"]): + area_array[row, col] = area elif gridtype == "vertex": area_array = util.get_da_from_da_ds(da, dims=("icell2d",), data=0) - cids = opp_cells.cellids - area = opp_cells.areas - area_array[cids.astype(int)] = area + area_array[df['cellid']] = df['areas'] return area_array @@ -1770,7 +1768,7 @@ def gdf_to_bool_da( # Rotate the polygon instead of the modelgrid if ix is None: modelgrid = modelgrid_from_ds(ds, rotated=False) - ix = GridIntersect(modelgrid, method="vertex") + ix = GridIntersect(modelgrid) grid_rotation = ds.attrs.get("angrot", 0.0) ix_rotation = ix.mfgrid.angrot @@ -1780,20 +1778,20 @@ def gdf_to_bool_da( multipolygon = affine_transform(multipolygon, affine) if kwargs or contains_centroid or min_area_fraction is not None: - r = ix.intersect( + df = ix.intersect( multipolygon, contains_centroid=contains_centroid, min_area_fraction=min_area_fraction, + geo_dataframe=True **kwargs, ) else: - r = ix.intersects(multipolygon) + df = ix.intersects(multipolygon, dataframe=True) - if r.size > 0 and ds.gridtype == "structured": - ixs, iys = zip(*r["cellids"], strict=True) - da.values[ixs, iys] = True - elif r.size > 0 and ds.gridtype == "vertex": - da[r["cellids"].astype(int)] = True + if df.shape[0] > 0 and ds.gridtype == "structured": + da.values[df['row'], df['col']] = True + elif df.shape[0] > 0 and ds.gridtype == "vertex": + da.values[df["cellid"]] = True return da @@ -1893,7 +1891,7 @@ def gdf_to_count_da(gdf, ds, ix=None, buffer=0.0, **kwargs): # build list of gridcells if ix is None: modelgrid = modelgrid_from_ds(ds, rotated=False) - ix = GridIntersect(modelgrid, method="vertex") + ix = GridIntersect(modelgrid) if ds.gridtype == "structured": da = util.get_da_from_da_ds(ds, dims=("y", "x"), data=0) @@ -1962,7 +1960,6 @@ def gdf_to_count_ds(gdf, ds, da_name, keep_coords=None, ix=None, buffer=0.0, **k def gdf_to_grid( gdf, ml=None, - method="vertex", ix=None, desc="Intersecting with grid", silent=False, @@ -1981,8 +1978,6 @@ def gdf_to_grid( The flopy model or xarray dataset that defines the grid. When a Dataset is supplied, and the grid is rotated, the geodataframe is transformed in model coordinates. The default is None. - method : string, optional - Method passed to the GridIntersect-class. The default is 'vertex'. ix : flopy.utils.GridIntersect, optional GridIntersect, if not provided the modelgrid in ml is used. desc : string, optional @@ -2019,19 +2014,19 @@ def gdf_to_grid( ) if ix is None: - ix = flopy.utils.GridIntersect(modelgrid, method=method) + ix = flopy.utils.GridIntersect(modelgrid) shps = [] geometry = gdf.geometry.name for _, shp in tqdm(gdf.iterrows(), total=gdf.shape[0], desc=desc, disable=silent): - r = ix.intersect(shp[geometry], **kwargs) - for i in range(r.shape[0]): + df = ix.intersect(shp[geometry], geo_dataframe=True, **kwargs) + for i in range(df.shape[0]): shpn = shp.copy() - shpn["cellid"] = r["cellids"][i] - shpn[geometry] = r["ixshapes"][i] + shpn["cellid"] = df["cellids"][i] + shpn[geometry] = df["geometry"][i] if shp[geometry].geom_type == ["LineString", "MultiLineString"]: - shpn["length"] = r["lengths"][i] + shpn["length"] = df["lengths"][i] elif shp[geometry].geom_type in ["Polygon", "MultiPolygon"]: - shpn["area"] = r["areas"][i] + shpn["area"] = df["areas"][i] shps.append(shpn) if len(shps) == 0: @@ -2115,7 +2110,7 @@ def gdf_area_to_da( affine = get_affine_world_to_mod(ds) gdf = affine_transform_gdf(gdf, affine) if ix is None: - ix = flopy.utils.GridIntersect(modelgrid, method="vertex") + ix = flopy.utils.GridIntersect(modelgrid) if index_name is None: index_name = gdf.index.name @@ -2135,12 +2130,12 @@ def gdf_area_to_da( for irow, index in tqdm( enumerate(gdf.index), total=gdf.shape[0], desc=desc, disable=silent ): - r = ix.intersect(gdf.at[index, geometry], **kwargs) + df = ix.intersect(gdf.at[index, geometry], geo_dataframe=True, **kwargs) if structured: - for i in range(r.shape[0]): - data[r["cellids"][i] + (irow,)] += r["areas"][i] + for i in range(df.shape[0]): + data[df["cellids"][i] + (irow,)] += df["areas"][i] else: - data[list(r["cellids"]), irow] = r["areas"] + data[list(df["cellids"]), irow] = df["areas"] if sparse: try: diff --git a/nlmod/dims/layers.py b/nlmod/dims/layers.py index b76c6fa9..7a43df6a 100644 --- a/nlmod/dims/layers.py +++ b/nlmod/dims/layers.py @@ -2151,7 +2151,7 @@ def get_modellayers_screens(ds, screen_top, screen_bottom, xy=None, icell2d=None """ if grid.is_vertex(ds): if icell2d is None: - gi = flopy.utils.GridIntersect(grid.modelgrid_from_ds(ds), method="vertex") + gi = flopy.utils.GridIntersect(grid.modelgrid_from_ds(ds)) icell2d = [grid.get_icell2d_from_xy(x, y, ds, gi=gi) for x, y in xy] # make dataset of observations ds_obs = ds.sel(icell2d=icell2d) diff --git a/nlmod/gwf/surface_water.py b/nlmod/gwf/surface_water.py index bdd452b6..75435ef0 100644 --- a/nlmod/gwf/surface_water.py +++ b/nlmod/gwf/surface_water.py @@ -1227,8 +1227,8 @@ def get_seaonal_timeseries( def rivdata_from_xylist(gwf, xylist, layer, stage, cond, rbot, aux=None): - gi = flopy.utils.GridIntersect(gwf.modelgrid, method="vertex") - cellids = gi.intersect(xylist, shapetype="linestring")["cellids"] + gi = flopy.utils.GridIntersect(gwf.modelgrid) + cellids = gi.intersect(xylist, shapetype="linestring", geo_dataframe=True)["cellids"] riv_data = [] for cid in cellids: if isinstance(cid, (list, tuple)) and len(cid) == 2: diff --git a/nlmod/mfoutput/mfoutput.py b/nlmod/mfoutput/mfoutput.py index 59102350..3049d67a 100644 --- a/nlmod/mfoutput/mfoutput.py +++ b/nlmod/mfoutput/mfoutput.py @@ -162,7 +162,7 @@ def _get_heads_da( flopy HeadFile object for binary heads modelgrid : flopy.discretization.Grid, optional flopy modelgrid object, default is None, in which case the modelgrid - is derived from `hobj.mg` + is derived from `hobj.modelgrid` Returns ------- @@ -175,7 +175,7 @@ def _get_heads_da( kstpkper = hobj.get_kstpkper() if modelgrid is None: - modelgrid = hobj.mg + modelgrid = hobj.modelgrid # shape is derived from hobj, not modelgrid as array read from # binary file always has 3 dimensions shape = (hobj.nlay, hobj.nrow, hobj.ncol) diff --git a/nlmod/plot/dcs.py b/nlmod/plot/dcs.py index 5f65614b..153d782b 100644 --- a/nlmod/plot/dcs.py +++ b/nlmod/plot/dcs.py @@ -66,11 +66,11 @@ def __init__( if self.icell2d in ds.dims: # determine the cells that are crossed modelgrid = modelgrid_from_ds(ds, rotated=False) - gi = flopy.utils.GridIntersect(modelgrid, method="vertex") - r = gi.intersect(line) + gi = flopy.utils.GridIntersect(modelgrid) + df = gi.intersect(line, geo_dataframe=True) s_cell = [] - for i, ic2d in enumerate(r["cellids"]): - intersection = r["ixshapes"][i] + for i, ic2d in enumerate(df["cellid"]): + intersection = df["geometry"][i] if intersection.length == 0: continue if isinstance(intersection, MultiLineString): diff --git a/pyproject.toml b/pyproject.toml index 0ddb2be7..9291b485 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -71,7 +71,7 @@ full = [ ] knmi = ["h5netcdf", "h5py", "nlmod[grib]"] grib = ["cfgrib", "ecmwflibs"] -test = ["pytest>=7", "pytest-cov", "pytest-dependency"] +test = ["pytest>=7", "pytest-cov", "lxml"] nbtest = ["nbformat", "nbconvert>6.4.5"] lint = ["ruff"] ci = ["nlmod[full,lint,test,nbtest]"] From 118ccf28e70398b0286b314b80735fce84d349a0 Mon Sep 17 00:00:00 2001 From: Onno Ebbens Date: Tue, 10 Feb 2026 09:35:12 +0100 Subject: [PATCH 2/3] remove cellids --- nlmod/dims/grid.py | 21 ++++++++++----------- nlmod/gwf/surface_water.py | 13 ++++++++----- 2 files changed, 18 insertions(+), 16 deletions(-) diff --git a/nlmod/dims/grid.py b/nlmod/dims/grid.py index f13fa543..631815e9 100644 --- a/nlmod/dims/grid.py +++ b/nlmod/dims/grid.py @@ -159,7 +159,7 @@ def get_icell2d_from_xy(x, y, ds, gi=None, rotated=True): gi = flopy.utils.GridIntersect( modelgrid_from_ds(ds, rotated=rotated) ) - cellids = gi.intersects(Point(x, y))["cellids"] + cellids = gi.intersects(Point(x, y), dataframe=True)["cellid"].values if len(cellids) < 1: raise (ValueError(f"Point ({x}, {y}) is outside of the model grid")) icell2d = cellids[0] @@ -227,7 +227,7 @@ def get_row_col_from_xy(x, y, ds, rotated=True, gi=None): msg = "get_row_col_from_xy can only be applied to a structured grid" assert ds.gridtype == "structured", msg if gi is not None: - cellids = gi.intersects(Point(x, y))["cellids"] + cellids = gi.intersects(Point(x, y), dataframe=True)[["row", "col"]].values if len(cellids) < 1: raise (ValueError(f"Point ({x}, {y}) is outside of the model grid")) row, col = cellids[0] @@ -1907,18 +1907,17 @@ def gdf_to_count_da(gdf, ds, ix=None, buffer=0.0, **kwargs): for geom in geoms: if buffer > 0.0: - cids = ix.intersects(geom.buffer(buffer), **kwargs)["cellids"] + df = ix.intersects(geom.buffer(buffer), dataframe=True, **kwargs) else: - cids = ix.intersects(geom, **kwargs)["cellids"] + df = ix.intersects(geom, dataframe=True, **kwargs) - if len(cids) == 0: + if len(df) == 0: continue if ds.gridtype == "structured": - ixs, iys = zip(*cids) - da.values[ixs, iys] += 1 + da.values[df['row'].values, df['col'].values] += 1 elif ds.gridtype == "vertex": - da[cids.astype(int)] += 1 + da[df['cellid'].values] += 1 return da @@ -2132,10 +2131,10 @@ def gdf_area_to_da( ): df = ix.intersect(gdf.at[index, geometry], geo_dataframe=True, **kwargs) if structured: - for i in range(df.shape[0]): - data[df["cellids"][i] + (irow,)] += df["areas"][i] + for row, col, area in zip(df['row'], df['col'], df['areas']): + data[(row, col, irow)] += area else: - data[list(df["cellids"]), irow] = df["areas"] + data[list(df["cellid"]), irow] = df["areas"] if sparse: try: diff --git a/nlmod/gwf/surface_water.py b/nlmod/gwf/surface_water.py index 75435ef0..ade38d71 100644 --- a/nlmod/gwf/surface_water.py +++ b/nlmod/gwf/surface_water.py @@ -1228,17 +1228,20 @@ def get_seaonal_timeseries( def rivdata_from_xylist(gwf, xylist, layer, stage, cond, rbot, aux=None): gi = flopy.utils.GridIntersect(gwf.modelgrid) - cellids = gi.intersect(xylist, shapetype="linestring", geo_dataframe=True)["cellids"] + df = gi.intersect(xylist, shapetype="linestring", geo_dataframe=True) riv_data = [] - for cid in cellids: - if isinstance(cid, (list, tuple)) and len(cid) == 2: - idata = [(layer, cid[0], cid[1]), stage, cond, rbot] + + if 'row' in df.columns: # structured grid + for row, col in zip(df['row'], df['col']): + idata = [(layer, row, col), stage, cond, rbot] if aux is not None: idata.append(aux) riv_data.append(idata) - else: + else: # unstructured grid + for cid in df['cellid']: idata = [(layer, cid), stage, cond, rbot] if aux is not None: idata.append(aux) riv_data.append(idata) + return riv_data From 56bfbd4820b32a41bb5eeec9c1528d16fca80c6e Mon Sep 17 00:00:00 2001 From: Onno Ebbens Date: Tue, 10 Feb 2026 09:39:45 +0100 Subject: [PATCH 3/3] ruff --- nlmod/dims/grid.py | 17 +++++++---------- nlmod/gwf/surface_water.py | 10 +++++----- 2 files changed, 12 insertions(+), 15 deletions(-) diff --git a/nlmod/dims/grid.py b/nlmod/dims/grid.py index 631815e9..bb5d4a63 100644 --- a/nlmod/dims/grid.py +++ b/nlmod/dims/grid.py @@ -156,9 +156,7 @@ def get_icell2d_from_xy(x, y, ds, gi=None, rotated=True): msg = "get_icell2d_from_xy can only be applied to a vertex grid" assert ds.gridtype == "vertex", msg if gi is None: - gi = flopy.utils.GridIntersect( - modelgrid_from_ds(ds, rotated=rotated) - ) + gi = flopy.utils.GridIntersect(modelgrid_from_ds(ds, rotated=rotated)) cellids = gi.intersects(Point(x, y), dataframe=True)["cellid"].values if len(cellids) < 1: raise (ValueError(f"Point ({x}, {y}) is outside of the model grid")) @@ -1354,7 +1352,7 @@ def polygon_to_area(modelgrid, polygon, da, gridtype="structured"): area_array[row, col] = area elif gridtype == "vertex": area_array = util.get_da_from_da_ds(da, dims=("icell2d",), data=0) - area_array[df['cellid']] = df['areas'] + area_array[df["cellid"]] = df["areas"] return area_array @@ -1782,14 +1780,13 @@ def gdf_to_bool_da( multipolygon, contains_centroid=contains_centroid, min_area_fraction=min_area_fraction, - geo_dataframe=True - **kwargs, + geo_dataframe=True**kwargs, ) else: df = ix.intersects(multipolygon, dataframe=True) if df.shape[0] > 0 and ds.gridtype == "structured": - da.values[df['row'], df['col']] = True + da.values[df["row"], df["col"]] = True elif df.shape[0] > 0 and ds.gridtype == "vertex": da.values[df["cellid"]] = True @@ -1915,9 +1912,9 @@ def gdf_to_count_da(gdf, ds, ix=None, buffer=0.0, **kwargs): continue if ds.gridtype == "structured": - da.values[df['row'].values, df['col'].values] += 1 + da.values[df["row"].values, df["col"].values] += 1 elif ds.gridtype == "vertex": - da[df['cellid'].values] += 1 + da[df["cellid"].values] += 1 return da @@ -2131,7 +2128,7 @@ def gdf_area_to_da( ): df = ix.intersect(gdf.at[index, geometry], geo_dataframe=True, **kwargs) if structured: - for row, col, area in zip(df['row'], df['col'], df['areas']): + for row, col, area in zip(df["row"], df["col"], df["areas"]): data[(row, col, irow)] += area else: data[list(df["cellid"]), irow] = df["areas"] diff --git a/nlmod/gwf/surface_water.py b/nlmod/gwf/surface_water.py index ade38d71..a15c74ca 100644 --- a/nlmod/gwf/surface_water.py +++ b/nlmod/gwf/surface_water.py @@ -1230,15 +1230,15 @@ def rivdata_from_xylist(gwf, xylist, layer, stage, cond, rbot, aux=None): gi = flopy.utils.GridIntersect(gwf.modelgrid) df = gi.intersect(xylist, shapetype="linestring", geo_dataframe=True) riv_data = [] - - if 'row' in df.columns: # structured grid - for row, col in zip(df['row'], df['col']): + + if "row" in df.columns: # structured grid + for row, col in zip(df["row"], df["col"]): idata = [(layer, row, col), stage, cond, rbot] if aux is not None: idata.append(aux) riv_data.append(idata) - else: # unstructured grid - for cid in df['cellid']: + else: # unstructured grid + for cid in df["cellid"]: idata = [(layer, cid), stage, cond, rbot] if aux is not None: idata.append(aux)