diff --git a/mikeio/dataset/_data_plot.py b/mikeio/dataset/_data_plot.py
index 6a2af5810..083aba85a 100644
--- a/mikeio/dataset/_data_plot.py
+++ b/mikeio/dataset/_data_plot.py
@@ -563,7 +563,7 @@ def _plot_FM_map(self, ax: Axes, **kwargs: Any) -> Axes:
             node_coordinates=geometry.node_coordinates,
             element_table=geometry.element_table,
             element_coordinates=geometry.element_coordinates,
-            boundary_polylines=geometry.boundary_polylines,
+            boundary_polylines=geometry.boundary_polylines.lines,
             projection=geometry.projection,
             z=values,
             ax=ax,
diff --git a/mikeio/spatial/_FM_geometry.py b/mikeio/spatial/_FM_geometry.py
index 6028f78d8..082d194da 100644
--- a/mikeio/spatial/_FM_geometry.py
+++ b/mikeio/spatial/_FM_geometry.py
@@ -1,5 +1,4 @@
 from __future__ import annotations
-from collections import namedtuple
 from functools import cached_property
 from pathlib import Path
 from typing import (
@@ -9,6 +8,7 @@
     Sized,
     TYPE_CHECKING,
 )
+import warnings
 
 
 import numpy as np
@@ -23,7 +23,8 @@
 from ._FM_utils import (
     _get_node_centered_data,
     _plot_map,
-    BoundaryPolylines,
+    BoundaryPolygons,
+    Polygon,
     _set_xy_label_by_projection,  # TODO remove
     _to_polygons,  # TODO remove
 )
@@ -113,7 +114,7 @@ def _plot_FM_map(self, ax: Axes, **kwargs: Any) -> Axes:
             node_coordinates=g.node_coordinates,
             element_table=g.element_table,
             element_coordinates=g.element_coordinates,
-            boundary_polylines=g.boundary_polylines,
+            boundary_polylines=g.boundary_polygons.lines,
             plot_type=plot_type,
             projection=g.projection,
             z=None,
@@ -158,9 +159,9 @@ def outline(
 
         linwid = 1.2
         out_col = "0.4"
-        for exterior in self.g.boundary_polylines.exteriors:
+        for exterior in self.g.boundary_polygons.exteriors:
             ax.plot(*exterior.xy.T, color=out_col, linewidth=linwid)
-        for interior in self.g.boundary_polylines.interiors:
+        for interior in self.g.boundary_polygons.interiors:
             ax.plot(*interior.xy.T, color=out_col, linewidth=linwid)
         if title is not None:
             ax.set_title(title)
@@ -843,9 +844,16 @@ def get_element_area(self) -> np.ndarray:
         return np.abs(area)
 
     @cached_property
-    def boundary_polylines(self) -> BoundaryPolylines:
-        """Lists of closed polylines defining domain outline."""
-        return self._get_boundary_polylines()
+    def boundary_polylines(self) -> BoundaryPolygons:
+        warnings.warn(
+            "boundary_polylines is renamed to boundary_polygons", FutureWarning
+        )
+        return self._get_boundary_polygons()
+
+    @cached_property
+    def boundary_polygons(self) -> BoundaryPolygons:
+        """Lists of polygons defining domain outline."""
+        return self._get_boundary_polygons()
 
     def contains(self, points: np.ndarray) -> np.ndarray:
         """Test if a list of points are contained by mesh.
@@ -861,34 +869,18 @@ def contains(self, points: np.ndarray) -> np.ndarray:
             True for points inside, False otherwise
 
         """
-        import matplotlib.path as mp  # type: ignore
-
         points = np.atleast_2d(points)
 
-        exterior = self.boundary_polylines.exteriors[0]
-        cnts = mp.Path(exterior.xy).contains_points(points)
-
-        if self.boundary_polylines.n_exteriors > 1:
-            # in case of several dis-joint outer domains
-            for exterior in self.boundary_polylines.exteriors[1:]:
-                in_domain = mp.Path(exterior.xy).contains_points(points)
-                cnts = np.logical_or(cnts, in_domain)
-
-        # subtract any holes
-        for interior in self.boundary_polylines.interiors:
-            in_hole = mp.Path(interior.xy).contains_points(points)
-            cnts = np.logical_and(cnts, ~in_hole)
-
-        return cnts
+        return self.boundary_polylines.contains(points)
 
     def __contains__(self, pt: np.ndarray) -> bool:
         return self.contains(pt)[0]
 
-    def _get_boundary_polylines_uncategorized(self) -> list[list[np.int64]]:
-        """Construct closed polylines for all boundary faces."""
+    def _get_boundary_polygons_uncategorized(self) -> list[list[np.int64]]:
+        """Construct closed polygons for all boundary faces."""
         boundary_faces = self._get_boundary_faces()
         face_remains = boundary_faces.copy()
-        polylines = []
+        polygons = []
         while face_remains.shape[0] > 1:
             n0 = face_remains[:, 0]
             n1 = face_remains[:, 1]
@@ -907,34 +899,26 @@ def _get_boundary_polylines_uncategorized(self) -> list[list[np.int64]]:
                     break
 
             face_remains = np.delete(face_remains, index_to_delete, axis=0)
-            polylines.append(polyline)
-        return polylines
+            polygons.append(polyline)
+        return polygons
 
-    def _get_boundary_polylines(self) -> BoundaryPolylines:
+    def _get_boundary_polygons(self) -> BoundaryPolygons:
         """Get boundary polylines and categorize as inner or outer by assessing the signed area."""
-        polylines = self._get_boundary_polylines_uncategorized()
-
-        poly_lines_int = []
-        poly_lines_ext = []
-        Polyline = namedtuple("Polyline", ["n_nodes", "nodes", "xy", "area"])
-
-        for polyline in polylines:
-            xy = self.node_coordinates[polyline, :2]
-            area = (
-                np.dot(xy[:, 1], np.roll(xy[:, 0], 1))
-                - np.dot(xy[:, 0], np.roll(xy[:, 1], 1))
-            ) * 0.5
-            poly_line = np.asarray(polyline)
-            xy = self.node_coordinates[poly_line, 0:2]
-            poly = Polyline(len(polyline), poly_line, xy, area)
-            if area > 0:
-                poly_lines_ext.append(poly)
+        polygons = self._get_boundary_polygons_uncategorized()
+
+        interiors = []
+        exteriors = []
+
+        for polygon in polygons:
+            polygon_np = np.asarray(polygon)
+            xy = self.node_coordinates[polygon_np, 0:2]
+            poly = Polygon(xy)
+            if poly.area > 0:
+                exteriors.append(poly)
             else:
-                poly_lines_int.append(poly)
+                interiors.append(poly)
 
-        n_ext = len(poly_lines_ext)
-        n_int = len(poly_lines_int)
-        return BoundaryPolylines(n_ext, poly_lines_ext, n_int, poly_lines_int)
+        return BoundaryPolygons(exteriors=exteriors, interiors=interiors)
 
     def _get_boundary_faces(self) -> np.ndarray:
         """Construct list of faces."""
diff --git a/mikeio/spatial/_FM_geometry_layered.py b/mikeio/spatial/_FM_geometry_layered.py
index d9f527372..9a7e7c866 100644
--- a/mikeio/spatial/_FM_geometry_layered.py
+++ b/mikeio/spatial/_FM_geometry_layered.py
@@ -3,6 +3,7 @@
 from pathlib import Path
 
 from typing import Any, Iterable, Literal, Sequence
+import warnings
 
 from matplotlib.axes import Axes
 import numpy as np
@@ -12,7 +13,7 @@
 from ._FM_geometry import GeometryFM2D, _GeometryFM, _GeometryFMPlotter
 from ._geometry import GeometryPoint3D
 
-from ._FM_utils import _plot_vertical_profile, BoundaryPolylines
+from ._FM_utils import _plot_vertical_profile, BoundaryPolygons
 
 from ._utils import _relative_cumulative_distance
 
@@ -653,7 +654,14 @@ def __init__(
         self.plot = _GeometryFMPlotter(self)
 
     @property
-    def boundary_polylines(self) -> BoundaryPolylines:
+    def boundary_polylines(self) -> BoundaryPolygons:
+        warnings.warn(
+            "boundary_polylines is renamed to boundary_polygons", FutureWarning
+        )
+        return self.geometry2d.boundary_polylines
+
+    @property
+    def boundary_polygons(self) -> BoundaryPolygons:
         return self.geometry2d.boundary_polylines
 
     def contains(self, points: np.ndarray) -> np.ndarray:
diff --git a/mikeio/spatial/_FM_utils.py b/mikeio/spatial/_FM_utils.py
index d6388f8f2..29a3ffbc3 100644
--- a/mikeio/spatial/_FM_utils.py
+++ b/mikeio/spatial/_FM_utils.py
@@ -1,5 +1,7 @@
-from __future__ import annotations
+from dataclasses import dataclass
 from typing import Any, Literal, Sequence
+
+from numpy.typing import NDArray
 from matplotlib.axes import Axes
 from matplotlib.cm import ScalarMappable
 from matplotlib.collections import PatchCollection
@@ -8,7 +10,6 @@
 from matplotlib.tri import Triangulation
 import numpy as np
 from scipy.sparse import csr_matrix
-from collections import namedtuple
 
 from ._utils import _relative_cumulative_distance
 
@@ -16,17 +17,66 @@
 MESH_COL = "0.95"
 MESH_COL_DARK = "0.6"
 
-BoundaryPolylines = namedtuple(
-    "BoundaryPolylines",
-    ["n_exteriors", "exteriors", "n_interiors", "interiors"],
-)
+
+@dataclass
+class Polygon:
+    xy: NDArray
+
+    @property
+    def area(self) -> float:
+        return (
+            np.dot(self.xy[:, 1], np.roll(self.xy[:, 0], 1))
+            - np.dot(self.xy[:, 0], np.roll(self.xy[:, 1], 1))
+        ) * 0.5
+
+
+@dataclass
+class BoundaryPolygons:
+    exteriors: list[Polygon]
+    interiors: list[Polygon]
+
+    @property
+    def lines(self) -> list[Polygon]:
+        return self.exteriors + self.interiors
+
+    def contains(self, points: np.ndarray) -> np.ndarray:
+        """Test if a list of points are contained by mesh.
+
+        Parameters
+        ----------
+        points : array-like n-by-2
+            x,y-coordinates of n points to be tested
+
+        Returns
+        -------
+        bool array
+            True for points inside, False otherwise
+
+        """
+        import matplotlib.path as mp  # type: ignore
+
+        exterior = self.exteriors[0]
+        cnts = mp.Path(exterior.xy).contains_points(points)
+
+        if len(self.exteriors) > 1:
+            # in case of several dis-joint outer domains
+            for exterior in self.exteriors[1:]:
+                in_domain = mp.Path(exterior.xy).contains_points(points)
+                cnts = np.logical_or(cnts, in_domain)
+
+        # subtract any holes
+        for interior in self.interiors:
+            in_hole = mp.Path(interior.xy).contains_points(points)
+            cnts = np.logical_and(cnts, ~in_hole)
+
+        return cnts
 
 
 def _plot_map(
     node_coordinates: np.ndarray,
     element_table: np.ndarray,
     element_coordinates: np.ndarray,
-    boundary_polylines: BoundaryPolylines,
+    boundary_polylines: list[Polygon],
     projection: str = "",
     z: np.ndarray | None = None,
     plot_type: Literal[
@@ -347,14 +397,14 @@ def __plot_mesh_only(ax: Axes, nc: np.ndarray, element_table: np.ndarray) -> Non
     ax.add_collection(fig_obj)
 
 
-def __plot_outline_only(ax: Axes, boundary_polylines: BoundaryPolylines) -> Axes:
+def __plot_outline_only(ax: Axes, boundary_polylines: list[Polygon]) -> Axes:
     """plot outline only (no data).
 
     Parameters
     ----------
     ax : matplotlib.axes.Axes
         axes object
-    boundary_polylines : BoundaryPolylines
+    boundary_polylines : list[PolyLine]
         boundary polylines
 
     Returns
@@ -583,14 +633,14 @@ def __add_non_tri_mesh(
     ax.add_collection(p)
 
 
-def __add_outline(ax: Axes, boundary_polylines: BoundaryPolylines) -> None:
+def __add_outline(ax: Axes, boundary_polylines: list[Polygon]) -> None:
     """add outline to axes.
 
     Parameters
     ----------
     ax : matplotlib.axes.Axes
         axes object
-    boundary_polylines: BoundaryPolylines
+    boundary_polylines: list[PolyLine]
         boundary polylines
 
     Returns
@@ -598,8 +648,7 @@ def __add_outline(ax: Axes, boundary_polylines: BoundaryPolylines) -> None:
     None
 
     """
-    lines = boundary_polylines.exteriors + boundary_polylines.interiors
-    for line in lines:
+    for line in boundary_polylines:
         ax.plot(*line.xy.T, color="0.4", linewidth=1.2)
 
 
diff --git a/mikeio/spatial/_geometry.py b/mikeio/spatial/_geometry.py
index 77af91e2f..05b47ad1e 100644
--- a/mikeio/spatial/_geometry.py
+++ b/mikeio/spatial/_geometry.py
@@ -1,11 +1,54 @@
 from abc import ABC, abstractmethod
 
-from collections import namedtuple
-from typing import Any
+from dataclasses import dataclass
+from typing import Any, Sequence
 
 from mikecore.Projections import MapProjection
 
-BoundingBox = namedtuple("BoundingBox", ["left", "bottom", "right", "top"])
+
+@dataclass
+class BoundingBox:
+    """Bounding box for spatial data."""
+
+    left: float
+    bottom: float
+    right: float
+    top: float
+
+    def __post_init__(self) -> None:
+        if self.left > self.right:
+            raise ValueError(
+                f"Invalid x axis, left: {self.left} must be smaller than right: {self.right}"
+            )
+
+        if self.bottom > self.top:
+            raise ValueError(
+                f"Invalid y axis, bottom: {self.bottom} must be smaller than top: {self.top}"
+            )
+
+    def overlaps(self, other: "BoundingBox") -> bool:
+        """Check if two bounding boxes overlap."""
+        return not (
+            self.left > other.right
+            or self.bottom > other.top
+            or self.right < other.left
+            or self.top < other.bottom
+        )
+
+    @staticmethod
+    def parse(
+        values: "BoundingBox" | Sequence[float],
+    ) -> "BoundingBox":
+        match values:
+            case BoundingBox():
+                bbox = values
+            case left, bottom, right, top:
+                bbox = BoundingBox(left, bottom, right, top)
+            case _:
+                raise ValueError(
+                    "values must be a bounding box of coordinates e.g. (-10.0, 10.0 20.0, 30.0)"
+                )
+        return bbox
 
 
 class _Geometry(ABC):
diff --git a/mikeio/spatial/_grid_geometry.py b/mikeio/spatial/_grid_geometry.py
index 2e9a8f124..19fa9ea13 100644
--- a/mikeio/spatial/_grid_geometry.py
+++ b/mikeio/spatial/_grid_geometry.py
@@ -433,7 +433,10 @@ def __init__(
         y0: float = 0.0,
         dy: float | None = None,
         ny: int | None = None,
-        bbox: tuple[float, float, float, float] | None = None,
+        bbox: BoundingBox
+        | Sequence[float]
+        | tuple[float, float, float, float]
+        | None = None,
         projection: str = "LONG/LAT",
         origin: tuple[float, float] | None = None,
         orientation: float = 0.0,
@@ -522,7 +525,7 @@ def _is_rotated(self) -> Any:
 
     def _create_in_bbox(
         self,
-        bbox: tuple[float, float, float, float],
+        bbox: BoundingBox | tuple[float, float, float, float] | Sequence[float],
         dx: float | tuple[float, float] | None = None,
         dy: float | None = None,
         nx: int | None = None,
@@ -549,10 +552,10 @@ def _create_in_bbox(
             in which case the value will be inferred
 
         """
-        left, bottom, right, top = self._parse_bbox(bbox)
+        box = BoundingBox.parse(bbox)
 
-        xr = right - left  # dx too large
-        yr = top - bottom  # dy too large
+        xr = box.right - box.left  # dx too large
+        yr = box.top - box.bottom  # dy too large
 
         if (dx is None and dy is None) and (nx is None and ny is None):
             if xr <= yr:
@@ -569,32 +572,12 @@ def _create_in_bbox(
             dx, dy = dx
         dy = dx if dy is None else dy
 
-        self._x0, self._dx, self._nx = self._create_in_bbox_1d("x", left, right, dx, nx)
-        self._y0, self._dy, self._ny = self._create_in_bbox_1d("y", bottom, top, dy, ny)
-
-    @staticmethod
-    def _parse_bbox(
-        bbox: tuple[float, float, float, float],
-    ) -> tuple[float, float, float, float]:
-        left = bbox[0]
-        bottom = bbox[1]
-        right = bbox[2]
-        top = bbox[3]
-
-        if left > right:
-            raise (
-                ValueError(
-                    f"Invalid x axis, left: {left} must be smaller than right: {right}"
-                )
-            )
-
-        if bottom > top:
-            raise (
-                ValueError(
-                    f"Invalid y axis, bottom: {bottom} must be smaller than top: {top}"
-                )
-            )
-        return left, bottom, right, top
+        self._x0, self._dx, self._nx = self._create_in_bbox_1d(
+            "x", box.left, box.right, dx, nx
+        )
+        self._y0, self._dy, self._ny = self._create_in_bbox_1d(
+            "y", box.bottom, box.top, dy, ny
+        )
 
     @staticmethod
     def _create_in_bbox_1d(
@@ -854,7 +837,8 @@ def find_index(
         if coords is not None:
             return self._xy_to_index(coords)
         elif area is not None:
-            return self._bbox_to_index(area)
+            bbox = BoundingBox.parse(area)
+            return self._bbox_to_index(bbox)
         else:
             raise ValueError("Provide x,y or coords")
 
@@ -877,22 +861,14 @@ def _xy_to_index(self, xy: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
 
         return ii, jj
 
-    def _bbox_to_index(
-        self, bbox: tuple[float, float, float, float] | BoundingBox
-    ) -> tuple[range, range]:
+    def _bbox_to_index(self, bbox: BoundingBox) -> tuple[range, range]:
         """Find subarea within this geometry."""
-        if not (len(bbox) == 4):
-            raise ValueError(
-                "area most be a bounding box of coordinates e.g. area=(-10.0, 10.0 20.0, 30.0)"
-            )
-
-        x0, y0, x1, y1 = bbox
-        if x0 > self.x[-1] or y0 > self.y[-1] or x1 < self.x[0] or y1 < self.y[0]:
+        if not bbox.overlaps(self.bbox):
             raise ValueError("area is outside grid")
 
-        mask = (self.x >= x0) & (self.x <= x1)
+        mask = (self.x >= bbox.left) & (self.x <= bbox.right)
         ii = np.where(mask)[0]
-        mask = (self.y >= y0) & (self.y <= y1)
+        mask = (self.y >= bbox.bottom) & (self.y <= bbox.top)
         jj = np.where(mask)[0]
 
         i = range(ii[0], ii[-1] + 1)