Repository: Deltares/xugrid (language: python)

Issue:
pcolormesh and line ignore levels
In plotting, levels are ignored.

The cause is to be found in these lines: 

```python
    norm = kwargs.pop("norm", None)
    vmin = kwargs.pop("vmin", None)
    vmax = kwargs.pop("vmax", None)

    collection = PolyCollection(vertices, **kwargs)
    collection.set_array(da.values.ravel())
    collection._scale_norm(norm, vmin, vmax)
```

As norm is popped from the kwargs, it's not passed during initialization. I have no idea what _scale_norm does with it, but it's not enough. Replacing pop by get does the trick.

See also: https://github.com/matplotlib/matplotlib/blob/29a86636a9c45ab5ac4d80ac76eaee497f460dce/lib/matplotlib/axes/_axes.py#L6094

(The matplotlib pcolormesh only plots quads, hence the reimplementation.)


Discussion / hints:

Fix patch:
diff --git a/.github/workflows/environment.yml b/.github/workflows/environment.yml
index af92f5e64..a45c998a9 100644
--- a/.github/workflows/environment.yml
+++ b/.github/workflows/environment.yml
@@ -7,6 +7,7 @@ channels:
 dependencies:
   - black
   - dask
+  - geopandas
   - flake8
   - isort
   - meshzoo
@@ -16,6 +17,8 @@ dependencies:
   - pip:
     - meshkernel
   - pydata-sphinx-theme
+  - pygeos
+  - pyproj
   - pytest
   - pytest-cov
   - python=3.9
diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index ddcbfd7db..907fe67c7 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -1,16 +1,16 @@
 repos:
   - repo: https://github.com/PyCQA/isort
-    rev: 5.9.1
+    rev: 5.10.1
     hooks:
       - id: isort
 
   - repo: https://github.com/python/black
-    rev: 21.6b0
+    rev: 22.6.0
     hooks:
       - id: black
         language_version: python3
 
   - repo: https://gitlab.com/pycqa/flake8
-    rev: 3.9.2
+    rev: 5.0.4
     hooks:
       - id: flake8
diff --git a/docs/_static/theme-deltares.css b/docs/_static/theme-deltares.css
index 7d102df78..a7b021dc5 100644
--- a/docs/_static/theme-deltares.css
+++ b/docs/_static/theme-deltares.css
@@ -1,34 +1,29 @@
 /* Override the default color set in the original theme */
 
-.navbar-light {
-    background: #080c80 !important;
-    /* color: #fff !important; */
+html[data-theme="light"] {
+    --pst-color-background: #fff !important;
+    --pst-color-on-background: #080c80 !important;
+    --pst-color-primary: #080c80 !important;
 }
-
-.bg-light {
-    background-color: #080c80 !important;
+html[data-theme="dark"] {
+    --pst-color-border: #080c80 !important;
+    --pst-color-background: #080c80 !important;
+    --pst-color-on-background: #080c80 !important;
 }
 
-.navbar-brand p {
-    color: #fff !important;
+div.admonition {
+    background-color:  #fff;
 }
 
-.navbar-nav>.active>.nav-link {
-    color: #fff !important;
+.nav-link {
+    color: #8486c0 !important;
 }
 
 .navbar-nav>:hover>.nav-link {
     color: #fff !important;
 }
 
-.nav-link {
-    color: #8486c0 !important;
-}
-
 img.icon-link-image {
     height: 2.5em !important;
 }
 
-.navbar-toggler-icon {
-    filter: invert(100%) ;
-}
\ No newline at end of file
diff --git a/docs/api.rst b/docs/api.rst
index 33fa3c2fc..5327f5d00 100644
--- a/docs/api.rst
+++ b/docs/api.rst
@@ -11,8 +11,8 @@ Top-level functions
 -------------------
 
 .. autosummary::
-    :toctree: api/
-    
+   :toctree: api/
+
     open_dataarray
     open_dataset
     open_mfdataset
@@ -20,35 +20,14 @@ Top-level functions
     full_like
     ones_like
     zeros_like
-
-UGRID Accessor
---------------
-
-.. autosummary::
-    :toctree: api/
-    
-    UgridAccessor
-    UgridAccessor.crs
-    UgridAccessor.isel
-    UgridAccessor.sel
-    UgridAccessor.sel_points
-    UgridAccessor.rasterize
-    UgridAccessor.rasterize_like
-    UgridAccessor.to_geodataframe
-    UgridAccessor.binary_dilation
-    UgridAccessor.binary_erosion
-    UgridAccessor.connected_components
-    UgridAccessor.reverse_cuthill_mckee
-    UgridAccessor.laplace_interpolate
-    UgridAccessor.to_dataset
-    UgridAccessor.to_netcdf
-    UgridAccessor.to_zarr
+    concat
+    merge
 
 UgridDataArray
 --------------
 
 .. autosummary::
-    :toctree: api/
+   :toctree: api/
 
     UgridDataArray
     UgridDataArray.ugrid
@@ -56,27 +35,88 @@ UgridDataArray
 
 UgridDataset
 ------------
-    
+
 .. autosummary::
-    :toctree: api/
-    
+   :toctree: api/
+
     UgridDataset
-    UgridDataArray.ugrid
+    UgridDataset.ugrid
     UgridDataset.from_geodataframe
 
+UGRID Accessor
+--------------
+
+These methods and attributes are available under the ``.ugrid`` attribute of a
+UgridDataArray or UgridDataset.
+
+.. autosummary::
+   :toctree: api/
+
+    UgridDatasetAccessor
+    UgridDatasetAccessor.assign_node_coords
+    UgridDatasetAccessor.set_node_coords
+    UgridDatasetAccessor.crs
+    UgridDatasetAccessor.isel
+    UgridDatasetAccessor.sel
+    UgridDatasetAccessor.sel_points
+    UgridDatasetAccessor.to_geodataframe
+    UgridDatasetAccessor.to_dataset
+    UgridDatasetAccessor.to_netcdf
+    UgridDatasetAccessor.to_zarr
+
+    UgridDataArrayAccessor
+    UgridDataArrayAccessor.assign_node_coords
+    UgridDataArrayAccessor.set_node_coords
+    UgridDataArrayAccessor.crs
+    UgridDataArrayAccessor.isel
+    UgridDataArrayAccessor.sel
+    UgridDataArrayAccessor.sel_points
+    UgridDataArrayAccessor.rasterize
+    UgridDataArrayAccessor.rasterize_like
+    UgridDataArrayAccessor.to_geodataframe
+    UgridDataArrayAccessor.binary_dilation
+    UgridDataArrayAccessor.binary_erosion
+    UgridDataArrayAccessor.connected_components
+    UgridDataArrayAccessor.reverse_cuthill_mckee
+    UgridDataArrayAccessor.laplace_interpolate
+    UgridDataArrayAccessor.to_dataset
+    UgridDataArrayAccessor.to_netcdf
+    UgridDataArrayAccessor.to_zarr
+
+Plotting
+--------
+
+These methods are also available under the ``.ugrid.plot`` attribute of a
+UgridDataArray.
+
+.. autosummary::
+   :toctree: api/
+
+   plot.contour
+   plot.contourf
+   plot.imshow
+   plot.line
+   plot.pcolormesh
+   plot.scatter
+   plot.surface
+   plot.tripcolor
+
 UGRID1D Topology
 ----------------
 
 .. autosummary::
-    :toctree: api/
-    
+   :toctree: api/
+
     Ugrid1d
 
     Ugrid1d.topology_dimension
+    Ugrid1d.dimensions
 
     Ugrid1d.n_node
     Ugrid1d.node_dimension
     Ugrid1d.node_coordinates
+    Ugrid1d.set_node_coords
+    Ugrid2d.assign_node_coords
 
     Ugrid1d.n_edge
     Ugrid1d.edge_dimension
@@ -85,25 +125,23 @@ UGRID1D Topology
     Ugrid1d.edge_y
 
     Ugrid1d.bounds
-    
+
     Ugrid1d.node_edge_connectivity
 
     Ugrid1d.copy
 
-    Ugrid1d.remove_topology
-    Ugrid1d.topology_coords
-    Ugrid1d.topology_dataset
-
+    Ugrid1d.isel
     Ugrid1d.sel
     Ugrid1d.topology_subset
 
     Ugrid1d.mesh
     Ugrid1d.meshkernel
-    
+
     Ugrid1d.set_crs
     Ugrid1d.to_crs
 
     Ugrid1d.from_dataset
+    Ugrid1d.to_dataset
     Ugrid1d.from_geodataframe
     Ugrid1d.to_pygeos
 
@@ -111,15 +149,18 @@ UGRID2D Topology
 ----------------
 
 .. autosummary::
-    :toctree: api/
+   :toctree: api/
 
     Ugrid2d
 
     Ugrid2d.topology_dimension
+    Ugrid2d.dimensions
 
     Ugrid2d.n_node
     Ugrid2d.node_dimension
     Ugrid2d.node_coordinates
+    Ugrid2d.set_node_coords
+    Ugrid2d.assign_node_coords
 
     Ugrid2d.n_edge
     Ugrid2d.edge_dimension
@@ -136,18 +177,18 @@ UGRID2D Topology
 
     Ugrid2d.bounds
     Ugrid2d.centroids
-    
+
     Ugrid2d.node_edge_connectivity
     Ugrid2d.node_face_connectivity
     Ugrid2d.edge_node_connectivity
     Ugrid2d.face_edge_connectivity
     Ugrid2d.face_face_connectivity
-    
+
     Ugrid2d.exterior_edges
     Ugrid2d.exterior_faces
-    
+
     Ugrid2d.copy
-    
+
     Ugrid2d.triangulate
     Ugrid2d.triangulation
     Ugrid2d.voronoi_topology
@@ -155,10 +196,7 @@ UGRID2D Topology
     Ugrid2d.tesselate_centroidal_voronoi
     Ugrid2d.reverse_cuthill_mckee
 
-    Ugrid2d.remove_topology
-    Ugrid2d.topology_coords
-    Ugrid2d.topology_dataset
-
+    Ugrid2d.isel
     Ugrid2d.sel
     Ugrid2d.sel_points
     Ugrid2d.celltree
@@ -171,11 +209,24 @@ UGRID2D Topology
 
     Ugrid2d.mesh
     Ugrid2d.meshkernel
-    
+
     Ugrid2d.set_crs
     Ugrid2d.to_crs
 
     Ugrid2d.from_dataset
+    Ugrid2d.to_dataset
     Ugrid2d.from_geodataframe
     Ugrid2d.from_structured
     Ugrid2d.to_pygeos
+
+UGRID Roles Accessor
+--------------------
+
+.. autosummary::
+   :toctree: api/
+
+    UgridRolesAccessor
+    UgridRolesAccessor.topology
+    UgridRolesAccessor.connectivity
+    UgridRolesAccessor.coordinates
+    UgridRolesAccessor.dimensions
diff --git a/docs/terminology.rst b/docs/terminology.rst
index 8cda19778..812fc99e2 100644
--- a/docs/terminology.rst
+++ b/docs/terminology.rst
@@ -53,7 +53,7 @@ recommended read.
         topology. "Vertex" is used interchangeably.
 
     Edge
-        A line or cuve bounded by two nodes.
+        A line or curve bounded by two nodes.
 
     Face
         A plane or surface enclosed by a set of edges. "Cell" is used somewhat
diff --git a/examples/connectivity.py b/examples/connectivity.py
index 7b7af756f..78d92f4eb 100644
--- a/examples/connectivity.py
+++ b/examples/connectivity.py
@@ -70,8 +70,8 @@
 
 ds = xugrid.data.disk()
 uda = xugrid.UgridDataArray(
-    xr.full_like(ds["face_z"].ugrid.obj, True, dtype=bool),
-    ds.grid,
+    xr.full_like(ds.obj["face_z"], True, dtype=bool),
+    ds.grids[0],
 )
 iter2 = uda.ugrid.binary_erosion(iterations=2)
 iter5 = uda.ugrid.binary_erosion(iterations=5)
@@ -87,7 +87,7 @@
 
 uda = xugrid.UgridDataArray(
     xr.full_like(ds["face_z"].ugrid.obj, False, dtype=bool),
-    ds.grid,
+    ds.grids[0],
 )
 uda.values[0] = True
 uda.ugrid.plot()
diff --git a/examples/plotting.py b/examples/plotting.py
index d19232ea0..57ade20d7 100644
--- a/examples/plotting.py
+++ b/examples/plotting.py
@@ -171,8 +171,8 @@
 # The plotting methods can also be called as a function, in which case they
 # take an xarray DataArray and a xugrid grid as arguments.
 
-grid = ds.ugrid.grid
-da = ds["face_z"].ugrid.obj
+grid = ds.ugrid.grids[0]
+da = ds.obj["face_z"]
 
 xugrid.plot.pcolormesh(grid, da)
 
diff --git a/examples/quick_overview.py b/examples/quick_overview.py
index f3e0d3edf..95750e3d8 100644
--- a/examples/quick_overview.py
+++ b/examples/quick_overview.py
@@ -121,7 +121,7 @@
 #
 #   ``ugrid.sel()`` currently only supports data on the faces for 2D
 #   topologies, and data on edges for 1D topologies. More flexibility
-#   will be added soon.
+#   may be added.
 #
 # Computation
 # -----------
@@ -135,9 +135,8 @@
 # ---------
 #
 # Xugrid objects provide a number of conversion functions from and to geopandas
-# GeoDataFrames using :py:meth:`xugrid.UgridDataArray.to_geodataframe`, and
-# :py:meth:`xugrid.UgridDataset.from_geodataframe`. Note that storing large
-# grids as GeoDataFrames can be very inefficient.
+# GeoDataFrames using :py:meth:`xugrid.UgridDataset.from_geodataframe`. Note
+# that storing large grids as GeoDataFrames can be very inefficient.
 
 gdf = uda.ugrid.to_geodataframe(name="test")
 gdf
@@ -162,7 +161,7 @@
 # A UgridDataset may be initialized without data variables, but this requires
 # a grid object:
 
-new_uds = xu.UgridDataset(grid=uds.ugrid.grid)
+new_uds = xu.UgridDataset(grids=uds.ugrid.grids)
 new_uds
 
 # %%
diff --git a/setup.cfg b/setup.cfg
index 2ad005752..539a4e6d9 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -1,6 +1,6 @@
 [metadata]
 name = xugrid
-version = 0.1.3
+version = 0.1.4
 author = Huite Bootsma
 author_email = huite.bootsma@deltares.nl
 description = Xarray extension for unstructured grids
diff --git a/xugrid/__init__.py b/xugrid/__init__.py
index 47d4a50e5..7e333456d 100644
--- a/xugrid/__init__.py
+++ b/xugrid/__init__.py
@@ -2,12 +2,15 @@
 
 from . import data
 from .plot import plot
-from .ugrid import Ugrid1d, Ugrid2d
+from .ugrid import Ugrid1d, Ugrid2d, UgridRolesAccessor
 from .ugrid_dataset import (
-    UgridAccessor,
     UgridDataArray,
+    UgridDataArrayAccessor,
     UgridDataset,
+    UgridDatasetAccessor,
+    concat,
     full_like,
+    merge,
     ones_like,
     open_dataarray,
     open_dataset,
diff --git a/xugrid/data/synthetic.py b/xugrid/data/synthetic.py
index 88c3e3d16..8c58b15bd 100644
--- a/xugrid/data/synthetic.py
+++ b/xugrid/data/synthetic.py
@@ -36,7 +36,7 @@ def function_z(x, y):
         z = -(
             2 * (np.exp((r1 / 10) ** 2) - 1) * 30.0 * np.cos(7.0 * theta1)
             + (np.exp((r2 / 10) ** 2) - 1) * 30.0 * np.cos(11.0 * theta2)
-            + 0.7 * (x ** 2 + y ** 2)
+            + 0.7 * (x**2 + y**2)
         )
         zmin = z.min()
         zmax = z.max()
diff --git a/xugrid/interpolate.py b/xugrid/interpolate.py
index bd714bd9a..498e4ce6a 100644
--- a/xugrid/interpolate.py
+++ b/xugrid/interpolate.py
@@ -125,7 +125,9 @@ def laplace_interpolate(
         )
         M = scipy.sparse.linalg.LinearOperator((n, n), ilu.solve)
         # Call conjugate gradient solver
-        x, info = scipy.sparse.linalg.cg(A, rhs, tol=tol, maxiter=maxiter, M=M)
+        x, info = scipy.sparse.linalg.cg(
+            A, rhs, tol=tol, maxiter=maxiter, M=M, atol="legacy"
+        )
         if info < 0:
             raise ValueError("scipy.sparse.linalg.cg: illegal input or breakdown")
         elif info > 0:
diff --git a/xugrid/plot/__init__.py b/xugrid/plot/__init__.py
index e69de29bb..3a7089e49 100644
--- a/xugrid/plot/__init__.py
+++ b/xugrid/plot/__init__.py
@@ -0,0 +1,10 @@
+from .plot import (
+    contour,
+    contourf,
+    imshow,
+    line,
+    pcolormesh,
+    scatter,
+    surface,
+    tripcolor,
+)
diff --git a/xugrid/plot/plot.py b/xugrid/plot/plot.py
index 54025af82..fcad4a634 100644
--- a/xugrid/plot/plot.py
+++ b/xugrid/plot/plot.py
@@ -4,7 +4,6 @@
 import functools
 
 import numpy as np
-import xarray as xr
 from xarray.plot.facetgrid import _easy_facetgrid
 from xarray.plot.utils import (
     _add_colorbar,
@@ -402,7 +401,8 @@ def line(grid, da, ax, **kwargs):
     edge_coords[:, 1, 0] = grid.node_x[node_1]
     edge_coords[:, 1, 1] = grid.node_y[node_1]
 
-    norm = kwargs.pop("norm", None)
+    # PolyCollection takes a norm, but not vmin, vmax.
+    norm = kwargs.get("norm", None)
     vmin = kwargs.pop("vmin", None)
     vmax = kwargs.pop("vmax", None)
 
@@ -503,7 +503,8 @@ def contourf(grid, da, ax, **kwargs):
 def pcolormesh(grid, da, ax, **kwargs):
     """
     Pseudocolor plot of 2D UgridDataArray.
-    Wraps :py:func:`matplotlib:matplotlib.pyplot.
+
+    Wraps matplotlib PolyCollection.
     """
     from matplotlib.collections import PolyCollection
 
@@ -517,7 +518,8 @@ def pcolormesh(grid, da, ax, **kwargs):
     # Replace fill value; PolyCollection ignores NaN.
     vertices[faces == -1] = np.nan
 
-    norm = kwargs.pop("norm", None)
+    # PolyCollection takes a norm, but not vmin, vmax.
+    norm = kwargs.get("norm", None)
     vmin = kwargs.pop("vmin", None)
     vmax = kwargs.pop("vmax", None)
 
@@ -584,14 +586,26 @@ def plot(
     """
     dim = darray.dims[0]
     kwargs["ax"] = ax
-    if dim == grid.face_dimension:
-        return pcolormesh(grid, darray, **kwargs)
-    elif dim == grid.node_dimension:
-        return tripcolor(grid, darray, **kwargs)
-    elif dim == grid.edge_dimension:
-        return line(grid, darray, **kwargs)
+    if grid.topology_dimension == 1:
+        if dim == grid.edge_dimension:
+            return line(grid, darray, **kwargs)
+        elif dim == grid.node_dimension:
+            return scatter(grid, darray, **kwargs)
+        else:
+            raise ValueError("Data dimensions is not one of node or edge dimension.")
+    elif grid.topology_dimension == 2:
+        if dim == grid.face_dimension:
+            return pcolormesh(grid, darray, **kwargs)
+        elif dim == grid.node_dimension:
+            return tripcolor(grid, darray, **kwargs)
+        elif dim == grid.edge_dimension:
+            return line(grid, darray, **kwargs)
+        else:
+            raise ValueError(
+                "Data dimensions is not one of face, node, or edge dimension."
+            )
     else:
-        raise ValueError("Data dimensions is not one of face, node, or edge dimension.")
+        raise ValueError("Topology dimension is not 1 or 2")
 
 
 class _PlotMethods:
@@ -606,14 +620,13 @@ def __init__(self, obj):
         darray = obj.obj
         grid = obj.grid
 
-        if isinstance(darray, xr.Dataset):
-            raise NotImplementedError("Cannot plot Datasets into a facetgrid (yet)")
-        if len(darray.dims) > 1:
-            msg = (
-                "Data contains more dimensions than just a topology dimension "
-                f"(face, node, edge): {', '.join(darray.dims)}"
+        invalid = set(darray.dims) - set(grid.dimensions)
+        if invalid:
+            raise ValueError(
+                f"UgridDataArray contains non-topology dimensions: {invalid}.\n"
+                f"Expected one of {grid.dimensions}"
             )
-            raise ValueError(msg)
+
         self.grid = grid
         self.darray = darray
 
diff --git a/xugrid/ugrid/__init__.py b/xugrid/ugrid/__init__.py
index 1af580436..f1d60e4fa 100644
--- a/xugrid/ugrid/__init__.py
+++ b/xugrid/ugrid/__init__.py
@@ -1,5 +1,9 @@
+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
@@ -28,3 +32,15 @@ def grid_from_geodataframe(geodataframe: "geopandas.GeoDataFrame"):  # type: ign
             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/ugrid/conventions.py b/xugrid/ugrid/conventions.py
new file mode 100644
index 000000000..0629a1ad0
--- /dev/null
+++ b/xugrid/ugrid/conventions.py
@@ -0,0 +1,449 @@
+"""
+This module deals with extracting the relevant data from the UGRID attributes.
+
+It takes some inspiration from: https://github.com/xarray-contrib/cf-xarray
+"""
+import warnings
+from collections import ChainMap
+from itertools import chain
+from typing import Tuple
+
+import xarray as xr
+
+
+class UgridDimensionError(Exception):
+    pass
+
+
+class UgridCoordinateError(Exception):
+    pass
+
+
+_DIM_NAMES = {
+    1: ("node_dimension", "edge_dimension"),
+    2: ("node_dimension", "face_dimension", "edge_dimension"),
+}
+
+_COORD_NAMES = {
+    1: ("node_coordinates", "edge_coordinates"),
+    2: ("node_coordinates", "face_coordinates", "edge_coordinates"),
+}
+_COORD_DIMS = {
+    "node_coordinates": "node_dimension",
+    "edge_coordinates": "node_dimension",
+    "face_coordinates": "face_dimension",
+}
+
+_CONNECTIVITY_NAMES = {
+    1: ("edge_node_connectivity",),
+    2: (
+        "face_node_connectivity",
+        "edge_node_connectivity",
+        "face_edge_connectivity",
+        # "face_face_connectivity",
+        "edge_face_connectivity",
+        #        "boundary_node_connectivity",
+    ),
+}
+
+_CONNECTIVITY_DIMS = {
+    "face_node_connectivity": ("face_dimension", None),
+    "edge_node_connectivity": ("edge_dimension", 2),
+    "face_edge_connectivity": ("face_dimension", None),
+    # "face_face_connectivity": ("face_dimension", None),
+    "edge_face_connectivity": ("edge_dimension", 2),
+    # "boundary_node_connectivity": ("???","???"),
+}
+
+X_STANDARD_NAMES = ("projection_x_coordinate", "longitude")
+Y_STANDARD_NAMES = ("projection_y_coordinate", "latitude")
+
+PROJECTED = True
+LATLON = False
+DEFAULT_ATTRS = {
+    "node_x": {
+        PROJECTED: {
+            "standard_name": "projection_x_coordinate",
+        },
+        LATLON: {
+            "standard_name": "longitude",
+        },
+    },
+    "node_y": {
+        PROJECTED: {
+            "standard_name": "projection_y_coordinate",
+        },
+        LATLON: {
+            "standard_name": "latitude",
+        },
+    },
+    "edge_x": {
+        PROJECTED: {
+            "standard_name": "projection_x_coordinate",
+        },
+        LATLON: {
+            "standard_name": "longitude",
+        },
+    },
+    "edge_y": {
+        PROJECTED: {
+            "standard_name": "projection_y_coordinate",
+        },
+        LATLON: {
+            "standard_name": "latitude",
+        },
+    },
+    "face_x": {
+        PROJECTED: {
+            "standard_name": "projection_x_coordinate",
+        },
+        LATLON: {
+            "standard_name": "longitude",
+        },
+    },
+    "face_y": {
+        PROJECTED: {
+            "standard_name": "projection_y_coordinate",
+        },
+        LATLON: {
+            "standard_name": "latitude",
+        },
+    },
+    "face_node_connectivity": {
+        "cf_role": "face_node_connectivity",
+        "start_index": 0,
+        "_FillValue": -1,
+    },
+    "edge_node_connectivity": {
+        "cf_role": "edge_node_connectivity",
+        "start_index": 0,
+        "_FillValue": -1,
+    },
+    "face_edge_connectivity": {
+        "cf_role": "face_edge_connectivity",
+        "start_index": 0,
+        "_FillValue": -1,
+    },
+    # "face_face_connectivity": {
+    #    "cf_role": "face_face_connectivity",
+    #    "start_index": 0,
+    #    "_FillValue": -1,
+    # },
+    "edge_face_connectivity": {
+        "cf_role": "edge_face_connectivity",
+        "start_index": 0,
+        "_FillValue": -1,
+    },
+}
+
+
+def default_topology_attrs(name: str, topology_dimension: int):
+    if topology_dimension == 1:
+        return {
+            "cf_role": "mesh_topology",
+            "long_name": "Topology data of 1D network",
+            "topology_dimension": 1,
+            "node_dimension": f"{name}_nNodes",
+            "edge_dimension": f"{name}_nEdges",
+            "edge_node_connectivity": f"{name}_edge_nodes",
+            "node_coordinates": f"{name}_node_x {name}_node_y",
+            "edge_coordinates": f"{name}_edge_x {name}_edge_y",
+        }
+    elif topology_dimension == 2:
+        return {
+            "cf_role": "mesh_topology",
+            "long_name": "Topology data of 2D mesh",
+            "topology_dimension": 2,
+            "node_dimension": f"{name}_nNodes",
+            "edge_dimension": f"{name}_nEdges",
+            "face_dimension": f"{name}_nFaces",
+            "max_face_nodes_dimension": f"{name}_nMax_face_nodes",
+            "edge_node_connectivity": f"{name}_edge_nodes",
+            "face_node_connectivity": f"{name}_face_nodes",
+            "node_coordinates": f"{name}_node_x {name}_node_y",
+            "edge_coordinates": f"{name}_edge_x {name}_edge_y",
+            "face_coordinates": f"{name}_face_x {name}_face_y",
+        }
+    else:
+        raise ValueError(
+            f"topology_dimensions should be 1 or 2, received {topology_dimension}"
+        )
+
+
+def _get_topology(ds: xr.Dataset) -> list[str]:
+    return [k for k in ds.data_vars if ds[k].attrs.get("cf_role") == "mesh_topology"]
+
+
+def _infer_xy_coords(
+    ds: xr.Dataset, candidates: list[str]
+) -> Tuple[list[str], list[str]]:
+    # TODO: add argument for latitude / longitude?
+    x = []
+    y = []
+    for candidate in candidates:
+        stdname = ds[candidate].attrs.get("standard_name")
+        if stdname in X_STANDARD_NAMES:
+            x.append(candidate)
+        elif stdname in Y_STANDARD_NAMES:
+            y.append(candidate)
+
+    if not x and not y:
+        first = candidates[0]
+        second = candidates[1]
+        warnings.warn(
+            f"No standard_name of {X_STANDARD_NAMES + Y_STANDARD_NAMES} in {candidates}.\n"
+            f"Using {first} and {second} as projected x and y coordinates."
+        )
+        x.append(first)
+        y.append(second)
+    elif not x:
+        raise UgridCoordinateError(
+            "No standard_name of {X_STANDARD_NAMES} in {candidates}"
+        )
+    elif not y:
+        raise UgridCoordinateError(
+            "No standard_name of {Y_STANDARD_NAMES} in {candidates}"
+        )
+
+    return x, y
+
+
+def _get_coordinates(
+    ds: xr.Dataset, topologies: list[str]
+) -> dict[str, dict[str, Tuple[list[str], list[str]]]]:
+    topology_dict = {}
+    for topology in topologies:
+        attrs = ds[topology].attrs
+        topodim = attrs["topology_dimension"]
+        vardict = {}
+        for name in _COORD_NAMES[topodim]:
+            if name in attrs:
+                candidates = [c for c in attrs[name].split(" ") if c in ds]
+                if len(candidates) < 2:
+                    raise UgridCoordinateError(
+                        f"{topology}: at least two values required for {name}"
+                    )
+                vardict[name] = _infer_xy_coords(ds, candidates)
+
+        topology_dict[topology] = vardict
+
+    return topology_dict
+
+
+def _infer_dims(
+    ds: xr.Dataset,
+    connectivities: dict[str, str],
+    coordinates: dict[str, dict[str, Tuple[list[str]]]],
+    vardict: dict[str, str],
+) -> dict[str, str]:
+    """
+    Infer dimensions based on connectivity and coordinates.
+    """
+    inferred = {}
+    for role, varname in connectivities.items():
+        expected_dims = _CONNECTIVITY_DIMS[role]
+        var_dims = ds[varname].dims
+        for key, dim in zip(expected_dims, var_dims):
+            if isinstance(key, str):  # skip None or integer
+                prev_dim = vardict.get(key)
+                # Not specified: default order can be used to infer dimensions.
+                if prev_dim is None:
+                    inferred[key] = dim
+                else:
+                    if prev_dim not in var_dims:
+                        raise UgridDimensionError(
+                            f"{key}: {prev_dim} not in {role}: {varname}"
+                            f" with dimensions: {var_dims}"
+                        )
+
+    for role, varnames in coordinates.items():
+        key = _COORD_DIMS[role]
+        for varname in chain.from_iterable(varnames):
+            var_dims = ds[varname].dims
+            if len(var_dims) != 1:
+                continue
+            var_dim = var_dims[0]
+
+            prev_dim = vardict.get(key) or inferred.get(key)
+            if prev_dim is None:
+                inferred[key] = var_dim
+            else:
+                if prev_dim != var_dim:
+                    raise UgridDimensionError(
+                        f"Conflicting names for {key}: {prev_dim} versus {var_dim}"
+                    )
+
+    return inferred
+
+
+def _get_dimensions(
+    ds: xr.Dataset,
+    topologies: list[str],
+    connectivity: dict[str, dict[str, str]],
+    coordinates: dict[str, dict[str, Tuple[list[str]]]],
+) -> dict[str, dict[str, str]]:
+    """
+    Get the dimensions from the topology attributes and infer them from
+    connectivity arrays or coordinates.
+    """
+    topology_dict = {}
+    for topology in topologies:
+        attrs = ds[topology].attrs
+        topodim = attrs["topology_dimension"]
+        # dimensions are optionally required: only if the dimension order is
+        # nonstandard in any of the connectivity variables.
+        vardict = {k: attrs[k] for k in _DIM_NAMES[topodim] if k in attrs}
+        inferred = _infer_dims(
+            ds, connectivity[topology], coordinates[topology], vardict
+        )
+        topology_dict[topology] = {**vardict, **inferred}
+
+    return topology_dict
+
+
+def _get_connectivity(
+    ds: xr.Dataset, topologies: list[str]
+) -> dict[str, dict[str, str]]:
+    topology_dict = {}
+    for topology in topologies:
+        attrs = ds[topology].attrs
+        topodim = attrs["topology_dimension"]
+        topology_dict[topology] = {
+            k: attrs[k]
+            for k in _CONNECTIVITY_NAMES[topodim]
+            if (k in attrs) and (attrs[k] in ds)
+        }
+    return topology_dict
+
+
+@xr.register_dataset_accessor("ugrid_roles")
+class UgridRolesAccessor:
+    """
+    Xarray Dataset "accessor" to retrieve the names of UGRID variables.
+
+    Examples
+    --------
+
+    To get a list of the UGRID dummy variables in the dataset:
+
+    >>> dataset.ugrid_roles.topology
+
+    To get the names of the connectivity variables in the dataset:
+
+    >>> dataset.ugrid_roles.connectivity
+
+    Names can also be accessed directly through the topology:
+
+    >>> dataset.ugrid_roles["mesh2d"]["node_dimension"]
+
+    """
+
+    def __init__(self, ds: xr.Dataset):
+        self._ds = ds
+
+    def __getitem__(self, key: str):
+        if key not in self.topology:
+            raise KeyError(key)
+        return ChainMap(
+            self.dimensions[key], self.coordinates[key], self.connectivity[key]
+        )
+
+    @property
+    def topology(self) -> list[str]:
+        """
+        Get the names of the topology dummy variables, marked by a CF-role of
+        ``mesh_topology``.
+
+        Returns
+        -------
+        topology: list[str]
+        """
+        return _get_topology(self._ds)
+
+    @property
+    def coordinates(self) -> dict[str, dict[str, Tuple[list[str], list[str]]]]:
+        """
+        Get the names of the coordinate variables from the topology attributes.
+
+        Returns a dictionary with the coordinates for the UGRID coordinates:
+
+            * node coordinates
+            * edge coordinates
+            * face coordinates
+
+        Multiple coordinates may be defined. The coordinates are grouped by
+        their role (x or y).
+
+        Returns
+        -------
+        coordinates: dict[str, dict[str, Tuple[list[str]]]]
+        """
+        return _get_coordinates(self._ds, self.topology)
+
+    @property
+    def dimensions(self) -> dict[str, dict[str, str]]:
+        """
+        Get the dimension names from the topology attributes and infer them
+        from connectivity arrays or coordinates.
+
+        Returns a dictionary with the UGRID dimensions per topology:
+
+            * node dimension
+            * edge dimension
+            * face dimension
+
+        Returns
+        -------
+        dimensions: dict[str, dict[str, str]]
+        """
+        return _get_dimensions(
+            self._ds, self.topology, self.connectivity, self.coordinates
+        )
+
+    @property
+    def connectivity(self) -> dict[str, dict[str, str]]:
+        """
+        Get the names of the variables containing the UGRID connectivity data.
+
+            * face_node_connectivity
+            * edge_node_connectivity
+            * face_edge_connectivity
+            * edge_face_connectivity
+
+        Returns
+        -------
+        connectivity: dict[str, dict[str, str]]
+        """
+        return _get_connectivity(self._ds, self.topology)
+
+    def __repr__(self):
+        dimensions = self.dimensions
+        coordinates = self.coordinates
+        connectivity = self.connectivity
+
+        def make_text_section(subtitle, entries, vardict):
+            tab = "    "
+            rows = [f"{tab}{subtitle}"]
+            for role in entries:
+                if role in vardict:
+                    rows += [f"{tab}{tab}{role}: {vardict[role]}"]
+                else:
+                    rows += [f"{tab}{tab}{role}: n/a"]
+            rows.append("")
+            return rows
+
+        rows = []
+        for topology in self.topology:
+            topodim = self._ds[topology].attrs["topology_dimension"]
+            rows += [f"UGRID {topodim}D Topology {topology}:"]
+            rows += make_text_section(
+                "Dimensions:", _DIM_NAMES[topodim], dimensions[topology]
+            )
+            rows += make_text_section(
+                "Connectivity:", _CONNECTIVITY_NAMES[topodim], connectivity[topology]
+            )
+            rows += make_text_section(
+                "Coordinates:", _COORD_NAMES[topodim], coordinates[topology]
+            )
+
+        return "\n".join(rows)
diff --git a/xugrid/ugrid/ugrid1d.py b/xugrid/ugrid/ugrid1d.py
index 89da8e14a..3c5f64262 100644
--- a/xugrid/ugrid/ugrid1d.py
+++ b/xugrid/ugrid/ugrid1d.py
@@ -1,3 +1,4 @@
+from itertools import chain
 from typing import Any, Tuple, Union
 
 import numpy as np
@@ -5,7 +6,7 @@
 
 from .. import conversion
 from ..typing import BoolArray, FloatArray, FloatDType, IntArray, IntDType
-from . import ugrid_io
+from . import conventions
 from .ugridbase import AbstractUgrid
 
 
@@ -35,18 +36,34 @@ def __init__(
         node_y: FloatArray,
         fill_value: int,
         edge_node_connectivity: IntArray = None,
+        name: str = "network1d",
         dataset: xr.Dataset = None,
-        name: str = None,
+        indexes: dict[str, str] = None,
+        projected: bool = True,
         crs: Any = None,
     ):
+        if dataset is not None and indexes is None:
+            raise ValueError("indexes must be provided for a dataset")
+
         self.node_x = np.ascontiguousarray(node_x)
         self.node_y = np.ascontiguousarray(node_y)
         self.fill_value = fill_value
         self.edge_node_connectivity = edge_node_connectivity
-        if name is None:
-            name = ugrid_io.UGRID1D_DEFAULT_NAME
         self.name = name
-        self.topology_attrs = None
+        self.projected = projected
+
+        # Store topology attributes
+        defaults = conventions.default_topology_attrs(name, self.topology_dimension)
+        if dataset is None:
+            self._attrs = defaults
+            x, y = defaults["node_coordinates"].split(" ")
+            self._indexes = {"node_x": x, "node_y": y}
+        else:
+            derived_dims = dataset.ugrid_roles.dimensions[name]
+            self._attrs = {**defaults, **derived_dims, **dataset[name].attrs}
+            self._indexes = indexes
+
+        self._dataset = dataset
 
         # Optional attributes, deferred initialization
         # Meshkernel
@@ -71,14 +88,9 @@ def __init__(
             import pyproj
 
             self.crs = pyproj.CRS.from_user_input(crs)
-        # Store dataset
-        if dataset is None:
-            dataset = self.topology_dataset()
-        self.dataset = dataset
-        self.topology_attrs = dataset[self.name].attrs
 
-    @staticmethod
-    def from_dataset(dataset: xr.Dataset):
+    @classmethod
+    def from_dataset(cls, dataset: xr.Dataset, topology: str = None):
         """
         Extract the 1D UGRID topology information from an xarray Dataset.
 
@@ -89,7 +101,7 @@ def from_dataset(dataset: xr.Dataset):
 
         Returns
         -------
-        grid: Ugrid1d
+        grid: Ugrid1dAdapter
         """
         ds = dataset
         if not isinstance(ds, xr.Dataset):
@@ -97,97 +109,74 @@ def from_dataset(dataset: xr.Dataset):
                 "Ugrid should be initialized with xarray.Dataset. "
                 f"Received instead: {type(ds)}"
             )
+        if topology is None:
+            topology = cls._single_topology(ds)
+
+        indexes = {}
 
         # Collect names
-        mesh_topology = ugrid_io.get_topology_variable(ds)
-        ugrid_roles = ugrid_io.get_ugrid2d_variables(dataset, mesh_topology)
-        # Rename the arrays to their standard UGRID roles
-        topology_ds = ds[set(ugrid_roles.keys())]
-        ugrid_ds = topology_ds.rename(ugrid_roles)
-        # Coerce type and fill value
-        ugrid_ds["node_x"] = ugrid_ds["node_x"].astype(FloatDType)
-        ugrid_ds["node_y"] = ugrid_ds["node_y"].astype(FloatDType)
-        fill_value = -1
-        ugrid_ds["edge_node_connectivity"] = ugrid_ds["edge_node_connectivity"].astype(
-            IntDType
+        connectivity = ds.ugrid_roles.connectivity[topology]
+        coordinates = ds.ugrid_roles.coordinates[topology]
+        ugrid_vars = (
+            [topology]
+            + list(connectivity.values())
+            + list(chain.from_iterable(chain.from_iterable(coordinates.values())))
         )
+        fill_value = -1
 
-        return Ugrid1d(
-            ugrid_ds["node_x"].values,
-            ugrid_ds["node_y"].values,
+        # Take the first coordinates by default.
+        # They can be reset with .set_node_coords()
+        x_index = coordinates["node_coordinates"][0][0]
+        y_index = coordinates["node_coordinates"][1][0]
+        node_x_coordinates = ds[x_index].astype(FloatDType).values
+        node_y_coordinates = ds[y_index].astype(FloatDType).values
+
+        edge_nodes = connectivity["edge_node_connectivity"]
+        edge_node_connectivity = cls._prepare_connectivity(
+            ds[edge_nodes], fill_value, dtype=IntDType
+        ).values
+
+        indexes["node_x"] = x_index
+        indexes["node_y"] = y_index
+        projected = False  # TODO
+
+        return cls(
+            node_x_coordinates,
+            node_y_coordinates,
             fill_value,
-            ugrid_ds["edge_node_connectivity"].values,
-            topology_ds,
-            mesh_topology.name,
+            edge_node_connectivity,
+            name=topology,
+            dataset=dataset[ugrid_vars],
+            indexes=indexes,
+            projected=projected,
+            crs=None,
         )
 
-    def remove_topology(
-        self, obj: Union[xr.DataArray, xr.Dataset]
-    ) -> Union[xr.DataArray, xr.Dataset]:
-        """
-        Remove UGRID specific variables from the object.
-
-        Parameters
-        ----------
-        obj: Union[xr.DataArray, xr.Dataset]
-
-        Returns
-        -------
-        sanitized: Union[xr.DataArray, xr.Dataset]
-        """
-        return self._remove_topology(obj, ugrid_io.UGRID1D_TOPOLOGY_VARIABLES)
-
-    def topology_coords(self, obj: Union[xr.DataArray, xr.Dataset]) -> dict:
-        """
-        Return a dictionary with the coordinates required for the dimension of
-        the object.
-
-        Parameters
-        ----------
-        obj: Union[xr.DataArray, xr.Dataset]
-
-        Returns
-        -------
-        coords: dict
-        """
-        coords = {}
-        dims = obj.dims
-        attrs = self.topology_attrs
-        edgedim = self.edge_dimension
-        nodedim = self.node_dimension
-        if edgedim in dims:
-            # It's not a required attribute.
-            if "edge_coordinates" in attrs:
-                name_x, name_y = attrs["edge_coordinates"].split()
-            else:
-                name_x = f"{self.name}_edge_x"
-                name_y = f"{self.name}_edge_y"
-            coords[name_x] = (edgedim, self.edge_x)
-            coords[name_y] = (edgedim, self.edge_y)
-            attrs["edge_coordinates"] = f"{name_x} {name_y}"
-        if nodedim in dims:
-            name_x, name_y = attrs["node_coordinates"].split()
-            coords[name_x] = (nodedim, self.node_x)
-            coords[name_y] = (nodedim, self.node_y)
-        return coords
-
-    def topology_dataset(self):
-        """
-        Store the UGRID topology information in an xarray Dataset according to
-        the UGRID conventions.
-
-        Returns
-        -------
-        ugrid_topology: xr.Dataset
-        """
-        return ugrid_io.ugrid1d_dataset(
-            self.node_x,
-            self.node_y,
-            self.fill_value,
-            self.edge_node_connectivity,
-            self.name,
-            self.topology_attrs,
-        )
+    def to_dataset(self, other: xr.Dataset = None) -> xr.Dataset:
+        node_x = self._indexes["node_x"]
+        node_y = self._indexes["node_y"]
+        edge_nodes = self._attrs["edge_node_connectivity"]
+        edge_nodes_attrs = conventions.DEFAULT_ATTRS["edge_node_connectivity"]
+
+        data_vars = {
+            self.name: 0,
+            edge_nodes: xr.DataArray(
+                data=self.edge_node_connectivity,
+                attrs=edge_nodes_attrs,
+                dims=(self.edge_dimension, "two"),
+            ),
+        }
+        dataset = xr.Dataset(data_vars, attrs={"Conventions": "CF-1.8 UGRID-1.0"})
+
+        if self._dataset:
+            dataset.update(self._dataset)
+        if other is not None:
+            dataset = dataset.merge(other)
+        if node_x not in dataset or node_y not in dataset:
+            dataset = self.assign_node_coords(dataset)
+
+        dataset[self.name].attrs = self._filtered_attrs(dataset)
+        return dataset
 
     def _clear_geometry_properties(self):
         """Clear all properties that may have been invalidated"""
@@ -205,20 +194,17 @@ def _clear_geometry_properties(self):
         self._edge_x = None
         self._edge_y = None
 
-    # These are all optional UGRID attributes. They are not computed by
-    # default, only when called upon.
     @property
     def topology_dimension(self):
         """Highest dimensionality of the geometric elements: 1"""
         return 1
 
-    def _get_dimension(self, dim):
-        key = f"{dim}_dimension"
-        if key not in self.topology_attrs:
-            self.topology_attrs[key] = ugrid_io.UGrid.mesh1d_get_attributes(self.name)[
-                key
-            ]
-        return self.topology_attrs[key]
+    @property
+    def dimensions(self):
+        return self.node_dimension, self.edge_dimension
+
+    # These are all optional UGRID attributes. They are not computed by
+    # default, only when called upon.
 
     @property
     def mesh(self) -> "mk.Mesh1d":  # type: ignore # noqa
@@ -309,6 +295,11 @@ 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 _validate_indexer(self, indexer) -> Tuple[float, float]:
         if isinstance(indexer, slice):
             if indexer.step is not None:
diff --git a/xugrid/ugrid/ugrid2d.py b/xugrid/ugrid/ugrid2d.py
index d25c26e45..a1f55a31f 100644
--- a/xugrid/ugrid/ugrid2d.py
+++ b/xugrid/ugrid/ugrid2d.py
@@ -1,3 +1,4 @@
+from itertools import chain
 from typing import Any, Tuple, Union
 
 import numpy as np
@@ -10,7 +11,7 @@
 from .. import meshkernel_utils as mku
 from ..typing import BoolArray, FloatArray, FloatDType, IntArray, IntDType, SparseMatrix
 from ..voronoi import voronoi_topology
-from . import ugrid_io
+from . import conventions
 from .ugridbase import AbstractUgrid
 
 
@@ -61,13 +62,21 @@ def __init__(
         node_y: FloatArray,
         fill_value: int,
         face_node_connectivity: Union[IntArray, SparseMatrix],
+        name: str = "mesh2d",
         edge_node_connectivity: IntArray = None,
         dataset: xr.Dataset = None,
-        name: str = None,
+        indexes: dict[str, str] = None,
+        projected: bool = True,
         crs: Any = None,
     ):
+        if dataset is not None and indexes is None:
+            raise ValueError("indexes must be provided for a dataset")
+
         self.node_x = np.ascontiguousarray(node_x)
         self.node_y = np.ascontiguousarray(node_y)
+        self.fill_value = fill_value
+        self.name = name
+        self.projected = projected
 
         if isinstance(face_node_connectivity, np.ndarray):
             face_node_connectivity = face_node_connectivity
@@ -80,14 +89,22 @@ def __init__(
                 "face_node_connectivity should be an array of integers or a sparse matrix"
             )
 
-        self.fill_value = fill_value
         self.face_node_connectivity = connectivity.counterclockwise(
             face_node_connectivity, self.fill_value, self.node_coordinates
         )
-        if name is None:
-            name = ugrid_io.UGRID2D_DEFAULT_NAME
-        self.name = name
-        self.topology_attrs = None
+
+        # Store topology attributes
+        defaults = conventions.default_topology_attrs(name, self.topology_dimension)
+        if dataset is None:
+            self._attrs = defaults
+            x, y = defaults["node_coordinates"].split(" ")
+            self._indexes = {"node_x": x, "node_y": y}
+        else:
+            derived_dims = dataset.ugrid_roles.dimensions[name]
+            self._attrs = {**defaults, **derived_dims, **dataset[name].attrs}
+            self._indexes = indexes
+
+        self._dataset = dataset
 
         # Optional attributes, deferred initialization
         # Meshkernel
@@ -123,11 +140,6 @@ def __init__(
             import pyproj
 
             self.crs = pyproj.CRS.from_user_input(crs)
-        # Store dataset
-        if dataset is None:
-            dataset = self.topology_dataset()
-        self.dataset = dataset
-        self.topology_attrs = dataset[self.name].attrs
 
     def _clear_geometry_properties(self):
         """Clear all properties that may have been invalidated"""
@@ -151,8 +163,8 @@ def _clear_geometry_properties(self):
         self._voronoi_topology = None
         self._centroid_triangulation = None
 
-    @staticmethod
-    def from_dataset(dataset: xr.Dataset):
+    @classmethod
+    def from_dataset(cls, dataset: xr.Dataset, topology: str = None):
         """
         Extract the 2D UGRID topology information from an xarray Dataset.
 
@@ -163,7 +175,7 @@ def from_dataset(dataset: xr.Dataset):
 
         Returns
         -------
-        grid: Ugrid2d
+        grid: Ugrid1dAdapter
         """
         ds = dataset
         if not isinstance(ds, xr.Dataset):
@@ -171,120 +183,91 @@ def from_dataset(dataset: xr.Dataset):
                 "Ugrid should be initialized with xarray.Dataset. "
                 f"Received instead: {type(ds)}"
             )
+        if topology is None:
+            topology = cls._single_topology(ds)
+
+        indexes = {}
 
         # Collect names
-        mesh_topology = ugrid_io.get_topology_variable(ds)
-        ugrid_roles = ugrid_io.get_ugrid2d_variables(dataset, mesh_topology)
-        # Rename the arrays to their standard UGRID roles
-        topology_ds = ds[set(ugrid_roles.keys())]
-        ugrid_ds = topology_ds.rename(ugrid_roles)
-        # Coerce type and fill value
-        ugrid_ds["node_x"] = ugrid_ds["node_x"].astype(FloatDType)
-        ugrid_ds["node_y"] = ugrid_ds["node_y"].astype(FloatDType)
-        fill_value = -1
-        ugrid_ds["face_node_connectivity"] = ugrid_io.cast(
-            ugrid_ds["face_node_connectivity"], fill_value, IntDType
+        connectivity = ds.ugrid_roles.connectivity[topology]
+        coordinates = ds.ugrid_roles.coordinates[topology]
+        ugrid_vars = (
+            [topology]
+            + list(connectivity.values())
+            + list(chain.from_iterable(chain.from_iterable(coordinates.values())))
         )
-        if "edge_node_connectivity" in ugrid_ds:
-            edge_node_connectivity = (
-                ugrid_ds["edge_node_connectivity"].astype(IntDType).values
-            )
+        fill_value = -1
+
+        x_index = coordinates["node_coordinates"][0][0]
+        y_index = coordinates["node_coordinates"][1][0]
+        node_x_coordinates = ds[x_index].astype(FloatDType).values
+        node_y_coordinates = ds[y_index].astype(FloatDType).values
+
+        face_nodes = connectivity["face_node_connectivity"]
+        face_node_connectivity = cls._prepare_connectivity(
+            ds[face_nodes], fill_value, dtype=IntDType
+        ).values
+
+        edge_nodes = connectivity.get("edge_node_connectivity")
+        if edge_nodes:
+            edge_node_connectivity = cls._prepare_connectivity(
+                ds[edge_nodes], fill_value, dtype=IntDType
+            ).values
         else:
             edge_node_connectivity = None
 
-        return Ugrid2d(
-            ugrid_ds["node_x"].values,
-            ugrid_ds["node_y"].values,
-            fill_value,
-            ugrid_ds["face_node_connectivity"].values,
-            edge_node_connectivity,
-            topology_ds,
-            mesh_topology.name,
-        )
-
-    def topology_dataset(self):
-        """
-        Store the UGRID topology information in an xarray Dataset according to
-        the UGRID conventions.
+        indexes["node_x"] = x_index
+        indexes["node_y"] = y_index
+        projected = False  # TODO
 
-        Returns
-        -------
-        ugrid_topology: xr.Dataset
-        """
-        return ugrid_io.ugrid2d_dataset(
-            self.node_x,
-            self.node_y,
-            self.fill_value,
-            self.face_node_connectivity,
-            self._edge_node_connectivity,
-            self.name,
-            self.topology_attrs,
+        return cls(
+            node_x_coordinates,
+            node_y_coordinates,
+            fill_value,
+            face_node_connectivity,
+            name=topology,
+            edge_node_connectivity=edge_node_connectivity,
+            dataset=ds[ugrid_vars],
+            indexes=indexes,
+            projected=projected,
+            crs=None,
         )
 
-    def remove_topology(
-        self, obj: Union[xr.DataArray, xr.Dataset]
-    ) -> Union[xr.DataArray, xr.Dataset]:
-        """
-        Remove UGRID specific variables from the object.
-
-        Parameters
-        ----------
-        obj: Union[xr.DataArray, xr.Dataset]
-
-        Returns
-        -------
-        sanitized: Union[xr.DataArray, xr.Dataset]
-        """
-        return self._remove_topology(obj, ugrid_io.UGRID2D_TOPOLOGY_VARIABLES)
-
-    def topology_coords(self, obj: Union[xr.DataArray, xr.Dataset]) -> dict:
-        """
-        Return a dictionary with the coordinates required for the dimension of
-        the object.
-
-        Parameters
-        ----------
-        obj: Union[xr.DataArray, xr.Dataset]
+    def to_dataset(self, other: xr.Dataset = None) -> xr.Dataset:
+        node_x = self._indexes["node_x"]
+        node_y = self._indexes["node_y"]
+        face_nodes = self._attrs["face_node_connectivity"]
+        face_nodes_attrs = conventions.DEFAULT_ATTRS["face_node_connectivity"]
+        nmax_node_dim = self._attrs["max_face_nodes_dimension"]
+        edge_nodes = self._attrs["edge_node_connectivity"]
+        edge_nodes_attrs = conventions.DEFAULT_ATTRS["edge_node_connectivity"]
+
+        data_vars = {
+            self.name: 0,
+            face_nodes: xr.DataArray(
+                data=self.face_node_connectivity,
+                attrs=face_nodes_attrs,
+                dims=(self.face_dimension, nmax_node_dim),
+            ),
+        }
+        if self.edge_node_connectivity is not None:
+            data_vars[edge_nodes] = xr.DataArray(
+                data=self.edge_node_connectivity,
+                attrs=edge_nodes_attrs,
+                dims=(self.edge_dimension, "two"),
+            )
 
-        Returns
-        -------
-        coords: dict
-        """
-        coords = {}
-        dims = obj.dims
-        facedim = self.face_dimension
-        edgedim = self.edge_dimension
-        nodedim = self.node_dimension
-        attrs = self.topology_attrs
-
-        if facedim in dims:
-            # It's not a required attribute.
-            if "face_coordinates" in attrs:
-                name_x, name_y = attrs["face_coordinates"].split()
-            else:
-                name_x = f"{self.name}_face_x"
-                name_y = f"{self.name}_face_y"
-            coords[name_x] = (facedim, self.face_x)
-            coords[name_y] = (facedim, self.face_y)
-            attrs["face-coordinates"] = f"{name_x} {name_y}"
-
-        if edgedim in dims:
-            # It's not a required attribute.
-            if "edge_coordinates" in attrs:
-                name_x, name_y = attrs["edge_coordinates"].split()
-            else:
-                name_x = f"{self.name}_edge_x"
-                name_y = f"{self.name}_edge_y"
-            coords[name_x] = (edgedim, self.edge_x)
-            coords[name_y] = (edgedim, self.edge_y)
-            attrs["edge_coordinates"] = f"{name_x} {name_y}"
+        dataset = xr.Dataset(data_vars, attrs={"Conventions": "CF-1.8 UGRID-1.0"})
 
-        if nodedim in dims:
-            name_x, name_y = attrs["node_coordinates"].split()
-            coords[name_x] = (nodedim, self.node_x)
-            coords[name_y] = (nodedim, self.node_y)
+        if self._dataset:
+            dataset.update(self._dataset)
+        if other is not None:
+            dataset = dataset.merge(other)
+        if node_x not in dataset or node_y not in dataset:
+            dataset = self.assign_node_coords(dataset)
 
-        return coords
+        dataset[self.name].attrs = self._filtered_attrs(dataset)
+        return dataset
 
     # These are all optional/derived UGRID attributes. They are not computed by
     # default, only when called upon.
@@ -295,25 +278,21 @@ def n_face(self):
         """
         return self.face_node_connectivity.shape[0]
 
+    @property
+    def dimensions(self):
+        return self.node_dimension, self.edge_dimension, self.face_dimension
+
     @property
     def topology_dimension(self):
         """Highest dimensionality of the geometric elements: 2"""
         return 2
 
-    def _get_dimension(self, dim):
-        key = f"{dim}_dimension"
-        if key not in self.topology_attrs:
-            self.topology_attrs[key] = ugrid_io.UGrid.mesh2d_get_attributes(self.name)[
-                key
-            ]
-        return self.topology_attrs[key]
-
     @property
     def face_dimension(self):
         """
         Return the name of the face dimension.
         """
-        return self._get_dimension("face")
+        return self._attrs["face_dimension"]
 
     def _edge_connectivity(self):
         (
@@ -690,6 +669,11 @@ def rasterize(self, resolution: float) -> Tuple[FloatArray, FloatArray, IntArray
         y = np.arange(ymax - 0.5 * d, ymin, -d)
         return self.rasterize_like(x, y)
 
+    def isel(self, dim, indexer):
+        if dim != self.face_dimension:
+            raise NotImplementedError("Can only index face_dimension for Ugrid2D")
+        return self.topology_subset(indexer)
+
     def _validate_indexer(self, indexer) -> Union[slice, np.ndarray]:
         if isinstance(indexer, slice):
             s = indexer
diff --git a/xugrid/ugrid/ugrid_io.py b/xugrid/ugrid/ugrid_io.py
deleted file mode 100644
index d8c2b3d7d..000000000
--- a/xugrid/ugrid/ugrid_io.py
+++ /dev/null
@@ -1,350 +0,0 @@
-"""
-Helper functions for parsing and composing UGRID files
-
-Most of these functions could be replaced by a centralized UGRID library with a
-C API.
-"""
-import warnings
-from typing import Dict, NamedTuple, Tuple, Union
-
-import numpy as np
-import xarray as xr
-
-from ..typing import FloatArray, IntArray
-
-# Ugrid convention data:
-
-
-class UgridTopologyAttributes(NamedTuple):
-    coordinates: Tuple[str]
-    dimensions: Tuple[str]
-    connectivity: Tuple[str]
-
-
-UGRID1D_DEFAULT_NAME = "network1d"
-UGRID1D_TOPOLOGY_VARIABLES = UgridTopologyAttributes(
-    coordinates=("node_coordinates", "edge_coordinates"),
-    dimensions=("node_dimension", "edge_dimension"),
-    connectivity=("edge_node_connectivity",),
-)
-
-UGRID2D_DEFAULT_NAME = "mesh2d"
-UGRID2D_TOPOLOGY_VARIABLES = UgridTopologyAttributes(
-    coordinates=("node_coordinates", "face_coordinates", "edge_coordinates"),
-    dimensions=("node_dimension", "face_dimension", "edge_dimension"),
-    connectivity=(
-        "face_node_connectivity",
-        "edge_node_connectivity",
-        "face_edge_connectivity",
-        "face_face_connectivity",
-        "edge_face_connectivity",
-        "boundary_node_connectivity",
-    ),
-)
-
-
-class UGrid:
-    # TODO: to be replaced by https://github.com/Deltares/UGridPy UGrid class
-    @staticmethod
-    def network1d_get_attributes(name: str):
-        """Get a dictionary of network1d default attribute names and the corresponding default values.
-        Args:
-            name (str): The network1d name.
-        Returns:
-            dict: A dictionary containing the attribute names and the corresponding default values.
-        """
-
-        return {
-            "cf_role": "mesh_topology",
-            "long_name": "Topology data of 1D network",
-            "topology_dimension": 1,
-            "node_dimension": f"{name}_nNodes",
-            "node_coordinates": f"{name}_node_x {name}_node_y",
-            "node_id": f"{name}_node_id",
-            "node_long_name": f"{name}_node_long_name",
-            "edge_dimension": f"{name}_nEdges",
-            "edge_node_connectivity": f"{name}_edge_nodes",
-            "edge_length": f"{name}_edge_length",
-            "edge_id": f"{name}_edge_id",
-            "edge_long_name": f"{name}_edge_long_name",
-            "edge_geometry": f"{name}_edge_geometry",
-        }
-
-    @staticmethod
-    def mesh2d_get_attributes(name: str):
-        """Get a dictionary of mesh2d default attribute names.
-        Args:
-            name (str): The mesh2d name.
-        Returns:
-            dict: A dictionary containing the attribute names and the corresponding default values.
-        """
-        return {
-            "cf_role": "mesh_topology",
-            "long_name": "Topology data of 2D mesh",
-            "topology_dimension": 2,
-            "node_dimension": f"{name}_nNodes",
-            "node_coordinates": f"{name}_node_x {name}_node_y",
-            "edge_dimension": f"{name}_nEdges",
-            "edge_node_connectivity": f"{name}_edge_nodes",
-            "face_dimension": f"{name}_nFaces",
-            "face_node_connectivity": f"{name}_face_nodes",
-            "max_face_nodes_dimension": f"{name}_nMax_face_nodes",
-            "face_coordinates": f"{name}_face_x {name}_face_y",
-        }
-
-
-# Reading / parsing:
-
-
-def get_topology_variable(dataset):
-    variables = [
-        da for da in dataset.values() if da.attrs.get("cf_role") == "mesh_topology"
-    ]
-    n = len(variables)
-    if n == 0:
-        raise ValueError("dataset does not contain a mesh topology variable")
-    elif n != 1:
-        names = [v.name for v in variables]
-        raise NotImplementedError(
-            f"dataset should contain a single mesh topology variable, contains: {names}"
-        )
-    return variables[0]
-
-
-def check_dim(attrs, dimrole, dimname, varname):
-    if dimrole in attrs and attrs[dimrole] != dimname:
-        warnings.warn(
-            f"{dimrole} of attributes does not match with dimension of "
-            f"variable: {attrs[dimrole]} in attrs versus {dimname} in "
-            f"variable {varname}"
-        )
-        return False
-    return True
-
-
-def _extract_topology_variables(
-    dataset, mesh_topology, ugrid_attrs: UgridTopologyAttributes
-) -> Dict[str, str]:
-    """
-    This standardizes the names of the UGRID attributes as variable names in
-    the dataset.
-
-    Their original names are preserved in the original attrs associated with
-    the dataset.
-    """
-
-    attrs = mesh_topology.attrs
-    variables = {mesh_topology.name: "mesh_topology"}
-    for role in ugrid_attrs.connectivity:
-        name = attrs.get(role)
-        if name and name in dataset:
-            variables[name] = role
-            # Also set the dimension names here.
-            dim = role.split("_")[0]
-            dimname = dataset[name].dims[0]
-            dimrole = f"{dim}_dimension"
-            variables[dimname] = dimrole
-            check_dim(attrs, dimrole, dimname, name)
-            attrs[dimrole] = dimname
-
-    for role in ugrid_attrs.coordinates:
-        name = attrs.get(role)
-        if name:
-            name_x, name_y = name.split()
-            dim = role.split("_")[0]
-            dimrole = f"{dim}_dimension"
-
-            has_x = name_x in dataset
-            has_y = name_y in dataset
-
-            if has_x:
-                variables[name_x] = f"{dim}_x"
-            if has_y:
-                variables[name_y] = f"{dim}_y"
-
-            if has_x and has_y:
-                dimname_x = dataset[name_x].dims[0]
-                dimname_y = dataset[name_y].dims[0]
-                if dimname_x != dimname_y:
-                    raise ValueError(
-                        f"dimensions of {name_x} and {name_y} do not match:"
-                        f"{dimname_x} versus {dimname_y}"
-                    )
-
-                if dimrole in attrs:
-                    if attrs[dimrole] != dimname_x:
-                        raise ValueError(
-                            f"{dimrole} from connectivity does not match dimension "
-                            f"{name_x}: {attrs[dimrole]} versus {dimname_x}"
-                        )
-                else:
-                    variables[dimname_x] = dimrole
-                    attrs[dimrole] = dimname_x
-
-    return variables
-
-
-def get_ugrid1d_variables(dataset, mesh_topology) -> Dict[str, str]:
-    return _extract_topology_variables(
-        dataset, mesh_topology, UGRID1D_TOPOLOGY_VARIABLES
-    )
-
-
-def get_ugrid2d_variables(dataset, mesh_topology) -> Dict[str, str]:
-    return _extract_topology_variables(
-        dataset, mesh_topology, UGRID2D_TOPOLOGY_VARIABLES
-    )
-
-
-def cast(da: xr.DataArray, fill_value: Union[float, int], dtype: type) -> xr.DataArray:
-    """
-    Set the appropriate fill value and cast to dtype.
-    """
-    data = da.values
-    # Note: Xarray always returns an array as floating type with NaNs for the
-    # fill value if _FillValue is set!
-    is_fill = np.isnan(data)
-    data[is_fill] = fill_value
-    # This returns a copy.
-    data = data.astype(dtype, copy=True)
-    return da.copy(data=data)
-
-
-# Writing / composing:
-
-
-def ugrid1d_dataset(
-    node_x: FloatArray,
-    node_y: FloatArray,
-    fill_value: int,
-    edge_node_connectivity: IntArray,
-    name: str = None,
-    attrs: Dict[str, str] = None,
-) -> xr.Dataset:
-    """
-    This creates a new dataset to hold all the UGRID variables for a network
-    ("1D").
-
-    If an attrs dictionary is provided, its names will be used to label the
-    variables; if an attrs dict is not provide, default names are provided.
-    """
-    if name is None:
-        name = UGRID1D_DEFAULT_NAME
-
-    if attrs is None:
-        attrs = UGrid.network1d_get_attributes(name)
-
-    node_dimension = attrs["node_dimension"]
-    edge_dimension = attrs["edge_dimension"]
-    node_x_name, node_y_name = attrs["node_coordinates"].split()
-    edge_node_name = attrs["edge_node_connectivity"]
-
-    ds = xr.Dataset()
-    ds[name] = xr.DataArray(
-        data=0,
-        attrs=attrs,
-    )
-    ds = ds.assign_coords(
-        {
-            node_x_name: xr.DataArray(
-                data=node_x,
-                dims=[node_dimension],
-            )
-        }
-    )
-    ds = ds.assign_coords(
-        {
-            node_y_name: xr.DataArray(
-                data=node_y,
-                dims=[node_dimension],
-            )
-        }
-    )
-    ds[edge_node_name] = xr.DataArray(
-        data=edge_node_connectivity,
-        dims=[edge_dimension, "two"],
-        attrs={
-            "cf_role": "edge_node_connectivity",
-            "long_name": "Vertex nodes of edges",
-            "start_index": 0,
-            "_FillValue": fill_value,
-        },
-    )
-    ds.attrs = {"Conventions": "CF-1.8 UGRID-1.0"}
-    return ds
-
-
-def ugrid2d_dataset(
-    node_x: FloatArray,
-    node_y: FloatArray,
-    fill_value: int,
-    face_node_connectivity: IntArray,
-    edge_node_connectivity: IntArray = None,
-    name: str = None,
-    attrs: Dict[str, str] = None,
-) -> xr.Dataset:
-    """
-    This creates a new dataset to hold all the UGRID variables for a 2D mesh.
-
-    If an attrs dictionary is provided, its names will be used to label the
-    variables; if an attrs dict is not provide, default names are provided.
-    """
-    if name is None:
-        name = UGRID2D_DEFAULT_NAME
-
-    if attrs is None:
-        attrs = UGrid.mesh2d_get_attributes(name)
-
-    node_dimension = attrs["node_dimension"]
-    face_dimension = attrs["face_dimension"]
-    edge_dimension = attrs["edge_dimension"]
-    node_x_name, node_y_name = attrs["node_coordinates"].split()
-    edge_node_name = attrs["edge_node_connectivity"]
-    face_node_name = attrs["face_node_connectivity"]
-    nmax_face_name = attrs["max_face_nodes_dimension"]
-
-    ds = xr.Dataset()
-    ds[name] = xr.DataArray(
-        data=0,
-        attrs=attrs,
-    )
-    ds = ds.assign_coords(
-        {
-            node_x_name: xr.DataArray(
-                data=node_x,
-                dims=[node_dimension],
-            )
-        }
-    )
-    ds = ds.assign_coords(
-        {
-            node_y_name: xr.DataArray(
-                data=node_y,
-                dims=[node_dimension],
-            )
-        }
-    )
-    ds[face_node_name] = xr.DataArray(
-        data=face_node_connectivity,
-        dims=[face_dimension, nmax_face_name],
-        attrs={
-            "cf_role": "face_node_connectivity",
-            "long_name": "Vertex nodes of mesh faces (counterclockwise)",
-            "start_index": 0,
-            "_FillValue": fill_value,
-        },
-    )
-    if edge_node_connectivity is not None:
-        ds[edge_node_name] = xr.DataArray(
-            data=edge_node_connectivity,
-            dims=[edge_dimension, "two"],
-            attrs={
-                "cf_role": "edge_node_connectivity",
-                "long_name": "Vertex nodes of edges",
-                "start_index": 0,
-                "_FillValue": fill_value,
-            },
-        )
-
-    ds.attrs = {"Conventions": "CF-1.8 UGRID-1.0"}
-    return ds
diff --git a/xugrid/ugrid/ugridbase.py b/xugrid/ugrid/ugridbase.py
index 81e2d149f..5a04f3b27 100644
--- a/xugrid/ugrid/ugridbase.py
+++ b/xugrid/ugrid/ugridbase.py
@@ -8,16 +8,16 @@
 
 from .. import connectivity
 from ..typing import BoolArray, FloatArray, IntArray
-from . import ugrid_io
+from . import conventions
 
 
 class AbstractUgrid(abc.ABC):
     @abc.abstractproperty
-    def topology_dimension(self):
+    def topology_dimension():
         """ """
 
-    @abc.abstractmethod
-    def _get_dimension(self):
+    @abc.abstractproperty
+    def dimensions():
         """ """
 
     @abc.abstractstaticmethod
@@ -25,28 +25,77 @@ def from_dataset():
         """ """
 
     @abc.abstractmethod
-    def sel(self):
+    def to_dataset():
         """ """
 
     @abc.abstractmethod
-    def topology_dataset(self):
+    def isel():
         """ """
 
     @abc.abstractmethod
-    def topology_subset(self):
+    def sel():
         """ """
 
     @abc.abstractmethod
-    def remove_topology(self):
+    def topology_subset():
         """ """
 
     @abc.abstractmethod
-    def topology_coords(self):
+    def _clear_geometry_properties():
         """ """
 
-    @abc.abstractmethod
-    def _clear_geometry_properties(self):
-        """ """
+    @staticmethod
+    def _single_topology(dataset: xr.Dataset):
+        topologies = dataset.ugrid_roles.topology
+        n_topology = len(topologies)
+        if n_topology == 0:
+            raise ValueError("Dataset contains no UGRID topology variable.")
+        elif n_topology > 1:
+            raise ValueError(
+                "Dataset contains {n_topology} topology variables, "
+                "please specify the topology variable name to use."
+            )
+        return topologies[0]
+
+    def _filtered_attrs(self, dataset: xr.Dataset):
+        """
+        Removes names that are not present in the dataset.
+        """
+        topodim = self.topology_dimension
+        attrs = self._attrs.copy()
+
+        for key in conventions._DIM_NAMES[topodim]:
+            if key in attrs:
+                if attrs[key] not in dataset.dims:
+                    attrs.pop(key)
+
+        for key in conventions._CONNECTIVITY_NAMES[topodim]:
+            if key in attrs:
+                if attrs[key] not in dataset:
+                    attrs.pop(key)
+
+        for coord in conventions._COORD_NAMES[topodim]:
+            if coord in attrs:
+                names = attrs[coord].split(" ")
+                present = [name for name in names if name in dataset]
+                if present:
+                    attrs[coord] = " ".join(present)
+                else:
+                    attrs.pop(coord)
+
+        return attrs
+
+    def __repr__(self):
+        if self._dataset:
+            return self._dataset.__repr__()
+        else:
+            return self.to_dataset.__repr__()
+
+    def equals(self, other):
+        if not isinstance(other, type(self)):
+            return False
+        # TODO: check values, etc.
+        return True
 
     def copy(self):
         """Creates deepcopy"""
@@ -55,12 +104,12 @@ def copy(self):
     @property
     def node_dimension(self):
         """Name of node dimension"""
-        return self._get_dimension("node")
+        return self._attrs["node_dimension"]
 
     @property
     def edge_dimension(self):
         """Name of edge dimension"""
-        return self._get_dimension("edge")
+        return self._attrs["edge_dimension"]
 
     @property
     def node_coordinates(self) -> FloatArray:
@@ -118,6 +167,28 @@ def bounds(self) -> Tuple[float, float, float, float]:
             self._ymax,
         )
 
+    @staticmethod
+    def _prepare_connectivity(
+        da: xr.DataArray, fill_value: Union[float, int], dtype: type
+    ) -> xr.DataArray:
+        start_index = da.attrs.get("start_index", 0)
+        if start_index not in (0, 1):
+            raise ValueError(f"start_index should be 0 or 1, received: {start_index}")
+
+        data = da.values
+        if "_FillValue" in da.attrs:
+            is_fill = data == da.attrs["_FillValue"]
+        else:
+            is_fill = np.isnan(data)
+
+        cast = data.astype(dtype, copy=True)
+        if start_index:
+            cast -= start_index
+        cast[is_fill] = fill_value
+        if (cast[~is_fill] < 0).any():
+            raise ValueError("connectivity contains negative values")
+        return da.copy(data=cast)
+
     def _topology_subset(
         self, indices: Union[BoolArray, IntArray], node_connectivity: IntArray
     ):
@@ -139,28 +210,92 @@ def _topology_subset(
             node_y = self.node_y[node_indices]
             return self.__class__(node_x, node_y, self.fill_value, new_connectivity)
 
-    def _remove_topology(
+    def set_node_coords(
         self,
-        obj: Union[xr.Dataset, xr.DataArray],
-        topology_variables: ugrid_io.UgridTopologyAttributes,
-    ) -> Union[xr.Dataset, xr.DataArray]:
+        node_x: str,
+        node_y: str,
+        obj: Union[xr.DataArray, xr.Dataset],
+        projected: bool = True,
+    ):
         """
-        Removes the grid topology data from a dataset. Use after creating an
-        Ugrid object from the dataset.
+        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.
+        """
+        if " " in node_x or " " in node_y:
+            raise ValueError("coordinate names may not contain spaces")
+
+        x = obj[node_x].values
+        y = obj[node_y].values
+
+        if (x.ndim != 1) or (x.size != self.n_node):
+            raise ValueError(
+                "shape of node_x does not match n_node of grid: "
+                f"{x.shape} versus {self.n_node}"
+            )
+        if (y.ndim != 1) or (y.size != self.n_node):
+            raise ValueError(
+                "shape of node_y does not match n_node of grid: "
+                f"{y.shape} versus {self.n_node}"
+            )
+
+        # Remove them, then append at the end.
+        node_coords = [
+            coord
+            for coord in self._attrs["node_coordinates"].split(" ")
+            if coord not in (node_x, node_y)
+        ]
+        node_coords.extend((node_x, node_y))
+
+        self._clear_geometry_properties()
+        self.node_x = np.ascontiguousarray(x)
+        self.node_y = np.ascontiguousarray(y)
+        self._attrs["node_coordinates"] = " ".join(node_coords)
+        self._indexes["node_x"] = node_x
+        self._indexes["node_y"] = node_y
+        self.projected = projected
+
+    def assign_node_coords(
+        self,
+        obj: Union[xr.DataArray, xr.Dataset],
+    ) -> Union[xr.DataArray, xr.Dataset]:
+        """
+        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.
+
+        Parameters
+        ----------
+        obj: xr.DataArray or xr.Dataset
+
+        Returns
+        -------
+        assigned (same type as obj)
         """
-        attrs = self.topology_attrs
-        names = []
-        if self.name in obj:
-            names.append(self.name)
-        for topology_attr in topology_variables.coordinates:
-            varname = attrs.get(topology_attr)
-            if varname and varname in obj:
-                names.extend(varname.split())
-        for topology_attr in topology_variables.connectivity:
-            varname = attrs.get(topology_attr)
-            if varname and varname in obj:
-                names.append(varname)
-        return obj.drop_vars(names)
+        xname = self._indexes["node_x"]
+        yname = self._indexes["node_y"]
+        x_attrs = conventions.DEFAULT_ATTRS["node_x"][self.projected]
+        y_attrs = conventions.DEFAULT_ATTRS["node_y"][self.projected]
+        coords = {
+            xname: xr.DataArray(
+                data=self.node_x,
+                dims=(self.node_dimension,),
+                attrs=x_attrs,
+            ),
+            yname: xr.DataArray(
+                data=self.node_y,
+                dims=(self.node_dimension,),
+                attrs=y_attrs,
+            ),
+        }
+        return obj.assign_coords(coords)
 
     @property
     def node_edge_connectivity(self) -> csr_matrix:
diff --git a/xugrid/ugrid_dataset.py b/xugrid/ugrid_dataset.py
index f6ccff7e4..ea09ab37b 100644
--- a/xugrid/ugrid_dataset.py
+++ b/xugrid/ugrid_dataset.py
@@ -1,20 +1,25 @@
+import abc
 import types
 from functools import wraps
-from typing import Any, Callable, Union
+from itertools import chain
+from typing import Any, Callable, Sequence, 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
+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 Ugrid1d, Ugrid2d, grid_from_geodataframe
+from .ugrid import AbstractUgrid, Ugrid2d, grid_from_dataset, grid_from_geodataframe
+
+UgridType = Type[AbstractUgrid]
 
 
 def xarray_wrapper(func, grid):
@@ -59,11 +64,12 @@ class UgridDataArray(DataArrayOpsMixin, DunderForwardMixin):
     Wraps an xarray DataArray, adding UGRID topology.
     """
 
-    def __init__(self, obj: xr.DataArray, grid: Union[Ugrid1d, Ugrid2d] = None):
-        if grid is None:
-            grid = Ugrid2d.from_dataset(obj)
+    def __repr__(self):
+        return self.obj.__repr__()
+
+    def __init__(self, obj: xr.DataArray, grid: UgridType):
         self.grid = grid
-        self.obj = obj.assign_coords(self.grid.topology_coords(obj))
+        self.obj = obj
 
     def __getitem__(self, key):
         """
@@ -85,6 +91,9 @@ 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)
@@ -120,7 +129,7 @@ def ugrid(self):
         UGRID Accessor. This "accessor" makes operations using the UGRID
         topology available.
         """
-        return UgridAccessor(self.obj, self.grid)
+        return UgridDataArrayAccessor(self.obj, self.grid)
 
     @staticmethod
     def from_structured(da: xr.DataArray):
@@ -152,42 +161,173 @@ def from_structured(da: xr.DataArray):
         return UgridDataArray(face_da, grid)
 
 
+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():
+        """ """
+
+    @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.
+        dim, ugrid, index, coords = grid.sel(x=x, y=y)
+        result = obj.isel({dim: index})
+
+        if not ugrid:
+            return result.assign_coords(coords)
+
+        grid = grid.topology_subset(index)
+        if isinstance(obj, xr.DataArray):
+            return UgridDataArray(result, grid)
+        elif isinstance(obj, xr.Dataset):
+            return UgridDataset(result, grid)
+        else:
+            raise TypeError(
+                f"Expected UgridDataArray or UgridDataset, got {type(result).__name__}"
+            )
+
+    @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 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 __init__(self, obj: xr.Dataset = None, grid: Ugrid2d = None):
-        if grid is None:
+    def __repr__(self):
+        return self.obj.__repr__()
+
+    def __init__(
+        self,
+        obj: xr.Dataset = None,
+        grids: Union[UgridType, Sequence[UgridType]] = None,
+    ):
+        if grids is None:
             if obj is None:
-                raise ValueError("At least either obj or grid is required")
-            # TODO: find out whether it's 1D or 2D topology
-            self.grid = Ugrid2d.from_dataset(obj)
+                raise ValueError("At least either obj or grids is required")
+            if not isinstance(obj, xr.Dataset):
+                raise TypeError(
+                    "UgridDataset should be initialized with xarray.Dataset. "
+                    f"Received instead: {type(obj)}"
+                )
+            topologies = obj.ugrid_roles.topology
+            grids = [grid_from_dataset(obj, topology) for topology in topologies]
         else:
-            self.grid = grid
+            # Make sure it's a new list
+            if isinstance(grids, (list, tuple, set)):
+                grids = [grid for grid in grids]
+            else:  # not iterable
+                grids = [grids]
 
         if obj is None:
             ds = xr.Dataset()
         else:
-            ds = self.grid.remove_topology(obj)
-        self.obj = ds.assign_coords(self.grid.topology_coords(ds))
+            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)
+
+        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):
-            return UgridDataArray(result, self.grid)
+            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, self.grid)
+            return UgridDataset(result, item_grids)
         else:
             return result
 
     def __setitem__(self, key, value):
         # TODO: check with topology
+        append = False
         if isinstance(value, UgridDataArray):
+            # Check if the dimensions occur in self.
+            # if they don't, the grid should be added.
+            matching_dims = [
+                dim for dim in value.grid.dimensions if dim in self.obj.dims
+            ]
+            if matching_dims:
+                # If they do match: the grids should match.
+                grids = {dim: grid for grid in self.grids for dim in grid.dimensions}
+                grid_to_match = grids[matching_dims[0]]
+                if not grid_to_match.equals(value.grid):
+                    raise ValueError("grids do not match")
+            else:
+                append = True
+
             self.obj[key] = value.obj
+            if append:
+                self.grids.append(value.grid)
         else:
             self.obj[key] = value
 
@@ -195,13 +335,23 @@ 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):
-            return UgridDataArray(result, self.grid)
+            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):
-            return UgridDataset(result, self.grid)
+            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.grid)
+            return xarray_wrapper(result, self.grids)
         else:
             return result
 
@@ -209,7 +359,7 @@ def _unary_op(
         self,
         f: Callable,
     ):
-        return UgridDataset(self.obj._unary_op(f), self.grid)
+        return UgridDataset(self.obj._unary_op(f), self.grids)
 
     def _binary_op(
         self,
@@ -219,12 +369,12 @@ def _binary_op(
     ):
         if isinstance(other, (UgridDataArray, UgridDataset)):
             other = other.obj
-        return UgridDataset(self.obj._binary_op(other, f, reflexive), self.grid)
+        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.grid)
+        return UgridDataset(self.obj._inplace_binary_op(other, f), self.grids)
 
     @property
     def ugrid(self):
@@ -232,7 +382,7 @@ def ugrid(self):
         UGRID Accessor. This "accessor" makes operations using the UGRID
         topology available.
         """
-        return UgridAccessor(self.obj, self.grid)
+        return UgridDatasetAccessor(self.obj, self.grids)
 
     @staticmethod
     def from_geodataframe(geodataframe: "geopandas.GeoDataFrame"):  # type: ignore # noqa
@@ -249,22 +399,236 @@ def from_geodataframe(geodataframe: "geopandas.GeoDataFrame"):  # type: ignore #
         """
         grid = grid_from_geodataframe(geodataframe)
         ds = xr.Dataset.from_dataframe(geodataframe.drop("geometry", axis=1))
-        return UgridDataset(ds, grid)
+        return UgridDataset(ds, [grid])
+
+
+class UgridDatasetAccessor(AbstractUgridAccessor):
+    def __init__(self, obj: xr.Dataset, grids: Sequence[UgridType]):
+        self.obj = obj
+        self.grids = grids
+
+    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 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 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.
 
-class UgridAccessor:
+        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: Union[xr.Dataset, xr.DataArray], grid: Ugrid2d):
+    def __init__(self, obj: xr.DataArray, grid: UgridType):
         self.obj = obj
         self.grid = grid
 
     plot = UncachedAccessor(_PlotMethods)
 
-    def isel(self, indexer):
+    def assign_node_coords(self) -> xr.DataArray:
+        """
+        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 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.
 
@@ -276,24 +640,24 @@ def isel(self, indexer):
         -------
         indexed: Union[UgridDataArray, UgridDataset]
         """
-        if self.grid.topology_dimension == 1:
-            dim = self.grid.edge_dimension
-        elif self.grid.topology_dimension == 2:
-            dim = self.grid.face_dimension
-        else:
-            raise NotImplementedError
+        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})
-        result = result.assign_coords({f"{dim}_index": (dim, indexer)})
-        grid = self.grid.topology_subset(indexer)
+        grid = self.grid.isel(dim, indexer)
         if isinstance(self.obj, xr.DataArray):
             return UgridDataArray(result, grid)
-        elif isinstance(self.obj, xr.Dataset):
-            return UgridDataset(result, grid)
         else:
-            raise TypeError(
-                f"Expected UgridDataArray or UgridDataset, got {type(result).__name__}"
-            )
+            raise TypeError(f"Expected UgridDataArray, got {type(result).__name__}")
 
     def sel(self, x=None, y=None):
         """
@@ -324,24 +688,7 @@ def sel(self, x=None, y=None):
         -------
         selection: Union[UgridDataArray, UgridDataset, xr.DataArray, xr.Dataset]
         """
-        # TODO: also do vectorized indexing like xarray?
-        # Might not be worth it, as orthogonal and vectorized indexing are
-        # quite confusing.
-        dim, ugrid, index, coords = self.grid.sel(x=x, y=y)
-        result = self.obj.isel({dim: index})
-
-        if not ugrid:
-            return result.assign_coords(coords)
-
-        grid = self.grid.topology_subset(index)
-        if isinstance(self.obj, xr.DataArray):
-            return UgridDataArray(result, grid)
-        elif isinstance(self.obj, xr.Dataset):
-            return UgridDataset(result, grid)
-        else:
-            raise TypeError(
-                f"Expected UgridDataArray or UgridDataset, got {type(result).__name__}"
-            )
+        return self._sel(self.obj, self.grid, x, y)
 
     def sel_points(self, x, y):
         """
@@ -356,11 +703,7 @@ def sel_points(self, x, y):
         -------
         points: Union[xr.DataArray, xr.Dataset]
         """
-        if self.grid.topology_dimension != 2:
-            raise NotImplementedError
-        dim, _, index, coords = self.grid.sel_points(x, y)
-        result = self.obj.isel({dim: index})
-        return result.assign_coords(coords)
+        return self._sel_points(self.obj, self.grid, x, y)
 
     def _raster(self, x, y, index) -> xr.DataArray:
         index = index.ravel()
@@ -417,11 +760,16 @@ 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.crs
+        return {self.grid.name: self.grid.crs}
 
     def to_geodataframe(
-        self, dim: str = None, name: str = None, dim_order=None
+        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
@@ -454,19 +802,17 @@ def to_geodataframe(
         """
         import geopandas as gpd
 
-        if isinstance(self.obj, xr.DataArray):
-            dim = self.obj.dims[-1]
-            if self.obj.name is None:
-                ds = self.obj.to_dataset(name=name)
-            else:
-                ds = self.obj.to_dataset()
+        dim = self.obj.dims[-1]
+        if name is not None:
+            ds = self.obj.to_dataset(name=name)
         else:
-            ds = self.obj
+            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 mult-index to generate them as regular columns.
+        # 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)
@@ -579,18 +925,10 @@ def reverse_cuthill_mckee(self):
         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.
-        if isinstance(self.obj, xr.DataArray):
-            return UgridDataArray(
-                reordered_data,
-                reordered_grid,
-            )
-        elif isinstance(self.obj, xr.Dataset):
-            return UgridDataset(
-                reordered_data,
-                reordered_grid,
-            )
-        else:
-            raise ValueError("object should be a xr.DataArray")
+        return UgridDataArray(
+            reordered_data,
+            reordered_grid,
+        )
 
     def laplace_interpolate(
         self,
@@ -645,8 +983,6 @@ def laplace_interpolate(
         """
         grid = self.grid
         da = self.obj
-        if isinstance(da, xr.Dataset):
-            raise NotImplementedError
         if grid.topology_dimension != 2:
             raise NotImplementedError
         if len(da.dims) > 1 or da.dims[0] != grid.face_dimension:
@@ -686,31 +1022,7 @@ def to_dataset(self):
         -------
         dataset: UgridDataset
         """
-        return self.grid.dataset.merge(self.obj)
-
-    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)
+        return self.grid.to_dataset(self.obj)
 
 
 # Wrapped IO methods
@@ -776,12 +1088,53 @@ def wrap_func_like(func):
     @wraps(func)
     def _like(other, *args, **kwargs):
         obj = func(other.obj, *args, **kwargs)
-        return type(other)(obj, other.grid)
+        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_conventions.py b/tests/test_conventions.py
new file mode 100644
index 000000000..4e71e744d
--- /dev/null
+++ b/tests/test_conventions.py
@@ -0,0 +1,145 @@
+from collections import ChainMap
+
+import pytest
+import xarray as xr
+
+import xugrid
+from xugrid.ugrid import conventions as cv
+
+
+def test_infer_xy_coords():
+    def assign_node_coords(ds: xr.Dataset, name: str, stdname: str = None):
+        if stdname:
+            attrs = {"standard_name": stdname}
+        else:
+            attrs = None
+
+        da = xr.DataArray([0.0, 1.0], dims=["node"], attrs=attrs)
+        return ds.assign_coords({name: da})
+
+    # All have standard names: take all
+    ds = xr.Dataset()
+    ds = assign_node_coords(ds, "x", "projection_x_coordinate")
+    ds = assign_node_coords(ds, "y", "projection_y_coordinate")
+    ds = assign_node_coords(ds, "lon", "longitude")
+    ds = assign_node_coords(ds, "lat", "latitude")
+    candidates = ["x", "y", "lon", "lat"]
+    x, y = cv._infer_xy_coords(ds, candidates)
+    assert x == ["x", "lon"]
+    assert y == ["y", "lat"]
+
+    # Only lon and lat have standard names: take them.
+    ds = xr.Dataset()
+    ds = assign_node_coords(ds, "x")
+    ds = assign_node_coords(ds, "y")
+    ds = assign_node_coords(ds, "lon", "longitude")
+    ds = assign_node_coords(ds, "lat", "latitude")
+    x, y = cv._infer_xy_coords(ds, candidates)
+    assert x == ["lon"]
+    assert y == ["lat"]
+
+    # Non have standard names: take the first two.
+    ds = xr.Dataset()
+    ds = assign_node_coords(ds, "x")
+    ds = assign_node_coords(ds, "y")
+    ds = assign_node_coords(ds, "lon")
+    ds = assign_node_coords(ds, "lat")
+    candidates = ["x", "y", "lon", "lat"]
+    with pytest.warns(UserWarning):
+        x, y = cv._infer_xy_coords(ds, candidates)
+    assert x == ["x"]
+    assert y == ["y"]
+
+    # Only one: error
+    ds = xr.Dataset()
+    ds = assign_node_coords(ds, "x", "projection_x_coordinate")
+    ds = assign_node_coords(ds, "y")
+    candidates = ["x", "y"]
+    with pytest.raises(cv.UgridCoordinateError):
+        x, y = cv._infer_xy_coords(ds, candidates)
+
+    # Only one: error
+    ds = xr.Dataset()
+    ds = assign_node_coords(ds, "x")
+    ds = assign_node_coords(ds, "y", "projection_y_coordinate")
+    candidates = ["x", "y"]
+    with pytest.raises(cv.UgridCoordinateError):
+        x, y = cv._infer_xy_coords(ds, candidates)
+
+
+class TestConventionsElevation:
+    @pytest.fixture(autouse=True)
+    def setup(self):
+        self.ds = xugrid.data.elevation_nl(xarray=True)
+        self.coordinates = {
+            "mesh2d": {
+                "node_coordinates": (
+                    ["mesh2d_node_x"],
+                    ["mesh2d_node_y"],
+                ),
+                "face_coordinates": (
+                    ["mesh2d_face_x"],
+                    ["mesh2d_face_y"],
+                ),
+            }
+        }
+        self.connectivity = {
+            "mesh2d": {
+                "face_node_connectivity": "mesh2d_face_nodes",
+            }
+        }
+        self.dimensions = {
+            "mesh2d": {
+                "edge_dimension": "mesh2d_nEdges",
+                "face_dimension": "mesh2d_nFaces",
+                "node_dimension": "mesh2d_nNodes",
+            },
+        }
+
+    def test_get_topology(self):
+        assert cv._get_topology(self.ds) == ["mesh2d"]
+
+    def test_get_coordinates(self):
+        ds = xugrid.data.elevation_nl(xarray=True)
+        actual = cv._get_coordinates(ds, ["mesh2d"])
+        assert actual == self.coordinates
+
+    def test_get_connectivity(self):
+        ds = xugrid.data.elevation_nl(xarray=True)
+        actual = cv._get_connectivity(ds, ["mesh2d"])
+        assert actual == self.connectivity
+
+    def test_get_dimensions(self):
+        ds = xugrid.data.elevation_nl(xarray=True)
+        connectivity = cv._get_connectivity(ds, ["mesh2d"])
+        coordinates = cv._get_coordinates(ds, ["mesh2d"])
+        actual = cv._get_dimensions(ds, ["mesh2d"], connectivity, coordinates)
+        assert actual == self.dimensions
+
+    def test_topology(self):
+        assert self.ds.ugrid_roles.topology == ["mesh2d"]
+
+    def test_coordinates(self):
+        assert self.ds.ugrid_roles.coordinates == self.coordinates
+
+    def test_dimensions(self):
+        assert self.ds.ugrid_roles.dimensions == self.dimensions
+
+    def test_connectivity(self):
+        assert self.ds.ugrid_roles.connectivity == self.connectivity
+
+    def test_getitem(self):
+        result = self.ds.ugrid_roles["mesh2d"]
+        assert isinstance(result, ChainMap)
+
+        with pytest.raises(KeyError):
+            self.ds.ugrid_roles["mesh1d"]
+
+        assert self.ds.ugrid_roles["mesh2d"]["node_coordinates"] == (
+            ["mesh2d_node_x"],
+            ["mesh2d_node_y"],
+        )
+
+    def test_repr(self):
+        result = self.ds.ugrid_roles.__repr__()
+        assert isinstance(result, str)
diff --git a/tests/test_plot.py b/tests/test_plot.py
index ffb5c635e..27a86058c 100644
--- a/tests/test_plot.py
+++ b/tests/test_plot.py
@@ -17,7 +17,7 @@ def setup(self):
         plt.close("all")
 
         self.ds = xugrid.data.disk()
-        self.grid = self.ds.ugrid.grid
+        self.grid = self.ds.ugrid.grids[0]
         self.node = self.ds["node_z"]
         self.edge = self.ds["edge_z"]
         self.face = self.ds["face_z"]
@@ -193,3 +193,16 @@ def test_default_size(self):
         h = fig.get_figheight()
         w = fig.get_figwidth()
         assert np.allclose(w / h, 1.26)
+
+    def test_error_dimension(self):
+        with pytest.raises(
+            ValueError, match="UgridDataArray contains non-topology dimensions"
+        ):
+            uda = xugrid.concat(
+                [
+                    self.face.assign_coords(time="2000-01-01"),
+                    self.face.assign_coords(time="2001-01-01"),
+                ],
+                dim="time",
+            )
+            uda.ugrid.plot()
diff --git a/tests/test_ugrid1d.py b/tests/test_ugrid1d.py
index 8c2047a2c..440bd78fe 100644
--- a/tests/test_ugrid1d.py
+++ b/tests/test_ugrid1d.py
@@ -15,10 +15,10 @@
 except ImportError:
     pass
 
-NAME = xugrid.ugrid.ugrid_io.UGRID1D_DEFAULT_NAME
+NAME = "network1d"
 
 
-def grid1d(dataset=None, name=None, crs=None):
+def grid1d(crs=None):
     xy = np.array(
         [
             [0.0, 0.0],
@@ -31,8 +31,6 @@ def grid1d(dataset=None, name=None, crs=None):
         node_y=xy[:, 1],
         fill_value=-1,
         edge_node_connectivity=np.array([[0, 1], [1, 2]]),
-        dataset=dataset,
-        name=name,
         crs=crs,
     )
     return grid
@@ -41,7 +39,7 @@ def grid1d(dataset=None, name=None, crs=None):
 def test_ugrid1d_init():
     grid = grid1d()
     assert grid.name == NAME
-    assert isinstance(grid.dataset, xr.Dataset)
+    assert grid._dataset is None
     assert grid.node_x.flags["C_CONTIGUOUS"]
     assert grid.node_y.flags["C_CONTIGUOUS"]
 
@@ -115,38 +113,9 @@ def test_to_crs():
     assert grid.crs == pyproj.CRS.from_epsg(4326)
 
 
-def test_ugrid1d_from_dataset():
+def test_to_dataset():
     grid = grid1d()
-    grid2 = xugrid.Ugrid1d.from_dataset(grid.dataset)
-    assert grid.dataset == grid2.dataset
-
-
-def test_remove_topology():
-
-    grid = grid1d()
-    ds = grid.dataset.copy()
-    ds["a"] = xr.DataArray(0)
-    actual = grid.remove_topology(ds)
-    print(actual)
-    assert set(actual.data_vars) == set(["a"])
-
-
-def test_topology_coords():
-    grid = grid1d()
-    ds = xr.Dataset()
-    ds["a"] = xr.DataArray([1, 2, 3], dims=[f"{NAME}_nNodes"])
-    ds["b"] = xr.DataArray([1, 2], dims=[f"{NAME}_nEdges"])
-    coords = grid.topology_coords(ds)
-    assert isinstance(coords, dict)
-    assert f"{NAME}_edge_x" in coords
-    assert f"{NAME}_edge_y" in coords
-    assert f"{NAME}_node_x" in coords
-    assert f"{NAME}_node_y" in coords
-
-
-def test_topology_dataset():
-    grid = grid1d()
-    ds = grid.topology_dataset()
+    ds = grid.to_dataset()
     assert isinstance(ds, xr.Dataset)
     assert f"{NAME}" in ds
     assert f"{NAME}_nNodes" in ds.dims
@@ -156,6 +125,14 @@ def test_topology_dataset():
     assert f"{NAME}_edge_nodes" in ds
 
 
+def test_ugrid1d_dataset_roundtrip():
+    grid = grid1d()
+    ds = grid.to_dataset()
+    grid2 = xugrid.Ugrid1d.from_dataset(grid.to_dataset())
+    assert isinstance(grid2._dataset, xr.Dataset)
+    assert grid2._dataset == ds
+
+
 def test_clear_geometry_properties():
     grid = grid1d()
     for attr in [
@@ -179,16 +156,11 @@ def test_topology_dimension():
     assert grid.topology_dimension == 1
 
 
-def test_get_dimension():
-    grid = grid1d()
-    assert grid._get_dimension("node") == f"{NAME}_nNodes"
-    assert grid._get_dimension("edge") == f"{NAME}_nEdges"
-
-
 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")
 
 
 @requires_meshkernel
diff --git a/tests/test_ugrid2d.py b/tests/test_ugrid2d.py
index cffb56a49..9dc814b77 100644
--- a/tests/test_ugrid2d.py
+++ b/tests/test_ugrid2d.py
@@ -16,7 +16,7 @@
 
 from . import requires_meshkernel
 
-NAME = xugrid.ugrid.ugrid_io.UGRID2D_DEFAULT_NAME
+NAME = "mesh2d"
 VERTICES = np.array(
     [
         [0.0, 0.0],  # 0
@@ -80,14 +80,13 @@
 NODE_FACE_CONNECTIVITY = sparse.coo_matrix((NFJ, (NFI, NFJ))).tocsr()
 
 
-def grid2d(dataset=None, name=None, crs=None):
+def grid2d(dataset=None, crs=None):
     grid = xugrid.Ugrid2d(
         node_x=VERTICES[:, 0],
         node_y=VERTICES[:, 1],
         fill_value=-1,
         face_node_connectivity=FACES,
         dataset=dataset,
-        name=name,
         crs=crs,
     )
     return grid
@@ -96,7 +95,7 @@ def grid2d(dataset=None, name=None, crs=None):
 def test_ugrid2d_init():
     grid = grid2d()
     assert grid.name == NAME
-    assert isinstance(grid.dataset, xr.Dataset)
+    assert grid._dataset is None
     assert grid.node_x.flags["C_CONTIGUOUS"]
     assert grid.node_y.flags["C_CONTIGUOUS"]
     assert grid._edge_node_connectivity is None
@@ -132,40 +131,9 @@ def test_to_crs():
     assert (~(grid.node_coordinates == reprojected.node_coordinates)).all()
 
 
-def test_ugrid2d_from_dataset():
+def test_to_dataset():
     grid = grid2d()
-    grid2 = xugrid.Ugrid2d.from_dataset(grid.dataset)
-    assert grid.dataset == grid2.dataset
-
-
-def test_remove_topology():
-    grid = grid2d()
-    ds = grid.dataset.copy()
-    ds["a"] = xr.DataArray(0)
-    actual = grid.remove_topology(ds)
-    print(actual)
-    assert set(actual.data_vars) == set(["a"])
-
-
-def test_topology_coords():
-    grid = grid2d()
-    ds = xr.Dataset()
-    ds["a"] = xr.DataArray([1, 2, 3], dims=[f"{NAME}_nNodes"])
-    ds["b"] = xr.DataArray([1, 2], dims=[f"{NAME}_nEdges"])
-    ds["c"] = xr.DataArray([1, 2], dims=[f"{NAME}_nFaces"])
-    coords = grid.topology_coords(ds)
-    assert isinstance(coords, dict)
-    assert f"{NAME}_edge_x" in coords
-    assert f"{NAME}_edge_y" in coords
-    assert f"{NAME}_node_x" in coords
-    assert f"{NAME}_node_y" in coords
-    assert f"{NAME}_face_x" in coords
-    assert f"{NAME}_face_y" in coords
-
-
-def test_topology_dataset():
-    grid = grid2d()
-    ds = grid.topology_dataset()
+    ds = grid.to_dataset()
     assert isinstance(ds, xr.Dataset)
     assert f"{NAME}" in ds
     assert f"{NAME}_nNodes" in ds.dims
@@ -175,6 +143,55 @@ def test_topology_dataset():
     assert f"{NAME}_face_nodes" in ds
 
 
+def test_ugrid2d_set_node_coords():
+    grid = grid2d()
+    ds = xr.Dataset()
+    lonvalues = VERTICES[:, 0] + 10.0
+    latvalues = VERTICES[:, 1] + 10.0
+    ds["lon"] = xr.DataArray(lonvalues, dims=[grid.node_dimension])
+    ds["lat"] = xr.DataArray(latvalues, dims=[grid.node_dimension])
+    ds["lon with space"] = ds["lon"]
+    ds["lat with space"] = ds["lat"]
+    ds["short_lon"] = xr.DataArray(np.arange(6.0), dims=["short_node"])
+    ds["long_lat"] = xr.DataArray(np.arange(8.0), dims=["long_node"])
+
+    with pytest.raises(ValueError, match="coordinate names may not contain spaces"):
+        grid.set_node_coords("lon with space", "lat with space", ds)
+    with pytest.raises(
+        ValueError, match="shape of node_x does not match n_node of grid: "
+    ):
+        grid.set_node_coords("short_lon", "lat", ds)
+    with pytest.raises(
+        ValueError, match="shape of node_y does not match n_node of grid: "
+    ):
+        grid.set_node_coords("lon", "long_lat", ds)
+
+    grid.set_node_coords("lon", "lat", ds, projected=False)
+    assert np.allclose(grid.node_x, lonvalues)
+    assert np.allclose(grid.node_y, latvalues)
+    assert grid._indexes["node_x"] == "lon"
+    assert grid._indexes["node_y"] == "lat"
+    assert not grid.projected
+
+
+def test_ugrid2d_dataset_roundtrip():
+    grid = grid2d()
+    ds = grid.to_dataset()
+    grid2 = xugrid.Ugrid2d.from_dataset(ds)
+    assert isinstance(grid2._dataset, xr.Dataset)
+    assert grid2._dataset == ds
+
+
+def test_assign_node_coords():
+    grid = grid2d()
+    ds = xr.Dataset()
+    # Place some data on the grid facets.
+    ds["a"] = xr.DataArray([1, 2, 3, 4, 5, 6, 7], dims=[f"{NAME}_nNodes"])
+    with_coords = grid.assign_node_coords(ds)
+    assert f"{NAME}_node_x" in with_coords
+    assert f"{NAME}_node_y" in with_coords
+
+
 def test_clear_geometry_properties():
     grid = grid2d()
     for attr in [
@@ -202,18 +219,12 @@ def test_topology_dimension():
     assert grid.topology_dimension == 2
 
 
-def test_get_dimension():
-    grid = grid2d()
-    assert grid._get_dimension("node") == f"{NAME}_nNodes"
-    assert grid._get_dimension("edge") == f"{NAME}_nEdges"
-    assert grid._get_dimension("face") == f"{NAME}_nFaces"
-
-
 def test_dimensions():
     grid = grid2d()
     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")
 
 
 def test_edge_node_connectivity():
diff --git a/tests/test_ugrid_dataset.py b/tests/test_ugrid_dataset.py
index 1fe644995..d14a211f8 100644
--- a/tests/test_ugrid_dataset.py
+++ b/tests/test_ugrid_dataset.py
@@ -36,11 +36,30 @@
     data=np.ones(GRID.n_face),
     dims=[GRID.face_dimension],
 )
-UGRID_DS = GRID.dataset.copy()
+UGRID_DS = GRID.to_dataset()
 UGRID_DS["a"] = DARRAY
 UGRID_DS["b"] = DARRAY * 2
 
 
+def ugrid1d_ds():
+    xy = np.array(
+        [
+            [0.0, 0.0],
+            [1.0, 1.0],
+            [2.0, 2.0],
+        ]
+    )
+    grid = xugrid.Ugrid1d(
+        node_x=xy[:, 0],
+        node_y=xy[:, 1],
+        fill_value=-1,
+        edge_node_connectivity=np.array([[0, 1], [1, 2]]),
+    )
+    ds = grid.to_dataset()
+    ds["a1d"] = xr.DataArray([1, 2, 3], dims=[grid.node_dimension])
+    return xugrid.UgridDataset(ds)
+
+
 class TestUgridDataArray:
     @pytest.fixture(autouse=True)
     def setup(self):
@@ -55,14 +74,8 @@ def test_dunder_forward(self):
         assert isinstance(int(self.uda[0]), int)
         assert isinstance(float(self.uda[0]), float)
 
-    def test_getitem(self):
-        assert isinstance(self.uda["mesh2d_face_x"], xugrid.UgridDataArray)
-
-    def test_setitem(self):
-        uda = self.uda.copy()
-        # This should be forwarded to the underlying xarray object
-        uda["mesh2d_face_x"] = 0.0
-        assert isinstance(uda["mesh2d_face_x"], xugrid.UgridDataArray)
+    def test_repr(self):
+        assert self.uda.__repr__() == self.uda.obj.__repr__()
 
     def test_getattr(self):
         # Get an attribute
@@ -71,13 +84,13 @@ def test_getattr(self):
         # should remain untouched.
         assert self.uda.dims == self.uda.ugrid.obj.dims
         assert isinstance(self.uda.data, np.ndarray)
-        # DataArrays are automatically wrapped
-        assert isinstance(self.uda.mesh2d_face_x, xugrid.UgridDataArray)
         # So are functions
         assert isinstance(self.uda.mean(), 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.UgridAccessor)
+        assert isinstance(self.uda.ugrid, xugrid.ugrid_dataset.UgridDataArrayAccessor)
 
     def test_from_structured(self):
         da = xr.DataArray([0.0, 1.0, 2.0], {"x": [5.0, 10.0, 15.0]}, ["x"])
@@ -127,11 +140,10 @@ def test_np_ops(self):
 
     # Accessor tests
     def test_isel(self):
-        actual = self.uda.ugrid.isel([0, 1])
+        actual = self.uda.ugrid.isel({GRID.face_dimension: [0, 1]})
         assert isinstance(actual, xugrid.UgridDataArray)
         assert actual.shape == (2,)
         assert actual.ugrid.grid.n_face == 2
-        assert "mesh2d_nFaces_index" in actual.coords
 
     def test_sel_points(self):
         with pytest.raises(ValueError, match="x and y must be 1d"):
@@ -182,11 +194,12 @@ def test_rasterize(self):
         assert np.allclose(actual["y"], y)
 
     def test_crs(self):
-        assert self.uda.ugrid.crs is None
+        crs = self.uda.ugrid.crs
+        assert crs == {"mesh2d": None}
 
     def test_to_geodataframe(self):
         with pytest.raises(ValueError, match="unable to convert unnamed"):
-            self.uda.ugrid.to_geodataframe("mesh2d_nFaces")
+            self.uda.ugrid.to_geodataframe()
         uda2 = self.uda.copy()
         uda2.ugrid.obj.name = "test"
         gdf = uda2.ugrid.to_geodataframe("mesh2d_nFaces")
@@ -250,10 +263,11 @@ def setup(self):
 
     def test_init(self):
         assert isinstance(self.uds.ugrid.obj, xr.Dataset)
-        assert isinstance(self.uds.ugrid.grid, xugrid.Ugrid2d)
-        assert "mesh2d_face_x" in self.uds.ugrid.obj
+        assert isinstance(self.uds.ugrid.grids[0], xugrid.Ugrid2d)
         # Try alternative initialization
-        uds = xugrid.UgridDataset(grid=GRID)
+        uds = xugrid.UgridDataset(grids=GRID)
+        assert isinstance(uds, xugrid.UgridDataset)
+        uds = xugrid.UgridDataset(grids=[GRID])
         assert isinstance(uds, xugrid.UgridDataset)
         uds["a"] = DARRAY
         assert "a" in uds.ugrid.obj
@@ -263,13 +277,15 @@ def test_init_from_dataset_only(self):
         assert isinstance(uds, xugrid.UgridDataset)
         assert "a" in uds.ugrid.obj
         assert "b" in uds.ugrid.obj
-        assert "mesh2d_face_nodes" in uds.ugrid.grid.dataset
+        assert "mesh2d_face_nodes" in uds.ugrid.grids[0].to_dataset()
         assert "mesh2d_face_nodes" not in uds.ugrid.obj
 
+    def test_repr(self):
+        assert self.uds.__repr__() == self.uds.obj.__repr__()
+
     def test_getitem(self):
         assert "a" in self.uds
         assert "b" in self.uds
-        assert "mesh2d_face_x" in self.uds
         assert isinstance(self.uds["a"], xugrid.UgridDataArray)
         assert isinstance(self.uds[["a", "b"]], xugrid.UgridDataset)
 
@@ -286,6 +302,8 @@ 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)
+        # obj should be accessible
+        assert isinstance(self.uds.obj, xr.Dataset)
 
     def test_unary_op(self):
         alltrue = self.uds.astype(bool)
@@ -303,7 +321,7 @@ def test_math(self):
         assert isinstance(actual, xugrid.UgridDataset)
 
     def test_ugrid_accessor(self):
-        assert isinstance(self.uds.ugrid, xugrid.ugrid_dataset.UgridAccessor)
+        assert isinstance(self.uds.ugrid, xugrid.ugrid_dataset.UgridDatasetAccessor)
 
     def test_from_geodataframe(self):
         xy = np.array(
@@ -324,10 +342,9 @@ def test_from_geodataframe(self):
 
     # Accessor tests
     def test_isel(self):
-        actual = self.uds.ugrid.isel([0, 1])
+        actual = self.uds.ugrid.isel({GRID.face_dimension: [0, 1]})
         assert isinstance(actual, xugrid.UgridDataset)
-        assert actual.ugrid.grid.n_face == 2
-        assert "mesh2d_nFaces_index" in actual.coords
+        assert actual.ugrid.grids[0].n_face == 2
         assert actual["a"].shape == (2,)
         assert actual["b"].shape == (2,)
 
@@ -360,35 +377,23 @@ def test_sel(self):
         assert isinstance(actual, xugrid.UgridDataset)
         assert actual["a"].shape == (2,)
         assert actual["b"].shape == (2,)
-        assert actual.ugrid.grid.n_face == 2
+        assert actual.ugrid.grids[0].n_face == 2
 
         actual = self.uds.ugrid.sel(x=slice(0, 1), y=slice(1, None))
         assert isinstance(actual, xugrid.UgridDataset)
         assert actual["a"].shape == (1,)
         assert actual["b"].shape == (1,)
-        assert actual.ugrid.grid.n_face == 1
+        assert actual.ugrid.grids[0].n_face == 1
 
     def test_crs(self):
-        assert self.uds.ugrid.crs is None
+        crs = self.uds.ugrid.crs
+        assert crs == {"mesh2d": None}
 
     def test_to_geodataframe(self):
-        gdf = self.uds.ugrid.to_geodataframe("mesh2d_nFaces")
+        gdf = self.uds.ugrid.to_geodataframe()
         assert isinstance(gdf, gpd.GeoDataFrame)
         assert (gdf.geometry.geom_type == "Polygon").all()
 
-    def test_connected_components(self):
-        actual = self.uds.ugrid.connected_components()
-        assert isinstance(actual, xugrid.UgridDataArray)
-        assert np.allclose(actual, 0)
-
-    def test_reverse_cuthill_mckee(self):
-        actual = self.uds.ugrid.reverse_cuthill_mckee()
-        assert isinstance(actual, xugrid.UgridDataset)
-
-    def test_laplace_interpolate(self):
-        with pytest.raises(NotImplementedError):
-            self.uds.ugrid.laplace_interpolate(direct_solve=True)
-
 
 def test_to_dataset():
     uds = xugrid.UgridDataset(UGRID_DS)
@@ -403,7 +408,7 @@ def test_open_dataset(tmp_path):
     back = xugrid.open_dataset(path)
     assert isinstance(back, xugrid.UgridDataset)
     assert "b" in back
-    assert "mesh2d_face_nodes" in back.ugrid.grid.dataset
+    assert "mesh2d_face_nodes" in back.ugrid.grids[0].to_dataset()
     assert "mesh2d_face_nodes" not in back.ugrid.obj
 
 
@@ -450,7 +455,7 @@ def test_zarr_roundtrip(tmp_path):
     assert isinstance(back, xugrid.UgridDataset)
     assert "a" in back
     assert "b" in back
-    assert "mesh2d_face_nodes" in back.ugrid.grid.dataset
+    assert "mesh2d_face_nodes" in back.ugrid.grids[0].to_dataset()
     assert "mesh2d_face_nodes" not in back.ugrid.obj
 
 
@@ -461,13 +466,13 @@ def test_func_like():
     assert isinstance(fullda, xugrid.UgridDataArray)
     assert (fullda == 2).all()
     # Topology should be untouched
-    assert fullda.ugrid.grid.dataset == uds.ugrid.grid.dataset
+    assert fullda.ugrid.grid.to_dataset() == uds.ugrid.grids[0].to_dataset()
 
     fullds = xugrid.full_like(uds, 2)
     assert isinstance(fullds, xugrid.UgridDataset)
     assert (fullds["a"] == 2).all()
     assert (fullds["b"] == 2).all()
-    assert fullds.ugrid.grid.dataset == uds.ugrid.grid.dataset
+    assert fullds.ugrid.grids[0].to_dataset == uds.ugrid.grids[0].to_dataset()
 
     fullda = xugrid.zeros_like(uds["a"])
     assert isinstance(fullda, xugrid.UgridDataArray)
@@ -476,3 +481,26 @@ def test_func_like():
     fullda = xugrid.ones_like(uds["a"])
     assert isinstance(fullda, xugrid.UgridDataArray)
     assert (fullda == 1).all()
+
+
+def test_concat():
+    uds = xugrid.UgridDataset(UGRID_DS)
+    uda = uds["a"]
+    uda1 = uda.assign_coords(layer=1)
+    uda2 = uda.assign_coords(layer=2)
+    result = xugrid.concat([uda1, uda2], dim="layer")
+    assert result.dims == ("layer", "mesh2d_nFaces")
+    assert np.array_equal(result["layer"], [1, 2])
+
+    uds1d = ugrid1d_ds()
+    uda3 = uds1d["a1d"].assign_coords(layer=2)
+    with pytest.raises(ValueError, match="All UgridDataArrays must have the same grid"):
+        xugrid.concat([uda1, uda3], dim="layer")
+
+
+def test_merge():
+    uds2d = xugrid.UgridDataset(UGRID_DS)
+    uds1d = ugrid1d_ds()
+    merged = xugrid.merge([uds2d, uds1d])
+    assert isinstance(merged, xugrid.UgridDataset)
+    assert len(merged.grids) == 2
diff --git a/tests/test_ugrid_io.py b/tests/test_ugrid_io.py
deleted file mode 100644
index a2eb7fd6a..000000000
--- a/tests/test_ugrid_io.py
+++ /dev/null
@@ -1,225 +0,0 @@
-import re
-
-import numpy as np
-import pytest
-import xarray as xr
-
-from xugrid.ugrid import ugrid_io
-
-
-def test_get_attributes():
-    d = ugrid_io.UGrid.network1d_get_attributes("network1d")
-    assert isinstance(d, dict)
-
-    d = ugrid_io.UGrid.mesh2d_get_attributes("network1d")
-    assert isinstance(d, dict)
-
-
-def test_get_topology_variable():
-    ds = xr.Dataset()
-    with pytest.raises(ValueError, match="dataset does not contain a mesh"):
-        ugrid_io.get_topology_variable(ds)
-
-    ds["mesh"] = xr.DataArray(0, attrs={"cf_role": "mesh_topology"})
-    assert ugrid_io.get_topology_variable(ds) == ds["mesh"]
-
-    ds["mesh2"] = xr.DataArray(0, attrs={"cf_role": "mesh_topology"})
-    with pytest.raises(
-        NotImplementedError,
-        match=re.escape(
-            "dataset should contain a single mesh topology variable, "
-            "contains: ['mesh', 'mesh2']"
-        ),
-    ):
-        ugrid_io.get_topology_variable(ds)
-
-
-def test_check_dim():
-    attrs = {"node_dimension": "node"}
-    dimrole = "node_dimension"
-    dimname = "node"
-    varname = "mesh"
-
-    assert ugrid_io.check_dim(attrs, dimrole, dimname, varname)
-    assert ugrid_io.check_dim({}, dimrole, dimname, varname)
-
-    with pytest.warns(UserWarning):
-        ok = ugrid_io.check_dim(attrs, dimrole, "n_nodes", varname)
-    assert not ok
-
-
-def test_cast():
-    da = xr.DataArray([1.0, 2.0, np.nan])
-    actual = ugrid_io.cast(da, fill_value=-1, dtype=int)
-    assert np.issubdtype(actual.dtype, np.integer)
-    assert np.array_equal(actual.values, [1, 2, -1])
-    # Check whether it's a copy, not a view
-    assert da.values is not actual.values
-
-    actual = ugrid_io.cast(da, fill_value=-1, dtype=float)
-    assert np.issubdtype(actual.dtype, np.floating)
-    assert np.array_equal(actual.values, [1.0, 2.0, -1.0])
-
-
-def test_ugrid1d_dataset__default_attrs():
-    for name in [None, "my-network"]:
-        ds = ugrid_io.ugrid1d_dataset(
-            node_x=np.array([0.0, 1.0]),
-            node_y=np.array([0.0, 1.0]),
-            fill_value=-1,
-            edge_node_connectivity=np.array([[0, 1]]),
-            name=name,
-            attrs=None,
-        )
-
-        if name is None:
-            name = ugrid_io.UGRID1D_DEFAULT_NAME
-
-        assert isinstance(ds, xr.Dataset)
-        assert f"{name}" in ds
-        assert f"{name}_nNodes" in ds.dims
-        assert f"{name}_nEdges" in ds.dims
-        assert f"{name}_node_x" in ds.coords
-        assert f"{name}_node_y" in ds.coords
-        assert f"{name}_edge_nodes" in ds
-        attrs = ds[f"{name}_edge_nodes"].attrs
-        assert attrs["cf_role"] == "edge_node_connectivity"
-        assert attrs["start_index"] == 0
-        assert attrs["_FillValue"] == -1
-
-        variables = ugrid_io.get_ugrid1d_variables(ds, ds[name])
-        assert isinstance(variables, dict)
-        assert variables[f"{name}_edge_nodes"] == "edge_node_connectivity"
-        assert variables[f"{name}_node_x"] == "node_x"
-        assert variables[f"{name}_node_y"] == "node_y"
-
-
-def test_ugrid1d_dataset():
-    # Test arbitrary names
-    attrs = {
-        "cf_role": "mesh_topology",
-        "long_name": "Topology data of 1D network",
-        "topology_dimension": 1,
-        "node_dimension": "vertex",
-        "node_coordinates": "vertex_x vertex_y",
-        "edge_dimension": "branch",
-        "edge_node_connectivity": "branch_vertices",
-    }
-    ds = ugrid_io.ugrid1d_dataset(
-        node_x=np.array([0.0, 1.0]),
-        node_y=np.array([0.0, 1.0]),
-        fill_value=-1,
-        edge_node_connectivity=np.array([[0, 1]]),
-        name="rivernetwork",
-        attrs=attrs,
-    )
-    assert isinstance(ds, xr.Dataset)
-    assert "rivernetwork" in ds
-    assert "vertex" in ds.dims
-    assert "branch" in ds.dims
-    assert "vertex_x" in ds.coords
-    assert "vertex_y" in ds.coords
-    assert "branch_vertices" in ds
-    attrs = ds["branch_vertices"].attrs
-    assert attrs["cf_role"] == "edge_node_connectivity"
-    assert attrs["start_index"] == 0
-    assert attrs["_FillValue"] == -1
-
-    variables = ugrid_io.get_ugrid1d_variables(ds, ds["rivernetwork"])
-    assert isinstance(variables, dict)
-    assert variables["branch_vertices"] == "edge_node_connectivity"
-    assert variables["vertex_x"] == "node_x"
-    assert variables["vertex_y"] == "node_y"
-
-
-def test_ugrid2d_dataset__default_attrs():
-    for name in [None, "my-mesh"]:
-        kwargs = dict(
-            node_x=np.array([0.0, 1.0, 0.0]),
-            node_y=np.array([0.0, 0.0, 1.0]),
-            fill_value=-1,
-            face_node_connectivity=np.array([[0, 1, 2]]),
-            name=name,
-            attrs=None,
-        )
-        ds = ugrid_io.ugrid2d_dataset(**kwargs)
-
-        if name is None:
-            name = ugrid_io.UGRID2D_DEFAULT_NAME
-
-        assert isinstance(ds, xr.Dataset)
-        assert f"{name}" in ds
-        assert f"{name}_nNodes" in ds.dims
-        assert f"{name}_nFaces" in ds.dims
-        assert f"{name}_nMax_face_nodes" in ds.dims
-        assert f"{name}_node_x" in ds.coords
-        assert f"{name}_node_y" in ds.coords
-        assert f"{name}_face_nodes" in ds
-        attrs = ds[f"{name}_face_nodes"].attrs
-        assert attrs["cf_role"] == "face_node_connectivity"
-        assert attrs["start_index"] == 0
-        assert attrs["_FillValue"] == -1
-
-        variables = ugrid_io.get_ugrid2d_variables(ds, ds[name])
-        assert isinstance(variables, dict)
-        assert variables[f"{name}_face_nodes"] == "face_node_connectivity"
-        assert variables[f"{name}_node_x"] == "node_x"
-        assert variables[f"{name}_node_y"] == "node_y"
-
-        kwargs["edge_node_connectivity"] = np.array([[0, 1], [1, 2], [2, 0]])
-        ds = ugrid_io.ugrid2d_dataset(**kwargs)
-        assert f"{name}_nEdges" in ds.dims
-        assert f"{name}_edge_nodes" in ds
-        attrs = ds[f"{name}_edge_nodes"].attrs
-        assert attrs["cf_role"] == "edge_node_connectivity"
-        assert attrs["start_index"] == 0
-        assert attrs["_FillValue"] == -1
-
-
-def test_ugrid2d_dataset():
-    # Test arbitrary names
-    attrs = {
-        "cf_role": "mesh_topology",
-        "long_name": "Topology data of 2D mesh",
-        "topology_dimension": 2,
-        "node_dimension": "vertex",
-        "node_coordinates": "vertex_x vertex_y",
-        "edge_dimension": "branch",
-        "edge_node_connectivity": "branch_vertices",
-        "face_dimension": "cell",
-        "face_node_connectivity": "cell_vertices",
-        "max_face_nodes_dimension": "node_nmax",
-    }
-    ds = ugrid_io.ugrid2d_dataset(
-        node_x=np.array([0.0, 1.0, 0.0]),
-        node_y=np.array([0.0, 0.0, 1.0]),
-        fill_value=-1,
-        face_node_connectivity=np.array([[0, 1, 2]]),
-        edge_node_connectivity=np.array([[0, 1], [1, 2], [2, 0]]),
-        name="unstructuredmesh",
-        attrs=attrs,
-    )
-    assert isinstance(ds, xr.Dataset)
-    assert "unstructuredmesh" in ds
-    assert "vertex" in ds.dims
-    assert "branch" in ds.dims
-    assert "cell" in ds.dims
-    assert "vertex_x" in ds.coords
-    assert "vertex_y" in ds.coords
-    assert "branch_vertices" in ds
-    assert "cell_vertices" in ds
-    attrs = ds["branch_vertices"].attrs
-    assert attrs["cf_role"] == "edge_node_connectivity"
-    assert attrs["start_index"] == 0
-    assert attrs["_FillValue"] == -1
-    attrs = ds["cell_vertices"].attrs
-    assert attrs["cf_role"] == "face_node_connectivity"
-    assert attrs["start_index"] == 0
-    assert attrs["_FillValue"] == -1
-
-    variables = ugrid_io.get_ugrid2d_variables(ds, ds["unstructuredmesh"])
-    assert isinstance(variables, dict)
-    assert variables["branch_vertices"] == "edge_node_connectivity"
-    assert variables["cell_vertices"] == "face_node_connectivity"
-    assert variables["vertex_x"] == "node_x"
-    assert variables["vertex_y"] == "node_y"