Release 0.11.2

docs/advanced_stress_packages/01_lake.ipynb

@@ -181,7 +181,7 @@
     "mask = lak_grid.area > 0.5 * ds[\"area\"].sel(icell2d=lak_grid.index)\n",
     "lak_grid = lak_grid[mask]\n",
     "# set the geometry to the entire cell\n",
-    "gi = flopy.utils.GridIntersect(nlmod.grid.modelgrid_from_ds(ds), method=\"vertex\")\n",
+    "gi = flopy.utils.GridIntersect(nlmod.grid.modelgrid_from_ds(ds))\n",
     "lak_grid.geometry = gi.geoms[lak_grid.index]\n",
     "\n",
     "# remove drains that overlap with the lake\n",

docs/generate_logo_text.py

@@ -0,0 +1,122 @@
+# %%
+import os
+import geopandas as gpd
+from shapely.geometry import Polygon, MultiPolygon
+from matplotlib.textpath import TextPath
+from matplotlib.font_manager import FontProperties
+import matplotlib.pyplot as plt
+
+import nlmod
+import numpy as np
+
+
+# %%
+def string_to_gdf(text, font, size=500, x_offset=0, **kwargs):
+    tp = TextPath((x_offset, 0), text, size=size, prop=font)
+    polygons = []
+
+    # Matplotlib returns list of polygons (outer + inner contours)
+    for poly in tp.to_polygons():
+        if len(poly) >= 3:
+            polygons.append(Polygon(poly))
+
+    # if a polygon has holes, we need to reconstruct it
+    final_polygons = []
+    used = [False] * len(polygons)
+    for i, outer in enumerate(polygons):
+        if used[i]:
+            continue
+        holes = []
+        for j, inner in enumerate(polygons):
+            if i != j and not used[j]:
+                if outer.contains(inner):
+                    holes.append(inner.exterior.coords)
+                    used[j] = True
+        final_polygons.append(Polygon(outer.exterior.coords, holes))
+        used[i] = True
+
+    gdf = gpd.GeoDataFrame(geometry=final_polygons, **kwargs)
+    return gdf
+
+
+text = "NLMOD"
+# font_path = "C:/Windows/Fonts/consola.ttf"  # Adjust for your system
+font_path = "C:/Windows/Fonts/cour.ttf"  # Courier font
+font = FontProperties(fname=font_path)
+text_size = 100  # Larger size -> smoother curves
+
+gdf = string_to_gdf(text, font, size=text_size)
+
+if False:
+    gdf.plot()
+    plt.axis("equal")
+    plt.show()
+
+dx = 10
+extent = gdf.total_bounds[[0, 2, 1, 3]]
+
+# add a buffer of 10 % of width
+buffer = 0.1 * (extent[1] - extent[0])
+extent = np.array(
+    [
+        extent[0] - buffer,
+        extent[1] + buffer,
+        extent[2] - buffer,
+        extent[3] + buffer,
+    ]
+)
+if False:
+    # make sure vertical extent is at least half of the horizontal extent
+    height = extent[3] - extent[2]
+    width = extent[1] - extent[0]
+    center_y = 0.5 * (extent[2] + extent[3])
+    if height < 0.5 * width:
+        height = 0.5 * width
+    extent = np.array(
+        [
+            extent[0],
+            extent[1],
+            center_y - 0.5 * height,
+            center_y + 0.5 * height,
+        ]
+    )
+
+extent = (extent / dx).round() * dx
+ds = nlmod.get_ds(extent, dx)
+ds = nlmod.grid.refine(ds, "logo", [(gdf, 2)])
+
+nlmod.plot.modelgrid(ds)
+
+# %% plot the logo
+# add 1 percent margin
+margin = 0.01 * (extent[1] - extent[0])
+extent = (
+    extent[0] - margin,
+    extent[1] + margin,
+    extent[2] - margin,
+    extent[3] + margin,
+)
+figwidth = 5
+figheight = figwidth * (extent[3] - extent[2]) / (extent[1] - extent[0])
+f = plt.figure(figsize=(figwidth, figheight))
+ax = f.add_axes([0, 0, 1, 1])
+ax.axis("equal")
+ax.axis(extent)
+color = "k"
+nlmod.plot.modelgrid(ds, color=color, ax=ax, linewidth=0.5, clip_on=False, antialiased=False)
+
+ax.set_xlabel("")
+ax.set_ylabel("")
+ax.set_title("")
+ax.axis("off")
+
+# %% save logo
+fname = f"logo_text"
+dpi = 300
+if figwidth != 5:
+    fname = f"{fname}_{figwidth}"
+    dpi = None
+f.savefig(os.path.join("_static", f"{fname}.png"), dpi=dpi)
+f.savefig(os.path.join("_static", f"{fname}.svg"))
+
+# %%

docs/getting_started.rst

@@ -111,7 +111,6 @@ On top of that there are some optional dependecies:
 - scikit-image (used in nlmod.read.rws.calculate_sea_coverage)
 - py7zr (used in nlmod.read.bofek.download_bofek_gdf)
 - joblib (used in nlmod.cache)
-- colorama (used in nlmod.util.get_color_logger)
 - tqdm (used for showing progress in long-running methods)
 - hydropandas (used in nlmod.read.knmi and nlmod.read.bro)
 - owslib (used in nlmod.read.ahn.get_latest_ahn_from_wcs)

nlmod/dims/grid.py

@@ -156,10 +156,8 @@ 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), method="vertex"
-        )
-    cellids = gi.intersects(Point(x, y))["cellids"]
+        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"))
     icell2d = cellids[0]
@@ -227,7 +225,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]
@@ -622,7 +620,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 +1343,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
+        area_array = util.get_da_from_da_ds(da, dims=("y", "x"), data=0.0)
+        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 = util.get_da_from_da_ds(da, dims=("icell2d",), data=0.0)
+        area_array[df["cellid"].values] = df["areas"]
 
     return area_array
 
@@ -1770,7 +1766,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 +1776,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 +1889,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)
@@ -1909,18 +1905,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
 
@@ -1962,7 +1957,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 +1975,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 +2011,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 +2107,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 +2127,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 row, col, area in zip(df["row"], df["col"], df["areas"]):
+                data[(row, col, irow)] += area
         else:
-            data[list(r["cellids"]), irow] = r["areas"]
+            data[list(df["cellid"]), irow] = df["areas"]
 
     if sparse:
         try:
@@ -2451,6 +2443,22 @@ def get_extent_polygon(ds, rotated=True):
 
 
 def get_extent_gdf(ds, rotated=True, crs="EPSG:28992"):
+    """Get the model extent as a GeoDataFrame with a polygon.
+
+    Parameters
+    ----------
+    ds : xr.Dataset
+        model dataset.
+    rotated : bool, optional
+        if True, the extent is corrected for angrot. The default is True.
+    crs : str, optional
+        Coordinate reference system. The default is "EPSG:28992".
+
+    Returns
+    -------
+    gdf : geopandas.GeoDataFrame
+        GeoDataFrame containing the model extent as a polygon geometry.
+    """
     polygon = get_extent_polygon(ds, rotated=rotated)
     return gpd.GeoDataFrame(geometry=[polygon], crs=crs)
 

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)

nlmod/dims/time.py

@@ -394,7 +394,7 @@ def set_ds_time_numeric(
 
 
 def set_time_variables(ds, start, time, steady, steady_start, time_units, nstp, tsmult):
-    """Add data variables: steady, nstp and tsmult, set attributes: start, time_units
+    """Add data variables: steady, nstp and tsmult, set attributes: start, time_units.
 
     Parameters
     ----------

nlmod/gwf/lake.py

@@ -405,7 +405,7 @@ def _copy_da_from_ds(gdf, ds, variable, boundname_column=None, set_to_0_in_ds=Fa
         else:
             da_cells = ds[variable].loc[cellids].copy()
             if set_to_0_in_ds:
-                ds[variable][:, cellids] = 0.0
+                ds[variable][cellids] = 0.0
             # calculate thea area-weighted mean
             df[column] = float((da_cells * area).sum("icell2d") / area.sum())
     return df