Improve lake and other stuff

docs/examples/09_schoonhoven.ipynb

@@ -119,7 +119,7 @@
     "norm = matplotlib.colors.Normalize(vmin=-3, vmax=1)\n",
     "cmap = \"viridis\"\n",
     "bgt.plot(\"summer_stage\", ax=ax, norm=norm, cmap=cmap)\n",
-    "nlmod.plot.colorbar_inside(norm=norm, cmap=cmap)"
+    "nlmod.plot.colorbar_inside(norm=norm, cmap=cmap);"
    ]
   },
   {
@@ -139,7 +139,7 @@
    "source": [
     "f, ax = nlmod.plot.get_map(extent)\n",
     "bgt.plot(\"ahn_min\", ax=ax, norm=norm, cmap=cmap)\n",
-    "nlmod.plot.colorbar_inside(norm=norm, cmap=cmap)"
+    "nlmod.plot.colorbar_inside(norm=norm, cmap=cmap);"
    ]
   },
   {
@@ -308,9 +308,7 @@
    "source": [
     "bed_resistance = 1.0\n",
     "\n",
-    "mg = nlmod.grid.modelgrid_from_ds(ds)\n",
-    "gi = flopy.utils.GridIntersect(mg, method=\"vertex\")\n",
-    "bgt_grid = nlmod.grid.gdf_to_grid(bgt, ix=gi).set_index(\"cellid\")\n",
+    "bgt_grid = nlmod.grid.gdf_to_grid(bgt, ds).set_index(\"cellid\")\n",
     "\n",
     "bgt_grid[\"cond\"] = bgt_grid.area / bed_resistance\n",
     "\n",
@@ -332,14 +330,17 @@
     "    \"W0656.cb3c3a25de4141d18c573b561f02e84a\",\n",
     "]\n",
     "\n",
-    "mask = bgt_grid[\"identificatie\"].isin(ids_grote_gracht) | bgt_grid[\n",
-    "    \"identificatie\"\n",
-    "].isin(ids_oude_haven)\n",
+    "mask = bgt_grid[\"identificatie\"].isin(ids_grote_gracht + ids_oude_haven)\n",
     "lakes = bgt_grid[mask].copy()\n",
     "lakes[\"name\"] = \"\"\n",
     "lakes.loc[lakes[\"identificatie\"].isin(ids_grote_gracht), \"name\"] = \"grotegracht\"\n",
     "lakes.loc[lakes[\"identificatie\"].isin(ids_oude_haven), \"name\"] = \"oudehaven\"\n",
-    "bgt_grid = bgt_grid[~mask]"
+    "bgt_grid = bgt_grid[~mask]\n",
+    "\n",
+    "# cut rainfall and evaporation from model dataset\n",
+    "lak_rainfall, lak_evaporation = nlmod.gwf.lake.cut_meteorological_data_from_ds(\n",
+    "    lakes, ds, boundname_column=\"name\"\n",
+    ")"
    ]
   },
   {
@@ -400,12 +401,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "# add specific properties to the lake gdf\n",
-    "lakes.loc[lakes[\"identificatie\"].isin(ids_grote_gracht), \"lakeno\"] = 0\n",
-    "lakes.loc[lakes[\"identificatie\"].isin(ids_oude_haven), \"lakeno\"] = 1\n",
-    "\n",
     "# add general properties to the lake gdf\n",
-    "lakes[\"elev\"] = lakes[\"ahn_min\"] - 0.5\n",
     "summer_months = (4, 5, 6, 7, 8, 9)\n",
     "if pd.to_datetime(ds.time.start).month in summer_months:\n",
     "    lakes[\"strt\"] = lakes[\"summer_stage\"]\n",
@@ -418,13 +414,11 @@
     "# lakes.loc[lakes['identificatie'].isin(ids_oude_haven), 'INFLOW'] = 'inflow_lake'\n",
     "\n",
     "# add outlet to Oude Haven, water flows from Oude Haven to Grote Gracht.\n",
-    "lakes.loc[lakes[\"identificatie\"].isin(ids_oude_haven), \"lakeout\"] = 0\n",
-    "lakes.loc[lakes[\"identificatie\"].isin(ids_oude_haven), \"outlet_invert\"] = (\n",
-    "    1.0  # overstort hoogte\n",
-    ")\n",
+    "lakes.loc[lakes['name'] == \"oudehaven\", \"lakeout\"] = \"grotegracht\"\n",
+    "lakes.loc[lakes['name'] == \"oudehaven\", \"outlet_invert\"] = 1.0 # overstort hoogte\n",
     "\n",
     "# add lake to groundwaterflow model\n",
-    "lak = nlmod.gwf.lake_from_gdf(gwf, lakes, ds, boundname_column=\"name\")"
+    "lak = nlmod.gwf.lake_from_gdf(gwf, lakes, ds, boundname_column=\"name\", rainfall=lak_rainfall, evaporation=lak_evaporation)"
    ]
   },
   {
@@ -501,7 +495,7 @@
     "    colorbar_label=\"head [m NAP]\",\n",
     "    title=\"mean head\",\n",
     ")\n",
-    "bgt.dissolve().plot(ax=pc.axes, edgecolor=\"k\", facecolor=\"none\")"
+    "bgt.dissolve().plot(ax=pc.axes, edgecolor=\"k\", facecolor=\"none\");"
    ]
   },
   {
@@ -784,7 +778,7 @@
     "line = ax.add_collection(lc)\n",
     "nlmod.plot.colorbar_inside(line, label=\"Travel time (years)\")\n",
     "\n",
-    "bgt.dissolve().plot(ax=ax, edgecolor=\"k\", facecolor=\"none\")"
+    "bgt.dissolve().plot(ax=ax, edgecolor=\"k\", facecolor=\"none\");"
    ]
   },
   {

nlmod/dims/grid.py

@@ -100,13 +100,60 @@ def xy_to_icell2d(xy, ds):
     icell2d : int
         number of the icell2d value of a cell containing the xy point.
     """
+    logger.warning(
+        "nlmod.grid.xy_to_icell2d is deprecated. "
+        "Use nlmod.grid.get_icell2d_from_xy instead"
+    )
     msg = "xy_to_icell2d can only be applied to a vertex grid"
     assert ds.gridtype == "vertex", msg
     icell2d = (np.abs(ds.x.data - xy[0]) + np.abs(ds.y.data - xy[1])).argmin().item()
 
     return icell2d
 
 
+def get_icell2d_from_xy(x, y, ds, gi=None, rotated=True):
+    """Get the icell2d value of a point defined by its x and y coordinates.
+
+    Parameters
+    ----------
+    x : float
+        The x-coordinate of the point.
+    y : float
+        The y-coordinate of the point.
+    ds : xarary dataset
+        The model dataset.
+    gi : flopy.utils.GridIntersect, optional
+        Can be supplied to speed up the calculation, as the generation of the
+        GridIntersect-instance can take some time. The default is None.
+    rotated : bool, optional
+        If the model grid has a rotation, and rotated is False, x and y are in model
+        coordinates. Otherwise x and y are in real world coordinates. The default is
+        True.
+
+    Raises
+    ------
+
+        Raises a ValueError if the point is outside of the model grid.
+
+    Returns
+    -------
+    icell2d : int
+        The icell2d-number of the model cell containing the point (zero-based)
+
+    """
+    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"]
+    if len(cellids) < 1:
+        raise (ValueError(f"Point ({x}, {y}) is outside of the model grid"))
+    icell2d = cellids[0]
+    return icell2d
+
+
 def xy_to_row_col(xy, ds):
     """Get the row and column values of a point defined by its x and y coordinates.
 
@@ -124,13 +171,71 @@ def xy_to_row_col(xy, ds):
     col : int
         number of the column value of a cell containing the xy point.
     """
+
+    logger.warning(
+        "nlmod.grid.xy_to_row_col is deprecated. "
+        "Use nlmod.grid.get_row_col_from_xy instead"
+    )
     msg = "xy_to_row_col can only be applied to a structured grid"
     assert ds.gridtype == "structured", msg
     row = np.abs(ds.y.data - xy[1]).argmin()
     col = np.abs(ds.x.data - xy[0]).argmin()
     return row, col
 
 
+def get_row_col_from_xy(x, y, ds, rotated=True, gi=None):
+    """Get the row and column of a point defined by a x and y coordinate.
+
+    Parameters
+    ----------
+    x : float
+        The x-coordinate of the point.
+    y : float
+        The y-coordinate of the point.
+    ds : xarary dataset
+        The model dataset.
+    rotated : bool, optional
+        If the model grid has a rotation, and rotated is False, x and y are in model
+        coordinates. Otherwise x and y are in real world coordinates. The default is
+        True.
+    gi : flopy.utils.GridIntersect, optional
+        Can be supplied to use a GridIntersect-class instead of our own calculation
+
+    Raises
+    ------
+
+        Raises a ValueError if the point is outside of the model grid.
+
+    Returns
+    -------
+    row : int
+        The row-number of the model cell containing the point (zero-based)
+    col : int
+        The column-number of the model cell containing the point (zero-based)
+    """
+    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"]
+        if len(cellids) < 1:
+            raise (ValueError(f"Point ({x}, {y}) is outside of the model grid"))
+        row, col = cellids[0]
+        return row, col
+    if rotated and ("angrot" in ds.attrs) and (ds.attrs["angrot"] != 0.0):
+        # calculate the x and y in model coordinates
+        affine = get_affine_world_to_mod(ds)
+        x, y = affine * (x, y)
+    if x < ds.extent[0] or x > ds.extent[1] or y < ds.extent[2] or y > ds.extent[3]:
+        raise (ValueError(f"Point ({x}, {y}) is outside of the model grid"))
+
+    y_bot = ds.y - get_delc(ds) / 2
+    row = np.where(y >= y_bot)[0][0]
+
+    x_left = ds.x - get_delr(ds) / 2
+    col = np.where(x >= x_left)[0][-1]
+    return row, col
+
+
 def xyz_to_cid(xyz, ds=None, modelgrid=None):
     """Get the icell2d value of a point defined by its x and y coordinates.
 
@@ -754,7 +859,7 @@ def col_to_list(col_in, ds, cellids):
     Parameters
     ----------
     col_in : xarray.DatArray, str, int or float
-        if col_in is a str type it is the name of the column in ds (if it exists).
+        if col_in is a str type it is the name of the variable in ds (if it exists).
         if col_in is an int or a float it is a value that will be used for all
         cells in cellids.
     ds : xarray.Dataset