Repository: Deltares/xugrid (language: python)

Issue:
UgridDataArray seems to behave differently than DataArray with where
When doing operations with a limited set of labels belonging to a coordinate with xarray, xarray automatically returns a reduced DataArray back. 

This reproducible example explains it best, but unfortunately for me it doesn't work with xugrid 0.1.10 anymore.

```python
import xarray as xr
import xugrid as xu

da_active = xr.DataArray([True, True], coords={"layer" : [1,2]}, dims=("layer",))
da_data = xr.DataArray([1], coords={"layer": [2]}, dims=("layer",))

da_data.where(da_active)  # -> returns DataArray([1.]) with coords 'layer' = 2

uda_active = xu.UgridDataArray(da_active, grid=None)
uda_data = xu.UgridDataArray(da_data, grid=None)

uda_data.where(uda_active)  # -> Throws MergeError
```

To reproduce with 0.1.10, @Huite made the following example:

```python
import xarray as xr
import xugrid as xu
import numpy as np
# %%
uda = xu.data.disk()["face_z"]
layered = (uda * xr.DataArray([1, 1, 1], {"layer": [1, 2, 3]}, ["layer"])).transpose("layer", "mesh2d_nFaces")
# %%
uda_data = layered.isel(layer=slice(1, None))
# %%
new = layered.where(layered["layer"] > 1)
# %%
out = uda_data.obj.where(new.obj)
```


Discussion / hints:

Fix patch:
diff --git a/.github/workflows/requirements/environment.yml b/.github/workflows/requirements/environment.yml
index 7bf4900cd..908f62ebe 100644
--- a/.github/workflows/requirements/environment.yml
+++ b/.github/workflows/requirements/environment.yml
@@ -3,7 +3,7 @@ name: xugrid
 channels:
   - conda-forge
     
-# Installing geopandas, pygeos, pyproj via mamba seems to break CI currently.
+# Installing geopandas, pygeos, pyproj via pip seems to break CI currently.
 dependencies:
   - black
   - dask
diff --git a/docs/api.rst b/docs/api.rst
index 7af4645a3..1208e91fd 100644
--- a/docs/api.rst
+++ b/docs/api.rst
@@ -7,6 +7,9 @@ API Reference
 
 This page provides an auto-generated summary of xugrid's API.
 
+Xugrid also almost completely exposes the attributes and methods of xarray
+DataArrays and Datasets. Refer to the Xarray documentation.
+
 .. toctree::
    :maxdepth: 1
 
@@ -67,7 +70,6 @@ UgridDataArray or UgridDataset.
     UgridDatasetAccessor.crs
     UgridDatasetAccessor.set_crs
     UgridDatasetAccessor.to_crs
-    UgridDatasetAccessor.isel
     UgridDatasetAccessor.sel
     UgridDatasetAccessor.sel_points
     UgridDatasetAccessor.to_geodataframe
@@ -85,7 +87,6 @@ UgridDataArray or UgridDataset.
     UgridDataArrayAccessor.crs
     UgridDataArrayAccessor.set_crs
     UgridDataArrayAccessor.to_crs
-    UgridDataArrayAccessor.isel
     UgridDataArrayAccessor.sel
     UgridDataArrayAccessor.sel_points
     UgridDataArrayAccessor.rasterize
diff --git a/docs/changelog.rst b/docs/changelog.rst
index da15f26f3..514104870 100644
--- a/docs/changelog.rst
+++ b/docs/changelog.rst
@@ -21,12 +21,36 @@ Fixed
 
 Changed
 ~~~~~~~
+- Forwarding to the internal xarray object is now setup at class definition of
+  :class:`UgridDataArray` and :class:`UgridDataset` rather than at runtime.
+  This means tab completion and docstrings for the xarray methods should work.
+- The UGRID dimensions in :class:`UgridDataArray` and :class:`UgridDataset` are
+  labelled at initialization. This allows us to track necessary changes to the
+  UGRID topology for general xarray operations. Forwarded methods (such as
+  :meth:`UgridDataArray.isel`) will now create a subset topology if possible, or
+  error if an invalid topology is created by the selection.
+- This also means that selection on one facet of the grid (e.g. the face
+  dimension) will also result in a valid selection of the data on another facet
+  (such as the edge dimension).
+- :meth:`xugrid.Ugrid1d.sel` and :meth:`xugrid.Ugrid2d.sel` 
+- Consequently, :meth:`xugrid.UgridDataArrayAccessor.isel` and
+  :meth:`xugrid.UgridDatasetAccessor.isel` have been removed.
+- :property:`xugrid.Ugrid1d.dimensions` and
+  :property:`xugrid.Ugrid2d.dimensions` will now return a dictionary with the
+  keys the dimension names and as the values the sizes of the dimensions.
 
 Added
 ~~~~~
 - :class:`xugrid.Ugrid1d` and :class:`xugrid.Ugrid2d` can now be initialized
   with an ``attrs`` argument to setup non-default UGRID attributes such as
   alternative node, edge, or face dimensions.
+- :meth:`xugrid.Ugrid1d.topology_subset`,
+  :meth:`xugrid.Ugrid2d.topology_subset`, :meth:`xugrid.Ugrid1d.isel`, and
+  :meth:`xugrid.Ugrid2d.isel` now take a ``return_index`` argument and will
+  to return UGRID dimension indexes if set to True.
+- :meth:`xugrid.UgridDataArrayAccessor.clip_box` and
+  :meth:`xugrid.UgridDatasetAccessor.clip_box` have been added to more easily
+  select data in a bounding box.
 
 [0.1.10] 2022-12-13
 -------------------
diff --git a/examples-dev/voronoi.py b/examples-dev/voronoi.py
index d9c2813ff..1153dc333 100644
--- a/examples-dev/voronoi.py
+++ b/examples-dev/voronoi.py
@@ -28,7 +28,7 @@
 # From xugrid, we need only import the ``connectivity`` and ``voronoi``
 # modules. The functions in these modules depend only on ``numpy`` and
 # ``scipy.sparse``.
-from xugrid import connectivity, voronoi
+from xugrid.ugrid import connectivity, voronoi
 
 
 def generate_disk(partitions: int, depth: int):
diff --git a/xugrid/__init__.py b/xugrid/__init__.py
index 7e333456d..1fba496cb 100644
--- a/xugrid/__init__.py
+++ b/xugrid/__init__.py
@@ -1,13 +1,7 @@
 import pkg_resources
 
-from . import data
-from .plot import plot
-from .ugrid import Ugrid1d, Ugrid2d, UgridRolesAccessor
-from .ugrid_dataset import (
-    UgridDataArray,
-    UgridDataArrayAccessor,
-    UgridDataset,
-    UgridDatasetAccessor,
+from xugrid import data
+from xugrid.core.common import (
     concat,
     full_like,
     merge,
@@ -18,6 +12,13 @@
     open_zarr,
     zeros_like,
 )
+from xugrid.core.dataarray_accessor import UgridDataArrayAccessor
+from xugrid.core.dataset_accessor import UgridDatasetAccessor
+from xugrid.core.wrap import UgridDataArray, UgridDataset
+from xugrid.plot import plot
+from xugrid.ugrid.conventions import UgridRolesAccessor
+from xugrid.ugrid.ugrid1d import Ugrid1d
+from xugrid.ugrid.ugrid2d import Ugrid2d
 
 try:
     __version__ = pkg_resources.get_distribution(__name__).version
diff --git a/xugrid/constants.py b/xugrid/constants.py
new file mode 100644
index 000000000..19952b536
--- /dev/null
+++ b/xugrid/constants.py
@@ -0,0 +1,54 @@
+from typing import NamedTuple, Union
+
+import numpy as np
+from scipy.sparse import coo_matrix, csr_matrix
+
+# import numpy.typing as npt
+
+FloatDType = np.float64
+IntDType = np.intp
+
+# Requires numpy 1.21, not on conda yet...
+# FloatArray = np.ndarray[FloatDType]
+# IntArray = np.ndarray[IntDType]
+# BoolArray = np.ndarray[np.bool_]
+
+FloatArray = np.ndarray
+IntArray = np.ndarray
+BoolArray = np.ndarray
+# Pygeos collections:
+PointArray = np.ndarray
+LineArray = np.ndarray
+PolygonArray = np.ndarray
+SparseMatrix = Union[coo_matrix, csr_matrix]
+
+
+class Point(NamedTuple):
+    x: float
+    y: float
+
+
+class Vector(NamedTuple):
+    x: float
+    y: float
+
+
+# Spatial coordinate epsilon for floating point comparison
+# Assuming world coordinates in meters: 40 000 m along equator:
+# 40 000 000 = 4e7 mm
+# np.spacing(4e7) == 7.45E-9 ~= 1E-8
+X_EPSILON = 1.0e-8
+X_OFFSET = 1.0e-8
+T_OFFSET = 1.0e-6
+
+
+class MissingOptionalModule:
+    """
+    Presents a clear error for optional modules.
+    """
+
+    def __init__(self, name):
+        self.name = name
+
+    def __getattr__(self, name):
+        raise ImportError(f"{self.name} is required for this functionality")
diff --git a/xugrid/conversion.py b/xugrid/conversion.py
index 3857e9f8f..b34357f01 100644
--- a/xugrid/conversion.py
+++ b/xugrid/conversion.py
@@ -10,8 +10,7 @@
 import numpy as np
 import xarray as xr
 
-from .connectivity import ragged_index
-from .typing import (
+from xugrid.constants import (
     FloatArray,
     IntArray,
     IntDType,
@@ -20,6 +19,9 @@
     PointArray,
     PolygonArray,
 )
+from xugrid.ugrid.connectivity import ragged_index
+from xugrid.ugrid.ugrid1d import Ugrid1d
+from xugrid.ugrid.ugrid2d import Ugrid2d
 
 try:
     import pygeos
@@ -228,3 +230,43 @@ def bounds_to_vertices(bounds):
     else:
         raise ValueError("Bounds are not monotonic ascending or monotonic descending")
     return vertices
+
+
+def grid_from_geodataframe(geodataframe: "geopandas.GeoDataFrame"):  # type: ignore # noqa
+    import geopandas as gpd
+
+    gdf = geodataframe
+    if not isinstance(gdf, gpd.GeoDataFrame):
+        raise TypeError(
+            f"Cannot convert a {type(gdf).__name__}, expected a GeoDataFrame"
+        )
+
+    geom_types = gdf.geom_type.unique()
+    if len(geom_types) == 0:
+        raise ValueError("geodataframe contains no geometry")
+    elif len(geom_types) > 1:
+        message = ", ".join(geom_types)
+        raise ValueError(f"Multiple geometry types detected: {message}")
+
+    geom_type = geom_types[0]
+    if geom_type == "LineString":
+        grid = Ugrid1d.from_geodataframe(gdf)
+    elif geom_type == "Polygon":
+        grid = Ugrid2d.from_geodataframe(gdf)
+    else:
+        raise ValueError(
+            f"Invalid geometry type: {geom_type}. Expected Linestring or Polygon."
+        )
+    return grid
+
+
+def grid_from_dataset(dataset: xr.Dataset, topology: str):
+    topodim = dataset[topology].attrs["topology_dimension"]
+    if topodim == 1:
+        return Ugrid1d.from_dataset(dataset, topology)
+    elif topodim == 2:
+        return Ugrid2d.from_dataset(dataset, topology)
+    elif topodim == 3:
+        raise NotImplementedError
+    else:
+        raise ValueError(f"Invalid topology dimension: {topodim}")
diff --git a/xugrid/core/accessorbase.py b/xugrid/core/accessorbase.py
new file mode 100644
index 000000000..0958e18e2
--- /dev/null
+++ b/xugrid/core/accessorbase.py
@@ -0,0 +1,95 @@
+import abc
+
+
+class AbstractUgridAccessor(abc.ABC):
+    @abc.abstractmethod
+    def to_dataset():
+        """ """
+
+    @abc.abstractmethod
+    def assign_node_coords():
+        """ """
+
+    @abc.abstractmethod
+    def set_node_coords():
+        """ """
+
+    @abc.abstractproperty
+    def crs():
+        """ """
+
+    @abc.abstractmethod
+    def set_crs():
+        """ """
+
+    @abc.abstractmethod
+    def to_crs():
+        """ """
+
+    @abc.abstractproperty
+    def bounds():
+        """ """
+
+    @abc.abstractproperty
+    def total_bounds():
+        """ """
+
+    @staticmethod
+    def _sel(obj, grid, x, y):
+        from xugrid.core.wrap import maybe_xugrid
+
+        # TODO: also do vectorized indexing like xarray?
+        # Might not be worth it, as orthogonal and vectorized indexing are
+        # quite confusing.
+        result = grid.sel(obj, x, y)
+        if isinstance(result, tuple):
+            return maybe_xugrid(*result)
+        else:
+            return result
+
+    def clip_box(
+        self,
+        xmin: float,
+        ymin: float,
+        xmax: float,
+        ymax: float,
+    ):
+        """
+        Clip the DataArray or Dataset by a bounding box.
+
+        Parameters
+        ----------
+        xmin: float
+        ymin: float
+        xmax: float
+        ymax: float
+
+        -------
+        clipped:
+            xugrid.UgridDataArray or xugrid.UgridDataset
+        """
+        return self.sel(x=slice(xmin, xmax), y=slice(ymin, ymax))
+
+    def to_netcdf(self, *args, **kwargs):
+        """
+        Write dataset contents to a UGRID compliant netCDF file.
+
+        This function wraps :py:meth:`xr.Dataset.to_netcdf`; it adds the UGRID
+        variables and coordinates to a standard xarray Dataset, then writes the
+        result to a netCDF.
+
+        All arguments are forwarded to :py:meth:`xr.Dataset.to_netcdf`.
+        """
+        self.to_dataset().to_netcdf(*args, **kwargs)
+
+    def to_zarr(self, *args, **kwargs):
+        """
+        Write dataset contents to a UGRID compliant Zarr file.
+
+        This function wraps :py:meth:`xr.Dataset.to_zarr`; it adds the UGRID
+        variables and coordinates to a standard xarray Dataset, then writes the
+        result to a Zarr file.
+
+        All arguments are forwarded to :py:meth:`xr.Dataset.to_zarr`.
+        """
+        self.to_dataset().to_zarr(*args, **kwargs)
diff --git a/xugrid/core/common.py b/xugrid/core/common.py
new file mode 100644
index 000000000..33850b14b
--- /dev/null
+++ b/xugrid/core/common.py
@@ -0,0 +1,117 @@
+from functools import wraps
+
+import xarray as xr
+from xarray.backends.api import DATAARRAY_NAME, DATAARRAY_VARIABLE
+
+from xugrid.core.wrap import UgridDataArray, UgridDataset
+
+
+def open_dataset(*args, **kwargs):
+    ds = xr.open_dataset(*args, **kwargs)
+    return UgridDataset(ds)
+
+
+def open_dataarray(*args, **kwargs):
+    ds = xr.open_dataset(*args, **kwargs)
+    dataset = UgridDataset(ds)
+
+    if len(dataset.data_vars) != 1:
+        raise ValueError(
+            "Given file dataset contains more than one data "
+            "variable. Please read with xarray.open_dataset and "
+            "then select the variable you want."
+        )
+    else:
+        (data_array,) = dataset.data_vars.values()
+
+    data_array.set_close(dataset._close)
+
+    # Reset names if they were changed during saving
+    # to ensure that we can 'roundtrip' perfectly
+    if DATAARRAY_NAME in dataset.attrs:
+        data_array.name = dataset.attrs[DATAARRAY_NAME]
+        del dataset.attrs[DATAARRAY_NAME]
+
+    if data_array.name == DATAARRAY_VARIABLE:
+        data_array.name = None
+
+    return data_array
+
+
+def open_mfdataset(*args, **kwargs):
+    if "data_vars" in kwargs:
+        raise ValueError("data_vars kwargs is not supported in xugrid.open_mfdataset")
+    kwargs["data_vars"] = "minimal"
+    ds = xr.open_mfdataset(*args, **kwargs)
+    return UgridDataset(ds)
+
+
+def open_zarr(*args, **kwargs):
+    ds = xr.open_zarr(*args, **kwargs)
+    return UgridDataset(ds)
+
+
+open_dataset.__doc__ = xr.open_dataset.__doc__
+open_dataarray.__doc__ = xr.open_dataarray.__doc__
+open_mfdataset.__doc__ = xr.open_mfdataset.__doc__
+open_zarr.__doc__ = xr.open_zarr.__doc__
+
+
+# Other utilities
+# ---------------
+
+
+def wrap_func_like(func):
+    @wraps(func)
+    def _like(other, *args, **kwargs):
+        obj = func(other.obj, *args, **kwargs)
+        if isinstance(obj, xr.DataArray):
+            return type(other)(obj, other.grid)
+        elif isinstance(obj, xr.Dataset):
+            return type(other)(obj, other.grids)
+        else:
+            raise TypeError(
+                f"Expected Dataset or DataArray, received {type(other).__name__}"
+            )
+
+    _like.__doc__ = func.__doc__
+    return _like
+
+
+def wrap_func_objects(func):
+    @wraps(func)
+    def _f(objects, *args, **kwargs):
+        grids = []
+        bare_objs = []
+        for obj in objects:
+            if isinstance(obj, UgridDataArray):
+                grids.append(obj.grid)
+            elif isinstance(obj, UgridDataset):
+                grids.extend(obj.grids)
+            else:
+                raise TypeError(
+                    "Can only concatenate xugrid UgridDataset and UgridDataArray "
+                    f"objects, got {type(obj).__name__}"
+                )
+
+            bare_objs.append(obj.obj)
+
+        grids = set(grids)
+        result = func(bare_objs, *args, **kwargs)
+        if isinstance(result, xr.DataArray):
+            if len(grids) > 1:
+                raise ValueError("All UgridDataArrays must have the same grid")
+            return UgridDataArray(result, next(iter(grids)))
+        else:
+            return UgridDataset(result, grids)
+
+    _f.__doc__ = func.__doc__
+    return _f
+
+
+full_like = wrap_func_like(xr.full_like)
+zeros_like = wrap_func_like(xr.zeros_like)
+ones_like = wrap_func_like(xr.ones_like)
+
+concat = wrap_func_objects(xr.concat)
+merge = wrap_func_objects(xr.merge)
diff --git a/xugrid/core/dataarray_accessor.py b/xugrid/core/dataarray_accessor.py
new file mode 100644
index 000000000..d1b325d02
--- /dev/null
+++ b/xugrid/core/dataarray_accessor.py
@@ -0,0 +1,521 @@
+from typing import Dict, Tuple, Union
+
+import numpy as np
+import scipy.sparse
+import xarray as xr
+from xarray.core.utils import UncachedAccessor
+
+# from xugrid.plot.pyvista import to_pyvista_grid
+from xugrid.core.accessorbase import AbstractUgridAccessor
+from xugrid.core.wrap import UgridDataArray
+from xugrid.plot.plot import _PlotMethods
+from xugrid.ugrid import connectivity
+from xugrid.ugrid.interpolate import laplace_interpolate
+from xugrid.ugrid.ugridbase import UgridType
+
+
+class UgridDataArrayAccessor(AbstractUgridAccessor):
+    """
+    This "accessor" makes operations using the UGRID topology available via the
+    ``.ugrid`` attribute for UgridDataArrays and UgridDatasets.
+    """
+
+    def __init__(self, obj: xr.DataArray, grid: UgridType):
+        self.obj = obj
+        self.grid = grid
+
+    @property
+    def bounds(self) -> Dict[str, Tuple]:
+        """
+        Mapping from grid name to tuple containing ``minx, miny, maxx, maxy``
+        values of the grid's node coordinates.
+        """
+        return {self.grid.name: self.grid.bounds}
+
+    @property
+    def total_bounds(self) -> Tuple:
+        """
+        Returns a tuple containing ``minx, miny, maxx, maxy`` values of the grid's
+        node coordinates.
+        """
+        return next(iter(self.bounds.values()))
+
+    plot = UncachedAccessor(_PlotMethods)
+
+    def assign_node_coords(self) -> UgridDataArray:
+        """
+        Assign node coordinates from the grid to the object.
+
+        Returns a new object with all the original data in addition to the new
+        node coordinates of the grid.
+
+        Returns
+        -------
+        assigned: UgridDataset
+        """
+        return UgridDataArray(self.grid.assign_node_coords(self.obj), self.grid)
+
+    def assign_edge_coords(self) -> UgridDataArray:
+        """
+        Assign edge coordinates from the grid to the object.
+
+        Returns a new object with all the original data in addition to the new
+        node coordinates of the grid.
+
+        Returns
+        -------
+        assigned: UgridDataset
+        """
+        return UgridDataArray(self.grid.assign_edge_coords(self.obj), self.grid)
+
+    def assign_face_coords(self) -> UgridDataArray:
+        """
+        Assign face coordinates from the grid to the object.
+
+        Returns a new object with all the original data in addition to the new
+        node coordinates of the grid.
+
+        Returns
+        -------
+        assigned: UgridDataset
+        """
+        if self.grid.topology_dimension == 1:
+            raise TypeError("Cannot set face coords from a Ugrid1D topology")
+        return UgridDataArray(self.grid.assign_face_coords(self.obj), self.grid)
+
+    def set_node_coords(self, node_x: str, node_y: str):
+        """
+        Given names of x and y coordinates of the nodes of an object, set them
+        as the coordinates in the grid.
+
+        Parameters
+        ----------
+        node_x: str
+            Name of the x coordinate of the nodes in the object.
+        node_y: str
+            Name of the y coordinate of the nodes in the object.
+        """
+        self.grid.set_node_coords(node_x, node_y, self.obj)
+
+    def sel(self, x=None, y=None):
+        """
+        Returns a new object, a subselection in the UGRID x and y coordinates.
+
+        The indexing for x and y always occurs orthogonally, i.e.:
+        ``.sel(x=[0.0, 5.0], y=[10.0, 15.0])`` results in a four points. For
+        vectorized indexing (equal to ``zip``ing through x and y), see
+        ``.sel_points``.
+
+        Depending on the nature of the x and y indexers, a xugrid or xarray
+        object is returned:
+
+        * slice without step: ``x=slice(-100, 100)``: returns xugrid object, a
+          part of the unstructured grid.
+        * slice with step: ``x=slice(-100, 100, 10)``: returns xarray object, a
+          series of points (x=[-100, -90, -80, ..., 90, 100]).
+        * a scalar: ``x=5.0``: returns xarray object, a point.
+        * an array: ``x=[1.0, 15.0, 17.0]``: returns xarray object, a series of
+          points.
+
+        Parameters
+        ----------
+        x: float, 1d array, slice
+        y: float, 1d array, slice
+
+        Returns
+        -------
+        selection: Union[UgridDataArray, UgridDataset, xr.DataArray, xr.Dataset]
+        """
+        return self._sel(self.obj, self.grid, x, y)
+
+    def sel_points(self, x, y):
+        """
+        Select points in the unstructured grid.
+
+        Parameters
+        ----------
+        x: ndarray of floats with shape ``(n_points,)``
+        y: ndarray of floats with shape ``(n_points,)``
+
+        Returns
+        -------
+        points: Union[xr.DataArray, xr.Dataset]
+        """
+        return self.grid.sel_points(self.obj, x, y)
+
+    def _raster(self, x, y, index) -> xr.DataArray:
+        index = index.ravel()
+        # Cast to float for nodata values (NaN)
+        data = self.obj.isel({self.grid.face_dimension: index}).astype(float).values
+        data[index == -1] = np.nan
+        out = xr.DataArray(
+            data=data.reshape(y.size, x.size),
+            coords={"y": y, "x": x},
+            dims=["y", "x"],
+        )
+        return out
+
+    def rasterize(self, resolution: float):
+        """
+        Rasterize unstructured grid by sampling.
+
+        Parameters
+        ----------
+        resolution: float
+            Spacing in x and y.
+
+        Returns
+        -------
+        rasterized: Union[xr.DataArray, xr.Dataset]
+        """
+        x, y, index = self.grid.rasterize(resolution)
+        return self._raster(x, y, index)
+
+    def rasterize_like(self, other: Union[xr.DataArray, xr.Dataset]):
+        """
+        Rasterize unstructured grid by sampling on the x and y coordinates
+        of ``other``.
+
+        Parameters
+        ----------
+        resolution: float
+            Spacing in x and y.
+        other: Union[xr.DataArray, xr.Dataset]
+            Object to take x and y coordinates from.
+
+        Returns
+        -------
+        rasterized: Union[xr.DataArray, xr.Dataset]
+        """
+        x, y, index = self.grid.rasterize_like(
+            x=other["x"].values,
+            y=other["y"].values,
+        )
+        return self._raster(x, y, index)
+
+    @property
+    def crs(self):
+        """
+        The Coordinate Reference System (CRS) represented as a ``pyproj.CRS`` object.
+
+        Returns None if the CRS is not set.
+
+        Returns
+        -------
+        crs: dict
+            A dictionary containing the name of the grid and its CRS.
+        """
+        return {self.grid.name: self.grid.crs}
+
+    def set_crs(
+        self,
+        crs: Union["pyproj.CRS", str] = None,  # type: ignore # noqa
+        epsg: int = None,
+        allow_override: bool = False,
+    ):
+        """
+        Set the Coordinate Reference System (CRS) of a UGRID topology.
+
+        NOTE: The underlying geometries are not transformed to this CRS. To
+        transform the geometries to a new CRS, use the ``to_crs`` method.
+
+        Parameters
+        ----------
+        crs : pyproj.CRS, optional if `epsg` is specified
+            The value can be anything accepted
+            by :meth:`pyproj.CRS.from_user_input() <pyproj.crs.CRS.from_user_input>`,
+            such as an authority string (eg "EPSG:4326") or a WKT string.
+        epsg : int, optional if `crs` is specified
+            EPSG code specifying the projection.
+        allow_override : bool, default False
+            If the the UGRID topology already has a CRS, allow to replace the
+            existing CRS, even when both are not equal.
+        """
+        self.grid.set_crs(crs, epsg, allow_override)
+
+    def to_crs(
+        self,
+        crs: Union["pyproj.CRS", str] = None,  # type: ignore # noqa
+        epsg: int = None,
+    ) -> UgridDataArray:
+        """
+        Transform geometries to a new coordinate reference system.
+        Transform all geometries in an active geometry column to a different coordinate
+        reference system. The ``crs`` attribute on the current Ugrid must
+        be set. Either ``crs`` or ``epsg`` may be specified for output.
+
+        This method will transform all points in all objects. It has no notion
+        of projecting the cells. All segments joining points are assumed to be
+        lines in the current projection, not geodesics. Objects crossing the
+        dateline (or other projection boundary) will have undesirable behavior.
+
+        Parameters
+        ----------
+        crs : pyproj.CRS, optional if `epsg` is specified
+            The value can be anything accepted by
+            :meth:`pyproj.CRS.from_user_input() <pyproj.crs.CRS.from_user_input>`,
+            such as an authority string (eg "EPSG:4326") or a WKT string.
+        epsg : int, optional if `crs` is specified
+            EPSG code specifying output projection.
+        """
+        uda = UgridDataArray(self.obj, self.grid.to_crs(crs, epsg))
+        if self.grid.node_dimension in self.obj.dims:
+            return uda.ugrid.assign_node_coords()
+        else:
+            return uda
+
+    def to_geodataframe(
+        self, name: str = None, dim_order=None
+    ) -> "geopandas.GeoDataFrame":  # type: ignore # noqa
+        """
+        Convert data and topology of one facet (node, edge, face) of the grid
+        to a geopandas GeoDataFrame. This also determines the geometry type of
+        the geodataframe:
+
+        * node: point
+        * edge: line
+        * face: polygon
+
+        Parameters
+        ----------
+        dim: str
+            node, edge, or face dimension. Inferred for DataArray.
+        name: str
+            Name to give to the array (required if unnamed).
+        dim_order:
+            Hierarchical dimension order for the resulting dataframe. Array content is
+            transposed to this order and then written out as flat vectors in contiguous
+            order, so the last dimension in this list will be contiguous in the resulting
+            DataFrame. This has a major influence on which operations are efficient on the
+            resulting dataframe.
+
+            If provided, must include all dimensions of this DataArray. By default,
+            dimensions are sorted according to the DataArray dimensions order.
+
+        Returns
+        -------
+        geodataframe: gpd.GeoDataFrame
+        """
+        import geopandas as gpd
+
+        dim = self.obj.dims[-1]
+        if name is not None:
+            ds = self.obj.to_dataset(name=name)
+        else:
+            ds = self.obj.to_dataset()
+
+        variables = [var for var in ds.data_vars if dim in ds[var].dims]
+        # TODO deal with time-dependent data, etc.
+        # Basically requires checking which variables are static, which aren't.
+        # For non-static, requires repeating all geometries.
+        # Call reset_index on multi-index to generate them as regular columns.
+        df = ds[variables].to_dataframe(dim_order=dim_order)
+        geometry = self.grid.to_pygeos(dim)
+        return gpd.GeoDataFrame(df, geometry=geometry)
+
+    def _binary_iterate(self, iterations: int, mask, value, border_value):
+        if border_value == value:
+            exterior = self.grid.exterior_faces
+        else:
+            exterior = None
+        if mask is not None:
+            mask = mask.values
+
+        obj = self.obj
+        if isinstance(obj, xr.DataArray):
+            output = connectivity._binary_iterate(
+                self.grid.face_face_connectivity,
+                obj.values,
+                value,
+                iterations,
+                mask,
+                exterior,
+                border_value,
+            )
+            da = obj.copy(data=output)
+            return UgridDataArray(da, self.grid.copy())
+        elif isinstance(obj, xr.Dataset):
+            raise NotImplementedError
+        else:
+            raise ValueError("object should be a xr.DataArray")
+
+    def binary_dilation(
+        self,
+        iterations: int = 1,
+        mask=None,
+        border_value=False,
+    ):
+        """
+        Binary dilation can be used on a boolean array to expand the "shape" of
+        features.
+
+        Compare with :py:func:`scipy.ndimage.binary_dilation`.
+
+        Parameters
+        ----------
+        iterations: int, default: 1
+        mask: 1d array of bool, optional
+        border_value: bool, default value: False
+
+        Returns
+        -------
+        dilated: UgridDataArray
+        """
+        return self._binary_iterate(iterations, mask, True, border_value)
+
+    def binary_erosion(
+        self,
+        iterations: int = 1,
+        mask=None,
+        border_value=False,
+    ):
+        """
+        Binary erosion can be used on a boolean array to shrink the "shape" of
+        features.
+
+        Compare with :py:func:`scipy.ndimage.binary_erosion`.
+
+        Parameters
+        ----------
+        iterations: int, default: 1
+        mask: 1d array of bool, optional
+        border_value: bool, default value: False
+
+        Returns
+        -------
+        eroded: UgridDataArray
+        """
+        return self._binary_iterate(iterations, mask, False, border_value)
+
+    def connected_components(self):
+        """
+        Every edge or face is given a component number. If all are connected,
+        all will have the same number.
+
+        Wraps :py:func:`scipy.sparse.csgraph.connected_components``.
+
+        Returns
+        -------
+        labelled: UgridDataArray
+        """
+        _, labels = scipy.sparse.csgraph.connected_components(
+            self.grid.face_face_connectivity
+        )
+        return UgridDataArray(
+            xr.DataArray(labels, dims=[self.grid.face_dimension]),
+            self.grid,
+        )
+
+    def reverse_cuthill_mckee(self):
+        """
+        Reduces bandwith of the connectivity matrix.
+
+        Wraps :py:func:`scipy.sparse.csgraph.reverse_cuthill_mckee`.
+
+        Returns
+        -------
+        reordered: Union[UgridDataArray, UgridDataset]
+        """
+        grid = self.grid
+        reordered_grid, reordering = self.grid.reverse_cuthill_mckee()
+        reordered_data = self.obj.isel({grid.face_dimension: reordering})
+        # TODO: this might not work properly if e.g. centroids are stored in obj.
+        # Not all metadata would be reordered.
+        return UgridDataArray(
+            reordered_data,
+            reordered_grid,
+        )
+
+    def laplace_interpolate(
+        self,
+        xy_weights: bool = True,
+        direct_solve: bool = False,
+        drop_tol: float = None,
+        fill_factor: float = None,
+        drop_rule: str = None,
+        options: dict = None,
+        tol: float = 1.0e-5,
+        maxiter: int = 250,
+    ):
+        """
+        Fill gaps in ``data`` (``np.nan`` values) using Laplace interpolation.
+
+        This solves Laplace's equation where where there is no data, with data
+        values functioning as fixed potential boundary conditions.
+
+        Note that an iterative solver method will be required for large grids.
+        In this case, some experimentation with the solver settings may be
+        required to find a converging solution of sufficient accuracy. Refer to
+        the documentation of :py:func:`scipy.sparse.linalg.spilu` and
+        :py:func:`scipy.sparse.linalg.cg`.
+
+        Parameters
+        ----------
+        xy_weights: bool, default False.
+            Wether to use the inverse of the centroid to centroid distance in
+            the coefficient matrix. If ``False``, defaults to uniform
+            coefficients of 1 so that each face connection has equal weight.
+        direct_solve: bool, optional, default ``False``
+            Whether to use a direct or an iterative solver or a conjugate gradient
+            solver. Direct method provides an exact answer, but are unsuitable
+            for large problems.
+        drop_tol: float, optional, default None.
+            Drop tolerance for ``scipy.sparse.linalg.spilu`` which functions as a
+            preconditioner for the conjugate gradient solver.
+        fill_factor: float, optional, default None.
+            Fill factor for ``scipy.sparse.linalg.spilu``.
+        drop_rule: str, optional default None.
+            Drop rule for ``scipy.sparse.linalg.spilu``.
+        options: dict, optional, default None.
+            Remaining other options for ``scipy.sparse.linalg.spilu``.
+        tol: float, optional, default 1.0e-5.
+            Convergence tolerance for ``scipy.sparse.linalg.cg``.
+        maxiter: int, default 250.
+            Maximum number of iterations for ``scipy.sparse.linalg.cg``.
+
+        Returns
+        -------
+        filled: UgridDataArray of floats
+        """
+        grid = self.grid
+        da = self.obj
+        if grid.topology_dimension != 2:
+            raise NotImplementedError
+        if len(da.dims) > 1 or da.dims[0] != grid.face_dimension:
+            raise NotImplementedError
+
+        connectivity = grid.face_face_connectivity.copy()
+        if xy_weights:
+            xy = grid.centroids
+            coo = connectivity.tocoo()
+            i = coo.row
+            j = coo.col
+            connectivity.data = 1.0 / np.linalg.norm(xy[j] - xy[i], axis=1)
+
+        filled = laplace_interpolate(
+            connectivity=connectivity,
+            data=da.values,
+            use_weights=xy_weights,
+            direct_solve=direct_solve,
+            drop_tol=drop_tol,
+            fill_factor=fill_factor,
+            drop_rule=drop_rule,
+            options=options,
+            tol=tol,
+            maxiter=maxiter,
+        )
+        da_filled = da.copy(data=filled)
+        return UgridDataArray(da_filled, grid)
+
+    def to_dataset(self):
+        """
+        Converts this UgridDataArray or UgridDataset into a standard
+        xarray.Dataset.
+
+        The UGRID topology information is added as standard data variables.
+
+        Returns
+        -------
+        dataset: UgridDataset
+        """
+        return self.grid.to_dataset(self.obj)
diff --git a/xugrid/core/dataset_accessor.py b/xugrid/core/dataset_accessor.py
new file mode 100644
index 000000000..bd41f44df
--- /dev/null
+++ b/xugrid/core/dataset_accessor.py
@@ -0,0 +1,323 @@
+from typing import Dict, Sequence, Tuple, Union
+
+import numpy as np
+import pandas as pd
+import xarray as xr
+
+# from xugrid.plot.pyvista import to_pyvista_grid
+from xugrid.core.accessorbase import AbstractUgridAccessor
+from xugrid.core.wrap import UgridDataset
+from xugrid.ugrid.ugridbase import UgridType
+
+
+class UgridDatasetAccessor(AbstractUgridAccessor):
+    def __init__(self, obj: xr.Dataset, grids: Sequence[UgridType]):
+        self.obj = obj
+        self.grids = grids
+
+    @property
+    def grid(self) -> UgridType:
+        ngrid = len(self.grids)
+        if ngrid == 1:
+            return self.grids[0]
+        else:
+            raise AttributeError(
+                "Can only access grid topology via `.grid` if dataset contains "
+                f"exactly one grid. Dataset contains {ngrid} grids. Use "
+                "`.grids` instead."
+            )
+
+    @property
+    def bounds(self) -> Dict[str, Tuple]:
+        """
+        Mapping from grid name to tuple containing ``minx, miny, maxx, maxy``
+        values of the grid's node coordinates for every grid in the dataset.
+        """
+        return {grid.name: grid.bounds for grid in self.grids}
+
+    @property
+    def total_bounds(self) -> Tuple:
+        """
+        Returns a tuple containing ``minx, miny, maxx, maxy`` values for the
+        bounds of the dataset as a whole. Currently does not check whether the
+        coordinate reference systems (CRS) of the grids in the dataset match.
+        """
+        bounds = np.column_stack([bound for bound in self.bounds.values()])
+        return (
+            bounds[0].min(),
+            bounds[1].min(),
+            bounds[2].max(),
+            bounds[3].max(),
+        )
+
+    def assign_node_coords(self) -> UgridDataset:
+        """
+        Assign node coordinates from the grid to the object.
+
+        Returns a new object with all the original data in addition to the new
+        node coordinates of the grid.
+
+        Returns
+        -------
+        assigned: UgridDataset
+        """
+        result = self.obj
+        for grid in self.grids:
+            result = grid.assign_node_coords(result)
+        return UgridDataset(result, self.grids)
+
+    def assign_edge_coords(self) -> UgridDataset:
+        """
+        Assign edge coordinates from the grid to the object.
+
+        Returns a new object with all the original data in addition to the new
+        node coordinates of the grid.
+
+        Returns
+        -------
+        assigned: UgridDataset
+        """
+        result = self.obj
+        for grid in self.grids:
+            result = grid.assign_edge_coords(result)
+        return UgridDataset(result, self.grids)
+
+    def assign_face_coords(self) -> UgridDataset:
+        """
+        Assign face coordinates from the grid to the object.
+
+        Returns a new object with all the original data in addition to the new
+        node coordinates of the grid.
+
+        Returns
+        -------
+        assigned: UgridDataset
+        """
+        result = self.obj
+        for grid in self.grids:
+            if grid.topology_dimension > 1:
+                result = grid.assign_face_coords(result)
+        return UgridDataset(result, self.grids)
+
+    def set_node_coords(self, node_x: str, node_y: str, topology: str = None):
+        """
+        Given names of x and y coordinates of the nodes of an object, set them
+        as the coordinates in the grid.
+
+        Parameters
+        ----------
+        node_x: str
+            Name of the x coordinate of the nodes in the object.
+        node_y: str
+            Name of the y coordinate of the nodes in the object.
+        topology: str, optional
+            Name of the grid topology in which to set the node_x and node_y
+            coordinates. Can be omitted if the UgridDataset contains only a
+            single grid.
+        """
+        if topology is None:
+            if len(self.grids) == 1:
+                grid = self.grids[0]
+            else:
+                raise ValueError(
+                    "topology must be specified when dataset contains multiple datasets"
+                )
+        else:
+            grids = [grid for grid in self.grid if grid.name == topology]
+            if len(grids != 1):
+                raise ValueError(
+                    f"Expected one grid with topology {topology}, found: {len(grids)}"
+                )
+            grid = grids[0]
+
+        grid.set_node_coords(node_x, node_y, self.obj)
+
+    def sel(self, x=None, y=None):
+        result = self.obj
+        for grid in self.grids:
+            result = self._sel(result, grid, x, y)
+        return result
+
+    def sel_points(self, x, y):
+        """
+        Select points in the unstructured grid.
+
+        Parameters
+        ----------
+        x: ndarray of floats with shape ``(n_points,)``
+        y: ndarray of floats with shape ``(n_points,)``
+
+        Returns
+        -------
+        points: Union[xr.DataArray, xr.Dataset]
+        """
+        result = self.obj
+        for grid in self.grids:
+            result = grid.sel_points(result, x, y)
+        return result
+
+    def to_dataset(self):
+        """
+        Converts this UgridDataArray or UgridDataset into a standard
+        xarray.Dataset.
+
+        The UGRID topology information is added as standard data variables.
+
+        Returns
+        -------
+        dataset: UgridDataset
+        """
+        return xr.merge([grid.to_dataset(self.obj) for grid in self.grids])
+
+    @property
+    def crs(self):
+        """
+        The Coordinate Reference System (CRS) represented as a ``pyproj.CRS`` object.
+
+        Returns None if the CRS is not set.
+
+        Returns
+        -------
+        crs: dict
+            A dictionary containing the names of the grids and their CRS.
+        """
+        return {grid.name: grid.crs for grid in self.grids}
+
+    def set_crs(
+        self,
+        crs: Union["pyproj.CRS", str] = None,  # type: ignore # noqa
+        epsg: int = None,
+        allow_override: bool = False,
+        topology: str = None,
+    ):
+        """
+        Set the Coordinate Reference System (CRS) of a UGRID topology.
+
+        NOTE: The underlying geometries are not transformed to this CRS. To
+        transform the geometries to a new CRS, use the ``to_crs`` method.
+
+        Parameters
+        ----------
+        crs : pyproj.CRS, optional if `epsg` is specified
+            The value can be anything accepted
+            by :meth:`pyproj.CRS.from_user_input() <pyproj.crs.CRS.from_user_input>`,
+            such as an authority string (eg "EPSG:4326") or a WKT string.
+        epsg : int, optional if `crs` is specified
+            EPSG code specifying the projection.
+        allow_override : bool, default False
+            If the the UGRID topology already has a CRS, allow to replace the
+            existing CRS, even when both are not equal.
+        topology: str, optional
+            Name of the grid topology in which to set the CRS.
+            Sets the CRS for all grids if left unspecified.
+        """
+        if topology is None:
+            grids = self.grids
+        else:
+            names = [grid.name for grid in self.grids]
+            if topology not in names:
+                raise ValueError(f"{topology} not found. Expected one of: {names}")
+            grids = [grid for grid in self.grids if grid.name == topology]
+
+        for grid in grids:
+            grid.set_crs(crs, epsg, allow_override)
+
+    def to_crs(
+        self,
+        crs: Union["pyproj.CRS", str] = None,  # type: ignore # noqa
+        epsg: int = None,
+        topology: str = None,
+    ) -> UgridDataset:
+        """
+        Transform geometries to a new coordinate reference system.
+        Transform all geometries in an active geometry column to a different coordinate
+        reference system. The ``crs`` attribute on the current Ugrid must
+        be set. Either ``crs`` or ``epsg`` may be specified for output.
+
+        This method will transform all points in all objects. It has no notion
+        of projecting the cells. All segments joining points are assumed to be
+        lines in the current projection, not geodesics. Objects crossing the
+        dateline (or other projection boundary) will have undesirable behavior.
+
+        Parameters
+        ----------
+        crs : pyproj.CRS, optional if `epsg` is specified
+            The value can be anything accepted by
+            :meth:`pyproj.CRS.from_user_input() <pyproj.crs.CRS.from_user_input>`,
+            such as an authority string (eg "EPSG:4326") or a WKT string.
+        epsg : int, optional if `crs` is specified
+            EPSG code specifying output projection.
+        topology: str, optional
+            Name of the grid topology to reproject.
+            Reprojects all grids if left unspecified.
+        """
+        if topology is None:
+            grids = [grid.to_crs(crs, epsg) for grid in self.grids]
+        else:
+            names = [grid.name for grid in self.grids]
+            if topology not in names:
+                raise ValueError(f"{topology} not found. Expected one of: {names}")
+            grids = [
+                grid.to_crs(crs, epsg) if grid.name == topology else grid.copy()
+                for grid in self.grids
+            ]
+
+        uds = UgridDataset(self.obj, grids)
+        return uds.ugrid.assign_node_coords()
+
+    def to_geodataframe(
+        self, dim_order=None
+    ) -> "geopandas.GeoDataFrame":  # type: ignore # noqa
+        """
+        Convert data and topology of one facet (node, edge, face) of the grid
+        to a geopandas GeoDataFrame. This also determines the geometry type of
+        the geodataframe:
+
+        * node: point
+        * edge: line
+        * face: polygon
+
+        Parameters
+        ----------
+        name: str
+            Name to give to the array (required if unnamed).
+        dim_order:
+            Hierarchical dimension order for the resulting dataframe. Array content is
+            transposed to this order and then written out as flat vectors in contiguous
+            order, so the last dimension in this list will be contiguous in the resulting
+            DataFrame. This has a major influence on which operations are efficient on the
+            resulting dataframe.
+
+            If provided, must include all dimensions of this DataArray. By default,
+            dimensions are sorted according to the DataArray dimensions order.
+
+        Returns
+        -------
+        geodataframe: gpd.GeoDataFrame
+        """
+        import geopandas as gpd
+
+        ds = self.obj
+
+        gdfs = []
+        for grid in self.grids:
+            if grid.topology_dimension == 1:
+                dim = grid.edge_dimension
+            elif grid.topology_dimension == 2:
+                dim = grid.face_dimension
+            else:
+                raise ValueError("invalid topology dimension on grid")
+
+            variables = [var for var in ds.data_vars if dim in ds[var].dims]
+            # TODO deal with time-dependent data, etc.
+            # Basically requires checking which variables are static, which aren't.
+            # For non-static, requires repeating all geometries.
+            # Call reset_index on multi-index to generate them as regular columns.
+            gdfs.append(
+                gpd.GeoDataFrame(
+                    data=ds[variables].to_dataframe(dim_order=dim_order),
+                    geometry=grid.to_pygeos(dim),
+                )
+            )
+
+        return pd.concat(gdfs)
diff --git a/xugrid/core/wrap.py b/xugrid/core/wrap.py
new file mode 100644
index 000000000..b845dc9b9
--- /dev/null
+++ b/xugrid/core/wrap.py
@@ -0,0 +1,366 @@
+"""
+Wrap in advance instead of overloading __getattr__.
+
+This allows for tab completion and documentation.
+"""
+
+import types
+from collections import ChainMap
+from functools import wraps
+from itertools import chain
+from typing import Sequence, Union
+
+import xarray as xr
+from pandas import RangeIndex
+
+import xugrid
+from xugrid.conversion import grid_from_dataset, grid_from_geodataframe
+from xugrid.ugrid.ugrid2d import Ugrid2d
+from xugrid.ugrid.ugridbase import AbstractUgrid, UgridType, align
+
+# Import entire module here for circular import of UgridDatasetAccessor and
+# UgridDataArrayAccessor. Note: can only be used in functions (since that code
+# is run at runtime).
+
+
+def maybe_xugrid(obj, topology, old_indexes=None):
+    if not isinstance(obj, (xr.DataArray, xr.Dataset)):
+        return obj
+
+    # Topology can either be a sequence of grids or a grid.
+    if isinstance(topology, (list, set, tuple)):
+        grids = {dim: grid for grid in topology for dim in grid.dimensions}
+    else:
+        grids = {dim: topology for dim in topology.dimensions}
+
+    item_grids = list(set(grids[dim] for dim in obj.dims if dim in grids))
+
+    if len(item_grids) == 0:
+        return obj
+    else:
+        result, aligned = align(obj, item_grids, old_indexes)
+
+        if isinstance(result, xr.DataArray):
+            if len(aligned) > 1:
+                raise RuntimeError("This shouldn't happen. Please open an issue.")
+            return UgridDataArray(result, aligned[0])
+
+        elif isinstance(result, xr.Dataset):
+            return UgridDataset(result, aligned)
+
+
+def maybe_xarray(arg):
+    if isinstance(arg, (UgridDataArray, UgridDataset)):
+        return arg.obj
+    else:
+        return arg
+
+
+def wraps_xarray(method):
+    @wraps(method)
+    def wrapped(*args, **kwargs):
+        self = args[0]
+        args = [maybe_xarray(arg) for arg in args]
+        kwargs = {k: maybe_xarray(v) for k, v in kwargs.items()}
+        result = method(*args, **kwargs)
+
+        # Sidestep staticmethods, classmethods
+        if isinstance(self, UgridDataArray):
+            return maybe_xugrid(result, [self.grid], self.obj.indexes)
+        elif isinstance(self, UgridDataset):
+            return maybe_xugrid(result, self.grids, self.obj.indexes)
+        else:
+            return result
+
+    return wrapped
+
+
+def wrap_accessor(accessor):
+    # TODO: This will not add dynamic accessors, those most be included at
+    # runtime instead?
+
+    def wrapped(*args, **kwargs):
+        args = [maybe_xarray(arg) for arg in args]
+        kwargs = {k: maybe_xarray(v) for k, v in kwargs.items()}
+        result = accessor(*args, **kwargs)
+        return result
+
+    return wrapped
+
+
+def wrap(
+    target_class_dict,
+    source_class,
+):
+    FuncType = (types.FunctionType, types.MethodType)
+
+    # Set every method, property from the xarray object to the UgridDataArray,
+    # UgridDataset. Don't set everything, as this will break the objects.
+    #
+    # class Empty:
+    #     pass
+    #
+    # keep = {
+    #     "__eq__",
+    #     "__ge__",
+    #     "__gt__",
+    #     "__le__",
+    #     "__lt__",
+    #     "__ne__",
+    #     "__repr__",
+    #     "__str__",
+    # }
+    #
+    # remove = set(dir(Empty)) - keep
+
+    remove = {
+        # These members are shared by all objects:
+        "__class__",
+        "__delattr__",
+        "__dict__",
+        "__dir__",
+        "__doc__",
+        "__format__",
+        "__getattribute__",
+        "__hash__",
+        "__init__",
+        "__init_subclass__",
+        "__module__",
+        "__new__",
+        "__reduce__",
+        "__reduce_ex__",
+        "__setattr__",
+        "__sizeof__",
+        "__subclasshook__",
+        "__weakref__"
+        # These are additionally included in xarray:
+        "__getatrr__",
+        "__slots__",
+        "__annotations__",
+    }
+
+    attr_names = set(dir(source_class)) - remove
+    all_attrs = {k: getattr(source_class, k) for k in attr_names}
+
+    methods = {k: v for k, v in all_attrs.items() if isinstance(v, FuncType)}
+    for name, method in methods.items():
+        wrapped = wraps_xarray(method)
+        setattr(wrapped, "__doc__", method.__doc__)
+        target_class_dict[name] = wrapped
+
+    properties = {k: v for k, v in all_attrs.items() if isinstance(v, property)}
+    for name, prop in properties.items():
+        wrapped = property(wraps_xarray(prop.__get__))
+        setattr(wrapped, "__doc__", prop.__doc__)
+        target_class_dict[name] = wrapped
+
+    accessors = {k: v for k, v in all_attrs.items() if isinstance(v, type)}
+    for name, accessor in accessors.items():
+        wrapped = property(wrap_accessor(accessor))
+        setattr(wrapped, "__doc__", accessor.__doc__)
+        target_class_dict[name] = wrapped
+
+    return
+
+
+class DataArrayForwardMixin:
+
+    wrap(
+        target_class_dict=vars(),
+        source_class=xr.DataArray,
+    )
+
+
+class DatasetForwardMixin:
+
+    wrap(
+        target_class_dict=vars(),
+        source_class=xr.Dataset,
+    )
+
+
+def assign_ugrid_coords(obj, grids):
+    grid_dims = ChainMap(*(grid.dimensions for grid in grids))
+    ugrid_dims = set(grid_dims.keys()).intersection(obj.dims)
+    ugrid_coords = {dim: RangeIndex(0, grid_dims[dim]) for dim in ugrid_dims}
+    obj = obj.assign_coords(ugrid_coords)
+    return obj
+
+
+class UgridDataArray(DataArrayForwardMixin):
+    def __init__(self, obj: xr.DataArray, grid: UgridType):
+        if not isinstance(obj, xr.DataArray):
+            raise TypeError(
+                "obj must be xarray.DataArray. Received instead: "
+                f"{type(obj).__name__}"
+            )
+        if not isinstance(grid, AbstractUgrid):
+            raise TypeError(
+                "grid must be Ugrid1d or Ugrid2d. Received instead: "
+                f"{type(grid).__name__}"
+            )
+
+        self._grid = grid
+        self._obj = assign_ugrid_coords(obj, [grid])
+
+    def __getattr__(self, attr):
+        result = getattr(self.obj, attr)
+        return maybe_xugrid(result, [self.grid])
+
+    @property
+    def obj(self):
+        return self._obj
+
+    @property
+    def grid(self):
+        return self._grid
+
+    @property
+    def ugrid(self):
+        """
+        UGRID Accessor. This "accessor" makes operations using the UGRID
+        topology available.
+        """
+        return xugrid.core.dataarray_accessor.UgridDataArrayAccessor(
+            self.obj, self.grid
+        )
+
+    @staticmethod
+    def from_structured(da: xr.DataArray):
+        """
+        Create a UgridDataArray from a (structured) xarray DataArray.
+
+        The spatial dimensions are flattened into a single UGRID face dimension.
+
+        Parameters
+        ----------
+        da: xr.DataArray
+            Last two dimensions must be ``("y", "x")``.
+
+        Returns
+        -------
+        unstructured: UgridDataArray
+        """
+        if da.dims[-2:] != ("y", "x"):
+            raise ValueError('Last two dimensions of da must be ("y", "x")')
+        grid = Ugrid2d.from_structured(da)
+        dims = da.dims[:-2]
+        coords = {k: da.coords[k] for k in dims}
+        face_da = xr.DataArray(
+            da.data.reshape(*da.shape[:-2], -1),
+            coords=coords,
+            dims=[*dims, grid.face_dimension],
+            name=da.name,
+        )
+        return UgridDataArray(face_da, grid)
+
+
+class UgridDataset(DatasetForwardMixin):
+    def __init__(
+        self,
+        obj: xr.Dataset = None,
+        grids: Union[UgridType, Sequence[UgridType]] = None,
+    ):
+        if obj is None and grids is None:
+            raise ValueError("At least either obj or grids is required")
+
+        if obj is None:
+            ds = xr.Dataset()
+        else:
+            if not isinstance(obj, xr.Dataset):
+                raise TypeError(
+                    "obj must be xarray.Dataset. Received instead: "
+                    f"{type(obj).__name__}"
+                )
+            connectivity_vars = [
+                name
+                for v in obj.ugrid_roles.connectivity.values()
+                for name in v.values()
+            ]
+            ds = obj.drop_vars(obj.ugrid_roles.topology + connectivity_vars)
+
+        if grids is None:
+            topologies = obj.ugrid_roles.topology
+            grids = [grid_from_dataset(obj, topology) for topology in topologies]
+        else:
+            # Make sure it's a new list
+            if isinstance(grids, (list, tuple, set)):
+                grids = [grid for grid in grids]
+            else:  # not iterable
+                grids = [grids]
+            # Now typecheck
+            for grid in grids:
+                if not isinstance(grid, AbstractUgrid):
+                    raise TypeError(
+                        "grid must be Ugrid1d or Ugrid2d. "
+                        f"Received instead: {type(grid).__name__}"
+                    )
+
+        self._grids = grids
+        self._obj = assign_ugrid_coords(ds, grids)
+
+    @property
+    def obj(self):
+        return self._obj
+
+    @property
+    def grids(self):
+        return self._grids
+
+    @property
+    def ugrid(self):
+        """
+        UGRID Accessor. This "accessor" makes operations using the UGRID
+        topology available.
+        """
+        return xugrid.core.dataset_accessor.UgridDatasetAccessor(self.obj, self.grids)
+
+    def __getattr__(self, attr):
+        result = getattr(self.obj, attr)
+        return maybe_xugrid(result, self.grids)
+
+    def __setitem__(self, key, value):
+        # TODO: check with topology
+        if isinstance(value, UgridDataArray):
+            append = True
+            # Check if the dimensions occur in self.
+            # if they don't, the grid should be added.
+            if self.grids is not None:
+                alldims = set(
+                    chain.from_iterable([grid.dimensions for grid in self.grids])
+                )
+                matching_dims = set(value.grid.dimensions).intersection(alldims)
+                if matching_dims:
+                    append = False
+                    # If they do match: the grids should match.
+                    grids = {
+                        dim: grid for grid in self.grids for dim in grid.dimensions
+                    }
+                    firstdim = next(iter(matching_dims))
+                    grid_to_check = grids[firstdim]
+                    if not grid_to_check.equals(value.grid):
+                        raise ValueError(
+                            "Grids share dimension names but do not are not identical. "
+                            f"Matching dimensions: {matching_dims}"
+                        )
+
+            self.obj[key] = value.obj
+            if append:
+                self._grids.append(value.grid)
+        else:
+            self.obj[key] = value
+
+    @staticmethod
+    def from_geodataframe(geodataframe: "geopandas.GeoDataFrame"):  # type: ignore # noqa
+        """
+        Convert a geodataframe into the appropriate Ugrid topology and dataset.
+        Parameters
+        ----------
+        geodataframe: gpd.GeoDataFrame
+        Returns
+        -------
+        dataset: UGridDataset
+        """
+        grid = grid_from_geodataframe(geodataframe)
+        ds = xr.Dataset.from_dataframe(geodataframe.drop("geometry", axis=1))
+        return UgridDataset(ds, [grid])
diff --git a/xugrid/data/__init__.py b/xugrid/data/__init__.py
index f507dee9f..b26fef0f3 100644
--- a/xugrid/data/__init__.py
+++ b/xugrid/data/__init__.py
@@ -1,2 +1,2 @@
-from .sample_data import adh_san_diego, elevation_nl, xoxo
-from .synthetic import disk, generate_disk
+from xugrid.data.sample_data import adh_san_diego, elevation_nl, xoxo
+from xugrid.data.synthetic import disk, generate_disk
diff --git a/xugrid/meshkernel_utils.py b/xugrid/meshkernel_utils.py
index c6eac0890..0eadf4bf7 100644
--- a/xugrid/meshkernel_utils.py
+++ b/xugrid/meshkernel_utils.py
@@ -6,8 +6,6 @@
 
 import numpy as np
 
-from .conversion import _to_pygeos
-
 
 def either_string_or_enum(value: Union[str, IntEnum], enum_class: EnumMeta) -> IntEnum:
     """Convert to enum if needed, check value"""
@@ -34,6 +32,8 @@ def to_geometry_list(
     import meshkernel
     import pygeos
 
+    from xugrid.conversion import _to_pygeos
+
     polygon = _to_pygeos([polygon])[0]
     xy = pygeos.get_coordinates(pygeos.get_exterior_ring(polygon))
     return meshkernel.GeometryList(np.array(xy[:, 0]), np.array(xy[:, 1]))
diff --git a/xugrid/plot/__init__.py b/xugrid/plot/__init__.py
index 3a7089e49..43f1dc367 100644
--- a/xugrid/plot/__init__.py
+++ b/xugrid/plot/__init__.py
@@ -1,4 +1,4 @@
-from .plot import (
+from xugrid.plot.plot import (
     contour,
     contourf,
     imshow,
diff --git a/xugrid/plot/plot.py b/xugrid/plot/plot.py
index fcad4a634..9fdfbb809 100644
--- a/xugrid/plot/plot.py
+++ b/xugrid/plot/plot.py
@@ -15,7 +15,7 @@
     label_from_attrs,
 )
 
-from ..typing import FloatDType
+from xugrid.constants import FloatDType
 
 NODE = 0
 EDGE = 1
diff --git a/xugrid/typing.py b/xugrid/typing.py
index 19952b536..0f0dac688 100644
--- a/xugrid/typing.py
+++ b/xugrid/typing.py
@@ -1,8 +1,10 @@
-from typing import NamedTuple, Union
+from typing import NamedTuple, Type, Union
 
 import numpy as np
 from scipy.sparse import coo_matrix, csr_matrix
 
+from xugrid.ugrid.ugridbase import AbstractUgrid
+
 # import numpy.typing as npt
 
 FloatDType = np.float64
@@ -33,6 +35,9 @@ class Vector(NamedTuple):
     y: float
 
 
+UgridType = Type[AbstractUgrid]
+
+
 # Spatial coordinate epsilon for floating point comparison
 # Assuming world coordinates in meters: 40 000 m along equator:
 # 40 000 000 = 4e7 mm
diff --git a/xugrid/ugrid/__init__.py b/xugrid/ugrid/__init__.py
index f1d60e4fa..e69de29bb 100644
--- a/xugrid/ugrid/__init__.py
+++ b/xugrid/ugrid/__init__.py
@@ -1,46 +0,0 @@
-import xarray as xr
-
-from .conventions import UgridRolesAccessor
-from .ugrid1d import Ugrid1d
-from .ugrid2d import Ugrid2d
-from .ugridbase import AbstractUgrid
-
-
-def grid_from_geodataframe(geodataframe: "geopandas.GeoDataFrame"):  # type: ignore # noqa
-    import geopandas as gpd
-
-    gdf = geodataframe
-    if not isinstance(gdf, gpd.GeoDataFrame):
-        raise TypeError(
-            f"Cannot convert a {type(gdf).__name__}, expected a GeoDataFrame"
-        )
-
-    geom_types = gdf.geom_type.unique()
-    if len(geom_types) == 0:
-        raise ValueError("geodataframe contains no geometry")
-    elif len(geom_types) > 1:
-        message = ", ".join(geom_types)
-        raise ValueError(f"Multiple geometry types detected: {message}")
-
-    geom_type = geom_types[0]
-    if geom_type == "LineString":
-        grid = Ugrid1d.from_geodataframe(gdf)
-    elif geom_type == "Polygon":
-        grid = Ugrid2d.from_geodataframe(gdf)
-    else:
-        raise ValueError(
-            f"Invalid geometry type: {geom_type}. Expected Linestring or Polygon."
-        )
-    return grid
-
-
-def grid_from_dataset(dataset: xr.Dataset, topology: str):
-    topodim = dataset[topology].attrs["topology_dimension"]
-    if topodim == 1:
-        return Ugrid1d.from_dataset(dataset, topology)
-    elif topodim == 2:
-        return Ugrid2d.from_dataset(dataset, topology)
-    elif topodim == 3:
-        raise NotImplementedError
-    else:
-        raise ValueError(f"Invalid topology dimension: {topodim}")
diff --git a/xugrid/connectivity.py b/xugrid/ugrid/connectivity.py
similarity index 99%
rename from xugrid/connectivity.py
rename to xugrid/ugrid/connectivity.py
index c8834f59a..ec72b213e 100644
--- a/xugrid/connectivity.py
+++ b/xugrid/ugrid/connectivity.py
@@ -4,7 +4,14 @@
 import numpy as np
 from scipy import sparse
 
-from .typing import X_EPSILON, BoolArray, FloatArray, IntArray, IntDType, SparseMatrix
+from xugrid.constants import (
+    X_EPSILON,
+    BoolArray,
+    FloatArray,
+    IntArray,
+    IntDType,
+    SparseMatrix,
+)
 
 
 class AdjacencyMatrix(NamedTuple):
diff --git a/xugrid/interpolate.py b/xugrid/ugrid/interpolate.py
similarity index 99%
rename from xugrid/interpolate.py
rename to xugrid/ugrid/interpolate.py
index 498e4ce6a..7711bcfae 100644
--- a/xugrid/interpolate.py
+++ b/xugrid/ugrid/interpolate.py
@@ -4,7 +4,7 @@
 import scipy.sparse
 import scipy.sparse.linalg
 
-from .typing import FloatArray
+from xugrid.constants import FloatArray
 
 
 def laplace_interpolate(
diff --git a/xugrid/snapping.py b/xugrid/ugrid/snapping.py
similarity index 98%
rename from xugrid/snapping.py
rename to xugrid/ugrid/snapping.py
index 79b4589d2..cc760be89 100644
--- a/xugrid/snapping.py
+++ b/xugrid/ugrid/snapping.py
@@ -12,11 +12,7 @@
 from scipy.sparse.csgraph import connected_components
 from scipy.spatial import cKDTree
 
-from xugrid.ugrid import Ugrid2d
-
-from . import connectivity
-from .connectivity import AdjacencyMatrix
-from .typing import (
+from xugrid.constants import (
     FloatArray,
     IntArray,
     LineArray,
@@ -24,6 +20,9 @@
     Point,
     Vector,
 )
+from xugrid.ugrid import connectivity
+from xugrid.ugrid.connectivity import AdjacencyMatrix
+from xugrid.ugrid.ugrid2d import Ugrid2d
 
 try:
     import geopandas as gpd
diff --git a/xugrid/ugrid/ugrid1d.py b/xugrid/ugrid/ugrid1d.py
index 7b154d145..d2d65cea1 100644
--- a/xugrid/ugrid/ugrid1d.py
+++ b/xugrid/ugrid/ugrid1d.py
@@ -2,12 +2,14 @@
 from typing import Any, Dict, Tuple, Union
 
 import numpy as np
+import pandas as pd
 import xarray as xr
+from xarray.core.utils import either_dict_or_kwargs
 
-from .. import connectivity, conversion
-from ..typing import BoolArray, FloatArray, FloatDType, IntArray, IntDType
-from . import conventions
-from .ugridbase import AbstractUgrid
+from xugrid import conversion
+from xugrid.constants import BoolArray, FloatArray, FloatDType, IntArray, IntDType
+from xugrid.ugrid import connectivity, conventions
+from xugrid.ugrid.ugridbase import AbstractUgrid, as_pandas_index
 
 
 class Ugrid1d(AbstractUgrid):
@@ -151,6 +153,22 @@ def from_dataset(cls, dataset: xr.Dataset, topology: str = None):
             crs=None,
         )
 
+    def _clear_geometry_properties(self):
+        """Clear all properties that may have been invalidated"""
+        # Meshkernel
+        self._mesh = None
+        self._meshkernel = None
+        # Celltree
+        self._celltree = None
+        # Bounds
+        self._xmin = None
+        self._xmax = None
+        self._ymin = None
+        self._ymax = None
+        # Edges
+        self._edge_x = None
+        self._edge_y = None
+
     @classmethod
     def from_meshkernel(
         cls,
@@ -195,27 +213,15 @@ def to_dataset(self, other: xr.Dataset = None) -> xr.Dataset:
         dataset[self.name].attrs = self._filtered_attrs(dataset)
         return dataset
 
-    def _clear_geometry_properties(self):
-        """Clear all properties that may have been invalidated"""
-        # Meshkernel
-        self._mesh = None
-        self._meshkernel = None
-        # Celltree
-        self._celltree = None
-        # Bounds
-        self._xmin = None
-        self._xmax = None
-        self._ymin = None
-        self._ymax = None
-        # Edges
-        self._edge_x = None
-        self._edge_y = None
-
     @property
     def topology_dimension(self):
         """Highest dimensionality of the geometric elements: 1"""
         return 1
 
+    @property
+    def core_dimension(self):
+        return self.node_dimension
+
     @property
     def dimensions(self):
         return {self.node_dimension: self.n_node, self.edge_dimension: self.n_edge}
@@ -312,10 +318,64 @@ def to_pygeos(self, dim):
                 " Ugrid1d topology."
             )
 
-    def isel(self, dim, indexer):
-        if dim != self.edge_dimension:
-            raise NotImplementedError("Can only index edge_dimension for Ugrid1D")
-        return self.topology_subset(indexer)
+    def sel_points(obj, x, y):
+        raise NotImplementedError("Cannot select points in a UGRID 1D topology")
+
+    def isel(self, indexers=None, return_index=False, **indexers_kwargs):
+        """
+        Select based on node or edge.
+
+        Edge selection always results in a valid UGRID topology. Node selection
+        may result in invalid topologies (incomplete edges), and will error in
+        such a case.
+
+        Parameters
+        ----------
+        indexers: dict of str to np.ndarray of integers or bools
+        return_index: bool, optional
+            Whether to return node_index, edge_index.
+
+        Returns
+        -------
+        obj: xr.Dataset or xr.DataArray
+        grid: Ugrid2d
+        indexes: dict
+            Dictionary with keys node dimension, edge dimension and values
+            their respective index. Only returned if return_index is True.
+        """
+        indexers = either_dict_or_kwargs(indexers, indexers_kwargs, "isel")
+        alldims = set(self.dimensions)
+        invalid = indexers.keys() - alldims
+        if invalid:
+            raise ValueError(
+                f"Dimensions {invalid} do not exist. Expected one of {alldims}"
+            )
+
+        indexers = {
+            k: as_pandas_index(v, self.dimensions[k]) for k, v in indexers.items()
+        }
+        nodedim, edgedim = self.dimensions
+        edge_index = {}
+        if nodedim in indexers:
+            node_index = indexers[nodedim]
+            edge_index[nodedim] = np.unique(
+                self.node_edge_connectivity[node_index].data
+            )
+        if edgedim in indexers:
+            edge_index[edgedim] = indexers[edgedim]
+
+        # Convert all to pandas index.
+        edge_index = {k: as_pandas_index(v, self.n_edge) for k, v in edge_index.items()}
+
+        # Check the indexes against each other.
+        index = self._precheck(edge_index)
+        grid, finalized_indexers = self.topology_subset(index, return_index=True)
+        self._postcheck(indexers, finalized_indexers)
+
+        if return_index:
+            return grid, finalized_indexers
+        else:
+            return grid
 
     def _validate_indexer(self, indexer) -> Tuple[float, float]:
         if isinstance(indexer, slice):
@@ -327,7 +387,7 @@ def _validate_indexer(self, indexer) -> Tuple[float, float]:
             raise ValueError("Ugrid1d only supports slice indexing")
         return indexer.start, indexer.stop
 
-    def sel(self, x, y) -> Tuple[str, bool, IntArray, dict]:
+    def sel(self, obj, x, y) -> Tuple:
         """
         Select a selection of edges, based on edge centroids.
 
@@ -345,15 +405,20 @@ def sel(self, x, y) -> Tuple[str, bool, IntArray, dict]:
         """
         xmin, xmax = self._validate_indexer(x)
         ymin, ymax = self._validate_indexer(y)
-        index = np.nonzero(
+        edge_index = np.nonzero(
             (self.edge_x >= xmin)
             & (self.edge_x < xmax)
             & (self.edge_y >= ymin)
             & (self.edge_y < ymax)
         )[0]
-        return self.edge_dimension, True, index, {}
+        grid, indexes = self.topology_subset(edge_index, return_index=True)
+        indexes = {k: v for k, v in indexes.items() if k in obj.dims}
+        new_obj = obj.isel(indexes)
+        return new_obj, grid
 
-    def topology_subset(self, edge_index: Union[BoolArray, IntArray]):
+    def topology_subset(
+        self, edge_index: Union[BoolArray, IntArray], return_index: bool = False
+    ):
         """
         Create a new UGRID1D topology for a subset of this topology.
 
@@ -361,21 +426,35 @@ def topology_subset(self, edge_index: Union[BoolArray, IntArray]):
         ----------
         edge_index: 1d array of integers or bool
             Edges of the subset.
+        return_index: bool, optional
+            Whether to return node_index, edge_index.
 
         Returns
         -------
         subset: Ugrid1d
+        indexes: dict
+            Dictionary with keys node dimension and edge dimension and values
+            their respective index. Only returned if return_index is True.
         """
-        index = np.atleast_1d(edge_index)
-        if self._check_index_identity(index, len(self.edge_node_connectivity)):
-            return self
-
-        edge_subset = self.edge_node_connectivity[index]
+        if not isinstance(edge_index, pd.Index):
+            edge_index = as_pandas_index(edge_index, self.n_edge)
+
+        if edge_index.size == self.n_edge:
+            if return_index:
+                indexes = {
+                    self.node_dimension: pd.RangeIndex(0, self.n_node),
+                    self.edge_dimension: pd.RangeIndex(0, self.n_edge),
+                }
+                return self, indexes
+            else:
+                return self
+
+        edge_subset = self.edge_node_connectivity[edge_index]
         node_index = np.unique(edge_subset.ravel())
         new_edges = connectivity.renumber(edge_subset)
         node_x = self.node_x[node_index]
         node_y = self.node_y[node_index]
-        return self.__class__(
+        grid = self.__class__(
             node_x,
             node_y,
             self.fill_value,
@@ -385,6 +464,14 @@ def topology_subset(self, edge_index: Union[BoolArray, IntArray]):
             crs=self.crs,
             attrs=self._attrs,
         )
+        if return_index:
+            indexes = {
+                self.node_dimension: pd.Index(node_index),
+                self.edge_dimension: edge_index,
+            }
+            return grid, indexes
+        else:
+            return grid
 
     def clip_box(
         self,
diff --git a/xugrid/ugrid/ugrid2d.py b/xugrid/ugrid/ugrid2d.py
index 9a75a1aee..9b08e7f1f 100644
--- a/xugrid/ugrid/ugrid2d.py
+++ b/xugrid/ugrid/ugrid2d.py
@@ -2,17 +2,26 @@
 from typing import Any, Dict, Tuple, Union
 
 import numpy as np
+import pandas as pd
 import xarray as xr
 from numba_celltree import CellTree2d
 from scipy.sparse import coo_matrix, csr_matrix
 from scipy.sparse.csgraph import reverse_cuthill_mckee
-
-from .. import connectivity, conversion
-from .. import meshkernel_utils as mku
-from ..typing import BoolArray, FloatArray, FloatDType, IntArray, IntDType, SparseMatrix
-from ..voronoi import voronoi_topology
-from . import conventions
-from .ugridbase import AbstractUgrid
+from xarray.core.utils import either_dict_or_kwargs
+
+from xugrid import conversion
+from xugrid import meshkernel_utils as mku
+from xugrid.constants import (
+    BoolArray,
+    FloatArray,
+    FloatDType,
+    IntArray,
+    IntDType,
+    SparseMatrix,
+)
+from xugrid.ugrid import connectivity, conventions
+from xugrid.ugrid.ugridbase import AbstractUgrid, as_pandas_index
+from xugrid.ugrid.voronoi import voronoi_topology
 
 
 def section_coordinates(
@@ -305,6 +314,10 @@ def n_face(self):
         """
         return self.face_node_connectivity.shape[0]
 
+    @property
+    def core_dimension(self):
+        return self.face_dimension
+
     @property
     def dimensions(self):
         return {
@@ -782,7 +795,7 @@ def topology_subset(
         self, face_index: Union[BoolArray, IntArray], return_index: bool = False
     ):
         """
-        Create a new (valid) UGRID1D topology for a subset of this topology.
+        Create a new UGRID1D topology for a subset of this topology.
 
         Parameters
         ----------
@@ -794,11 +807,28 @@ def topology_subset(
         Returns
         -------
         subset: Ugrid2d
-        """
-        index = np.atleast_1d(face_index)
-        if self._check_index_identity(index, self.n_face):
-            return self
+        indexes: dict
+            Dictionary with keys node dimension, edge dimension, face dimension
+            and values their respective index. Only returned if return_index is
+            True.
+        """
+        if not isinstance(face_index, pd.Index):
+            face_index = as_pandas_index(face_index, self.n_face)
+
+        # The pandas index may only contain uniques. So if size matches, it's the identity.
+        if face_index.size == self.n_face:
+            # TODO: return self.copy instead?
+            if return_index:
+                indexes = {
+                    self.node_dimension: pd.RangeIndex(0, self.n_node),
+                    self.edge_dimension: pd.RangeIndex(0, self.n_edge),
+                    self.face_dimension: pd.RangeIndex(0, self.n_face),
+                }
+                return self, indexes
+            else:
+                return self
 
+        index = face_index.values
         face_subset = self.face_node_connectivity[index]
         node_index = np.unique(face_subset.ravel())
         new_faces = connectivity.renumber(face_subset)
@@ -808,7 +838,8 @@ def topology_subset(
         edge_index = None
         new_edges = None
         if self.edge_node_connectivity is not None:
-            edge_index = np.sort(self.face_edge_connectivity[index].ravel())
+            edge_index = np.unique(self.face_edge_connectivity[index].ravel())
+            edge_index = edge_index[edge_index != self.fill_value]
             edge_subset = self.edge_node_connectivity[edge_index]
             new_edges = connectivity.renumber(edge_subset)
 
@@ -824,7 +855,13 @@ def topology_subset(
             attrs=self._attrs,
         )
         if return_index:
-            return grid, (node_index, edge_index, face_index)
+            indexes = {
+                self.node_dimension: pd.Index(node_index),
+                self.face_dimension: face_index,
+            }
+            if edge_index is not None:
+                indexes[self.edge_dimension] = pd.Index(edge_index)
+            return grid, indexes
         else:
             return grid
 
@@ -840,19 +877,7 @@ def clip_box(
         face_index = self.locate_bounding_box(*bounds)
         return self.topology_subset(face_index)
 
-    def _isel_obj(
-        self,
-        obj,
-        indexes,
-    ):
-        indexes = {k: v for k, v in zip(self.dimensions, indexes)}
-        dims = set(self.dimensions).intersection(obj.dims)
-        new_obj = obj
-        for dim in dims:
-            new_obj = new_obj.isel({dim: indexes[dim]})
-        return new_obj
-
-    def isel(self, dim, indexer, obj):
+    def isel(self, indexers=None, return_index=False, **indexers_kwargs):
         """
         Select based on node, edge, or face.
 
@@ -862,43 +887,56 @@ def isel(self, dim, indexer, obj):
 
         Parameters
         ----------
-        indexer: str
-            One of the grid's node, edge, or face dimensions.
-        indexer: np.ndarray of integers or bools
-        obj: xr.Dataset or xr.DataArray
+        indexers: dict of str to np.ndarray of integers or bools
+        return_index: bool, optional
+            Whether to return node_index, edge_index, face_index.
 
         Returns
         -------
         obj: xr.Dataset or xr.DataArray
         grid: Ugrid2d
-        """
-        msg = (
-            "This {entity} selection results in an invalid Ugrid2d topology. "
-            "Use .ugrid.sel() instead."
-        )
+        indexes: dict
+            Dictionary with keys node dimension, edge dimension, face dimension
+            and values their respective index. Only returned if return_index is
+            True.
+        """
+        indexers = either_dict_or_kwargs(indexers, indexers_kwargs, "isel")
+        alldims = set(self.dimensions)
+        invalid = indexers.keys() - alldims
+        if invalid:
+            raise ValueError(
+                f"Dimensions {invalid} do not exist. Expected one of {alldims}"
+            )
 
-        if dim == self.face_dimension:
-            grid, indexes = self.topology_subset(indexer, True)
+        indexers = {
+            k: as_pandas_index(v, self.dimensions[k]) for k, v in indexers.items()
+        }
+        nodedim, edgedim, facedim = self.dimensions
+        face_index = {}
+        if nodedim in indexers:
+            node_index = indexers[nodedim]
+            face_index[nodedim] = np.unique(
+                self.node_face_connectivity[node_index].data
+            )
+        if edgedim in indexers:
+            edge_index = indexers[edgedim]
+            index = np.unique(self.edge_face_connectivity[edge_index])
+            face_index[edgedim] = index[index != self.fill_value]
+        if facedim in indexers:
+            face_index[facedim] = indexers[facedim]
 
-        elif dim == self.edge_dimension:
-            face_index = np.unique(self.node_face_connectivity[indexer].data)
-            grid, indexes = self.topology_subset(self, face_index, True)
-            other_index = indexes[0]
-            if not np.array_equal(np.sort(indexer), other_index):
-                raise ValueError(msg.format(entity="edge"))
+        # Convert all to pandas index.
+        face_index = {k: as_pandas_index(v, self.n_face) for k, v in face_index.items()}
 
-        elif dim == self.node_dimension:
-            face_index = np.unique(self.edge_face_connectivity[indexer])
-            grid, indexes = self.topology_subset(self, face_index, True)
-            other_index = indexes[1]
-            if not np.array_equal(np.sort(indexer), other_index):
-                raise ValueError(msg.format(entity="node"))
+        # Check the indexes against each other.
+        index = self._precheck(face_index)
+        grid, finalized_indexers = self.topology_subset(index, return_index=True)
+        self._postcheck(indexers, finalized_indexers)
 
+        if return_index:
+            return grid, finalized_indexers
         else:
-            raise ValueError(f"dim should be one of {self.dimensions}, received: {dim}")
-
-        result = self._isel_obj(obj, indexes)
-        return result, grid
+            return grid
 
     def _validate_indexer(self, indexer) -> Union[slice, np.ndarray]:
         if isinstance(indexer, slice):
@@ -934,9 +972,9 @@ def _validate_indexer(self, indexer) -> Union[slice, np.ndarray]:
 
     def _sel_box(
         self,
+        obj,
         x,
         y,
-        obj,
     ):
         xmin, ymin, xmax, ymax = self.bounds
         bounds = [
@@ -947,14 +985,15 @@ def _sel_box(
         ]
         face_index = self.locate_bounding_box(*bounds)
         grid, indexes = self.topology_subset(face_index, return_index=True)
-        new_obj = self._isel_obj(obj, indexes)
+        indexes = {k: v for k, v in indexes.items() if k in obj.dims}
+        new_obj = obj.isel(indexes)
         return new_obj, grid
 
     def _sel_yline(
         self,
+        obj,
         x,
         y,
-        obj,
     ):
         xmin, _, xmax, _ = self.bounds
         if y.size != 1:
@@ -971,9 +1010,9 @@ def _sel_yline(
 
     def _sel_xline(
         self,
+        obj,
         x,
         y,
-        obj,
     ):
         _, ymin, _, ymax = self.bounds
         if x.size != 1:
@@ -988,7 +1027,7 @@ def _sel_xline(
         coords, index = section_coordinates(edges, xy, self.face_dimension, index)
         return obj.isel({self.face_dimension: index}).assign_coords(coords)
 
-    def _sel_points(self, x: FloatArray, y: FloatArray, obj):
+    def sel_points(self, obj, x: FloatArray, y: FloatArray):
         """
         Select points in the unstructured grid.
 
@@ -996,13 +1035,11 @@ def _sel_points(self, x: FloatArray, y: FloatArray, obj):
         ----------
         x: 1d array of floats with shape ``(n_points,)``
         y: 1d array of floats with shape ``(n_points,)``
+        obj: xr.DataArray or xr.Dataset
 
         Returns
         -------
-        dimension: str
-        as_ugrid: bool
-        index: 1d array of integers
-        coords: dict
+        selection: xr.DataArray or xr.Dataset
         """
         x = np.atleast_1d(x)
         y = np.atleast_1d(y)
@@ -1052,19 +1089,20 @@ def sel(self, obj, x=None, y=None):
         x = self._validate_indexer(x)
         y = self._validate_indexer(y)
         if isinstance(x, slice) and isinstance(y, slice):
-            return self._sel_box(x, y, obj)
+            f = self._sel_box
         elif isinstance(x, slice) and isinstance(y, np.ndarray):
-            return self._sel_yline(x, y, obj)
+            f = self._sel_yline
         elif isinstance(x, np.ndarray) and isinstance(y, slice):
-            return self._sel_xline(x, y, obj)
+            f = self._sel_xline
         elif isinstance(x, np.ndarray) and isinstance(y, np.ndarray):
             # Orthogonal points
-            yy, xx = np.meshgrid(y, x, indexing="ij")
-            return self.sel_points(xx.ravel(), yy.ravel())
+            y, x = [a.ravel() for a in np.meshgrid(y, x, indexing="ij")]
+            f = self.sel_points
         else:
             raise TypeError(
                 f"Invalid indexer types: {type(x).__name__}, and {type(y).__name__}"
             )
+        return f(obj, x, y)
 
     def triangulate(self):
         """
diff --git a/xugrid/ugrid/ugridbase.py b/xugrid/ugrid/ugridbase.py
index e3647cfc2..e69ea34e3 100644
--- a/xugrid/ugrid/ugridbase.py
+++ b/xugrid/ugrid/ugridbase.py
@@ -1,14 +1,85 @@
 import abc
 import copy
-from typing import Tuple, Union
+from itertools import chain
+from typing import Tuple, Type, Union
 
 import numpy as np
+import pandas as pd
 import xarray as xr
 from scipy.sparse import csr_matrix
 
-from .. import connectivity
-from ..typing import BoolArray, FloatArray, IntArray
-from . import conventions
+from xugrid.constants import BoolArray, FloatArray, IntArray
+from xugrid.ugrid import connectivity, conventions
+
+
+def as_pandas_index(index: Union[BoolArray, IntArray, pd.Index], n: int):
+    if isinstance(index, np.ndarray):
+        if index.size > n:
+            raise ValueError(
+                f"index size {index.size} is larger than dimension size: {n}"
+            )
+        if np.issubdtype(index.dtype, np.bool_):
+            # Significantly quicker if all true.
+            if index.all():
+                pd_index = pd.RangeIndex(0, n)
+            else:
+                pd_index = pd.Index((np.arange(n)[index]))
+        elif np.issubdtype(index.dtype, np.integer):
+            pd_index = pd.Index(index)
+        else:
+            raise TypeError(f"index should be bool or integer. Received: {index.dtype}")
+
+    elif isinstance(index, pd.Index):
+        pd_index = index
+
+    else:
+        raise TypeError(
+            "index should be pandas Index or numpy array. Received: "
+            f"{type(index).__name__}"
+        )
+
+    if not pd_index.is_unique:
+        raise ValueError(
+            "index contains repeated values; only subsets will result "
+            "in valid UGRID topology."
+        )
+    if not pd_index.is_monotonic_increasing:
+        raise NotImplementedError("UGRID indexes must be sorted and unique.")
+
+    return pd_index
+
+
+def align(obj, grids, old_indexes):
+    """
+    Check which indexes have changed. Index on those new values.
+    If none are changed, return (obj, grids) as is.
+    """
+    if old_indexes is None:
+        return obj, grids
+
+    ugrid_dims = set(
+        chain.from_iterable(grid.dimensions for grid in grids)
+    ).intersection(old_indexes)
+    new_indexes = {
+        k: index
+        for k, index in obj.indexes.items()
+        if (k in ugrid_dims) and (not index.equals(old_indexes[k]))
+    }
+    if not new_indexes:
+        return obj, grids
+
+    # Group the indexers by grid
+    new_grids = []
+    for grid in grids:
+        ugrid_dims = set(grid.dimensions).intersection(new_indexes)
+        ugrid_indexes = {dim: new_indexes[dim] for dim in ugrid_dims}
+        newgrid, indexers = grid.isel(indexers=ugrid_indexes, return_index=True)
+        indexers = {
+            k: v for k, v in indexers.items() if k in obj.dims and k not in new_indexes
+        }
+        obj = obj.isel(indexers)
+        new_grids.append(newgrid)
+    return obj, new_grids
 
 
 class AbstractUgrid(abc.ABC):
@@ -16,6 +87,10 @@ class AbstractUgrid(abc.ABC):
     def topology_dimension():
         """ """
 
+    @abc.abstractproperty
+    def core_dimension():
+        """ """
+
     @abc.abstractproperty
     def dimensions():
         """ """
@@ -37,7 +112,7 @@ def clip_box():
         """ """
 
     @abc.abstractmethod
-    def isel():
+    def sel_points():
         """ """
 
     @abc.abstractmethod
@@ -234,17 +309,24 @@ def _prepare_connectivity(
             raise ValueError("connectivity contains negative values")
         return da.copy(data=cast)
 
-    @staticmethod
-    def _check_index_identity(index: Union[BoolArray, IntArray], n: int) -> bool:
-        """
-        Check whether index would return all values.
-        """
-        if np.issubdtype(index.dtype, np.bool_):
-            return index.all()
-        elif np.issubdtype(index.dtype, np.integer):
-            return np.array_equal(index, np.arange(n))
-        else:
-            raise TypeError(f"index should be bool or integer. Received: {index.dtype}")
+    def _precheck(self, multi_index):
+        dim, index = multi_index.popitem()
+        for check_dim, check_index in multi_index.items():
+            if not index.equals(check_index):
+                raise ValueError(
+                    f"UGRID dimensions do not align: {dim} versus {check_dim}"
+                )
+        return index
+
+    def _postcheck(self, indexers, finalized_indexers):
+        for dim, indexer in indexers.items():
+            if dim != self.core_dimension:
+                if not indexer.equals(finalized_indexers[dim]):
+                    raise ValueError(
+                        f"This subset selection of UGRID dimension {dim} results"
+                        "in an invalid topology "
+                    )
+        return
 
     def set_node_coords(
         self,
@@ -490,6 +572,9 @@ def plot(self, **kwargs):
         **kwargs : optional
             Additional keyword arguments to ``matplotlib.pyplot.line``.
         """
-        from ..plot import line
+        from xugrid.plot import line
 
         return line(self, **kwargs)
+
+
+UgridType = Type[AbstractUgrid]
diff --git a/xugrid/voronoi.py b/xugrid/ugrid/voronoi.py
similarity index 99%
rename from xugrid/voronoi.py
rename to xugrid/ugrid/voronoi.py
index 0ef6f9e8a..0ffa05c4f 100644
--- a/xugrid/voronoi.py
+++ b/xugrid/ugrid/voronoi.py
@@ -20,8 +20,8 @@
 import pandas as pd
 from scipy import sparse
 
-from .connectivity import renumber
-from .typing import FloatArray, IntArray
+from xugrid.constants import FloatArray, IntArray
+from xugrid.ugrid.connectivity import renumber
 
 
 def dot_product2d(U: FloatArray, V: FloatArray):
diff --git a/xugrid/ugrid_dataset.py b/xugrid/ugrid_dataset.py
deleted file mode 100644
index 13688d9fd..000000000
--- a/xugrid/ugrid_dataset.py
+++ /dev/null
@@ -1,1568 +0,0 @@
-import abc
-import types
-from functools import wraps
-from itertools import chain
-from typing import Any, Callable, Dict, Sequence, Tuple, Type, Union
-
-import numpy as np
-import pandas as pd
-import scipy.sparse
-import xarray as xr
-from xarray.backends.api import DATAARRAY_NAME, DATAARRAY_VARIABLE
-from xarray.core._typed_ops import DataArrayOpsMixin, DatasetOpsMixin
-from xarray.core.utils import UncachedAccessor, either_dict_or_kwargs
-
-from . import connectivity
-from .interpolate import laplace_interpolate
-from .plot.plot import _PlotMethods
-
-# from .plot.pyvista import to_pyvista_grid
-from .ugrid import AbstractUgrid, Ugrid2d, grid_from_dataset, grid_from_geodataframe
-
-UgridType = Type[AbstractUgrid]
-
-
-def maybe_xugrid(obj, grid):
-    if isinstance(obj, xr.DataArray):
-        return UgridDataArray(obj, grid)
-    elif isinstance(obj, xr.Dataset):
-        return UgridDataset(obj, grid)
-    else:
-        return obj
-        # raise TypeError(f"Expected Dataset or DataArray, received {type(obj).__name__}")
-
-
-def check_ugrid_alignment(obj, grid):
-    """
-    Check whether the xarray object dimensions still align with the grid.
-    """
-    griddims = grid.dimensions
-    shared_dims = set(griddims).intersection(obj.dims)
-    if shared_dims:
-        for dim in shared_dims:
-            ugridsize = griddims[dim]
-            objsize = obj[dim].size
-            if ugridsize != objsize:
-                raise ValueError(
-                    f"conflicting sizes for dimension '{dim}': "
-                    f"length {ugridsize} in UGRID topology and "
-                    f"length {objsize} in xarray dimension"
-                )
-        return maybe_xugrid(obj, grid)
-    else:
-        return obj
-
-
-def xarray_wrapper(func, grid):
-    """
-    Runs a function, and if the result is an xarray Dataset or an xarray
-    DataArray, it creates an UgridDataset or an UgridDataArray around it
-    if UGRID dimensions are present.
-    """
-
-    @wraps(func)
-    def wrapped(*args, **kwargs):
-        result = func(*args, **kwargs)
-        if isinstance(result, (xr.Dataset, xr.DataArray)):
-            return check_ugrid_alignment(result, grid)
-        else:
-            return result
-
-    return wrapped
-
-
-def filter_indexers(indexers, grids):
-    indexers = indexers.copy()
-    ugrid_indexers = []
-    for grid in grids:
-        ugrid_dims = set(grid.dimensions).intersection(indexers)
-        if ugrid_dims:
-            if len(ugrid_dims) > 1:
-                raise ValueError(
-                    "Can only index along one UGRID dimension at a time. "
-                    f"Received for topology {grid.name}: {ugrid_dims}"
-                )
-            dim = ugrid_dims.pop()
-            ugrid_indexers.append((grid, (dim, indexers.pop(dim))))
-        else:
-            ugrid_indexers.append((grid, None))
-    return indexers, ugrid_indexers
-
-
-def ugrid_sel(obj, ugrid_indexers):
-    grids = []
-    for (grid, indexer_args) in ugrid_indexers:
-        if indexer_args is None:
-            grids.append(grid)
-        else:
-            obj, new_grid = grid.isel(*indexer_args, obj)
-            grids.append(new_grid)
-    return obj, grids
-
-
-class DunderForwardMixin:
-    """
-    These methods are not present in the xarray OpsMixin classes.
-    """
-
-    def __bool__(self):
-        return bool(self.obj)
-
-    def __contains__(self, key: Any):
-        return key in self.obj
-
-    def __int__(self):
-        return int(self.obj)
-
-    def __float__(self):
-        return float(self.obj)
-
-
-class UgridDataArray(DataArrayOpsMixin, DunderForwardMixin):
-    """
-    Wraps an xarray DataArray, adding UGRID topology.
-    """
-
-    def __repr__(self):
-        return self.obj.__repr__()
-
-    def __init__(self, obj: xr.DataArray, grid: UgridType):
-        if not isinstance(obj, xr.DataArray):
-            raise TypeError(
-                "obj must be xarray.DataArray. Received instead: "
-                f"{type(obj).__name__}"
-            )
-        if not isinstance(grid, AbstractUgrid):
-            raise TypeError(
-                "grid must be Ugrid1d or Ugrid2d. Received instead: "
-                f"{type(grid).__name__}"
-            )
-        self.grid = grid
-        self.obj = obj
-
-    def __getitem__(self, key):
-        """
-        forward getters to xr.DataArray. Wrap result if necessary
-        """
-        result = self.obj[key]
-        return maybe_xugrid(result, self.grid)
-
-    def __setitem__(self, key, value):
-        """
-        forward setters to  xr.DataArray
-        """
-        self.obj[key] = value
-
-    def __getattr__(self, attr):
-        """
-        Appropriately wrap result if necessary.
-        """
-        if attr == "obj":
-            return self.obj
-
-        result = getattr(self.obj, attr)
-        if isinstance(result, xr.DataArray):
-            return UgridDataArray(result, self.grid)
-        elif isinstance(result, types.MethodType):
-            return xarray_wrapper(result, self.grid)
-        else:
-            return result
-
-    def _unary_op(
-        self,
-        f: Callable,
-    ):
-        return UgridDataArray(self.obj._unary_op(f), self.grid)
-
-    def _binary_op(
-        self,
-        other,
-        f: Callable,
-        reflexive: bool = False,
-    ):
-        if isinstance(other, (UgridDataArray, UgridDataset)):
-            other = other.obj
-        return UgridDataArray(self.obj._binary_op(other, f, reflexive), self.grid)
-
-    def _inplace_binary_op(self, other, f: Callable):
-        if isinstance(other, (UgridDataArray, UgridDataset)):
-            other = other.obj
-        return UgridDataArray(self.obj._inplace_binary_op(other, f), self.grid)
-
-    def isel(
-        self,
-        indexers=None,
-        drop: bool = False,
-        missing_dims="raise",
-        **indexers_kwargs,
-    ):
-        indexers = either_dict_or_kwargs(indexers, indexers_kwargs, "sel")
-        obj_indexers, ugrid_indexers = filter_indexers(indexers, self.grids)
-        result = self.obj.isel(obj_indexers, drop=drop, missing_dims=missing_dims)
-        result, grids = ugrid_sel(result, ugrid_indexers)
-        return UgridDataset(result, grids[0])
-
-    def sel(
-        self,
-        indexers=None,
-        method=None,
-        tolerance=None,
-        drop=False,
-        **indexers_kwargs,
-    ):
-        indexers = either_dict_or_kwargs(indexers, indexers_kwargs, "sel")
-        obj_indexers, ugrid_indexers = filter_indexers(indexers, [self.grid])
-        result = self.obj.sel(
-            obj_indexers, method=method, tolerance=tolerance, drop=drop
-        )
-        result, grids = ugrid_sel(result, ugrid_indexers)
-        return UgridDataset(result, grids[0])
-
-    @property
-    def ugrid(self):
-        """
-        UGRID Accessor. This "accessor" makes operations using the UGRID
-        topology available.
-        """
-        return UgridDataArrayAccessor(self.obj, self.grid)
-
-    @staticmethod
-    def from_structured(da: xr.DataArray):
-        """
-        Create a UgridDataArray from a (structured) xarray DataArray.
-
-        The spatial dimensions are flattened into a single UGRID face dimension.
-
-        Parameters
-        ----------
-        da: xr.DataArray
-            Last two dimensions must be ``("y", "x")``.
-
-        Returns
-        -------
-        unstructured: UgridDataArray
-        """
-        if da.dims[-2:] != ("y", "x"):
-            raise ValueError('Last two dimensions of da must be ("y", "x")')
-        grid = Ugrid2d.from_structured(da)
-        dims = da.dims[:-2]
-        coords = {k: da.coords[k] for k in dims}
-        face_da = xr.DataArray(
-            da.data.reshape(*da.shape[:-2], -1),
-            coords=coords,
-            dims=[*dims, grid.face_dimension],
-            name=da.name,
-        )
-        return UgridDataArray(face_da, grid)
-
-
-UgridDataArray.isel.__doc__ = xr.DataArray.isel.__doc__
-UgridDataArray.sel.__doc__ = xr.DataArray.sel.__doc__
-
-
-class AbstractUgridAccessor(abc.ABC):
-    @abc.abstractmethod
-    def to_dataset():
-        """ """
-
-    @abc.abstractmethod
-    def assign_node_coords():
-        """ """
-
-    @abc.abstractmethod
-    def set_node_coords():
-        """ """
-
-    @abc.abstractproperty
-    def crs():
-        """ """
-
-    @abc.abstractmethod
-    def set_crs():
-        """ """
-
-    @abc.abstractmethod
-    def to_crs():
-        """ """
-
-    @abc.abstractproperty
-    def bounds():
-        """ """
-
-    @abc.abstractproperty
-    def total_bounds():
-        """ """
-
-    @staticmethod
-    def _sel(obj, grid, x, y):
-        # TODO: also do vectorized indexing like xarray?
-        # Might not be worth it, as orthogonal and vectorized indexing are
-        # quite confusing.
-        result = grid.sel(obj, x, y)
-        if isinstance(result, tuple):
-            return maybe_xugrid(*result)
-        else:
-            return result
-
-    @staticmethod
-    def _sel_points(obj, grid, x, y):
-        """
-        Select points in the unstructured grid.
-
-        Parameters
-        ----------
-        x: ndarray of floats with shape ``(n_points,)``
-        y: ndarray of floats with shape ``(n_points,)``
-
-        Returns
-        -------
-        points: Union[xr.DataArray, xr.Dataset]
-        """
-        if grid.topology_dimension != 2:
-            raise NotImplementedError
-        dim, _, index, coords = grid.sel_points(x, y)
-        result = obj.isel({dim: index})
-        return result.assign_coords(coords)
-
-    def clip_box(
-        self,
-        xmin: float,
-        ymin: float,
-        xmax: float,
-        ymax: float,
-    ):
-        """
-        Clip the DataArray or Dataset by a bounding box.
-
-        Parameters
-        ----------
-        xmin: float
-        ymin: float
-        xmax: float
-        ymax: float
-
-        -------
-        clipped:
-            xugrid.UgridDataArray or xugrid.UgridDataset
-        """
-        return self.sel(x=slice(xmin, xmax), y=slice(ymin, ymax))
-
-    def to_netcdf(self, *args, **kwargs):
-        """
-        Write dataset contents to a UGRID compliant netCDF file.
-
-        This function wraps :py:meth:`xr.Dataset.to_netcdf`; it adds the UGRID
-        variables and coordinates to a standard xarray Dataset, then writes the
-        result to a netCDF.
-
-        All arguments are forwarded to :py:meth:`xr.Dataset.to_netcdf`.
-        """
-        self.to_dataset().to_netcdf(*args, **kwargs)
-
-    def to_zarr(self, *args, **kwargs):
-        """
-        Write dataset contents to a UGRID compliant Zarr file.
-
-        This function wraps :py:meth:`xr.Dataset.to_zarr`; it adds the UGRID
-        variables and coordinates to a standard xarray Dataset, then writes the
-        result to a Zarr file.
-
-        All arguments are forwarded to :py:meth:`xr.Dataset.to_zarr`.
-        """
-        self.to_dataset().to_zarr(*args, **kwargs)
-
-
-class UgridDataset(DatasetOpsMixin, DunderForwardMixin):
-    """
-    Wraps an xarray Dataset, adding UGRID topology.
-    """
-
-    def __repr__(self):
-        return self.obj.__repr__()
-
-    def __init__(
-        self,
-        obj: xr.Dataset = None,
-        grids: Union[UgridType, Sequence[UgridType]] = None,
-    ):
-        if obj is None and grids is None:
-            raise ValueError("At least either obj or grids is required")
-
-        if obj is None:
-            ds = xr.Dataset()
-        else:
-            if not isinstance(obj, xr.Dataset):
-                raise TypeError(
-                    "obj must be xarray.Dataset. Received instead: "
-                    f"{type(obj).__name__}"
-                )
-            connectivity_vars = [
-                name
-                for v in obj.ugrid_roles.connectivity.values()
-                for name in v.values()
-            ]
-            ds = obj.drop_vars(obj.ugrid_roles.topology + connectivity_vars)
-
-        if grids is None:
-            topologies = obj.ugrid_roles.topology
-            grids = [grid_from_dataset(obj, topology) for topology in topologies]
-        else:
-            # Make sure it's a new list
-            if isinstance(grids, (list, tuple, set)):
-                grids = [grid for grid in grids]
-            else:  # not iterable
-                grids = [grids]
-            # Now typecheck
-            for grid in grids:
-                if not isinstance(grid, AbstractUgrid):
-                    raise TypeError(
-                        "grid must be Ugrid1d or Ugrid2d. "
-                        f"Received instead: {type(grid).__name__}"
-                    )
-
-        self.grids = grids
-        self.obj = ds
-
-    def __contains__(self, key: Any):
-        return key in self.obj
-
-    def __getitem__(self, key):
-        result = self.obj[key]
-        grids = {dim: grid for grid in self.grids for dim in grid.dimensions}
-        item_grids = list(set(grids[dim] for dim in result.dims if dim in grids))
-
-        # Now find out which grid to return for this key
-        if isinstance(result, xr.DataArray):
-            if len(item_grids) == 0:
-                return result
-            elif len(item_grids) > 1:
-                raise RuntimeError("This shouldn't happen. Please open an issue.")
-            return UgridDataArray(result, item_grids[0])
-        elif isinstance(result, xr.Dataset):
-            return UgridDataset(result, item_grids)
-        else:
-            return result
-
-    def __setitem__(self, key, value):
-        # TODO: check with topology
-        if isinstance(value, UgridDataArray):
-            append = True
-            # Check if the dimensions occur in self.
-            # if they don't, the grid should be added.
-            if self.grids is not None:
-                alldims = set(
-                    chain.from_iterable([grid.dimensions for grid in self.grids])
-                )
-                matching_dims = set(value.grid.dimensions).intersection(alldims)
-                if matching_dims:
-                    # If they do match: the grids should match.
-                    grids = {
-                        dim: grid for grid in self.grids for dim in grid.dimensions
-                    }
-                    firstdim = next(iter(matching_dims))
-                    if not grids[firstdim].equals(value.grid):
-                        raise ValueError(
-                            "Grids share dimension names but do not are not identical. "
-                            f"Matching dimensions: {matching_dims}"
-                        )
-                    append = False
-
-            self.obj[key] = value.obj
-            if append:
-                self.grids.append(value.grid)
-        else:
-            self.obj[key] = value
-
-    def __getattr__(self, attr):
-        """
-        Appropriately wrap result if necessary.
-        """
-        if attr == "obj":
-            return self.obj
-
-        result = getattr(self.obj, attr)
-        grids = {dim: grid for grid in self.grids for dim in grid.dimensions}
-        if isinstance(result, xr.DataArray):
-            item_grids = list(set(grids[dim] for dim in result.dims if dim in grids))
-            if len(item_grids) == 0:
-                return result
-            if len(item_grids) > 1:
-                raise RuntimeError("This shouldn't happen. Please open an issue.")
-            return UgridDataArray(result, item_grids[0])
-        elif isinstance(result, xr.Dataset):
-            item_grids = list(set(grids[dim] for dim in result.dims if dim in grids))
-            return UgridDataset(result, item_grids)
-        elif isinstance(result, types.MethodType):
-            return xarray_wrapper(result, self.grids)
-        else:
-            return result
-
-    def _unary_op(
-        self,
-        f: Callable,
-    ):
-        return UgridDataset(self.obj._unary_op(f), self.grids)
-
-    def _binary_op(
-        self,
-        other,
-        f: Callable,
-        reflexive: bool = False,
-    ):
-        if isinstance(other, (UgridDataArray, UgridDataset)):
-            other = other.obj
-        return UgridDataset(self.obj._binary_op(other, f, reflexive), self.grids)
-
-    def _inplace_binary_op(self, other, f: Callable):
-        if isinstance(other, (UgridDataArray, UgridDataset)):
-            other = other.obj
-        return UgridDataset(self.obj._inplace_binary_op(other, f), self.grids)
-
-    def isel(
-        self,
-        indexers=None,
-        drop: bool = False,
-        missing_dims="warn",
-        **indexers_kwargs: Any,
-    ):
-        indexers = either_dict_or_kwargs(indexers, indexers_kwargs, "isel")
-        obj_indexers, ugrid_indexers = filter_indexers(indexers, self.grids)
-        result = self.obj.isel(obj_indexers, drop=drop, missing_dims=missing_dims)
-        result, grids = ugrid_sel(result, ugrid_indexers)
-        return UgridDataset(result, grids)
-
-    def sel(
-        self,
-        indexers=None,
-        method=None,
-        tolerance=None,
-        drop: bool = False,
-        **indexers_kwargs: Any,
-    ):
-        indexers = either_dict_or_kwargs(indexers, indexers_kwargs, "sel")
-        obj_indexers, ugrid_indexers = filter_indexers(indexers, self.grids)
-        result = self.obj.sel(
-            obj_indexers, method=method, tolerance=tolerance, drop=drop
-        )
-        result, grids = ugrid_sel(result, ugrid_indexers)
-        return UgridDataset(result, grids)
-
-    @property
-    def ugrid(self):
-        """
-        UGRID Accessor. This "accessor" makes operations using the UGRID
-        topology available.
-        """
-        return UgridDatasetAccessor(self.obj, self.grids)
-
-    @staticmethod
-    def from_geodataframe(geodataframe: "geopandas.GeoDataFrame"):  # type: ignore # noqa
-        """
-        Convert a geodataframe into the appropriate Ugrid topology and dataset.
-
-        Parameters
-        ----------
-        geodataframe: gpd.GeoDataFrame
-
-        Returns
-        -------
-        dataset: UGridDataset
-        """
-        grid = grid_from_geodataframe(geodataframe)
-        ds = xr.Dataset.from_dataframe(geodataframe.drop("geometry", axis=1))
-        return UgridDataset(ds, [grid])
-
-
-# Set docstrings
-UgridDataset.sel.__doc__ = xr.Dataset.sel.__doc__
-UgridDataset.isel.__doc__ = xr.Dataset.isel.__doc__
-
-
-class UgridDatasetAccessor(AbstractUgridAccessor):
-    def __init__(self, obj: xr.Dataset, grids: Sequence[UgridType]):
-        self.obj = obj
-        self.grids = grids
-
-    @property
-    def grid(self) -> UgridType:
-        ngrid = len(self.grids)
-        if ngrid == 1:
-            return self.grids[0]
-        else:
-            raise AttributeError(
-                "Can only access grid topology via `.grid` if dataset contains "
-                f"exactly one grid. Dataset contains {ngrid} grids. Use "
-                "`.grids` instead."
-            )
-
-    @property
-    def bounds(self) -> Dict[str, Tuple]:
-        """
-        Mapping from grid name to tuple containing ``minx, miny, maxx, maxy``
-        values of the grid's node coordinates for every grid in the dataset.
-        """
-        return {grid.name: grid.bounds for grid in self.grids}
-
-    @property
-    def total_bounds(self) -> Tuple:
-        """
-        Returns a tuple containing ``minx, miny, maxx, maxy`` values for the
-        bounds of the dataset as a whole. Currently does not check whether the
-        coordinate reference systems (CRS) of the grids in the dataset match.
-        """
-        bounds = np.column_stack([bound for bound in self.bounds.values()])
-        return (
-            bounds[0].min(),
-            bounds[1].min(),
-            bounds[2].max(),
-            bounds[3].max(),
-        )
-
-    def assign_node_coords(self) -> UgridDataset:
-        """
-        Assign node coordinates from the grid to the object.
-
-        Returns a new object with all the original data in addition to the new
-        node coordinates of the grid.
-
-        Returns
-        -------
-        assigned: UgridDataset
-        """
-        result = self.obj
-        for grid in self.grids:
-            result = grid.assign_node_coords(result)
-        return UgridDataset(result, self.grids)
-
-    def assign_edge_coords(self) -> UgridDataset:
-        """
-        Assign edge coordinates from the grid to the object.
-
-        Returns a new object with all the original data in addition to the new
-        node coordinates of the grid.
-
-        Returns
-        -------
-        assigned: UgridDataset
-        """
-        result = self.obj
-        for grid in self.grids:
-            result = grid.assign_edge_coords(result)
-        return UgridDataset(result, self.grids)
-
-    def assign_face_coords(self) -> UgridDataset:
-        """
-        Assign face coordinates from the grid to the object.
-
-        Returns a new object with all the original data in addition to the new
-        node coordinates of the grid.
-
-        Returns
-        -------
-        assigned: UgridDataset
-        """
-        result = self.obj
-        for grid in self.grids:
-            if grid.topology_dimension > 1:
-                result = grid.assign_face_coords(result)
-        return UgridDataset(result, self.grids)
-
-    def set_node_coords(self, node_x: str, node_y: str, topology: str = None):
-        """
-        Given names of x and y coordinates of the nodes of an object, set them
-        as the coordinates in the grid.
-
-        Parameters
-        ----------
-        node_x: str
-            Name of the x coordinate of the nodes in the object.
-        node_y: str
-            Name of the y coordinate of the nodes in the object.
-        topology: str, optional
-            Name of the grid topology in which to set the node_x and node_y
-            coordinates. Can be omitted if the UgridDataset contains only a
-            single grid.
-        """
-        if topology is None:
-            if len(self.grids) == 1:
-                grid = self.grids[0]
-            else:
-                raise ValueError(
-                    "topology must be specified when dataset contains multiple datasets"
-                )
-        else:
-            grids = [grid for grid in self.grid if grid.name == topology]
-            if len(grids != 1):
-                raise ValueError(
-                    f"Expected one grid with topology {topology}, found: {len(grids)}"
-                )
-            grid = grids[0]
-
-        grid.set_node_coords(node_x, node_y, self.obj)
-
-    def isel(self, indexers, **indexers_kwargs):
-        """
-        Returns a new object with arrays indexed along edges or faces.
-
-        Parameters
-        ----------
-        indexer: 1d array of integer or bool
-
-        Returns
-        -------
-        indexed: Union[UgridDataArray, UgridDataset]
-        """
-        indexers = either_dict_or_kwargs(indexers, indexers_kwargs, "isel")
-        alldims = set(chain.from_iterable([grid.dimensions for grid in self.grids]))
-        invalid = indexers.keys() - alldims
-        if invalid:
-            raise ValueError(
-                f"Dimensions {invalid} do not exist. Expected one of {alldims}"
-            )
-
-        if len(indexers) > 1:
-            raise NotImplementedError("Can only index a single dimension at a time")
-        dim, indexer = next(iter(indexers.items()))
-
-        # Find grid to use for selection
-        grids = [
-            grid.isel(dim, indexer) if dim in grid.dimensions else grid
-            for grid in self.grids
-        ]
-        result = self.obj.isel({dim: indexer})
-
-        if isinstance(self.obj, xr.Dataset):
-            return UgridDataset(result, grids)
-        else:
-            raise TypeError(f"Expected UgridDataset, got {type(result).__name__}")
-
-    def sel(self, x=None, y=None):
-        result = self.obj
-        for grid in self.grids:
-            result = self._sel(result, grid, x, y)
-        return result
-
-    def sel_points(self, x, y):
-        result = self.obj
-        for grid in self.grids:
-            result = self._sel_points(result, grid, x, y)
-        return result
-
-    def to_dataset(self):
-        """
-        Converts this UgridDataArray or UgridDataset into a standard
-        xarray.Dataset.
-
-        The UGRID topology information is added as standard data variables.
-
-        Returns
-        -------
-        dataset: UgridDataset
-        """
-        return xr.merge([grid.to_dataset(self.obj) for grid in self.grids])
-
-    @property
-    def crs(self):
-        """
-        The Coordinate Reference System (CRS) represented as a ``pyproj.CRS`` object.
-
-        Returns None if the CRS is not set.
-
-        Returns
-        -------
-        crs: dict
-            A dictionary containing the names of the grids and their CRS.
-        """
-        return {grid.name: grid.crs for grid in self.grids}
-
-    def set_crs(
-        self,
-        crs: Union["pyproj.CRS", str] = None,  # type: ignore # noqa
-        epsg: int = None,
-        allow_override: bool = False,
-        topology: str = None,
-    ):
-        """
-        Set the Coordinate Reference System (CRS) of a UGRID topology.
-
-        NOTE: The underlying geometries are not transformed to this CRS. To
-        transform the geometries to a new CRS, use the ``to_crs`` method.
-
-        Parameters
-        ----------
-        crs : pyproj.CRS, optional if `epsg` is specified
-            The value can be anything accepted
-            by :meth:`pyproj.CRS.from_user_input() <pyproj.crs.CRS.from_user_input>`,
-            such as an authority string (eg "EPSG:4326") or a WKT string.
-        epsg : int, optional if `crs` is specified
-            EPSG code specifying the projection.
-        allow_override : bool, default False
-            If the the UGRID topology already has a CRS, allow to replace the
-            existing CRS, even when both are not equal.
-        topology: str, optional
-            Name of the grid topology in which to set the CRS.
-            Sets the CRS for all grids if left unspecified.
-        """
-        if topology is None:
-            grids = self.grids
-        else:
-            names = [grid.name for grid in self.grids]
-            if topology not in names:
-                raise ValueError(f"{topology} not found. Expected one of: {names}")
-            grids = [grid for grid in self.grids if grid.name == topology]
-
-        for grid in grids:
-            grid.set_crs(crs, epsg, allow_override)
-
-    def to_crs(
-        self,
-        crs: Union["pyproj.CRS", str] = None,  # type: ignore # noqa
-        epsg: int = None,
-        topology: str = None,
-    ) -> UgridDataset:
-        """
-        Transform geometries to a new coordinate reference system.
-        Transform all geometries in an active geometry column to a different coordinate
-        reference system. The ``crs`` attribute on the current Ugrid must
-        be set. Either ``crs`` or ``epsg`` may be specified for output.
-
-        This method will transform all points in all objects. It has no notion
-        of projecting the cells. All segments joining points are assumed to be
-        lines in the current projection, not geodesics. Objects crossing the
-        dateline (or other projection boundary) will have undesirable behavior.
-
-        Parameters
-        ----------
-        crs : pyproj.CRS, optional if `epsg` is specified
-            The value can be anything accepted by
-            :meth:`pyproj.CRS.from_user_input() <pyproj.crs.CRS.from_user_input>`,
-            such as an authority string (eg "EPSG:4326") or a WKT string.
-        epsg : int, optional if `crs` is specified
-            EPSG code specifying output projection.
-        topology: str, optional
-            Name of the grid topology to reproject.
-            Reprojects all grids if left unspecified.
-        """
-        if topology is None:
-            grids = [grid.to_crs(crs, epsg) for grid in self.grids]
-        else:
-            names = [grid.name for grid in self.grids]
-            if topology not in names:
-                raise ValueError(f"{topology} not found. Expected one of: {names}")
-            grids = [
-                grid.to_crs(crs, epsg) if grid.name == topology else grid.copy()
-                for grid in self.grids
-            ]
-
-        uds = UgridDataset(self.obj, grids)
-        return uds.ugrid.assign_node_coords()
-
-    def to_geodataframe(
-        self, dim_order=None
-    ) -> "geopandas.GeoDataFrame":  # type: ignore # noqa
-        """
-        Convert data and topology of one facet (node, edge, face) of the grid
-        to a geopandas GeoDataFrame. This also determines the geometry type of
-        the geodataframe:
-
-        * node: point
-        * edge: line
-        * face: polygon
-
-        Parameters
-        ----------
-        name: str
-            Name to give to the array (required if unnamed).
-        dim_order:
-            Hierarchical dimension order for the resulting dataframe. Array content is
-            transposed to this order and then written out as flat vectors in contiguous
-            order, so the last dimension in this list will be contiguous in the resulting
-            DataFrame. This has a major influence on which operations are efficient on the
-            resulting dataframe.
-
-            If provided, must include all dimensions of this DataArray. By default,
-            dimensions are sorted according to the DataArray dimensions order.
-
-        Returns
-        -------
-        geodataframe: gpd.GeoDataFrame
-        """
-        import geopandas as gpd
-
-        ds = self.obj
-
-        gdfs = []
-        for grid in self.grids:
-            if grid.topology_dimension == 1:
-                dim = grid.edge_dimension
-            elif grid.topology_dimension == 2:
-                dim = grid.face_dimension
-            else:
-                raise ValueError("invalid topology dimension on grid")
-
-            variables = [var for var in ds.data_vars if dim in ds[var].dims]
-            # TODO deal with time-dependent data, etc.
-            # Basically requires checking which variables are static, which aren't.
-            # For non-static, requires repeating all geometries.
-            # Call reset_index on multi-index to generate them as regular columns.
-            gdfs.append(
-                gpd.GeoDataFrame(
-                    data=ds[variables].to_dataframe(dim_order=dim_order),
-                    geometry=grid.to_pygeos(dim),
-                )
-            )
-
-        return pd.concat(gdfs)
-
-
-class UgridDataArrayAccessor(AbstractUgridAccessor):
-    """
-    This "accessor" makes operations using the UGRID topology available via the
-    ``.ugrid`` attribute for UgridDataArrays and UgridDatasets.
-    """
-
-    def __init__(self, obj: xr.DataArray, grid: UgridType):
-        self.obj = obj
-        self.grid = grid
-
-    @property
-    def bounds(self) -> Dict[str, Tuple]:
-        """
-        Mapping from grid name to tuple containing ``minx, miny, maxx, maxy``
-        values of the grid's node coordinates.
-        """
-        return {self.grid.name: self.grid.bounds}
-
-    @property
-    def total_bounds(self) -> Tuple:
-        """
-        Returns a tuple containing ``minx, miny, maxx, maxy`` values of the grid's
-        node coordinates.
-        """
-        return next(iter(self.bounds.values()))
-
-    plot = UncachedAccessor(_PlotMethods)
-
-    def assign_node_coords(self) -> UgridDataArray:
-        """
-        Assign node coordinates from the grid to the object.
-
-        Returns a new object with all the original data in addition to the new
-        node coordinates of the grid.
-
-        Returns
-        -------
-        assigned: UgridDataset
-        """
-        return UgridDataArray(self.grid.assign_node_coords(self.obj), self.grid)
-
-    def assign_edge_coords(self) -> UgridDataArray:
-        """
-        Assign edge coordinates from the grid to the object.
-
-        Returns a new object with all the original data in addition to the new
-        node coordinates of the grid.
-
-        Returns
-        -------
-        assigned: UgridDataset
-        """
-        return UgridDataArray(self.grid.assign_edge_coords(self.obj), self.grid)
-
-    def assign_face_coords(self) -> UgridDataArray:
-        """
-        Assign face coordinates from the grid to the object.
-
-        Returns a new object with all the original data in addition to the new
-        node coordinates of the grid.
-
-        Returns
-        -------
-        assigned: UgridDataset
-        """
-        if self.grid.topology_dimension == 1:
-            raise TypeError("Cannot set face coords from a Ugrid1D topology")
-        return UgridDataArray(self.grid.assign_face_coords(self.obj), self.grid)
-
-    def set_node_coords(self, node_x: str, node_y: str):
-        """
-        Given names of x and y coordinates of the nodes of an object, set them
-        as the coordinates in the grid.
-
-        Parameters
-        ----------
-        node_x: str
-            Name of the x coordinate of the nodes in the object.
-        node_y: str
-            Name of the y coordinate of the nodes in the object.
-        """
-        self.grid.set_node_coords(node_x, node_y, self.obj)
-
-    def isel(self, indexers, **indexers_kwargs):
-        """
-        Returns a new object with arrays indexed along edges or faces.
-
-        Parameters
-        ----------
-        indexer: 1d array of integer or bool
-
-        Returns
-        -------
-        indexed: Union[UgridDataArray, UgridDataset]
-        """
-        indexers = either_dict_or_kwargs(indexers, indexers_kwargs, "isel")
-        dims = self.grid.dimensions
-        invalid = indexers.keys() - set(dims)
-        if invalid:
-            raise ValueError(
-                f"Dimensions {invalid} do not exist. Expected one of {dims}"
-            )
-
-        if len(indexers) > 1:
-            raise NotImplementedError("Can only index a single dimension at a time")
-        dim, indexer = next(iter(indexers.items()))
-
-        result = self.obj.isel({dim: indexer})
-        grid = self.grid.isel(dim, indexer)
-        if isinstance(self.obj, xr.DataArray):
-            return UgridDataArray(result, grid)
-        else:
-            raise TypeError(f"Expected UgridDataArray, got {type(result).__name__}")
-
-    def sel(self, x=None, y=None):
-        """
-        Returns a new object, a subselection in the UGRID x and y coordinates.
-
-        The indexing for x and y always occurs orthogonally, i.e.:
-        ``.sel(x=[0.0, 5.0], y=[10.0, 15.0])`` results in a four points. For
-        vectorized indexing (equal to ``zip``ing through x and y), see
-        ``.sel_points``.
-
-        Depending on the nature of the x and y indexers, a xugrid or xarray
-        object is returned:
-
-        * slice without step: ``x=slice(-100, 100)``: returns xugrid object, a
-          part of the unstructured grid.
-        * slice with step: ``x=slice(-100, 100, 10)``: returns xarray object, a
-          series of points (x=[-100, -90, -80, ..., 90, 100]).
-        * a scalar: ``x=5.0``: returns xarray object, a point.
-        * an array: ``x=[1.0, 15.0, 17.0]``: returns xarray object, a series of
-          points.
-
-        Parameters
-        ----------
-        x: float, 1d array, slice
-        y: float, 1d array, slice
-
-        Returns
-        -------
-        selection: Union[UgridDataArray, UgridDataset, xr.DataArray, xr.Dataset]
-        """
-        return self._sel(self.obj, self.grid, x, y)
-
-    def sel_points(self, x, y):
-        """
-        Select points in the unstructured grid.
-
-        Parameters
-        ----------
-        x: ndarray of floats with shape ``(n_points,)``
-        y: ndarray of floats with shape ``(n_points,)``
-
-        Returns
-        -------
-        points: Union[xr.DataArray, xr.Dataset]
-        """
-        return self._sel_points(self.obj, self.grid, x, y)
-
-    def _raster(self, x, y, index) -> xr.DataArray:
-        index = index.ravel()
-        # Cast to float for nodata values (NaN)
-        data = self.obj.isel({self.grid.face_dimension: index}).astype(float).values
-        data[index == -1] = np.nan
-        out = xr.DataArray(
-            data=data.reshape(y.size, x.size),
-            coords={"y": y, "x": x},
-            dims=["y", "x"],
-        )
-        return out
-
-    def rasterize(self, resolution: float):
-        """
-        Rasterize unstructured grid by sampling.
-
-        Parameters
-        ----------
-        resolution: float
-            Spacing in x and y.
-
-        Returns
-        -------
-        rasterized: Union[xr.DataArray, xr.Dataset]
-        """
-        x, y, index = self.grid.rasterize(resolution)
-        return self._raster(x, y, index)
-
-    def rasterize_like(self, other: Union[xr.DataArray, xr.Dataset]):
-        """
-        Rasterize unstructured grid by sampling on the x and y coordinates
-        of ``other``.
-
-        Parameters
-        ----------
-        resolution: float
-            Spacing in x and y.
-        other: Union[xr.DataArray, xr.Dataset]
-            Object to take x and y coordinates from.
-
-        Returns
-        -------
-        rasterized: Union[xr.DataArray, xr.Dataset]
-        """
-        x, y, index = self.grid.rasterize_like(
-            x=other["x"].values,
-            y=other["y"].values,
-        )
-        return self._raster(x, y, index)
-
-    @property
-    def crs(self):
-        """
-        The Coordinate Reference System (CRS) represented as a ``pyproj.CRS`` object.
-
-        Returns None if the CRS is not set.
-
-        Returns
-        -------
-        crs: dict
-            A dictionary containing the name of the grid and its CRS.
-        """
-        return {self.grid.name: self.grid.crs}
-
-    def set_crs(
-        self,
-        crs: Union["pyproj.CRS", str] = None,  # type: ignore # noqa
-        epsg: int = None,
-        allow_override: bool = False,
-    ):
-        """
-        Set the Coordinate Reference System (CRS) of a UGRID topology.
-
-        NOTE: The underlying geometries are not transformed to this CRS. To
-        transform the geometries to a new CRS, use the ``to_crs`` method.
-
-        Parameters
-        ----------
-        crs : pyproj.CRS, optional if `epsg` is specified
-            The value can be anything accepted
-            by :meth:`pyproj.CRS.from_user_input() <pyproj.crs.CRS.from_user_input>`,
-            such as an authority string (eg "EPSG:4326") or a WKT string.
-        epsg : int, optional if `crs` is specified
-            EPSG code specifying the projection.
-        allow_override : bool, default False
-            If the the UGRID topology already has a CRS, allow to replace the
-            existing CRS, even when both are not equal.
-        """
-        self.grid.set_crs(crs, epsg, allow_override)
-
-    def to_crs(
-        self,
-        crs: Union["pyproj.CRS", str] = None,  # type: ignore # noqa
-        epsg: int = None,
-    ) -> UgridDataArray:
-        """
-        Transform geometries to a new coordinate reference system.
-        Transform all geometries in an active geometry column to a different coordinate
-        reference system. The ``crs`` attribute on the current Ugrid must
-        be set. Either ``crs`` or ``epsg`` may be specified for output.
-
-        This method will transform all points in all objects. It has no notion
-        of projecting the cells. All segments joining points are assumed to be
-        lines in the current projection, not geodesics. Objects crossing the
-        dateline (or other projection boundary) will have undesirable behavior.
-
-        Parameters
-        ----------
-        crs : pyproj.CRS, optional if `epsg` is specified
-            The value can be anything accepted by
-            :meth:`pyproj.CRS.from_user_input() <pyproj.crs.CRS.from_user_input>`,
-            such as an authority string (eg "EPSG:4326") or a WKT string.
-        epsg : int, optional if `crs` is specified
-            EPSG code specifying output projection.
-        """
-        uda = UgridDataArray(self.obj, self.grid.to_crs(crs, epsg))
-        if self.grid.node_dimension in self.obj.dims:
-            return uda.ugrid.assign_node_coords()
-        else:
-            return uda
-
-    def to_geodataframe(
-        self, name: str = None, dim_order=None
-    ) -> "geopandas.GeoDataFrame":  # type: ignore # noqa
-        """
-        Convert data and topology of one facet (node, edge, face) of the grid
-        to a geopandas GeoDataFrame. This also determines the geometry type of
-        the geodataframe:
-
-        * node: point
-        * edge: line
-        * face: polygon
-
-        Parameters
-        ----------
-        dim: str
-            node, edge, or face dimension. Inferred for DataArray.
-        name: str
-            Name to give to the array (required if unnamed).
-        dim_order:
-            Hierarchical dimension order for the resulting dataframe. Array content is
-            transposed to this order and then written out as flat vectors in contiguous
-            order, so the last dimension in this list will be contiguous in the resulting
-            DataFrame. This has a major influence on which operations are efficient on the
-            resulting dataframe.
-
-            If provided, must include all dimensions of this DataArray. By default,
-            dimensions are sorted according to the DataArray dimensions order.
-
-        Returns
-        -------
-        geodataframe: gpd.GeoDataFrame
-        """
-        import geopandas as gpd
-
-        dim = self.obj.dims[-1]
-        if name is not None:
-            ds = self.obj.to_dataset(name=name)
-        else:
-            ds = self.obj.to_dataset()
-
-        variables = [var for var in ds.data_vars if dim in ds[var].dims]
-        # TODO deal with time-dependent data, etc.
-        # Basically requires checking which variables are static, which aren't.
-        # For non-static, requires repeating all geometries.
-        # Call reset_index on multi-index to generate them as regular columns.
-        df = ds[variables].to_dataframe(dim_order=dim_order)
-        geometry = self.grid.to_pygeos(dim)
-        return gpd.GeoDataFrame(df, geometry=geometry)
-
-    def _binary_iterate(self, iterations: int, mask, value, border_value):
-        if border_value == value:
-            exterior = self.grid.exterior_faces
-        else:
-            exterior = None
-        if mask is not None:
-            mask = mask.values
-
-        obj = self.obj
-        if isinstance(obj, xr.DataArray):
-            output = connectivity._binary_iterate(
-                self.grid.face_face_connectivity,
-                obj.values,
-                value,
-                iterations,
-                mask,
-                exterior,
-                border_value,
-            )
-            da = obj.copy(data=output)
-            return UgridDataArray(da, self.grid.copy())
-        elif isinstance(obj, xr.Dataset):
-            raise NotImplementedError
-        else:
-            raise ValueError("object should be a xr.DataArray")
-
-    def binary_dilation(
-        self,
-        iterations: int = 1,
-        mask=None,
-        border_value=False,
-    ):
-        """
-        Binary dilation can be used on a boolean array to expand the "shape" of
-        features.
-
-        Compare with :py:func:`scipy.ndimage.binary_dilation`.
-
-        Parameters
-        ----------
-        iterations: int, default: 1
-        mask: 1d array of bool, optional
-        border_value: bool, default value: False
-
-        Returns
-        -------
-        dilated: UgridDataArray
-        """
-        return self._binary_iterate(iterations, mask, True, border_value)
-
-    def binary_erosion(
-        self,
-        iterations: int = 1,
-        mask=None,
-        border_value=False,
-    ):
-        """
-        Binary erosion can be used on a boolean array to shrink the "shape" of
-        features.
-
-        Compare with :py:func:`scipy.ndimage.binary_erosion`.
-
-        Parameters
-        ----------
-        iterations: int, default: 1
-        mask: 1d array of bool, optional
-        border_value: bool, default value: False
-
-        Returns
-        -------
-        eroded: UgridDataArray
-        """
-        return self._binary_iterate(iterations, mask, False, border_value)
-
-    def connected_components(self):
-        """
-        Every edge or face is given a component number. If all are connected,
-        all will have the same number.
-
-        Wraps :py:func:`scipy.sparse.csgraph.connected_components``.
-
-        Returns
-        -------
-        labelled: UgridDataArray
-        """
-        _, labels = scipy.sparse.csgraph.connected_components(
-            self.grid.face_face_connectivity
-        )
-        return UgridDataArray(
-            xr.DataArray(labels, dims=[self.grid.face_dimension]),
-            self.grid,
-        )
-
-    def reverse_cuthill_mckee(self):
-        """
-        Reduces bandwith of the connectivity matrix.
-
-        Wraps :py:func:`scipy.sparse.csgraph.reverse_cuthill_mckee`.
-
-        Returns
-        -------
-        reordered: Union[UgridDataArray, UgridDataset]
-        """
-        grid = self.grid
-        reordered_grid, reordering = self.grid.reverse_cuthill_mckee()
-        reordered_data = self.obj.isel({grid.face_dimension: reordering})
-        # TODO: this might not work properly if e.g. centroids are stored in obj.
-        # Not all metadata would be reordered.
-        return UgridDataArray(
-            reordered_data,
-            reordered_grid,
-        )
-
-    def laplace_interpolate(
-        self,
-        xy_weights: bool = True,
-        direct_solve: bool = False,
-        drop_tol: float = None,
-        fill_factor: float = None,
-        drop_rule: str = None,
-        options: dict = None,
-        tol: float = 1.0e-5,
-        maxiter: int = 250,
-    ):
-        """
-        Fill gaps in ``data`` (``np.nan`` values) using Laplace interpolation.
-
-        This solves Laplace's equation where where there is no data, with data
-        values functioning as fixed potential boundary conditions.
-
-        Note that an iterative solver method will be required for large grids.
-        In this case, some experimentation with the solver settings may be
-        required to find a converging solution of sufficient accuracy. Refer to
-        the documentation of :py:func:`scipy.sparse.linalg.spilu` and
-        :py:func:`scipy.sparse.linalg.cg`.
-
-        Parameters
-        ----------
-        xy_weights: bool, default False.
-            Wether to use the inverse of the centroid to centroid distance in
-            the coefficient matrix. If ``False``, defaults to uniform
-            coefficients of 1 so that each face connection has equal weight.
-        direct_solve: bool, optional, default ``False``
-            Whether to use a direct or an iterative solver or a conjugate gradient
-            solver. Direct method provides an exact answer, but are unsuitable
-            for large problems.
-        drop_tol: float, optional, default None.
-            Drop tolerance for ``scipy.sparse.linalg.spilu`` which functions as a
-            preconditioner for the conjugate gradient solver.
-        fill_factor: float, optional, default None.
-            Fill factor for ``scipy.sparse.linalg.spilu``.
-        drop_rule: str, optional default None.
-            Drop rule for ``scipy.sparse.linalg.spilu``.
-        options: dict, optional, default None.
-            Remaining other options for ``scipy.sparse.linalg.spilu``.
-        tol: float, optional, default 1.0e-5.
-            Convergence tolerance for ``scipy.sparse.linalg.cg``.
-        maxiter: int, default 250.
-            Maximum number of iterations for ``scipy.sparse.linalg.cg``.
-
-        Returns
-        -------
-        filled: UgridDataArray of floats
-        """
-        grid = self.grid
-        da = self.obj
-        if grid.topology_dimension != 2:
-            raise NotImplementedError
-        if len(da.dims) > 1 or da.dims[0] != grid.face_dimension:
-            raise NotImplementedError
-
-        connectivity = grid.face_face_connectivity.copy()
-        if xy_weights:
-            xy = grid.centroids
-            coo = connectivity.tocoo()
-            i = coo.row
-            j = coo.col
-            connectivity.data = 1.0 / np.linalg.norm(xy[j] - xy[i], axis=1)
-
-        filled = laplace_interpolate(
-            connectivity=connectivity,
-            data=da.values,
-            use_weights=xy_weights,
-            direct_solve=direct_solve,
-            drop_tol=drop_tol,
-            fill_factor=fill_factor,
-            drop_rule=drop_rule,
-            options=options,
-            tol=tol,
-            maxiter=maxiter,
-        )
-        da_filled = da.copy(data=filled)
-        return UgridDataArray(da_filled, grid)
-
-    def to_dataset(self):
-        """
-        Converts this UgridDataArray or UgridDataset into a standard
-        xarray.Dataset.
-
-        The UGRID topology information is added as standard data variables.
-
-        Returns
-        -------
-        dataset: UgridDataset
-        """
-        return self.grid.to_dataset(self.obj)
-
-
-# Wrapped IO methods
-# ------------------
-
-
-def open_dataset(*args, **kwargs):
-    ds = xr.open_dataset(*args, **kwargs)
-    return UgridDataset(ds)
-
-
-def open_dataarray(*args, **kwargs):
-    ds = xr.open_dataset(*args, **kwargs)
-    dataset = UgridDataset(ds)
-
-    if len(dataset.data_vars) != 1:
-        raise ValueError(
-            "Given file dataset contains more than one data "
-            "variable. Please read with xarray.open_dataset and "
-            "then select the variable you want."
-        )
-    else:
-        (data_array,) = dataset.data_vars.values()
-
-    data_array.set_close(dataset._close)
-
-    # Reset names if they were changed during saving
-    # to ensure that we can 'roundtrip' perfectly
-    if DATAARRAY_NAME in dataset.attrs:
-        data_array.name = dataset.attrs[DATAARRAY_NAME]
-        del dataset.attrs[DATAARRAY_NAME]
-
-    if data_array.name == DATAARRAY_VARIABLE:
-        data_array.name = None
-
-    return data_array
-
-
-def open_mfdataset(*args, **kwargs):
-    if "data_vars" in kwargs:
-        raise ValueError("data_vars kwargs is not supported in xugrid.open_mfdataset")
-    kwargs["data_vars"] = "minimal"
-    ds = xr.open_mfdataset(*args, **kwargs)
-    return UgridDataset(ds)
-
-
-def open_zarr(*args, **kwargs):
-    ds = xr.open_zarr(*args, **kwargs)
-    return UgridDataset(ds)
-
-
-open_dataset.__doc__ = xr.open_dataset.__doc__
-open_dataarray.__doc__ = xr.open_dataarray.__doc__
-open_mfdataset.__doc__ = xr.open_mfdataset.__doc__
-open_zarr.__doc__ = xr.open_zarr.__doc__
-
-
-# Other utilities
-# ---------------
-
-
-def wrap_func_like(func):
-    @wraps(func)
-    def _like(other, *args, **kwargs):
-        obj = func(other.obj, *args, **kwargs)
-        if isinstance(obj, xr.DataArray):
-            return type(other)(obj, other.grid)
-        elif isinstance(obj, xr.Dataset):
-            return type(other)(obj, other.grids)
-        else:
-            raise TypeError(
-                f"Expected Dataset or DataArray, received {type(other).__name__}"
-            )
-
-    _like.__doc__ = func.__doc__
-    return _like
-
-
-def wrap_func_objects(func):
-    @wraps(func)
-    def _f(objects, *args, **kwargs):
-        grids = []
-        bare_objs = []
-        for obj in objects:
-            if isinstance(obj, UgridDataArray):
-                grids.append(obj.grid)
-            elif isinstance(obj, UgridDataset):
-                grids.extend(obj.grids)
-            else:
-                raise TypeError(
-                    "Can only concatenate xugrid UgridDataset and UgridDataArray "
-                    f"objects, got {type(obj).__name__}"
-                )
-
-            bare_objs.append(obj.obj)
-
-        grids = set(grids)
-        result = func(bare_objs, *args, **kwargs)
-        if isinstance(result, xr.DataArray):
-            if len(grids) > 1:
-                raise ValueError("All UgridDataArrays must have the same grid")
-            return UgridDataArray(result, next(iter(grids)))
-        else:
-            return UgridDataset(result, grids)
-
-    _f.__doc__ = func.__doc__
-    return _f
-
-
-full_like = wrap_func_like(xr.full_like)
-zeros_like = wrap_func_like(xr.zeros_like)
-ones_like = wrap_func_like(xr.ones_like)
-
-concat = wrap_func_objects(xr.concat)
-merge = wrap_func_objects(xr.merge)


Test patch:
diff --git a/tests/test_connectivity.py b/tests/test_connectivity.py
index c3ae9ef6c..0cbbd5b98 100644
--- a/tests/test_connectivity.py
+++ b/tests/test_connectivity.py
@@ -2,7 +2,7 @@
 import pytest
 from scipy import sparse
 
-from xugrid import connectivity
+from xugrid.ugrid import connectivity
 
 
 @pytest.fixture(scope="function")
diff --git a/tests/test_interpolate.py b/tests/test_interpolate.py
index 7a09e3a28..78cc98adb 100644
--- a/tests/test_interpolate.py
+++ b/tests/test_interpolate.py
@@ -2,7 +2,7 @@
 import pytest
 from scipy import sparse
 
-from xugrid import interpolate
+from xugrid.ugrid import interpolate
 
 
 def test_laplace_interpolate():
diff --git a/tests/test_snap.py b/tests/test_snap.py
index 35c7114e3..a06e549f4 100644
--- a/tests/test_snap.py
+++ b/tests/test_snap.py
@@ -9,7 +9,7 @@
 import shapely.geometry as sg
 import xarray as xr
 
-from xugrid.snapping import snap_nodes, snap_to_nodes, snap_to_structured_grid
+from xugrid.ugrid.snapping import snap_nodes, snap_to_nodes, snap_to_structured_grid
 
 
 def test_snap__three_points():
diff --git a/tests/test_ugrid1d.py b/tests/test_ugrid1d.py
index 41168caeb..5e80ac9b1 100644
--- a/tests/test_ugrid1d.py
+++ b/tests/test_ugrid1d.py
@@ -227,7 +227,10 @@ def test_dimensions():
     grid = grid1d()
     assert grid.node_dimension == f"{NAME}_nNodes"
     assert grid.edge_dimension == f"{NAME}_nEdges"
-    assert grid.dimensions == (f"{NAME}_nNodes", f"{NAME}_nEdges")
+    assert grid.dimensions == {
+        f"{NAME}_nNodes": 3,
+        f"{NAME}_nEdges": 2,
+    }
 
 
 @requires_meshkernel
@@ -262,17 +265,26 @@ def test_to_pygeos():
 
 def test_sel():
     grid = grid1d()
+    obj = xr.DataArray(
+        data=[0, 1],
+        dims=[grid.edge_dimension],
+    )
     with pytest.raises(ValueError, match="Ugrid1d only supports slice indexing"):
-        grid.sel(x=1.0, y=1.0)
+        grid.sel(obj=obj, x=1.0, y=1.0)
     with pytest.raises(ValueError, match="Ugrid1d does not support steps"):
-        grid.sel(x=slice(0, 2, 1), y=slice(0, 2, 1))
+        grid.sel(obj=obj, x=slice(0, 2, 1), y=slice(0, 2, 1))
     with pytest.raises(ValueError, match="slice start should be smaller"):
-        grid.sel(x=slice(2, 0), y=slice(0, 2))
-    dim, as_ugrid, index, coords = grid.sel(x=slice(0, 1), y=slice(0, 1))
-    assert dim == f"{NAME}_nEdges"
-    assert as_ugrid
-    assert np.allclose(index, [0])
-    assert coords == {}
+        grid.sel(obj, x=slice(2, 0), y=slice(0, 2))
+
+    actual = grid.sel(obj=obj, x=slice(0, 1), y=slice(0, 1))
+    assert isinstance(actual, tuple)
+    new_obj, new_grid = actual
+
+    assert isinstance(new_obj, xr.DataArray)
+    assert isinstance(new_grid, xugrid.Ugrid1d)
+    assert new_obj.dims[0] == f"{NAME}_nEdges"
+    assert new_grid.edge_dimension == f"{NAME}_nEdges"
+    assert np.array_equal(new_obj.values, [0])
 
 
 def test_topology_subset():
diff --git a/tests/test_ugrid2d.py b/tests/test_ugrid2d.py
index 45873b8a1..17afd472f 100644
--- a/tests/test_ugrid2d.py
+++ b/tests/test_ugrid2d.py
@@ -154,8 +154,6 @@ def test_ugrid2d_edge_bounds():
         ]
     )
     actual = grid.edge_bounds
-    print(actual)
-    print(expected)
     assert actual.shape == (10, 4)
     assert np.allclose(actual, expected)
 
@@ -332,7 +330,11 @@ def test_dimensions():
     assert grid.node_dimension == f"{NAME}_nNodes"
     assert grid.edge_dimension == f"{NAME}_nEdges"
     assert grid.face_dimension == f"{NAME}_nFaces"
-    assert grid.dimensions == (f"{NAME}_nNodes", f"{NAME}_nEdges", f"{NAME}_nFaces")
+    assert grid.dimensions == {
+        f"{NAME}_nNodes": 7,
+        f"{NAME}_nEdges": 10,
+        f"{NAME}_nFaces": 4,
+    }
 
 
 def test_edge_node_connectivity():
@@ -482,197 +484,217 @@ def test_rasterize():
     assert np.array_equal(index, expected_index)
 
 
-def test_sel_points():
-    grid = grid2d()
-    x = [0.5, 1.5]
-    y = [0.5, 1.25]
-
-    with pytest.raises(ValueError, match="shape of x does not match shape of y"):
-        grid.sel_points(x=[0.5, 1.5], y=[0.5])
-    with pytest.raises(ValueError, match="x and y must be 1d"):
-        grid.sel_points(x=[x], y=[y])
-
-    dim, as_ugrid, index, coords = grid.sel_points(x=x, y=y)
-    assert dim == f"{NAME}_nFaces"
-    assert not as_ugrid
-    assert np.array_equal(index, [0, 3])
-    assert coords["x"][0] == dim
-    assert coords["y"][0] == dim
-    assert np.array_equal(coords["x"][1], x)
-    assert np.array_equal(coords["y"][1], y)
-
-
-def test_validate_indexer():
-    grid = grid2d()
-    with pytest.raises(ValueError, match="slice stop should be larger than"):
-        grid._validate_indexer(slice(2, 0))
-    with pytest.raises(ValueError, match="step should be None"):
-        grid._validate_indexer(slice(None, 2, 1))
-    with pytest.raises(ValueError, match="step should be None"):
-        grid._validate_indexer(slice(0, None, 1))
-
-    expected = np.arange(0.0, 2.0, 0.5)
-    assert np.allclose(grid._validate_indexer(slice(0, 2, 0.5)), expected)
-    assert grid._validate_indexer(slice(None, 2)) == slice(None, 2)
-    assert grid._validate_indexer(slice(0, None)) == slice(0, None)
-
-    with pytest.raises(TypeError, match="Invalid indexer type"):
-        grid._validate_indexer((0, 1, 2))
-
-    # list
-    actual = grid._validate_indexer([0.0, 1.0, 2.0])
-    assert isinstance(actual, np.ndarray)
-    assert np.allclose(actual, [0.0, 1.0, 2.0])
-
-    # numpy array
-    actual = grid._validate_indexer(np.array([0.0, 1.0, 2.0]))
-    assert isinstance(actual, np.ndarray)
-    assert np.allclose(actual, [0.0, 1.0, 2.0])
-
-    # xarray DataArray
-    indexer = xr.DataArray([0.0, 1.0, 2.0], {"x": [0, 1, 2]}, ["x"])
-    actual = grid._validate_indexer(indexer)
-    assert isinstance(actual, np.ndarray)
-    assert np.allclose(actual, [0.0, 1.0, 2.0])
-
-    # float
-    actual = grid._validate_indexer(1.0)
-    assert isinstance(actual, np.ndarray)
-    assert np.allclose(actual, [1.0])
-
-    # int
-    actual = grid._validate_indexer(1)
-    assert isinstance(actual, np.ndarray)
-    assert np.allclose(actual, [1])
-
-
-def test_sel__bounding_box():
-    grid = grid2d()
-
-    def check_output(dim, as_ugrid, index, coords, expected_index):
-        assert dim == f"{NAME}_nFaces"
-        assert as_ugrid
-        assert isinstance(index, np.ndarray)
-        assert np.allclose(index, expected_index)
-        assert coords == {}
-
-    dim, as_ugrid, index, coords = grid.sel(x=slice(0.0, 2.0), y=slice(0.0, 1.0))
-    check_output(dim, as_ugrid, index, coords, [0, 1])
-
-    dim, as_ugrid, index, coords = grid.sel(x=slice(None, None), y=slice(None, 1.0))
-    check_output(dim, as_ugrid, index, coords, [0, 1])
-
-    dim, as_ugrid, index, coords = grid.sel(x=slice(0.0, 1.0), y=slice(0.0, 2.0))
-    check_output(dim, as_ugrid, index, coords, [0, 2])
-
-    dim, as_ugrid, index, coords = grid.sel(x=slice(None, 1.0), y=slice(None, None))
-    check_output(dim, as_ugrid, index, coords, [0, 2])
-
-    for x, y in zip([None, None, slice(0, 2)], [None, slice(0, 2), None]):
-        dim, as_ugrid, index, coords = grid.sel(x=x, y=y)
-        check_output(dim, as_ugrid, index, coords, [0, 1, 2, 3])
-
-    # Check default arguments, should return entire grid
-    dim, as_ugrid, index, coords = grid.sel()
-    check_output(dim, as_ugrid, index, coords, [0, 1, 2, 3])
-
-
-def test_sel__points_from_scalar():
-    grid = grid2d()
-
-    def check_output(dim, as_ugrid, index, coords):
-        assert dim == f"{NAME}_nFaces"
-        assert not as_ugrid
-        assert np.allclose(index, [0])
-        assert coords["x"][0] == dim
-        assert coords["y"][0] == dim
-        assert np.allclose(coords["x"][1], [0.5])
-        assert np.allclose(coords["y"][1], [0.5])
-
-    dim, as_ugrid, index, coords = grid.sel(x=0.5, y=0.5)
-    check_output(dim, as_ugrid, index, coords)
-
-    dim, as_ugrid, index, coords = grid.sel(x=[0.5], y=[0.5])
-    check_output(dim, as_ugrid, index, coords)
-
-    with pytest.raises(TypeError, match="Invalid indexer type"):
-        grid.sel(x=(0.5,), y=[0.5])
-
-
-def test_sel__points_from_arrays_and_slice():
-    grid = grid2d()
-
-    def check_output(dim, as_ugrid, index, coords):
-        assert dim == f"{NAME}_nFaces"
-        assert not as_ugrid
-        assert coords["x"][0] == dim
-        assert coords["y"][0] == dim
-        assert np.allclose(index, [0, 0, 1, 2, 2, 3])
-        assert np.allclose(coords["x"][1], [0.4, 0.8, 1.2, 0.4, 0.8, 1.2])
-        assert np.allclose(coords["y"][1], [0.5, 0.5, 0.5, 1.1, 1.1, 1.1])
-
-    x = [0.4, 0.8, 1.2]
-    y = [0.5, 1.1]
-    dim, as_ugrid, index, coords = grid.sel(x=x, y=y)
-    check_output(dim, as_ugrid, index, coords)
-
-    x = slice(0.4, 1.5, 0.4)  # Evaluates to: [0.4, 0.8, 1.2]
-    dim, as_ugrid, index, coords = grid.sel(x=x, y=y)
-    check_output(dim, as_ugrid, index, coords)
-
-
-def test_sel__edges_from_slice():
-    grid = grid2d()
-
-    with pytest.raises(ValueError, match="If x is a slice without steps"):
-        grid.sel(x=slice(None, None), y=[0.25, 0.75])
-    with pytest.raises(ValueError, match="If x is a slice without steps"):
-        grid.sel(x=slice(None, None), y=slice(0.25, 1.0, 0.25))
-    with pytest.raises(ValueError, match="If y is a slice without steps"):
-        grid.sel(x=[0.25, 0.75], y=slice(None, None))
-
-    dim, as_ugrid, index, coords = grid.sel(x=slice(None, None), y=0.5)
-    assert dim == f"{NAME}_nFaces"
-    assert not as_ugrid
-    assert coords["x"][0] == dim
-    assert coords["y"][0] == dim
-    assert np.allclose(index, [0, 1])
-    assert np.allclose(coords["x"][1], [0.5, 1.5])
-    assert np.allclose(coords["y"][1], [0.5, 0.5])
-
-    dim, as_ugrid, index, coords = grid.sel(x=0.5, y=slice(None, None))
-    assert dim == f"{NAME}_nFaces"
-    assert not as_ugrid
-    assert coords["x"][0] == dim
-    assert coords["y"][0] == dim
-    assert np.allclose(index, [0, 2])
-    assert np.allclose(coords["x"][1], [0.5, 0.5])
-    assert np.allclose(coords["y"][1], [0.5, 1.25])
+class TestUgrid2dSelection:
+    @pytest.fixture(autouse=True)
+    def setup(self):
+        self.grid = grid2d()
+        self.obj = xr.DataArray([0, 1, 2, 3], dims=[self.grid.face_dimension])
+
+    def test_sel_points(self):
+        x = [0.5, 1.5]
+        y = [0.5, 1.25]
+
+        with pytest.raises(ValueError, match="shape of x does not match shape of y"):
+            self.grid.sel_points(obj=self.obj, x=[0.5, 1.5], y=[0.5])
+        with pytest.raises(ValueError, match="x and y must be 1d"):
+            self.grid.sel_points(obj=self.obj, x=[x], y=[y])
+
+        actual = self.grid.sel_points(obj=self.obj, x=x, y=y)
+        assert isinstance(actual, xr.DataArray)
+
+        dim = f"{NAME}_nFaces"
+        expected = xr.DataArray(
+            data=[0, 3],
+            coords={
+                "x": (dim, x),
+                "y": (dim, y),
+            },
+            dims=[dim],
+        )
+        assert expected.equals(actual)
+
+    def test_validate_indexer(self):
+        with pytest.raises(ValueError, match="slice stop should be larger than"):
+            self.grid._validate_indexer(slice(2, 0))
+        with pytest.raises(ValueError, match="step should be None"):
+            self.grid._validate_indexer(slice(None, 2, 1))
+        with pytest.raises(ValueError, match="step should be None"):
+            self.grid._validate_indexer(slice(0, None, 1))
+
+        expected = np.arange(0.0, 2.0, 0.5)
+        assert np.allclose(self.grid._validate_indexer(slice(0, 2, 0.5)), expected)
+        assert self.grid._validate_indexer(slice(None, 2)) == slice(None, 2)
+        assert self.grid._validate_indexer(slice(0, None)) == slice(0, None)
+
+        with pytest.raises(TypeError, match="Invalid indexer type"):
+            self.grid._validate_indexer((0, 1, 2))
+
+        # list
+        actual = self.grid._validate_indexer([0.0, 1.0, 2.0])
+        assert isinstance(actual, np.ndarray)
+        assert np.allclose(actual, [0.0, 1.0, 2.0])
+
+        # numpy array
+        actual = self.grid._validate_indexer(np.array([0.0, 1.0, 2.0]))
+        assert isinstance(actual, np.ndarray)
+        assert np.allclose(actual, [0.0, 1.0, 2.0])
+
+        # xarray DataArray
+        indexer = xr.DataArray([0.0, 1.0, 2.0], {"x": [0, 1, 2]}, ["x"])
+        actual = self.grid._validate_indexer(indexer)
+        assert isinstance(actual, np.ndarray)
+        assert np.allclose(actual, [0.0, 1.0, 2.0])
+
+        # float
+        actual = self.grid._validate_indexer(1.0)
+        assert isinstance(actual, np.ndarray)
+        assert np.allclose(actual, [1.0])
+
+        # int
+        actual = self.grid._validate_indexer(1)
+        assert isinstance(actual, np.ndarray)
+        assert np.allclose(actual, [1])
+
+    def test_sel__bounding_box(self):
+        def check_output(actual, expected):
+            assert isinstance(actual, tuple)
+            new_obj, new_grid = actual
+            assert isinstance(new_obj, xr.DataArray)
+            assert isinstance(new_grid, xugrid.Ugrid2d)
+            assert new_obj.dims[0] == f"{NAME}_nFaces"
+            assert new_grid.face_dimension == f"{NAME}_nFaces"
+            assert np.array_equal(new_obj.values, expected)
+
+        actual = self.grid.sel(obj=self.obj, x=slice(0.0, 2.0), y=slice(0.0, 1.0))
+        check_output(actual, [0, 1])
+
+        actual = self.grid.sel(obj=self.obj, x=slice(None, None), y=slice(None, 1.0))
+        check_output(actual, [0, 1])
+
+        actual = self.grid.sel(obj=self.obj, x=slice(0.0, 1.0), y=slice(0.0, 2.0))
+        check_output(actual, [0, 2])
+
+        actual = self.grid.sel(obj=self.obj, x=slice(None, 1.0), y=slice(None, None))
+        check_output(actual, [0, 2])
+
+        for x, y in zip([None, None, slice(0, 2)], [None, slice(0, 2), None]):
+            actual = self.grid.sel(obj=self.obj, x=x, y=y)
+            check_output(actual, [0, 1, 2, 3])
+
+        # Check default arguments, should return entire grid
+        actual = self.grid.sel(obj=self.obj)
+        check_output(actual, [0, 1, 2, 3])
+
+    def test_sel__points_from_scalar(self):
+        def check_output(actual):
+            assert isinstance(actual, xr.DataArray)
+            dim = f"{NAME}_nFaces"
+            expected = xr.DataArray(
+                data=[0],
+                coords={
+                    "x": (dim, [0.5]),
+                    "y": (dim, [0.5]),
+                },
+                dims=[dim],
+            )
+            assert expected.equals(actual)
+
+        actual = self.grid.sel(obj=self.obj, x=0.5, y=0.5)
+        check_output(actual)
+
+        actual = self.grid.sel(obj=self.obj, x=[0.5], y=[0.5])
+        check_output(actual)
+
+        with pytest.raises(TypeError, match="Invalid indexer type"):
+            self.grid.sel(obj=self.obj, x=(0.5,), y=[0.5])
+
+    def test_sel__points_from_arrays_and_slice(self):
+        def check_output(actual):
+            assert isinstance(actual, xr.DataArray)
+            dim = f"{NAME}_nFaces"
+            expected = xr.DataArray(
+                data=[0, 0, 1, 2, 2, 3],
+                coords={
+                    "x": (dim, [0.4, 0.8, 1.2, 0.4, 0.8, 1.2]),
+                    "y": (dim, [0.5, 0.5, 0.5, 1.1, 1.1, 1.1]),
+                },
+                dims=[dim],
+            )
+            # This fails for some reason:
+            # assert expected.equals(actual)
+            assert np.array_equal(expected.values, actual.values)
+            assert expected.dims == actual.dims
+            assert np.allclose(expected["y"].values, actual["y"].values)
+            assert np.allclose(expected["x"].values, actual["x"].values)
+
+        x = [0.4, 0.8, 1.2]
+        y = [0.5, 1.1]
+        actual = self.grid.sel(obj=self.obj, x=x, y=y)
+        check_output(actual)
+
+        x = slice(0.4, 1.5, 0.4)  # Evaluates to: [0.4, 0.8, 1.2]
+        actual = self.grid.sel(obj=self.obj, x=x, y=y)
+        check_output(actual)
+
+    def test_sel__edges_from_slice(self):
+        with pytest.raises(ValueError, match="If x is a slice without steps"):
+            self.grid.sel(obj=self.obj, x=slice(None, None), y=[0.25, 0.75])
+        with pytest.raises(ValueError, match="If x is a slice without steps"):
+            self.grid.sel(obj=self.obj, x=slice(None, None), y=slice(0.25, 1.0, 0.25))
+        with pytest.raises(ValueError, match="If y is a slice without steps"):
+            self.grid.sel(obj=self.obj, x=[0.25, 0.75], y=slice(None, None))
+
+        actual = self.grid.sel(obj=self.obj, x=slice(None, None), y=0.5)
+        assert isinstance(actual, xr.DataArray)
+        dim = f"{NAME}_nFaces"
+        expected = xr.DataArray(
+            data=[0, 1],
+            coords={
+                "x": (dim, [0.5, 1.5]),
+                "y": (dim, [0.5, 0.5]),
+                "s": (dim, [0.5, 1.5]),
+            },
+            dims=[dim],
+        )
+        assert expected.equals(actual)
+
+        actual = self.grid.sel(obj=self.obj, x=0.5, y=slice(None, None))
+        assert isinstance(actual, xr.DataArray)
+        expected = xr.DataArray(
+            data=[0, 2],
+            coords={
+                "x": (dim, [0.5, 0.5]),
+                "y": (dim, [0.5, 1.25]),
+                "s": (dim, [0.5, 1.25]),
+            },
+            dims=[dim],
+        )
+        assert expected.equals(actual)
 
 
 def test_topology_subset():
     grid = grid2d()
-    edge_indices = np.array([1])
-    actual = grid.topology_subset(edge_indices)
+    face_index = np.array([1])
+    actual = grid.topology_subset(face_index)
     assert np.array_equal(actual.face_node_connectivity, [[0, 1, 3, 2]])
     assert np.array_equal(actual.node_x, [1.0, 2.0, 1.0, 2.0])
     assert np.array_equal(actual.node_y, [0.0, 0.0, 1.0, 1.0])
 
-    edge_indices = np.array([False, True, False, False])
-    actual = grid.topology_subset(edge_indices)
+    face_index = np.array([False, True, False, False])
+    actual = grid.topology_subset(face_index)
     assert np.array_equal(actual.face_node_connectivity, [[0, 1, 3, 2]])
     assert np.array_equal(actual.node_x, [1.0, 2.0, 1.0, 2.0])
     assert np.array_equal(actual.node_y, [0.0, 0.0, 1.0, 1.0])
 
     # Entire mesh
-    edge_indices = np.array([0, 1, 2, 3])
-    actual = grid.topology_subset(edge_indices)
+    face_index = np.array([0, 1, 2, 3])
+    actual = grid.topology_subset(face_index)
     assert actual is grid
 
     # Check that alternative attrs are preserved.
     grid = grid2d(attrs={"node_dimension": "nNetNode"})
-    edge_indices = np.array([1])
-    actual = grid.topology_subset(edge_indices)
+    face_index = np.array([1])
+    actual = grid.topology_subset(face_index)
     assert actual.node_dimension == "nNetNode"
 
 
@@ -746,10 +768,10 @@ def test_to_pygeos():
 
 def test_grid_from_geodataframe():
     with pytest.raises(TypeError, match="Cannot convert a list"):
-        xugrid.ugrid.grid_from_geodataframe([])
+        xugrid.conversion.grid_from_geodataframe([])
 
     with pytest.raises(ValueError, match="geodataframe contains no geometry"):
-        xugrid.ugrid.grid_from_geodataframe(gpd.GeoDataFrame(geometry=[]))
+        xugrid.conversion.grid_from_geodataframe(gpd.GeoDataFrame(geometry=[]))
 
     x = np.array([0.0, 1.0, 2.0])
     y = np.array([0.0, 0.0, 0.0])
@@ -766,13 +788,17 @@ def test_grid_from_geodataframe():
     points = pygeos.creation.points(x, y)
 
     with pytest.raises(ValueError, match="Multiple geometry types detected"):
-        xugrid.ugrid.grid_from_geodataframe(gpd.GeoDataFrame(geometry=[line, polygon]))
+        xugrid.conversion.grid_from_geodataframe(
+            gpd.GeoDataFrame(geometry=[line, polygon])
+        )
     with pytest.raises(ValueError, match="Invalid geometry type"):
-        xugrid.ugrid.grid_from_geodataframe(gpd.GeoDataFrame(geometry=points))
+        xugrid.conversion.grid_from_geodataframe(gpd.GeoDataFrame(geometry=points))
 
-    grid = xugrid.ugrid.grid_from_geodataframe(gpd.GeoDataFrame(geometry=[line]))
+    grid = xugrid.conversion.grid_from_geodataframe(gpd.GeoDataFrame(geometry=[line]))
     assert isinstance(grid, xugrid.Ugrid1d)
-    grid = xugrid.ugrid.grid_from_geodataframe(gpd.GeoDataFrame(geometry=[polygon]))
+    grid = xugrid.conversion.grid_from_geodataframe(
+        gpd.GeoDataFrame(geometry=[polygon])
+    )
     assert isinstance(grid, xugrid.Ugrid2d)
 
 
diff --git a/tests/test_ugrid_dataset.py b/tests/test_ugrid_dataset.py
index 056823494..aeeabfed4 100644
--- a/tests/test_ugrid_dataset.py
+++ b/tests/test_ugrid_dataset.py
@@ -92,6 +92,7 @@ def setup(self):
     def test_init(self):
         assert isinstance(self.uda.ugrid.obj, xr.DataArray)
         assert isinstance(self.uda.ugrid.grid, xugrid.Ugrid2d)
+        assert self.uda.grid.face_dimension in self.uda.coords
 
     def test_reinit_error(self):
         # Should not be able to initialize using a UgridDataArray.
@@ -114,12 +115,12 @@ def test_getattr(self):
         assert self.uda.dims == self.uda.ugrid.obj.dims
         assert isinstance(self.uda.data, np.ndarray)
         # So are functions
-        assert isinstance(self.uda.mean(), xugrid.UgridDataArray)
+        assert isinstance(self.uda.isnull(), xugrid.UgridDataArray)
         # obj should be accessible
         assert isinstance(self.uda.obj, xr.DataArray)
 
     def test_ugrid_accessor(self):
-        assert isinstance(self.uda.ugrid, xugrid.ugrid_dataset.UgridDataArrayAccessor)
+        assert isinstance(self.uda.ugrid, xugrid.UgridDataArrayAccessor)
 
     def test_from_structured(self):
         da = xr.DataArray([0.0, 1.0, 2.0], {"x": [5.0, 10.0, 15.0]}, ["x"])
@@ -137,7 +138,6 @@ def test_from_structured(self):
         assert uda.name == "grid"
         assert uda.dims == ("layer", "mesh2d_nFaces")
         assert uda.shape == (2, 12)
-        assert uda.ugrid.grid.face_dimension not in uda.coords
         assert np.allclose(uda.ugrid.sel(x=2.0, y=5.0), [[0], [12]])
         # Check whether flipping the y-axis doesn't cause any problems
         flipped = da.isel(y=slice(None, None, -1))
@@ -174,7 +174,7 @@ def test_np_ops(self):
 
     # Accessor tests
     def test_isel(self):
-        actual = self.uda.ugrid.isel({GRID().face_dimension: [0, 1]})
+        actual = self.uda.isel({GRID().face_dimension: [0, 1]})
         assert isinstance(actual, xugrid.UgridDataArray)
         assert actual.shape == (2,)
         assert actual.ugrid.grid.n_face == 2
@@ -248,7 +248,7 @@ def test_to_geodataframe(self):
             self.uda.ugrid.to_geodataframe()
         uda2 = self.uda.copy()
         uda2.ugrid.obj.name = "test"
-        gdf = uda2.ugrid.to_geodataframe("mesh2d_nFaces")
+        gdf = uda2.ugrid.to_geodataframe("mesh2d")
         assert isinstance(gdf, gpd.GeoDataFrame)
         assert (gdf.geometry.geom_type == "Polygon").all()
 
@@ -373,7 +373,7 @@ def test_setitem(self):
     def test_getattr(self):
         assert tuple(self.uds.dims) == ("mesh2d_nFaces",)
         assert isinstance(self.uds.a, xugrid.UgridDataArray)
-        assert isinstance(self.uds.mean(), xugrid.UgridDataset)
+        assert isinstance(self.uds.notnull(), xugrid.UgridDataset)
         # obj should be accessible
         assert isinstance(self.uds.obj, xr.Dataset)
 
@@ -393,7 +393,7 @@ def test_math(self):
         assert isinstance(actual, xugrid.UgridDataset)
 
     def test_ugrid_accessor(self):
-        assert isinstance(self.uds.ugrid, xugrid.ugrid_dataset.UgridDatasetAccessor)
+        assert isinstance(self.uds.ugrid, xugrid.UgridDatasetAccessor)
 
     def test_from_geodataframe(self):
         xy = np.array(
@@ -414,7 +414,7 @@ def test_from_geodataframe(self):
 
     # Accessor tests
     def test_isel(self):
-        actual = self.uds.ugrid.isel({GRID().face_dimension: [0, 1]})
+        actual = self.uds.isel({GRID().face_dimension: [0, 1]})
         assert isinstance(actual, xugrid.UgridDataset)
         assert actual.ugrid.grids[0].n_face == 2
         assert actual["a"].shape == (2,)
diff --git a/tests/test_voronoi.py b/tests/test_voronoi.py
index 481ed0681..6924ce359 100644
--- a/tests/test_voronoi.py
+++ b/tests/test_voronoi.py
@@ -1,7 +1,7 @@
 import numpy as np
 import pytest
 
-from xugrid import connectivity, voronoi
+from xugrid.ugrid import connectivity, voronoi
 
 
 def rowsort(a):