diff --git a/docs/examples/02_surface_water.ipynb b/docs/examples/02_surface_water.ipynb
index e5a42897..25ef0574 100644
--- a/docs/examples/02_surface_water.ipynb
+++ b/docs/examples/02_surface_water.ipynb
@@ -47,7 +47,9 @@
     "import geopandas as gpd\n",
     "import matplotlib.pyplot as plt\n",
     "import nlmod\n",
-    "import numpy as np\n"
+    "import numpy as np\n",
+    "import warnings\n",
+    "from shapely.errors import ShapelyDeprecationWarning"
    ]
   },
   {
@@ -57,10 +59,12 @@
    "metadata": {},
    "outputs": [],
    "source": [
+    "warnings.filterwarnings(\"ignore\", category=ShapelyDeprecationWarning)\n",
+    "\n",
     "print(f\"nlmod version: {nlmod.__version__}\")\n",
     "\n",
     "# toon informatie bij het aanroepen van functies\n",
-    "logging.basicConfig(level=logging.INFO)\n"
+    "logging.basicConfig(level=logging.INFO)"
    ]
   },
   {
@@ -106,7 +110,7 @@
     "# set the stage of the Lek to 0.0 m NAP and the botm to -3 m NAP\n",
     "mask = sfw[\"bronhouder\"] == \"L0002\"\n",
     "sfw.loc[mask, \"stage\"] = 0.0\n",
-    "sfw.loc[mask, \"botm\"] = -3.0\n"
+    "sfw.loc[mask, \"botm\"] = -3.0"
    ]
   },
   {
@@ -147,7 +151,7 @@
     "ax.set_xlabel(\"X (m RD)\")\n",
     "ax.set_ylabel(\"Y (m RD)\")\n",
     "plt.yticks(rotation=90, va=\"center\")\n",
-    "fig.tight_layout()\n"
+    "fig.tight_layout()"
    ]
   },
   {
@@ -172,7 +176,7 @@
     "ax.set_xlabel(\"X (m RD)\")\n",
     "ax.set_ylabel(\"Y (m RD)\")\n",
     "plt.yticks(rotation=90, va=\"center\")\n",
-    "fig.tight_layout()\n"
+    "fig.tight_layout()"
    ]
   },
   {
@@ -206,7 +210,7 @@
     "figdir, cachedir = nlmod.util.get_model_dirs(model_ws)\n",
     "\n",
     "delr = delc = 50.0\n",
-    "start_time = \"2021-01-01\"\n"
+    "start_time = \"2021-01-01\""
    ]
   },
   {
@@ -224,7 +228,7 @@
     "    cachedir=cachedir,\n",
     "    cachename=\"combined_layer_ds.nc\",\n",
     ")\n",
-    "layer_model\n"
+    "layer_model"
    ]
   },
   {
@@ -242,7 +246,7 @@
     "# create model time dataset\n",
     "ds = nlmod.mdims.set_ds_time(ds, start_time=start_time, steady_state=True)\n",
     "\n",
-    "ds\n"
+    "ds"
    ]
   },
   {
@@ -274,7 +278,7 @@
     "ic = nlmod.gwf.ic(ds, gwf, starting_head=1.0)\n",
     "\n",
     "# Create the output control package\n",
-    "oc = nlmod.gwf.oc(ds, gwf)\n"
+    "oc = nlmod.gwf.oc(ds, gwf)"
    ]
   },
   {
@@ -311,7 +315,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "sfw_grid = nlmod.mdims.gdf_to_grid(sfw, gwf)\n"
+    "sfw_grid = nlmod.mdims.gdf_to_grid(sfw, gwf)"
    ]
   },
   {
@@ -359,7 +363,7 @@
     "cid = (107, 6)  # for 50 x 50 m grid\n",
     "# cid = (5, 45)  # for 100 x 100 m grid\n",
     "mask = sfw_grid.cellid == cid\n",
-    "sfw_grid.loc[mask]\n"
+    "sfw_grid.loc[mask]"
    ]
   },
   {
@@ -421,7 +425,7 @@
     "try:\n",
     "    nlmod.gwf.surface_water.aggregate(sfw_grid, \"area_weighted\")\n",
     "except ValueError as e:\n",
-    "    print(e)\n"
+    "    print(e)"
    ]
   },
   {
@@ -439,7 +443,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "sfw_grid[\"c0\"] = 1.0  # days\n"
+    "sfw_grid[\"c0\"] = 1.0  # days"
    ]
   },
   {
@@ -459,7 +463,7 @@
    "source": [
     "celldata = nlmod.gwf.surface_water.aggregate(\n",
     "    sfw_grid, \"area_weighted\"\n",
-    ")\n"
+    ")"
    ]
   },
   {
@@ -512,7 +516,7 @@
    "outputs": [],
    "source": [
     "new_celldata = celldata.loc[~celldata.rbot.isna()]\n",
-    "print(f\"removed {len(celldata)-len(new_celldata)} reaches because rbot is nan\")\n"
+    "print(f\"removed {len(celldata)-len(new_celldata)} reaches because rbot is nan\")"
    ]
   },
   {
@@ -522,7 +526,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "riv_spd = nlmod.gwf.surface_water.build_spd(new_celldata, \"RIV\", ds)\n"
+    "riv_spd = nlmod.gwf.surface_water.build_spd(new_celldata, \"RIV\", ds)"
    ]
   },
   {
@@ -560,7 +564,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "riv = flopy.mf6.ModflowGwfriv(gwf, stress_period_data=riv_spd)\n"
+    "riv = flopy.mf6.ModflowGwfriv(gwf, stress_period_data=riv_spd)"
    ]
   },
   {
@@ -580,7 +584,7 @@
    "source": [
     "fig, ax = plt.subplots(1, 1, figsize=(10, 8), constrained_layout=True)\n",
     "mv = flopy.plot.PlotMapView(model=gwf, ax=ax, layer=0)\n",
-    "mv.plot_bc(\"RIV\")\n"
+    "mv.plot_bc(\"RIV\")"
    ]
   },
   {
@@ -602,7 +606,7 @@
    "source": [
     "nlmod.sim.write_and_run(\n",
     "    sim, ds, write_ds=True, nb_path=\"02_surface_water.ipynb\"\n",
-    ")\n"
+    ")"
    ]
   },
   {
@@ -694,14 +698,6 @@
     "col = gwf.modelgrid.ncol // 2\n",
     "ax.set_title(f\"Cross-section along column {col}\");\n"
    ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "b4478d66",
-   "metadata": {},
-   "outputs": [],
-   "source": []
   }
  ],
  "metadata": {
@@ -720,7 +716,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.9.4"
+   "version": "3.9.7"
   },
   "vscode": {
    "interpreter": {
diff --git a/docs/examples/03_local_grid_refinement.ipynb b/docs/examples/03_local_grid_refinement.ipynb
index b3abdb37..0cc13fcd 100644
--- a/docs/examples/03_local_grid_refinement.ipynb
+++ b/docs/examples/03_local_grid_refinement.ipynb
@@ -37,7 +37,9 @@
     "import flopy\n",
     "import geopandas as gpd\n",
     "import matplotlib.pyplot as plt\n",
-    "import nlmod\n"
+    "import nlmod\n",
+    "import warnings\n",
+    "from shapely.errors import ShapelyDeprecationWarning"
    ]
   },
   {
@@ -46,6 +48,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
+    "warnings.filterwarnings(\"ignore\", category=ShapelyDeprecationWarning)\n",
+    "\n",
     "print(f'nlmod version: {nlmod.__version__}')\n",
     "\n",
     "# toon informatie bij het aanroepen van functies\n",
@@ -85,7 +89,7 @@
     "regis_botm_layer = \"MSz1\"\n",
     "use_geotop = True\n",
     "add_northsea = True\n",
-    "starting_head = 1.0\n"
+    "starting_head = 1.0"
    ]
   },
   {
@@ -116,7 +120,7 @@
     "    steady_start=steady_start,\n",
     "    transient_timesteps=transient_timesteps,\n",
     "    perlen=perlen,\n",
-    ")\n"
+    ")"
    ]
   },
   {
@@ -140,7 +144,7 @@
     ")\n",
     "\n",
     "if add_northsea:\n",
-    "    ds = nlmod.read.rws.add_northsea(ds, cachedir=cachedir)\n"
+    "    ds = nlmod.read.rws.add_northsea(ds, cachedir=cachedir)"
    ]
   },
   {
@@ -171,7 +175,7 @@
     "ic = nlmod.gwf.ic(ds, gwf, starting_head=starting_head)\n",
     "\n",
     "# Create the output control package\n",
-    "oc = nlmod.gwf.oc(ds, gwf)\n"
+    "oc = nlmod.gwf.oc(ds, gwf)"
    ]
   },
   {
@@ -197,7 +201,7 @@
     "\n",
     "# add constant head cells at model boundaries\n",
     "ds.update(nlmod.mgrid.mask_model_edge(ds, ds[\"idomain\"]))\n",
-    "chd = nlmod.gwf.chd(ds, gwf, chd=\"edge_mask\", head=\"starting_head\")\n"
+    "chd = nlmod.gwf.chd(ds, gwf, chd=\"edge_mask\", head=\"starting_head\")"
    ]
   },
   {
@@ -213,7 +217,7 @@
     "ds.update(knmi_ds)\n",
     "\n",
     "# create recharge package\n",
-    "rch = nlmod.gwf.rch(ds, gwf)\n"
+    "rch = nlmod.gwf.rch(ds, gwf)"
    ]
   },
   {
@@ -240,7 +244,7 @@
    "source": [
     "nlmod.sim.write_and_run(\n",
     "    sim, ds, write_ds=True, nb_path=\"03_local_grid_refinement.ipynb\"\n",
-    ")\n"
+    ")"
    ]
   },
   {
@@ -318,7 +322,7 @@
     "nlmod.plot.data_array(ds[\"bathymetry\"], ds, ax=axes[0][0])\n",
     "nlmod.plot.data_array(ds[\"northsea\"], ds, ax=axes[0][1])\n",
     "nlmod.plot.data_array(ds[\"kh\"][1], ds, ax=axes[1][0])\n",
-    "nlmod.plot.data_array(ds[\"recharge\"][:, 0], ds, ax=axes[1][1])\n"
+    "nlmod.plot.data_array(ds[\"recharge\"][:, 0], ds, ax=axes[1][1])"
    ]
   },
   {
diff --git a/docs/examples/07_resampling.ipynb b/docs/examples/07_resampling.ipynb
index cbe633e9..2d281b13 100644
--- a/docs/examples/07_resampling.ipynb
+++ b/docs/examples/07_resampling.ipynb
@@ -9,7 +9,7 @@
     "\n",
     "# Resampling \n",
     "\n",
-    "Resampling data is a very common operation when building a Modflow model. Usually it is used to project data from one grid onto the other. There are many different ways to do this. This notebook shows some examples of resampling methods that are incorporated in the `nlmod` package. These methods rely heavily on resampling methods in other packages such as `scipy.interpolate` and `xarray`."
+    "Resampling data is a very common operation when building a Modflow model. Usually it is used to project data from one grid onto the other. There are many different ways to do this. This notebook shows some examples of resampling methods that are incorporated in the `nlmod` package. These methods rely heavily on resampling methods in packages such as `rioxarray` and `scipy.interpolate`."
    ]
   },
   {
@@ -198,7 +198,7 @@
     "# create vertextured dataarray\n",
     "coords = dict(x=xr.DataArray(xvc, dims=['icell2d',]), y=xr.DataArray(yvc, dims=['icell2d',]))\n",
     "vertex1 = xr.DataArray(values, dims=('icell2d'), coords=coords)\n",
-    "nlmod.visualise.plots.plot_vertex_array(vertex1, vertices, gridkwargs={'edgecolor': 'k'});"
+    "nlmod.plot.plot_vertex_array(vertex1, vertices, gridkwargs={'edgecolor': 'k'});"
    ]
   },
   {
@@ -217,7 +217,7 @@
     "vertex1_nan = vertex1.copy().astype(float)\n",
     "vertex1_nan.values[7] = np.nan\n",
     "\n",
-    "nlmod.visualise.plots.plot_vertex_array(vertex1_nan, vertices, gridkwargs={'edgecolor': 'k'});"
+    "nlmod.plot.plot_vertex_array(vertex1_nan, vertices, gridkwargs={'edgecolor': 'k'});"
    ]
   },
   {
@@ -247,10 +247,10 @@
    "source": [
     "def compare_structured_data_arrays(da1, da2, method, edgecolor='k'):\n",
     "    fig, axes = plt.subplots(ncols=2, figsize=(12,6))\n",
-    "    da1.plot(ax=axes[0], edgecolor=edgecolor)\n",
+    "    da1.plot(ax=axes[0], edgecolor=edgecolor, vmin=0, vmax=9)\n",
     "    axes[0].set_aspect('equal')\n",
     "    axes[0].set_title('original grid')\n",
-    "    da2.plot(ax=axes[1], edgecolor=edgecolor)\n",
+    "    da2.plot(ax=axes[1], edgecolor=edgecolor, vmin=0, vmax=9)\n",
     "    axes[1].set_aspect('equal')\n",
     "    axes[1].set_title(f'resampled grid, method {method}')"
    ]
@@ -367,10 +367,10 @@
    "source": [
     "def compare_struct_to_vertex(struc2d, res_vertex2d_n, vertices, method):\n",
     "    fig, axes = plt.subplots(ncols=2, figsize=(12,6))\n",
-    "    struc2d.plot(ax=axes[0], edgecolor='k')\n",
+    "    struc2d.plot(ax=axes[0], edgecolor='k', vmin=0, vmax=9)\n",
     "    axes[0].set_aspect('equal')\n",
     "    axes[0].set_title('structured grid')\n",
-    "    nlmod.visualise.plots.plot_vertex_array(res_vertex2d_n, vertices, ax=axes[1], gridkwargs={'edgecolor': 'k'})\n",
+    "    nlmod.plot.plot_vertex_array(res_vertex2d_n, vertices, ax=axes[1], gridkwargs={'edgecolor': 'k'}, vmin=0, vmax=9)\n",
     "    axes[1].set_title(f'locally refined grid, method {method}')"
    ]
   },
@@ -419,9 +419,9 @@
    "source": [
     "def compare_vertex_to_struct(vertex1, struc_out_n, method):\n",
     "    fig, axes = plt.subplots(ncols=2, figsize=(12,6))\n",
-    "    nlmod.visualise.plots.plot_vertex_array(vertex1, vertices, ax=axes[0], gridkwargs={'edgecolor': 'k'})\n",
+    "    nlmod.plot.plot_vertex_array(vertex1, vertices, ax=axes[0], gridkwargs={'edgecolor': 'k'}, vmin=0, vmax=9)\n",
     "    axes[0].set_title('original')\n",
-    "    struc_out_n.plot(ax=axes[1], edgecolor='k')\n",
+    "    struc_out_n.plot(ax=axes[1], edgecolor='k', vmin=0, vmax=9)\n",
     "    axes[1].set_title(f'resampled, method {method}')\n",
     "    axes[1].set_aspect('equal')"
    ]
@@ -502,9 +502,9 @@
    "source": [
     "def compare_vertex_arrays(vertex1, vertex2, method):\n",
     "    fig, axes = plt.subplots(ncols=2, figsize=(12,6))\n",
-    "    nlmod.visualise.plots.plot_vertex_array(vertex1, vertices, ax=axes[0], gridkwargs={'edgecolor': 'k'})\n",
+    "    nlmod.plot.plot_vertex_array(vertex1, vertices, ax=axes[0], gridkwargs={'edgecolor': 'k'}, vmin=0, vmax=9)\n",
     "    axes[0].set_title('original')\n",
-    "    nlmod.visualise.plots.plot_vertex_array(vertex2, vertices, ax=axes[1], gridkwargs={'edgecolor': 'k'})\n",
+    "    nlmod.plot.plot_vertex_array(vertex2, vertices, ax=axes[1], gridkwargs={'edgecolor': 'k'}, vmin=0, vmax=9)\n",
     "    axes[1].set_title(f'resampled, method {method}')\n",
     "    axes[1].set_aspect('equal')"
    ]
@@ -549,7 +549,9 @@
    "source": [
     "point_geom = [Point(x,y) for x, y in zip([1000, 1200, 1225, 1300],[20200, 20175, 20175, 20425])]\n",
     "point_gdf = gpd.GeoDataFrame({'values':[1,52,66,24]}, geometry=point_geom)\n",
-    "line_geom = [LineString([point_geom[0], point_geom[1]]), LineString([point_geom[2], point_geom[3]]), LineString([point_geom[0], point_geom[3]])]\n",
+    "line_geom = [LineString([point_geom[0], point_geom[1]]),\n",
+    "             LineString([point_geom[2], point_geom[3]]),\n",
+    "             LineString([point_geom[0], point_geom[3]])]\n",
     "line_gdf = gpd.GeoDataFrame({'values':[1,52,66]}, geometry=line_geom)\n",
     "pol_geom = [shp_polygon([point_geom[0], point_geom[1], point_geom[2], point_geom[3], point_geom[0]]),\n",
     "            shp_polygon([point_geom[0], point_geom[1], point_geom[2], Point(1200,20300), point_geom[0]])]\n",
@@ -1100,13 +1102,6 @@
     "    axes[1].set_aspect('equal')\n",
     "    axes[1].set_title(f'resampled grid, method {method}')"
    ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": []
   }
  ],
  "metadata": {
diff --git a/docs/examples/09_schoonhoven.ipynb b/docs/examples/09_schoonhoven.ipynb
index 36243bc8..12e29a20 100644
--- a/docs/examples/09_schoonhoven.ipynb
+++ b/docs/examples/09_schoonhoven.ipynb
@@ -522,8 +522,6 @@
     "    head.sel(layer=\"HLc\").mean(\"time\"), ds=ds, edgecolor=\"k\", norm=norm\n",
     ")\n",
     "cbar = nlmod.plot.colorbar_inside(pc)\n",
-    "for label in cbar.ax.yaxis.get_ticklabels():\n",
-    "    label.set_bbox(dict(facecolor=\"w\", alpha=0.5))\n",
     "bgt.plot(ax=ax, edgecolor=\"k\", facecolor=\"none\")"
    ]
   },
@@ -550,8 +548,6 @@
     "pc = dcs.plot_array(head.mean(\"time\"), norm=norm, head=head.mean(\"time\"))\n",
     "# add labels with layer names\n",
     "cbar = nlmod.plot.colorbar_inside(pc, bounds=[0.05, 0.05, 0.02, 0.9])\n",
-    "for label in cbar.ax.yaxis.get_ticklabels():\n",
-    "    label.set_bbox(dict(facecolor=\"w\", alpha=0.5))\n",
     "dcs.plot_grid()\n",
     "dcs.plot_layers(alpha=0.0, min_label_area=1000)\n",
     "f.tight_layout(pad=0.0)"
diff --git a/docs/examples/11_grid_rotation.ipynb b/docs/examples/11_grid_rotation.ipynb
index 40799f57..68839849 100644
--- a/docs/examples/11_grid_rotation.ipynb
+++ b/docs/examples/11_grid_rotation.ipynb
@@ -34,10 +34,14 @@
     "import matplotlib\n",
     "import nlmod\n",
     "import pandas as pd\n",
+    "import warnings\n",
+    "from shapely.errors import ShapelyDeprecationWarning\n",
+    "\n",
+    "warnings.filterwarnings(\"ignore\", category=ShapelyDeprecationWarning)\n",
     "\n",
     "# set the log-level to INFO, so more information is shown\n",
     "# (compared to the default setting of WARNING)\n",
-    "logging.basicConfig(level=logging.INFO)\n"
+    "logging.basicConfig(level=logging.INFO)"
    ]
   },
   {
@@ -66,7 +70,7 @@
     "    model_ws=\"model11\",\n",
     ")\n",
     "\n",
-    "ds = nlmod.mdims.set_ds_time(ds, time=pd.date_range(\"2015\", \"2022\", freq=\"MS\"))\n"
+    "ds = nlmod.mdims.set_ds_time(ds, time=pd.date_range(\"2015\", \"2022\", freq=\"MS\"))"
    ]
   },
   {
@@ -119,7 +123,7 @@
     "pc = nlmod.plot.data_array(\n",
     "    ds[\"ahn\"], ds=ds, ax=axes[1], rotated=True, norm=norm, edgecolor=\"face\"\n",
     ")\n",
-    "nlmod.plot.colorbar_inside(pc, ax=axes[1])\n"
+    "nlmod.plot.colorbar_inside(pc, ax=axes[1])"
    ]
   },
   {
@@ -215,7 +219,7 @@
     "rch = nlmod.gwf.rch(ds, gwf)\n",
     "\n",
     "# create storage package\n",
-    "sto = nlmod.gwf.sto(ds, gwf)\n"
+    "sto = nlmod.gwf.sto(ds, gwf)"
    ]
   },
   {
@@ -280,9 +284,7 @@
     "    head.sel(layer=1).mean(\"time\"), ds=ds, edgecolor=\"k\", rotated=True\n",
     ")\n",
     "cbar = nlmod.plot.colorbar_inside(pc)\n",
-    "for label in cbar.ax.yaxis.get_ticklabels():\n",
-    "    label.set_bbox(dict(facecolor=\"w\", alpha=0.5))\n",
-    "bgt.plot(ax=ax, edgecolor=\"k\", facecolor=\"none\")\n"
+    "bgt.plot(ax=ax, edgecolor=\"k\", facecolor=\"none\")"
    ]
   },
   {
@@ -303,14 +305,6 @@
     "fname = os.path.join(ds.model_ws, 'ugrid_ds.nc')\n",
     "nlmod.visualise.gis.model_dataset_to_ugrid_nc_file(ds, fname)"
    ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "2daf659e-7330-48ec-9b38-7ae135feedfd",
-   "metadata": {},
-   "outputs": [],
-   "source": []
   }
  ],
  "metadata": {
diff --git a/docs/examples/13_plot_methods.ipynb b/docs/examples/13_plot_methods.ipynb
index 87db3c6b..832127c7 100644
--- a/docs/examples/13_plot_methods.ipynb
+++ b/docs/examples/13_plot_methods.ipynb
@@ -6,7 +6,9 @@
    "metadata": {},
    "source": [
     "# Plot methods\n",
-    "This notebook shows the plot methods that are available in nlmod. Most plot methods use a model Dataset as input, which is an xarray Dataset with some required variables and attributes. There are some plot methods in flopy as well, whcih require a grounwwater flow model, or a modelgrid?"
+    "This notebook shows the plot methods that are available in nlmod. Most plot methods use a model Dataset as input, which is an xarray Dataset with some required variables and attributes.\n",
+    "\n",
+    "There are some plot methods in flopy as well, which require a flopy modelgrid. This notebook shows the plot methods in nlmod and flopy create similar plots, so the user can choose the methods that best fit their needs."
    ]
   },
   {
@@ -23,6 +25,10 @@
     "import flopy\n",
     "import nlmod\n",
     "from nlmod.visualise.netcdf import DatasetCrossSection\n",
+    "import warnings\n",
+    "from shapely.errors import ShapelyDeprecationWarning\n",
+    "\n",
+    "warnings.filterwarnings(\"ignore\", category=ShapelyDeprecationWarning)\n",
     "\n",
     "# set the log-level to INFO, so more information is shown (compared to the default setting of WARNING)\n",
     "logging.basicConfig(level=logging.INFO)"
diff --git a/nlmod/mdims/mgrid.py b/nlmod/mdims/mgrid.py
index 483bb80d..1281ac0f 100644
--- a/nlmod/mdims/mgrid.py
+++ b/nlmod/mdims/mgrid.py
@@ -58,7 +58,7 @@ def xy_to_icell2d(xy, ds):
         number of the icell2d value of a cell containing the xy point.
     """
 
-    icell2d = (np.abs(ds.x - xy[0]) + np.abs(ds.y - xy[1])).argmin().item()
+    icell2d = (np.abs(ds.x.data - xy[0]) + np.abs(ds.y.data - xy[1])).argmin().item()
 
     return icell2d
 
@@ -292,7 +292,7 @@ def refine(
                         )
                     mask = geom_types == geom_type
                     # features = [gdf[mask].unary_union]
-                    features = list(gdf[mask].geometry.explode())
+                    features = list(gdf[mask].geometry.explode(index_parts=True))
                     g.add_refinement_features(features, geom_type, level, layers=[0])
     g.build()
     gridprops = g.get_gridprops_disv()
diff --git a/nlmod/visualise/plots.py b/nlmod/visualise/plots.py
index e78b7c65..e6c3f131 100644
--- a/nlmod/visualise/plots.py
+++ b/nlmod/visualise/plots.py
@@ -594,7 +594,9 @@ def fmt_rd_ticks(x, _):
     ax.yaxis.set_major_formatter(FuncFormatter(fmt_rd_ticks))
 
 
-def colorbar_inside(mappable=None, ax=None, norm=None, cmap=None, bounds=None, **kw):
+def colorbar_inside(
+    mappable=None, ax=None, norm=None, cmap=None, bounds=None, bbox_labels=True, **kw
+):
     """Place a colorbar inside an axes."""
     if ax is None:
         ax = plt.gca()
@@ -609,6 +611,12 @@ def colorbar_inside(mappable=None, ax=None, norm=None, cmap=None, bounds=None, *
     if bounds[0] > 0.5:
         cax.yaxis.tick_left()
         cax.yaxis.set_label_position("left")
+    if isinstance(bbox_labels, bool) and bbox_labels is True:
+        bbox_labels = dict(facecolor="w", alpha=0.5)
+    if isinstance(bbox_labels, dict):
+        for label in cb.ax.yaxis.get_ticklabels():
+            label.set_bbox(bbox_labels)
+
     return cb