diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml
index 63251fa..10aa933 100644
--- a/.github/workflows/tests.yml
+++ b/.github/workflows/tests.yml
@@ -11,7 +11,7 @@ jobs:
     strategy:
       matrix:
         python-version: ["3.8", "3.9", "3.10", "3.11"]
-        os: [windows-latest, ubuntu-latest, macos-latest]
+        os: [ubuntu-latest]
       fail-fast: false
 
     steps:
diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index a88636b..d17d214 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -28,7 +28,7 @@ repos:
     - id: blackdoc
 
 - repo: https://github.com/charliermarsh/ruff-pre-commit
-  rev: v0.0.253
+  rev: v0.0.254
   hooks:
     - id: ruff
 
@@ -39,7 +39,7 @@ repos:
     language_version: python3
 
 - repo: https://github.com/codespell-project/codespell
-  rev: v2.2.2
+  rev: v2.2.4
   hooks:
     - id: codespell
       args:
diff --git a/notebooks/fetch_datasets.ipynb b/notebooks/fetch_datasets.ipynb
deleted file mode 100644
index ff689d4..0000000
--- a/notebooks/fetch_datasets.ipynb
+++ /dev/null
@@ -1,180 +0,0 @@
-{
- "cells": [
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "# Bathymetry"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "import pandas as pd\n",
-    "from oceans.datasets import get_depth\n",
-    "\n",
-    "\n",
-    "df = pd.DataFrame(\n",
-    "    data=[[-40, -20], [-32, -20]],\n",
-    "    columns=[\"lon\", \"lat\"],\n",
-    ")\n",
-    "\n",
-    "df[\"depth\"] = get_depth(df[\"lon\"], df[\"lat\"])\n",
-    "df"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "from oceans.datasets import etopo_subset, get_isobath\n",
-    "\n",
-    "\n",
-    "bbox = [-43, -30, -22, -17]\n",
-    "\n",
-    "segments = get_isobath(bbox=bbox, iso=-200, smoo=True)\n",
-    "lon, lat, bathy = etopo_subset(bbox=bbox, smoo=True)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "import matplotlib.pyplot as plt\n",
-    "import cartopy.crs as ccrs\n",
-    "\n",
-    "\n",
-    "fig, ax = plt.subplots(\n",
-    "    figsize=(12, 12),\n",
-    "    subplot_kw={\"projection\": ccrs.PlateCarree()}\n",
-    ")\n",
-    "\n",
-    "ax.plot(df[\"lon\"], df[\"lat\"], \"ro\")\n",
-    "for k, row in df.iterrows():\n",
-    "    ax.text(row[\"lon\"], row[\"lat\"], f\"{row['depth']:.2f} m\", color=\"white\")\n",
-    "\n",
-    "cs = ax.pcolormesh(lon, lat, bathy)\n",
-    "for segment in segments[0]:\n",
-    "    lines = ax.plot(segment[:, 0], segment[:, -1], \"darkorange\", linewidth=2)\n",
-    "\n",
-    "ax.coastlines(\"10m\")\n",
-    "ax.gridlines(draw_labels=True);"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "from oceans.datasets import woa_profile\n",
-    "\n",
-    "\n",
-    "cube = woa_profile(-143, 10, variable=\"temperature\", time_period=\"annual\", resolution=\"5\")\n",
-    "\n",
-    "fig, ax = plt.subplots(figsize=(4.25, 7))\n",
-    "z = cube.coord(axis=\"Z\").points\n",
-    "l = ax.plot(cube[0, :].data, z)\n",
-    "ax.grid(True)\n",
-    "ax.invert_yaxis()"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "# Extract a 2D surface -- Annual temperature climatology"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "import iris.plot as iplt\n",
-    "from oceans.datasets import woa_subset\n",
-    "from oceans.colormaps import cm\n",
-    "\n",
-    "\n",
-    "bbox = [2.5, 357.5, -87.5, 87.5]\n",
-    "cube = woa_subset(bbox, variable=\"temperature\", time_period=\"annual\", resolution=\"5\")\n",
-    "c = cube[0, 0, ...]  # Slice singleton time and first level.\n",
-    "\n",
-    "fig, ax = plt.subplots(figsize=(12, 12))\n",
-    "\n",
-    "cs = iplt.pcolormesh(c, cmap=cm.avhrr)\n",
-    "cbar = plt.colorbar(cs)"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "# Extract a square around the Mariana Trench averaging into a profile"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "scrolled": false
-   },
-   "outputs": [],
-   "source": [
-    "import warnings\n",
-    "import iris\n",
-    "from oceans.colormaps import get_color\n",
-    "\n",
-    "\n",
-    "months = \"Jan Feb Apr Mar May Jun Jul Aug Sep Oct Nov Dec\".split()\n",
-    "colors = get_color(len(months))\n",
-    "bbox = [-143, -141, 10, 12]\n",
-    "\n",
-    "fig, ax = plt.subplots(figsize=(5, 5))\n",
-    "\n",
-    "for month in months:\n",
-    "    cube = woa_subset(bbox, time_period=month, variable=\"temperature\", resolution=\"1\")\n",
-    "    with warnings.catch_warnings():\n",
-    "        warnings.simplefilter(\"ignore\")\n",
-    "        grid_areas = iris.analysis.cartography.area_weights(cube)\n",
-    "        c = cube.collapsed([\"longitude\", \"latitude\"], iris.analysis.MEAN, weights=grid_areas)\n",
-    "    z = c.coord(axis=\"Z\").points\n",
-    "    l = ax.plot(c[0, :].data, z, label=month, color=next(colors))\n",
-    "\n",
-    "ax.grid(True)\n",
-    "ax.invert_yaxis()\n",
-    "leg = ax.legend(loc=\"lower left\")\n",
-    "ax.set_ylim(200, 0);"
-   ]
-  }
- ],
- "metadata": {
-  "kernelspec": {
-   "display_name": "Python 3 (ipykernel)",
-   "language": "python",
-   "name": "python3"
-  },
-  "language_info": {
-   "codemirror_mode": {
-    "name": "ipython",
-    "version": 3
-   },
-   "file_extension": ".py",
-   "mimetype": "text/x-python",
-   "name": "python",
-   "nbconvert_exporter": "python",
-   "pygments_lexer": "ipython3",
-   "version": "3.11.0"
-  }
- },
- "nbformat": 4,
- "nbformat_minor": 2
-}
diff --git a/oceans/datasets.py b/oceans/datasets.py
index 63ef84e..c52ce90 100644
--- a/oceans/datasets.py
+++ b/oceans/datasets.py
@@ -1,3 +1,4 @@
+import functools
 import warnings
 
 import numpy as np
@@ -101,9 +102,10 @@ def _woa_url(variable, time_period, resolution):
     return url
 
 
+@functools.lru_cache(maxsize=256)
 def woa_profile(lon, lat, variable="temperature", time_period="annual", resolution="1"):
     """
-    Return an iris.cube instance from a World Ocean Atlas variable at a
+    Return a xarray DAtaset instance from a World Ocean Atlas variable at a
     given lon, lat point.
 
     Parameters
@@ -123,95 +125,88 @@ def woa_profile(lon, lat, variable="temperature", time_period="annual", resoluti
 
     Returns
     -------
-    Iris.cube instance with the climatology.
+    xr.Dataset instance with the climatology.
 
     Examples
     --------
     >>> import matplotlib.pyplot as plt
     >>> from oceans.datasets import woa_profile
-    >>> cube = woa_profile(
+    >>> woa = woa_profile(
     ...     -143, 10, variable="temperature", time_period="annual", resolution="5"
     ... )
     >>> fig, ax = plt.subplots(figsize=(2.25, 5))
-    >>> z = cube.coord(axis="Z").points
-    >>> l = ax.plot(cube[0, :].data, z)
+    >>> l = woa.plot(ax=ax, y="depth")
     >>> ax.grid(True)
     >>> ax.invert_yaxis()
 
     """
-    import iris
+    import cf_xarray  # noqa
+    import xarray as xr
 
     url = _woa_url(variable=variable, time_period=time_period, resolution=resolution)
-
-    with warnings.catch_warnings():
-        warnings.simplefilter("ignore")
-        cubes = iris.load_raw(url)
-
-    # TODO: should we be using `an` instead of `mn`?
     v = _woa_variable(variable)
-    cube = [c for c in cubes if c.var_name == f"{v}_mn"][0]
-    scheme = iris.analysis.Nearest()
-    sample_points = [("longitude", lon), ("latitude", lat)]
-    kw = {
-        "sample_points": sample_points,
-        "scheme": scheme,
-        "collapse_scalar": True,
-    }
-    return cube.interpolate(**kw)
+
+    ds = xr.open_dataset(url, decode_times=False)
+    ds = ds[f"{v}_mn"]
+    return ds.cf.sel({"X": lon, "Y": lat}, method="nearest")
 
 
+@functools.lru_cache(maxsize=256)
 def woa_subset(
-    bbox,
+    min_lon,
+    max_lon,
+    min_lat,
+    max_lat,
     variable="temperature",
     time_period="annual",
     resolution="5",
     full=False,
 ):
     """
-    Return an iris.cube instance from a World Ocean Atlas variable at a
+    Return an xarray Dataset instance from a World Ocean Atlas variable at a
     given lon, lat bounding box.
 
     Parameters
     ----------
-    bbox: list, tuple
-          minx, maxx, miny, maxy positions to extract.
+    min_lon, max_lon, min_lat, max_lat: positions to extract.
     See `woa_profile` for the other options.
 
     Returns
     -------
-    `iris.Cube` instance with the climatology.
+    `xr.Dataset` instance with the climatology.
 
     Examples
     --------
     >>> # Extract a 2D surface -- Annual temperature climatology:
-    >>> import iris.plot as iplt
     >>> import matplotlib.pyplot as plt
-    >>> from oceans.colormaps import cm
-    >>> bbox = [2.5, 357.5, -87.5, 87.5]
-    >>> cube = woa_subset(
-    ...     bbox, variable="temperature", time_period="annual", resolution="5"
+    >>> from cmcrameri import cm
+    >>> from oceans.datasets import woa_subset
+    >>> bbox = [-177.5, 177.5, -87.5, 87.5]
+    >>> woa = woa_subset(
+    ...     *bbox, variable="temperature", time_period="annual", resolution="5"
     ... )
-    >>> c = cube[0, 0, ...]  # Slice singleton time and first level.
-    >>> cs = iplt.pcolormesh(c, cmap=cm.avhrr)
-    >>> cbar = plt.colorbar(cs)
+    >>> cs = woa["t_mn"].sel(depth=0).plot(cmap=cm.lajolla)
 
     >>> # Extract a square around the Mariana Trench averaging into a profile.
-    >>> import iris
+    >>> import matplotlib.pyplot as plt
+    >>> import numpy as np
     >>> from oceans.colormaps import get_color
     >>> colors = get_color(12)
     >>> months = "Jan Feb Apr Mar May Jun Jul Aug Sep Oct Nov Dec".split()
+    >>> def area_weights_avg(woa):
+    ...     woa = woa["t_mn"].squeeze()
+    ...     weights = np.cos(np.deg2rad(woa["lat"])).where(~woa.isnull())
+    ...     weights /= weights.mean()
+    ...     return (woa * weights).mean(dim=["lon", "lat"])
+    ...
     >>> bbox = [-143, -141, 10, 12]
     >>> fig, ax = plt.subplots(figsize=(5, 5))
     >>> for month in months:
-    ...     cube = woa_subset(
-    ...         bbox, time_period=month, variable="temperature", resolution="1"
-    ...     )
-    ...     grid_areas = iris.analysis.cartography.area_weights(cube)
-    ...     c = cube.collapsed(
-    ...         ["longitude", "latitude"], iris.analysis.MEAN, weights=grid_areas
+    ...     woa = woa_subset(
+    ...         *bbox, time_period=month, variable="temperature", resolution="1"
     ...     )
-    ...     z = c.coord(axis="Z").points
-    ...     l = ax.plot(c[0, :].data, z, label=month, color=next(colors))
+    ...     profile = area_weights_avg(woa)
+    ...     l = profile.plot(ax=ax, y="depth", label=month, color=next(colors))
     ...
     >>> ax.grid(True)
     >>> ax.invert_yaxis()
@@ -219,27 +214,20 @@ def woa_subset(
     >>> _ = ax.set_ylim(200, 0)
 
     """
-    import iris
+    import cf_xarray  # noqa
+    import xarray as xr
 
-    v = _woa_variable(variable)
     url = _woa_url(variable, time_period, resolution)
-    cubes = iris.load_raw(url)
-    cubes = [
-        cube.intersection(longitude=(bbox[0], bbox[1]), latitude=(bbox[2], bbox[3]))
-        for cube in cubes
-    ]
-
-    with warnings.catch_warnings():
-        warnings.simplefilter("ignore")
-        cubes = iris.cube.CubeList(cubes)
-
+    ds = xr.open_dataset(url, decode_times=False)
+    ds = ds.cf.sel({"X": slice(min_lon, max_lon), "Y": slice(min_lat, max_lat)})
+    v = _woa_variable(variable)
     if full:
-        return cubes
-    else:
-        return [c for c in cubes if c.var_name == f"{v}_mn"][0]
+        return ds
+    return ds[[f"{v}_mn"]]  # always return a dataset
 
 
-def etopo_subset(bbox, tfile=None, smoo=False):
+@functools.lru_cache(maxsize=256)
+def etopo_subset(min_lon, max_lon, min_lat, max_lat, tfile=None, smoo=False):
     """
     Get a etopo subset.
     Should work on any netCDF with x, y, data
@@ -249,7 +237,7 @@ def etopo_subset(bbox, tfile=None, smoo=False):
     >>> from oceans.datasets import etopo_subset
     >>> import matplotlib.pyplot as plt
     >>> bbox = [-43, -30, -22, -17]
-    >>> lon, lat, bathy = etopo_subset(bbox=bbox, smoo=True)
+    >>> lon, lat, bathy = etopo_subset(*bbox, smoo=True)
     >>> fig, ax = plt.subplots()
     >>> cs = ax.pcolormesh(lon, lat, bathy)
 
@@ -262,6 +250,7 @@ def etopo_subset(bbox, tfile=None, smoo=False):
         lons = etopo.variables["x"][:]
         lats = etopo.variables["y"][:]
 
+        bbox = min_lon, max_lon, min_lat, max_lat
         imin, imax, jmin, jmax = _get_indices(bbox, lons, lats)
         lon, lat = np.meshgrid(lons[imin:imax], lats[jmin:jmax])
 
@@ -298,7 +287,7 @@ def get_depth(lon, lat, tfile=None):
         lat.min() - offset,
         lat.max() + offset,
     ]
-    lons, lats, bathy = etopo_subset(bbox, tfile=tfile, smoo=False)
+    lons, lats, bathy = etopo_subset(*bbox, tfile=tfile, smoo=False)
 
     return get_profile(lons, lats, bathy, lon, lat, mode="nearest", order=3)
 
@@ -314,17 +303,17 @@ def get_isobath(bbox, iso=-200, tfile=None, smoo=False):
     >>> import matplotlib.pyplot as plt
     >>> bbox = [-43, -30, -22, -17]
     >>> segments = get_isobath(bbox=bbox, iso=-200, smoo=True)
-    >>> lon, lat, bathy = etopo_subset(bbox=bbox, smoo=True)
+    >>> lon, lat, bathy = etopo_subset(*bbox, smoo=True)
     >>> fig, ax = plt.subplots()
     >>> cs = ax.pcolormesh(lon, lat, bathy)
-    >>> for segment in segments[0]:
+    >>> for segment in segments:
     ...     lines = ax.plot(segment[:, 0], segment[:, -1], "k", linewidth=2)
     ...
 
     """
     import contourpy
 
-    lon, lat, topo = etopo_subset(bbox, tfile=tfile, smoo=smoo)
+    lon, lat, topo = etopo_subset(*bbox, tfile=tfile, smoo=smoo)
 
     c = contourpy.contour_generator(
         lon,
@@ -339,6 +328,8 @@ def get_isobath(bbox, iso=-200, tfile=None, smoo=False):
     res = c.create_contour(iso)
     nseg = len(res) // 2
     segments = res[:nseg]
+    if len(segments) == 1:
+        segments = segments[0]
     return segments
 
 
diff --git a/requirements-dev.txt b/requirements-dev.txt
index 4a68520..baf312e 100644
--- a/requirements-dev.txt
+++ b/requirements-dev.txt
@@ -1,7 +1,7 @@
 black
 cartopy
+cf_xarray
 check-manifest
-iris
 nbsphinx
 netcdf4
 pandas
@@ -11,3 +11,4 @@ pytest-xdist
 scipy
 sphinx
 twine
+xarray
diff --git a/requirements.txt b/requirements.txt
index 9675995..d3a8bbc 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,3 +1,4 @@
+cmcrameri
 contourpy
 gsw
 matplotlib>=3.6.0