Repository: kmnhan/erlabpy (language: python)

Issue:
Make accessor namespace opt-in
### Description

It is recommended to provide a single accessor under the package namespace. (See https://github.com/pydata/xarray/issues/1080)

Currently, we have multiple callable accessors in the top level namespace of `xarray.DataArray`, but this is not good practice!

### Possible solutions

The best way is to move all these accessors to a single accessor namespace:

| Current           | New                  |
| ----------------- | -------------------- |
| `da.qplot`        | `da.er.qplot`        |
| `da.qshow`        | `da.er.qshow`        |
| `da.qsel`         | `da.er.qsel`         |
| `da.kspace`       | `da.er.kspace`       |
| `da.modelfit`     | `da.er.modelfit`     |
| `da.parallel_fit` | `da.er.parallel_fit` |
| `ds.modelfit`     | `ds.er.modelfit`     |

However, we lose the convenience of calling these functions directly. As a workaround, we should add an opt-in environment variable to enable the top-level accessors. This way, users can choose to use the old syntax if they prefer!
Integrate DataTree
### Description

Integrate our `io` modules and plugins with `xarray.DataTree` once it is merged with `xarray` core. This will allow us to represent any data in a single `DataTree` object, or even all data from a single experiment. 

The merge of `DataTree` with `xarray` core is currently in progress. It seems to be tracked [here](https://github.com/pydata/xarray/issues/8572).


Discussion / hints:

Fix patch:
diff --git a/docs/environment.yml b/docs/environment.yml
index f536e225..fba41872 100644
--- a/docs/environment.yml
+++ b/docs/environment.yml
@@ -32,7 +32,6 @@ dependencies:
   - nbsphinx
   - furo
   - sphinx-design
-  - pyperclip
   - ipywidgets
   - pip:
     - pyqt6>=6.2.2
diff --git a/docs/requirements.txt b/docs/requirements.txt
index 2860dd18..ed62c0c3 100644
--- a/docs/requirements.txt
+++ b/docs/requirements.txt
@@ -1,36 +1,4 @@
 dask>=2024.4.1
-h5netcdf>=1.2.0
-igor2>=0.5.6
-iminuit>=2.25.2
-joblib>=1.3.2
-lmfit>=1.2.0,!=1.3.0
-matplotlib>=3.8.0
-numba>=0.59.0
-numbagg>=0.8.1
-numpy>=1.26.0
-pyperclip>=1.8.2
 pyqt6>=6.2.2
-PyQt6-Qt6<6.7.0
-pyqtgraph>=0.13.1
-qtawesome>=1.3.1
-qtpy>=2.4.1
-scipy>=1.12.0
-tqdm>=4.66.2
-uncertainties>=3.1.4
-varname>=0.13.0
-xarray>=2024.02.0
-sphinx
-sphinx-autodoc-typehints
-setuptools
-nbsphinx
-furo
-sphinx-design
-sphinx-qt-documentation
-sphinx-copybutton
-sphinxcontrib-bibtex
-pybtex
-hvplot
-csaps
 ipykernel
-ipywidgets
--e .
+-e '.[viz,perf,misc,docs]'
diff --git a/docs/source/conf.py b/docs/source/conf.py
index 777e6dbb..80543380 100644
--- a/docs/source/conf.py
+++ b/docs/source/conf.py
@@ -71,9 +71,7 @@
 
 # based on numpy doc/source/conf.py
 def linkcode_resolve(domain, info):
-    """
-    Determine the URL corresponding to Python object
-    """
+    """Determine the URL corresponding to Python object."""
     if domain != "py":
         return None
 
@@ -123,23 +121,20 @@ def linkcode_resolve(domain, info):
 autosummary_generate = True
 autosummary_generate_overwrite = True
 
-autodoc_inherit_docstrings = False
-autodoc_typehints = "description"
-autodoc_typehints_format = "short"
 autodoc_class_signature = "mixed"
-autodoc_member_order = "bysource"
+autodoc_member_order = "groupwise"
 autodoc_default_options = {
     "members": True,
     "undoc-members": True,
-    # "exclude-members":("sigDataChanged",),
     "show-inheritance": False,
 }
+autodoc_typehints = "description"
 autodoc_typehints_description_target = "documented"
+# autodoc_type_aliases = {}
+autodoc_typehints_format = "short"
+autodoc_preserve_defaults = True
+autodoc_inherit_docstrings = False
 
-# autodoc_type_aliases = {
-# "numpy.float64": "float",
-# "float64": "float",
-# }
 
 # -- Napoleon settings -------------------------------------------------------
 
@@ -149,17 +144,26 @@ def linkcode_resolve(domain, info):
 # napoleon_include_private_with_doc = False
 napoleon_include_special_with_doc = True
 # napoleon_use_admonition_for_examples = False
-# napoleon_use_admonition_for_notes = False
+napoleon_use_admonition_for_notes = True
 # napoleon_use_admonition_for_references = False
-# napoleon_use_ivar = False
+napoleon_use_ivar = True
 napoleon_use_param = True
+napoleon_use_keyword = False
 # napoleon_use_rtype = True
 napoleon_preprocess_types = True
-# napoleon_type_aliases = {
-#     "array-like": "ndarray <numpy.ndarray>",
-#     "numpy.float64": "float",
-#     "float64": "float",
-# }
+
+napoleon_type_aliases = {
+    "np.float32": "float32",
+    "numpy.float32": "float32",
+    "np.float64": "float64",
+    "numpy.float64": "float64",
+    "xr.DataArray": "xarray.DataArray",
+    "array-like": "`array-like <numpy.typing.ArrayLike>`",
+    "array_like": "`array-like <numpy.typing.ArrayLike>`",
+    "ColorType": "`ColorType <matplotlib.typing.ColorType>`",
+    "RGBColorType": "`RGBColorType <matplotlib.typing.RGBColorType>`",
+    "RGBAColorType": "`RGBAColorType <matplotlib.typing.RGBAColorType>`",
+}
 napoleon_attr_annotations = True
 napoleon_custom_sections = [("Signals", "params_style")]
 
@@ -214,8 +218,9 @@ def top_level(children, data):
 
 
 class APSStyle(pybtex.style.formatting.unsrt.Style):
-    """
-    APS style for BibTeX formatting, adapted from the conf.py file of the `mitiq
+    """APS style for BibTeX formatting.
+
+    Adapted from the conf.py file of the `mitiq
     library<https://github.com/unitaryfund/mitiq>`_.
     """
 
@@ -465,14 +470,14 @@ def get_inproceedings_template(self, e):
             "class": "",
         },
     ],
-    "light_css_variables": {
-        "color-brand-primary": "#6d50bf",
-        "color-brand-content": "#6d50bf",
-    },
-    "dark_css_variables": {
-        "color-brand-primary": "#a180ff",
-        "color-brand-content": "#a180ff",
-    },
+    # "light_css_variables": {
+    #     "color-brand-primary": "#6d50bf",
+    #     "color-brand-content": "#6d50bf",
+    # },
+    # "dark_css_variables": {
+    #     "color-brand-primary": "#a180ff",
+    #     "color-brand-content": "#a180ff",
+    # },
 }
 
 # -- LaTeX options -----------------------------------------------------------
diff --git a/docs/source/contributing.rst b/docs/source/contributing.rst
index d41b985c..fa09197e 100644
--- a/docs/source/contributing.rst
+++ b/docs/source/contributing.rst
@@ -138,7 +138,7 @@ reflected in the package without having to reinstall it. Before installing:
 
    .. code-block:: sh
 
-     pip install -e . --config-settings editable_mode=compat
+     pip install -e ".[dev]" --config-settings editable_mode=compat
 
 
 Updating the editable installation
diff --git a/docs/source/getting-started.rst b/docs/source/getting-started.rst
index f15672b3..7cb475e1 100644
--- a/docs/source/getting-started.rst
+++ b/docs/source/getting-started.rst
@@ -11,6 +11,8 @@ conda environment, you can install ERLabPy with the conda command line tool: ::
 
   conda install -c conda-forge erlab
 
+Add any `optional dependencies`_ you want to install to the command above.
+
 .. hint::
 
   If you are using macOS, you might experience degraded performance with the
@@ -31,7 +33,14 @@ conda environment, you can install ERLabPy with the conda command line tool: ::
 
 If you don’t use conda, you can install ERLabPy with pip: ::
 
-  pip install erlab
+  python -m pip install erlab
+
+Optional dependency groups can be installed with the following commands: ::
+
+  python -m pip install erlab[viz]       # Install optional dependencies for visualization
+  python -m pip install erlab[perf]      # Install optional dependencies for performance
+  python -m pip install erlab[misc]      # Install miscellaneous optional dependencies
+  python -m pip install erlab[complete]  # Install all optional dependencies
 
 If you wish to install ERLabPy from source, see the :doc:`contributing`.
 
@@ -70,26 +79,14 @@ select the appropriate library based on what is installed.
 
 See the :doc:`user-guide/index` to start using ERLabPy!
 
-Notes on compatibility
-----------------------
-
-- ERLabPy is tested on Python 3.11 and 3.12. It is not guaranteed to work on older
-  versions of Python.
-- There are some `known compatibility issues
-  <https://github.com/kmnhan/erlabpy/issues/17>`_ with PyQt5 and PySide2, so it is
-  recommended to use the newer PyQt6 or PySide6 if possible.
-- If you meet any unexpected behaviour while using IPython's `autoreload extension
-  <https://ipython.readthedocs.io/en/stable/config/extensions/autoreload.html>`_, try
-  excluding the following modules: ::
-
-    %aimport -erlab.io.dataloader -erlab.accessors
+.. _optional dependencies:
 
 Optional dependencies
 ---------------------
 
-The following packages are optional dependencies that are not installed by
-default. They are used in only specific functions, or is not used at all but is
-listed just for convenience.
+The following packages are optional dependencies that are not installed by default. They
+are only used in specific functions, or is not used at all but is listed just for
+convenience.
 
 .. list-table::
     :header-rows: 1
@@ -100,6 +97,8 @@ listed just for convenience.
       - Description
     * - `csaps <https://github.com/espdev/csaps>`_
       - Multidimensional smoothing splines
+    * - `ipywidgets <https://github.com/jupyter-widgets/ipywidgets>`_
+      - Interactive widgets
     * - `hvplot <https://github.com/holoviz/hvplot>`_ and `bokeh
         <https://github.com/bokeh/bokeh>`_
       - Interactive plotting
@@ -107,11 +106,30 @@ listed just for convenience.
         `cmocean <https://matplotlib.org/cmocean/>`_, and
         `colorcet <https://colorcet.holoviz.org>`_
       - More colormaps!
+    * - `numbagg <https://github.com/numbagg/numbagg>`_ and `bottleneck
+        <https://github.com/pydata/bottleneck>`_
+      - Fast multidimensional aggregation, accelerates xarray
 
-For a full list of dependencies and optional dependencies, take a look at the ``[project]`` and ``[project.optional-dependencies]`` section in
-`pyproject.toml <https://github.com/kmnhan/erlabpy/blob/main/pyproject.toml>`_:
+For a full list of dependencies and optional dependencies, take a look at the
+``[project]`` and ``[project.optional-dependencies]`` section in `pyproject.toml
+<https://github.com/kmnhan/erlabpy/blob/main/pyproject.toml>`_:
 
 .. literalinclude:: ../../pyproject.toml
    :language: toml
    :start-at: dependencies = [
    :end-before: [project.urls]
+
+
+Notes on compatibility
+----------------------
+
+- ERLabPy is tested on Python 3.11 and 3.12. It is not guaranteed to work on older
+  versions of Python.
+- There are some `known compatibility issues
+  <https://github.com/kmnhan/erlabpy/issues/17>`_ with PyQt5 and PySide2, so it is
+  recommended to use the newer PyQt6 or PySide6 if possible.
+- If you meet any unexpected behaviour while using IPython's `autoreload extension
+  <https://ipython.readthedocs.io/en/stable/config/extensions/autoreload.html>`_, try
+  excluding the following modules: ::
+
+    %aimport -erlab.io.dataloader -erlab.accessors
diff --git a/docs/source/user-guide/curve-fitting.ipynb b/docs/source/user-guide/curve-fitting.ipynb
index dd9fac6c..807e6b71 100644
--- a/docs/source/user-guide/curve-fitting.ipynb
+++ b/docs/source/user-guide/curve-fitting.ipynb
@@ -45,7 +45,6 @@
                 "import matplotlib.pyplot as plt\n",
                 "import numpy as np\n",
                 "import xarray as xr\n",
-                "\n",
                 "from erlab.io.exampledata import generate_data"
             ]
         },
@@ -423,7 +422,7 @@
                 "model = MultiPeakModel(npeaks=4, peak_shapes=[\"lorentzian\"], fd=False, convolve=True)\n",
                 "\n",
                 "params = model.make_params(\n",
-                "    p0_center=-0.59,\n",
+                "    p0_center=-0.6,\n",
                 "    p1_center=-0.45,\n",
                 "    p2_center=0.45,\n",
                 "    p3_center=0.6,\n",
@@ -995,7 +994,7 @@
                 "    model,\n",
                 "    mdc,\n",
                 "    mdc.kx,\n",
-                "    p0_center=-0.59,\n",
+                "    p0_center=-0.6,\n",
                 "    p1_center=-0.45,\n",
                 "    p2_center=0.45,\n",
                 "    p3_center=0.6,\n",
diff --git a/docs/source/user-guide/indexing.ipynb b/docs/source/user-guide/indexing.ipynb
index 5ae0fd74..16ff86ed 100644
--- a/docs/source/user-guide/indexing.ipynb
+++ b/docs/source/user-guide/indexing.ipynb
@@ -47,7 +47,7 @@
             "source": [
                 "from erlab.io.exampledata import generate_data\n",
                 "\n",
-                "dat = generate_data(seed=1)"
+                "dat = generate_data(seed=1).T"
             ]
         },
         {
@@ -228,13 +228,93 @@
                 "Interpolation\n",
                 "-------------\n",
                 "\n",
-                ".. admonition:: Work in Progress\n",
-                "   :class: warning\n",
+                "In addition to the `powerful interpolation methods\n",
+                "<https://docs.xarray.dev/en/latest/user-guide/interpolation.html>`_ provided by\n",
+                ":mod:`xarray`, ERLabPy provides a convenient way to interpolate data along an arbitrary\n",
+                "path.\n",
+                "\n",
+                "Consider a Γ-M-K-Γ high symmetry path given as a list of kx and ky coordinates:"
+            ]
+        },
+        {
+            "cell_type": "code",
+            "execution_count": null,
+            "id": "4ab13824",
+            "metadata": {},
+            "outputs": [],
+            "source": [
+                "import erlab.plotting.erplot as eplt\n",
+                "import matplotlib.pyplot as plt\n",
+                "import numpy as np\n",
+                "\n",
+                "a = 6.97\n",
+                "kx = [0, 2 * np.pi / (a * np.sqrt(3)), 2 * np.pi / (a * np.sqrt(3)), 0]\n",
+                "ky = [0, 0, 2 * np.pi / (a * 3), 0]\n",
+                "\n",
                 "\n",
-                "   This part of the user guide is still under construction. For now, please refer to the\n",
-                "   `interpolation section of the xarray user guide\n",
-                "   <https://docs.xarray.dev/en/latest/user-guide/interpolation.html>`_. For the full\n",
-                "   list of packages and modules provided by ERLabPy, see :doc:`../reference`."
+                "dat.qsel(eV=-0.2).qplot(aspect=\"equal\", cmap=\"Greys\")\n",
+                "plt.plot(kx, ky, \"o-\")"
+            ]
+        },
+        {
+            "cell_type": "markdown",
+            "id": "be5e2a2d",
+            "metadata": {},
+            "source": [
+                "To interpolate the data along this path with a step of 0.01 Å$^{-1}$, we can use the following code:"
+            ]
+        },
+        {
+            "cell_type": "code",
+            "execution_count": null,
+            "id": "bf4ec3c2",
+            "metadata": {},
+            "outputs": [],
+            "source": [
+                "import erlab.analysis as era\n",
+                "\n",
+                "dat_sliced = era.interpolate.slice_along_path(\n",
+                "    dat, vertices={\"kx\": kx, \"ky\": ky}, step_size=0.01\n",
+                ")\n",
+                "dat_sliced"
+            ]
+        },
+        {
+            "cell_type": "markdown",
+            "id": "84418499",
+            "metadata": {},
+            "source": [
+                "We can see that the data has been interpolated along the path. The new coordinate `path` contains the distance along the path, and the dimensions `kx` and `ky` are now expressed in terms of `path`.\n",
+                "\n",
+                "The distance along the path can be calculated as the sum of the distances between consecutive points in the path."
+            ]
+        },
+        {
+            "cell_type": "code",
+            "execution_count": null,
+            "id": "c00bb21c",
+            "metadata": {},
+            "outputs": [],
+            "source": [
+                "dat_sliced.plot(cmap=\"Greys\")\n",
+                "eplt.fermiline()\n",
+                "\n",
+                "# Distance between each pair of consecutive points\n",
+                "distances = np.linalg.norm(np.diff(np.vstack([kx, ky]), axis=-1), axis=0)\n",
+                "seg_coords = np.concatenate(([0], np.cumsum(distances)))\n",
+                "\n",
+                "plt.xticks(seg_coords, labels=[\"Γ\", \"M\", \"K\", \"Γ\"])\n",
+                "plt.xlim(0, seg_coords[-1])\n",
+                "for seg in seg_coords[1:-1]:\n",
+                "    plt.axvline(seg, ls=\"--\", c=\"k\", lw=1)"
+            ]
+        },
+        {
+            "cell_type": "markdown",
+            "id": "5a9036d9",
+            "metadata": {},
+            "source": [
+                "You will learn more about plotting in the next section."
             ]
         }
     ],
@@ -254,7 +334,7 @@
             "name": "python",
             "nbconvert_exporter": "python",
             "pygments_lexer": "ipython3",
-            "version": "3.11.8"
+            "version": "3.11.9"
         }
     },
     "nbformat": 4,
diff --git a/docs/source/user-guide/io.ipynb b/docs/source/user-guide/io.ipynb
index 3dae1c84..d767c6f3 100644
--- a/docs/source/user-guide/io.ipynb
+++ b/docs/source/user-guide/io.ipynb
@@ -677,7 +677,14 @@
                 "    additional_attrs = {\n",
                 "        \"configuration\": 1,  # Experimental geometry. Required for momentum conversion\n",
                 "        \"sample_workfunction\": 4.3,\n",
-                "    }  # Any additional metadata you want to add to the data\n",
+                "    }\n",
+                "    # Any additional metadata you want to add to the data. Note that attributes defined\n",
+                "    # here will not be transformed into coordinates. If you wish to promote some fixed\n",
+                "    # attributes to coordinates, add them to additional_coords.\n",
+                "\n",
+                "    additional_coords = {}\n",
+                "    # Additional non-dimension coordinates to be added to the data, for instance the\n",
+                "    # photon energy for lab-based ARPES.\n",
                 "\n",
                 "    skip_validate = False\n",
                 "\n",
diff --git a/docs/source/user-guide/plotting.ipynb b/docs/source/user-guide/plotting.ipynb
index b0ed0369..85090826 100644
--- a/docs/source/user-guide/plotting.ipynb
+++ b/docs/source/user-guide/plotting.ipynb
@@ -277,48 +277,63 @@
             "source": [
                 "What if we want to plot multiple slices at once? We should create subplots to place the slices. *plt.subplots* is very useful in managing multiple axes and figures. If you are unfamiliar with the syntax, visit the [relevant matplotlib documentation](https://matplotlib.org/stable/gallery/subplots_axes_and_figures/subplots_demo.html).\n",
                 "\n",
-                "Suppose we want to plot constant energy surfaces at specific binding energies, say, at `[-0.4, -0.2, 0.0]`. We could create three subplots and iterate over the list, plotting each slice on each subplot. We would also want to label each plot with their binding energy. "
+                "Suppose we want to plot constant energy surfaces at specific binding energies, say, at `[-0.4, -0.2, 0.0]`. We could create three subplots and iterate over the axes."
             ]
         },
         {
-            "cell_type": "raw",
+            "cell_type": "code",
+            "execution_count": null,
             "metadata": {
                 "editable": true,
-                "raw_mimetype": "text/restructuredtext",
                 "slideshow": {
                     "slide_type": ""
                 },
                 "tags": []
             },
+            "outputs": [],
             "source": [
-                "We use three functions from ERLabPy to achieve this:\n",
-                ":func:`plot_array <erlab.plotting.general.plot_array>`\\ ,\n",
-                ":func:`label_subplot_properties <erlab.plotting.annotations.label_subplot_properties>`\\ , and\n",
-                ":func:`clean_labels <erlab.plotting.erplot.clean_labels>`\\ ."
+                "energies = [-0.4, -0.2, 0.0]\n",
+                "\n",
+                "fig, axs = plt.subplots(1, 3, layout=\"compressed\", sharey=True)\n",
+                "for energy, ax in zip(energies, axs):\n",
+                "    const_energy_surface = dat.qsel(eV=energy)\n",
+                "    eplt.plot_array(const_energy_surface, ax=ax, gamma=0.5, aspect=\"equal\")"
             ]
         },
         {
-            "cell_type": "code",
-            "execution_count": null,
+            "cell_type": "raw",
             "metadata": {
                 "editable": true,
+                "raw_mimetype": "text/restructuredtext",
                 "slideshow": {
                     "slide_type": ""
                 },
-                "tags": []
+                "tags": [],
+                "vscode": {
+                    "languageId": "raw"
+                }
             },
+            "source": [
+                "Here, we plotted each constant energy surface with :func:`plot_array\n",
+                "<erlab.plotting.general.plot_array>`. To remove the duplicated y axis labels and add\n",
+                "some annotations, we can use :func:`clean_labels <erlab.plotting.erplot.clean_labels>`\n",
+                "and :func:`label_subplot_properties\n",
+                "<erlab.plotting.annotations.label_subplot_properties>`:"
+            ]
+        },
+        {
+            "cell_type": "code",
+            "execution_count": null,
+            "metadata": {},
             "outputs": [],
             "source": [
-                "energies = [-0.4, -0.2, 0.0]\n",
-                "\n",
                 "fig, axs = plt.subplots(1, 3, layout=\"compressed\", sharey=True)\n",
-                "\n",
                 "for energy, ax in zip(energies, axs):\n",
                 "    const_energy_surface = dat.qsel(eV=energy)\n",
                 "    eplt.plot_array(const_energy_surface, ax=ax, gamma=0.5, aspect=\"equal\")\n",
                 "\n",
                 "eplt.clean_labels(axs)  # removes shared y labels\n",
-                "eplt.label_subplot_properties(axs, values=dict(Eb=energies))  # annotates energy"
+                "eplt.label_subplot_properties(axs, values={\"Eb\": energies})  # annotates energy"
             ]
         },
         {
@@ -399,7 +414,9 @@
             "cell_type": "markdown",
             "metadata": {},
             "source": [
-                "Here, we notice that the first two plots slices through regions with less spectral weight, so the color across the three subplots are not on the same scale. This may be misleading in some occasions where intensity across different slices are important. Luckily, we have a function that can unify the color limits across multiple axes. The same effect can be achieved by passing on `same_limits=True` to `plot_slices`."
+                "Here, we notice that the first two plots slices through regions with less spectral weight, so the color across the three subplots are not on the same scale. This may be misleading in some occasions where intensity across different slices are important. Luckily, we have a function that can unify the color limits across multiple axes.\n",
+                "\n",
+                "The same effect can be achieved by passing on `same_limits=True` to `plot_slices`."
             ]
         },
         {
@@ -428,7 +445,7 @@
             "outputs": [],
             "source": [
                 "fig, axs = eplt.plot_slices([dat], ky=[0.0, 0.1, 0.3], gamma=0.5, figsize=(6, 2))\n",
-                "eplt.unify_clim(axs, target=axs[0, 1])"
+                "eplt.unify_clim(axs, target=axs.flat[1])"
             ]
         },
         {
diff --git a/environment.yml b/environment.yml
index 40ae75cd..1da942fd 100644
--- a/environment.yml
+++ b/environment.yml
@@ -22,5 +22,4 @@ dependencies:
   - uncertainties>=3.1.4
   - varname>=0.13.0
   - xarray>=2024.02.0
-  - pip:
-    - pyqt6>=6.2.2
+  - pyside6>=6.7.0
diff --git a/pyproject.toml b/pyproject.toml
index 5c93ebab..131d89b6 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -28,7 +28,6 @@ dependencies = [
     "lmfit>=1.2.0,!=1.3.0",
     "matplotlib>=3.8.0",
     "numba>=0.59.0",
-    "numbagg>=0.8.1",
     "numpy>=1.26.0",
     "pyperclip>=1.8.2",
     "pyqtgraph>=0.13.1",
@@ -42,24 +41,30 @@ dependencies = [
 ]
 
 [project.optional-dependencies]
+complete = ["erlab[viz,perf,misc,dev]"]
+viz = ["cmasher", "cmocean", "colorcet", "hvplot", "ipywidgets"]
+perf = ["numbagg>=0.8.1", "bottleneck>=1.3.8"]
+misc = ["iminuit>=2.25.2", "csaps>=1.1.0"]
+dev = [
+    "mypy",
+    "pre-commit",
+    "pytest-xdist",
+    "pytest",
+    "python-semantic-release",
+    "ruff",
+]
 docs = [
     "sphinx",
     "sphinx-autodoc-typehints",
     "sphinx-copybutton",
     "sphinxcontrib-bibtex",
+    "sphinx-qt-documentation",
     "pybtex",
     "nbsphinx",
     "furo",
     "sphinx-design",
 ]
-dev = [
-    "pytest",
-    "pytest-xdist",
-    "python-semantic-release",
-    "ruff",
-    "pre-commit",
-]
-viz = ["cmasher", "cmocean", "colorcet", "hvplot", "ipywidgets"]
+
 
 [project.urls]
 Documentation = "https://erlabpy.readthedocs.io"
@@ -132,10 +137,14 @@ indent-width = 4
 
 [tool.ruff.lint]
 select = [
+    "F",
     "E4",
     "E7",
     "E9",
-    "F",
+    "W",
+    "I",
+    "D",
+    "UP",
     "B",
     "C4",
     "FA",
@@ -146,21 +155,28 @@ select = [
     "RSE",
     "TID",
     "TCH",
+    "INT",
+    # "ARG",
     # "PTH",
     # "FIX",
     "TRY",
     "FLY",
     "NPY",
     "PERF",
+    # "FURB",
     "RUF",
 ]
 ignore = [
     "ICN001", # Import conventions
     "TRY003", # Long exception messages
+    "D100",   # Missing docstring in public module
+    "D101",   # Missing docstring in public class
+    "D102",   # Missing docstring in public method
+    "D103",   # Missing docstring in public function
+    "D104",   # Missing docstring in public package
+    "D105",   # Missing docstring in magic method
 ]
-extend-select = [
-    "UP", # pyupgrade
-]
+extend-select = []
 allowed-confusables = [
     "×",
     "−",
@@ -179,6 +195,9 @@ allowed-confusables = [
     "α",
 ]
 
+[tool.ruff.lint.pydocstyle]
+convention = "numpy"
+
 [tool.ruff.format]
 quote-style = "double"
 indent-style = "space"
@@ -201,12 +220,15 @@ warn_unused_configs = true
 warn_redundant_casts = true
 warn_unused_ignores = true
 allow_redefinition = true
+check_untyped_defs = false
 exclude = [
     "^docs/",
     "^tests/",
     "_deprecated/",
     "interactive/fermiedge.py",
-    "io/",
+    "interactive/bzplot.py",
+    "interactive/curvefittingtool.py",
+    "io/plugins/",
 ]
 
 [[tool.mypy.overrides]]
diff --git a/src/erlab/accessors/__init__.py b/src/erlab/accessors/__init__.py
index f8912bcf..8fb8cb89 100644
--- a/src/erlab/accessors/__init__.py
+++ b/src/erlab/accessors/__init__.py
@@ -15,15 +15,15 @@
    kspace
    fit
 
-"""
+"""  # noqa: D205
 
 __all__ = [
+    "ImageToolAccessor",
     "ModelFitDataArrayAccessor",
     "ModelFitDatasetAccessor",
-    "ParallelFitDataArrayAccessor",
     "MomentumAccessor",
     "OffsetView",
-    "ImageToolAccessor",
+    "ParallelFitDataArrayAccessor",
     "PlotAccessor",
     "SelectionAccessor",
 ]
diff --git a/src/erlab/accessors/fit.py b/src/erlab/accessors/fit.py
index 094aefcd..6cd8df01 100644
--- a/src/erlab/accessors/fit.py
+++ b/src/erlab/accessors/fit.py
@@ -10,7 +10,7 @@
 import itertools
 import warnings
 from collections.abc import Hashable, Iterable, Mapping, Sequence
-from typing import TYPE_CHECKING, Any, Literal
+from typing import TYPE_CHECKING, Any, Literal, cast
 
 import joblib
 import lmfit
@@ -45,11 +45,10 @@ def _broadcast_dict_values(d: dict[str, Any]) -> dict[str, xr.DataArray]:
         else:
             to_broadcast[k] = xr.DataArray(v)
 
-    for k, v in zip(
-        to_broadcast.keys(), xr.broadcast(*to_broadcast.values()), strict=True
-    ):
-        d[k] = v
-    return d
+    d = dict(
+        zip(to_broadcast.keys(), xr.broadcast(*to_broadcast.values()), strict=True)
+    )
+    return cast(dict[str, xr.DataArray], d)
 
 
 def _concat_along_keys(d: dict[str, xr.DataArray], dim_name: str) -> xr.DataArray:
@@ -375,7 +374,9 @@ def _wrapper(Y, *args, **kwargs):
                 x = x[:, mask]
                 y = y[mask]
                 if not len(y):
-                    modres = lmfit.model.ModelResult(model, model.make_params(), data=y)
+                    modres: lmfit.model.ModelResult = lmfit.model.ModelResult(
+                        model, model.make_params(), data=y
+                    )
                     modres.success = False
                     return popt, perr, pcov, stats, data, best, modres
             else:
@@ -417,7 +418,7 @@ def _wrapper(Y, *args, **kwargs):
                         )
                         initial_params = model.make_params().update(initial_params)
             try:
-                modres: lmfit.model.ModelResult = model.fit(
+                modres = model.fit(
                     y, **indep_var_kwargs, params=initial_params, **kwargs
                 )
             except ValueError:
@@ -428,16 +429,20 @@ def _wrapper(Y, *args, **kwargs):
                 return popt, perr, pcov, stats, data, best, modres
             else:
                 if modres.success:
-                    popt, perr = [], []
+                    popt_list, perr_list = [], []
                     for name in param_names:
                         p = modres.params[name]
-                        popt.append(p.value if p.value is not None else np.nan)
-                        perr.append(p.stderr if p.stderr is not None else np.nan)
+                        popt_list.append(p.value if p.value is not None else np.nan)
+                        perr_list.append(p.stderr if p.stderr is not None else np.nan)
 
-                    popt, perr = np.array(popt), np.array(perr)
+                    popt, perr = np.array(popt_list), np.array(perr_list)
 
-                    stats = [getattr(modres, s) for s in stat_names]
-                    stats = np.array([s if s is not None else np.nan for s in stats])
+                    stats = np.array(
+                        [
+                            s if s is not None else np.nan
+                            for s in [getattr(modres, s) for s in stat_names]
+                        ]
+                    )
 
                     if modres.covar is not None:
                         var_names = modres.var_names
@@ -550,19 +555,21 @@ def _output_wrapper(name, da, out=None) -> dict:
             if parallel_obj.return_generator:
                 out_dicts = tqdm.auto.tqdm(  # type: ignore[call-overload]
                     parallel_obj(
-                        joblib.delayed(_output_wrapper)(name, da)
-                        for name, da in self._obj.data_vars.items()
+                        itertools.starmap(
+                            joblib.delayed(_output_wrapper), self._obj.data_vars.items()
+                        )
                     ),
                     **tqdm_kw,
                 )
             else:
                 with joblib_progress(**tqdm_kw) as _:
                     out_dicts = parallel_obj(
-                        joblib.delayed(_output_wrapper)(name, da)
-                        for name, da in self._obj.data_vars.items()
+                        itertools.starmap(
+                            joblib.delayed(_output_wrapper), self._obj.data_vars.items()
+                        )
                     )
             result = type(self._obj)(
-                dict(itertools.chain.from_iterable(d.items() for d in out_dicts))  # type: ignore[call-overload]
+                dict(itertools.chain.from_iterable(d.items() for d in out_dicts))
             )
             del out_dicts
 
@@ -602,11 +609,7 @@ def __call__(self, *args, **kwargs) -> xr.Dataset:
 
 @xr.register_dataarray_accessor("parallel_fit")
 class ParallelFitDataArrayAccessor(ERLabDataArrayAccessor):
-    """
-    `xarray.DataArray.parallel_fit` accessor for fitting lmfit models in parallel along
-    a single dimension.
-
-    """
+    """`xarray.DataArray.parallel_fit` accessor for fitting lmfit models in parallel."""
 
     _VAR_KEYS: tuple[str, ...] = (
         "modelfit_results",
diff --git a/src/erlab/accessors/kspace.py b/src/erlab/accessors/kspace.py
index ab950519..e02b351f 100644
--- a/src/erlab/accessors/kspace.py
+++ b/src/erlab/accessors/kspace.py
@@ -16,13 +16,12 @@
 from erlab.analysis.interpolate import interpn
 from erlab.analysis.kspace import AxesConfiguration, get_kconv_func, kz_func
 from erlab.constants import rel_kconv, rel_kzconv
-from erlab.interactive.kspace import KspaceTool, ktool
 
 
 def only_angles(method=None):
-    """
-    A decorator that ensures the data is in angle space before executing the decorated
-    method.
+    """Decorate methods that require data to be in angle space.
+
+    Ensures the data is in angle space before executing the decorated method.
 
     If the data is not in angle space (i.e., if "kx" or "ky" dimensions are present), a
     `ValueError` is raised.
@@ -45,9 +44,9 @@ def _impl(self, *args, **kwargs):
 
 
 def only_momentum(method=None):
-    """
-    A decorator that ensures the data is in momentum space before executing the
-    decorated method.
+    """Decorate methods that require data to be in momentum space.
+
+    Ensure the data is in momentum space before executing the decorated method.
 
     If the data is not in momentum space (i.e., if "kx" nor "ky" dimensions are
     present), a `ValueError` is raised.
@@ -155,7 +154,7 @@ def update(
         other: dict[str, float] | Iterable[tuple[str, float]] | None = None,
         **kwargs,
     ) -> "OffsetView":
-        """Updates the offset view with the provided key-value pairs."""
+        """Update the offset view with the provided key-value pairs."""
         if other is not None:
             for k, v in other.items() if isinstance(other, dict) else other:
                 self[str(k)] = v
@@ -164,7 +163,7 @@ def update(
         return self
 
     def items(self) -> ItemsView[str, float]:
-        """Returns a view of the offset view as a list of (key, value) pairs."""
+        """Return a view of the offset view as a list of (key, value) pairs."""
         return dict(self).items()
 
     def reset(self) -> "OffsetView":
@@ -187,7 +186,7 @@ class MomentumAccessor(ERLabDataArrayAccessor):
 
     @property
     def configuration(self) -> AxesConfiguration:
-        """Returns the experimental configuration.
+        """Return the experimental configuration.
 
         For a properly implemented data loader, the configuration attribute must be set
         on data import. See :class:`erlab.analysis.kspace.AxesConfiguration` for
@@ -223,7 +222,6 @@ def inner_potential(self) -> float:
         >>> data.kspace.inner_potential
         13.0
         """
-
         if "inner_potential" in self._obj.attrs:
             return float(self._obj.attrs["inner_potential"])
         else:
@@ -256,7 +254,6 @@ def work_function(self) -> float:
         >>> data.kspace.work_function
         4.5
         """
-
         if "sample_workfunction" in self._obj.attrs:
             return float(self._obj.attrs["sample_workfunction"])
         else:
@@ -280,7 +277,6 @@ def angle_resolution(self) -> float:
         This property is used in `best_kp_resolution` upon estimating momentum step
         sizes through `estimate_resolution`.
         """
-
         try:
             return float(self._obj.attrs["angle_resolution"])
         except KeyError:
@@ -295,7 +291,7 @@ def angle_resolution(self, value: float):
 
     @property
     def slit_axis(self) -> Literal["kx", "ky"]:
-        """Returns the momentum axis parallel to the slit.
+        """Return the momentum axis parallel to the slit.
 
         Returns
         -------
@@ -310,7 +306,7 @@ def slit_axis(self) -> Literal["kx", "ky"]:
 
     @property
     def other_axis(self) -> Literal["kx", "ky"]:
-        """Returns the momentum axis perpendicular to the slit.
+        """Return the momentum axis perpendicular to the slit.
 
         Returns
         -------
@@ -326,7 +322,7 @@ def other_axis(self) -> Literal["kx", "ky"]:
     @property
     @only_angles
     def momentum_axes(self) -> tuple[Literal["kx", "ky", "kz"], ...]:
-        """Returns the momentum axes of the data after conversion.
+        """Return the momentum axes of the data after conversion.
 
         Returns
         -------
@@ -388,13 +384,13 @@ def kinetic_energy(self) -> xr.DataArray:
 
     @property
     def has_eV(self) -> bool:
-        """Returns `True` if object has an energy axis."""
+        """Return `True` if object has an energy axis."""
         return "eV" in self._obj.dims
 
     @property
     @only_angles
     def has_hv(self) -> bool:
-        """Returns `True` for photon energy dependent data."""
+        """Return `True` for photon energy dependent data."""
         return self._obj["hv"].size > 1
 
     @property
@@ -446,7 +442,6 @@ def offsets(self) -> OffsetView:
 
         Examples
         --------
-
         - View all offsets
 
           >>> data.kspace.offsets
@@ -500,9 +495,9 @@ def offsets(self, offset_dict: dict[str, float]):
     @property
     @only_angles
     def best_kp_resolution(self) -> float:
-        r"""
-        Estimates the minimum in-plane momentum resolution based on the kinetic energy
-        and angular resolution:
+        r"""Estimate the minimum in-plane momentum resolution.
+
+        The resolution is estimated with the kinetic energy and angular resolution:
 
         .. math:: \Delta k_{\parallel} \sim \sqrt{2 m_e E_k/\hbar^2} \cos(\alpha) \Delta\alpha
 
@@ -519,9 +514,10 @@ def best_kp_resolution(self) -> float:
     @property
     @only_angles
     def best_kz_resolution(self) -> float:
-        r"""
-        Estimates the minimum out-of-plane momentum resolution based on the mean free
-        path :cite:p:`seah1979imfp` and the kinetic energy.
+        r"""Estimate the minimum out-of-plane momentum resolution.
+
+        The resolution is estimated based on the mean free path :cite:p:`seah1979imfp`
+        and the kinetic energy.
 
         .. math:: \Delta k_z \sim 1/\lambda
 
@@ -541,8 +537,7 @@ def _get_transformed_coords(self) -> dict[Literal["kx", "ky", "kz"], xr.DataArra
 
     def estimate_bounds(self) -> dict[Literal["kx", "ky", "kz"], tuple[float, float]]:
         """
-        Estimates the bounds of the data in momentum space based on the available
-        parameters.
+        Estimate the bounds of the data in momentum space.
 
         Returns
         -------
@@ -552,7 +547,6 @@ def estimate_bounds(self) -> dict[Literal["kx", "ky", "kz"], tuple[float, float]
             values are tuples representing the minimum and maximum values.
 
         """
-
         return {
             k: (v.values.min(), v.values.max())
             for k, v in self._get_transformed_coords().items()
@@ -672,9 +666,9 @@ def convert_coords(self) -> xr.DataArray:
 
     @only_angles
     def _get_coord_for_conversion(self, name: Hashable) -> xr.DataArray:
-        """
-        Get the coordinte array for given dimension name. This just ensures that the
-        energy coordinates are given as binding energy.
+        """Get the coordinte array for given dimension name.
+
+        This just ensures that the energy coordinates are given as binding energy.
         """
         if name == "eV":
             return self.binding_energy
@@ -683,9 +677,7 @@ def _get_coord_for_conversion(self, name: Hashable) -> xr.DataArray:
 
     @only_angles
     def _data_ensure_binding(self) -> xr.DataArray:
-        """
-        Returns the data while ensuring that the energy axis is in binding energy.
-        """
+        """Return the data while ensuring that the energy axis is in binding energy."""
         return self._obj.assign_coords(eV=self.binding_energy)
 
     @only_angles
@@ -757,7 +749,6 @@ def convert(
             converted_data = data.kspace.convert(bounds, resolution)
 
         """
-
         if bounds is None:
             bounds = {}
 
@@ -870,15 +861,17 @@ def _wrap_interpn(arr, *args):
 
         return out
 
-    def interactive(self, **kwargs) -> KspaceTool:
+    def interactive(self, **kwargs):
         """Open the interactive momentum space conversion tool."""
+        from erlab.interactive.kspace import ktool
+
         if self._obj.ndim < 3:
             raise ValueError("Interactive tool requires three-dimensional data.")
         return ktool(self._obj, **kwargs)
 
     @only_momentum
     def hv_to_kz(self, hv: float) -> xr.DataArray:
-        """Returns :math:`k_z` for a given photon energy.
+        """Return :math:`k_z` for a given photon energy.
 
         Useful when creating overlays on :math:`hν`-dependent data.
 
diff --git a/src/erlab/accessors/utils.py b/src/erlab/accessors/utils.py
index 607f211b..3bfdfab1 100644
--- a/src/erlab/accessors/utils.py
+++ b/src/erlab/accessors/utils.py
@@ -12,7 +12,6 @@
 import xarray as xr
 
 import erlab.plotting.erplot as eplt
-from erlab.interactive.imagetool import ImageTool, itool
 
 T = TypeVar("T")
 
@@ -63,10 +62,11 @@ class PlotAccessor(ERLabDataArrayAccessor):
     """`xarray.DataArray.qplot` accessor for plotting data."""
 
     def __call__(self, *args, **kwargs):
-        """
-        Plot the data. If a 2D data array is provided, it is plotted using
-        :func:`plot_array <erlab.plotting.general.plot_array>`. Otherwise, it is
-        equivalent to calling :meth:`xarray.DataArray.plot`.
+        """Plot the data.
+
+        If a 2D data array is provided, it is plotted using :func:`plot_array
+        <erlab.plotting.general.plot_array>`. Otherwise, it is equivalent to calling
+        :meth:`xarray.DataArray.plot`.
 
         Parameters
         ----------
@@ -93,7 +93,9 @@ def __call__(self, *args, **kwargs):
 class ImageToolAccessor(ERLabDataArrayAccessor):
     """`xarray.DataArray.qshow` accessor for interactive visualization."""
 
-    def __call__(self, *args, **kwargs) -> ImageTool | list[ImageTool] | None:
+    def __call__(self, *args, **kwargs):
+        from erlab.interactive.imagetool import itool
+
         if len(self._obj.dims) >= 2:
             return itool(self._obj, *args, **kwargs)
         else:
@@ -102,10 +104,7 @@ def __call__(self, *args, **kwargs) -> ImageTool | list[ImageTool] | None:
 
 @xr.register_dataarray_accessor("qsel")
 class SelectionAccessor(ERLabDataArrayAccessor):
-    """
-    `xarray.DataArray.qsel` accessor for conveniently selecting and averaging
-    data.
-    """
+    """`xarray.DataArray.qsel` accessor for convenient selection and averaging."""
 
     def __call__(
         self,
@@ -154,7 +153,6 @@ def __call__(
         ValueError
             If a specified dimension is not present in the data.
         """
-
         indexers = either_dict_or_kwargs(indexers, indexers_kwargs, "qsel")
 
         # Bin widths for each dimension, zero if width not specified
diff --git a/src/erlab/analysis/__init__.py b/src/erlab/analysis/__init__.py
index 3740d7f6..8a0de0cd 100644
--- a/src/erlab/analysis/__init__.py
+++ b/src/erlab/analysis/__init__.py
@@ -30,9 +30,12 @@
     "rotateinplane",
     "rotatestackinplane",
     "shift",
+    "slice_along_path",
 ]
 
+from erlab.analysis import fit, gold, image, interpolate, mask  # noqa: F401
 from erlab.analysis.gold import correct_with_edge
+from erlab.analysis.interpolate import slice_along_path
 from erlab.analysis.mask import (
     mask_with_hex_bz,
     mask_with_polygon,
diff --git a/src/erlab/analysis/correlation.py b/src/erlab/analysis/correlation.py
index 30f46be9..7b7ba4b2 100644
--- a/src/erlab/analysis/correlation.py
+++ b/src/erlab/analysis/correlation.py
@@ -1,6 +1,6 @@
 """Macros for correlation analysis."""
 
-__all__ = ["acf2", "acf2stack", "match_dims", "xcorr1d"]
+__all__ = ["acf2", "acf2stack", "xcorr1d"]
 
 
 import itertools
@@ -48,9 +48,7 @@ def nanacf(arr, *args, **kwargs):
 
 
 def acf2(arr, mode: str = "full", method: str = "fft"):
-    """
-    Calculate the autocorrelation function (ACF) of a two-dimensional array including
-    nan values.
+    """Calculate the autocorrelation function (ACF) of a 2D array including nans.
 
     Parameters
     ----------
@@ -152,18 +150,9 @@ def acf2stack(arr, stack_dims=("eV",), mode: str = "full", method: str = "fft"):
     return out
 
 
-def match_dims(da1: xr.DataArray, da2: xr.DataArray):
-    """
-    Returns the second array interpolated with the coordinates of the first array,
-    making them the same size.
-
-    """
-    return da2.interp({dim: da1[dim] for dim in da2.dims})
-
-
 def xcorr1d(in1: xr.DataArray, in2: xr.DataArray, method="direct"):
-    """Performs 1-dimensional correlation analysis on `xarray.DataArray` s."""
-    in2 = match_dims(in1, in2)
+    """Perform 1-dimensional correlation analysis on `xarray.DataArray` s."""
+    in2 = in2.interp_like(in1)
     out = in1.copy(deep=True)
     xind = scipy.signal.correlation_lags(in1.values.size, in2.values.size, mode="same")
     xzero = np.flatnonzero(xind == 0)[0]
diff --git a/src/erlab/analysis/fit/functions/__init__.py b/src/erlab/analysis/fit/functions/__init__.py
index 5014be1a..a93e50a4 100644
--- a/src/erlab/analysis/fit/functions/__init__.py
+++ b/src/erlab/analysis/fit/functions/__init__.py
@@ -25,10 +25,10 @@
     "fermi_dirac",
     "fermi_dirac_linbkg",
     "fermi_dirac_linbkg_broad",
-    "gaussian_wh",
     "gaussian",
-    "lorentzian_wh",
+    "gaussian_wh",
     "lorentzian",
+    "lorentzian_wh",
     "step_broad",
     "step_linbkg_broad",
 ]
@@ -47,10 +47,10 @@
     fermi_dirac,
     fermi_dirac_linbkg,
     fermi_dirac_linbkg_broad,
-    gaussian_wh,
     gaussian,
-    lorentzian_wh,
+    gaussian_wh,
     lorentzian,
+    lorentzian_wh,
     step_broad,
     step_linbkg_broad,
 )
diff --git a/src/erlab/analysis/fit/functions/dynamic.py b/src/erlab/analysis/fit/functions/dynamic.py
index 909b0fe8..f50f12ad 100644
--- a/src/erlab/analysis/fit/functions/dynamic.py
+++ b/src/erlab/analysis/fit/functions/dynamic.py
@@ -14,7 +14,7 @@
 import functools
 import inspect
 from collections.abc import Callable, Sequence
-from typing import Any, TypedDict, no_type_check, ClassVar
+from typing import Any, ClassVar, TypedDict, no_type_check
 
 import numpy as np
 import numpy.typing as npt
@@ -125,16 +125,16 @@ def __init__(self, degree: int = 1) -> None:
     def argnames(self) -> list[str]:
         return ["x"] + [f"c{i}" for i in range(self.degree + 1)]
 
-    def __call__(self, x, *coeffs, **params):
+    def __call__(self, x, *coeffs: float, **params):
         if len(coeffs) != self.degree + 1:
-            coeffs = [params[f"c{d}"] for d in range(self.degree + 1)]
+            coeffs = tuple(params[f"c{d}"] for d in range(self.degree + 1))
         if isinstance(x, np.ndarray):
             return np.polynomial.polynomial.polyval(x, coeffs)
         else:
-            coeffs = xr.DataArray(
+            coeffs_xr = xr.DataArray(
                 np.asarray(coeffs), coords={"degree": np.arange(self.degree + 1)}
             )
-            return xr.polyval(x, coeffs)
+            return xr.polyval(x, coeffs_xr)
 
 
 class MultiPeakFunction(DynamicFunction):
@@ -231,14 +231,14 @@ def argnames(self) -> list[str]:
     @property
     def kwargs(self):
         kws = [
-            ("lin_bkg", 0),
-            ("const_bkg", 0),
+            ("lin_bkg", 0.0),
+            ("const_bkg", 0.0),
         ]
         if self.fd:
             kws += [
                 ("efermi", 0.0),  # fermi level
-                ("temp", 30),  # temperature
-                ("offset", 0),
+                ("temp", 30.0),  # temperature
+                ("offset", 0.0),
             ]
         if self.convolve:
             kws += [("resolution", 0.02)]
diff --git a/src/erlab/analysis/fit/functions/general.py b/src/erlab/analysis/fit/functions/general.py
index 7a6b2335..6fdc0b1e 100644
--- a/src/erlab/analysis/fit/functions/general.py
+++ b/src/erlab/analysis/fit/functions/general.py
@@ -12,15 +12,16 @@
     "fermi_dirac",
     "fermi_dirac_linbkg",
     "fermi_dirac_linbkg_broad",
-    "gaussian_wh",
     "gaussian",
-    "lorentzian_wh",
+    "gaussian_wh",
     "lorentzian",
+    "lorentzian_wh",
     "step_broad",
     "step_linbkg_broad",
 ]
 
 from collections.abc import Callable
+
 import numba
 import numpy as np
 import numpy.typing as npt
@@ -83,7 +84,6 @@ def _gen_kernel(
         The gaussian kernel defined on `extended`.
 
     """
-
     delta_x = x[1] - x[0]
 
     sigma = abs(resolution) / np.sqrt(8 * np.log(2))  # resolution given in FWHM
@@ -263,10 +263,7 @@ def fermi_dirac_linbkg_broad(
     dos0: float,
     dos1: float,
 ) -> npt.NDArray[np.float64]:
-    """
-    Resolution-broadened Fermi-dirac edge with linear backgrounds above and below the
-    fermi level.
-    """
+    """Resolution-broadened Fermi edge with linear backgrounds above and below EF."""
     return do_convolve(
         x,
         fermi_dirac_linbkg,
@@ -303,10 +300,7 @@ def step_linbkg_broad(
     dos0: float,
     dos1: float,
 ):
-    """
-    A linear density of states multiplied with a resolution broadened step function with
-    a linear background.
-    """
+    """Resolution broadened step function with linear backgrounds."""
     return (back0 + back1 * x) + (dos0 - back0 + (dos1 - back1) * x) * (
         step_broad(x, center, sigma, 1.0)
     )
@@ -325,7 +319,7 @@ def bcs_gap(
 
     Parameters
     ----------
-    x : array_like
+    x : array-like
         The temperature values in kelvins at which to calculate the BCS gap.
     a
         Proportionality constant. Default is 1.76.
@@ -345,7 +339,7 @@ def bcs_gap(
 
 
 def dynes(x, n0=1.0, gamma=0.003, delta=0.01):
-    r"""The Dynes formula for superconducting density of states.
+    r"""Dynes formula for superconducting density of states.
 
     The formula is given by :cite:p:`dynes1978dynes`:
 
@@ -360,7 +354,7 @@ def dynes(x, n0=1.0, gamma=0.003, delta=0.01):
 
     Parameters
     ----------
-    x : array_like
+    x : array-like
         The input array of energy in eV.
     n0
         :math:`N_0`, by default 1.0.
diff --git a/src/erlab/analysis/fit/minuit.py b/src/erlab/analysis/fit/minuit.py
index 3cdae183..0df8199f 100644
--- a/src/erlab/analysis/fit/minuit.py
+++ b/src/erlab/analysis/fit/minuit.py
@@ -71,8 +71,7 @@ def __exit__(self, *args: object) -> None:
 
 
 class Minuit(iminuit.Minuit):
-    """A subclass of `iminuit.Minuit` that provides additional functionality for fitting
-    `lmfit` library models.
+    """`iminuit.Minuit` with additional functionality.
 
     This class extends the functionality of the `iminuit.Minuit` class by providing a
     convenient method `from_lmfit` to initialize the `Minuit` object from an
@@ -82,7 +81,6 @@ class Minuit(iminuit.Minuit):
 
     Examples
     --------
-
     >>> import lmfit.models
     >>> import numpy as np
     >>> from erlab.analysis.fit.minuit import Minuit
diff --git a/src/erlab/analysis/fit/models.py b/src/erlab/analysis/fit/models.py
index 85e6d654..0c5900c6 100644
--- a/src/erlab/analysis/fit/models.py
+++ b/src/erlab/analysis/fit/models.py
@@ -3,8 +3,8 @@
 __all__ = [
     "BCSGapModel",
     "DynesModel",
-    "FermiEdgeModel",
     "FermiEdge2dModel",
+    "FermiEdgeModel",
     "MultiPeakModel",
     "PolynomialModel",
     "StepEdgeModel",
@@ -112,10 +112,10 @@ def fit_edges_linear(x, data, len_fit) -> tuple[float, float, float, float]:
 
 
 class FermiEdgeModel(lmfit.Model):
-    """
-    Fermi-dirac function with linear background above and below the fermi level,
-    convolved with a gaussian kernel.
+    """Model for fitting a Fermi edge with a linear background.
 
+    The model function is a Fermi-dirac function with linear background above and below
+    the fermi level, convolved with a gaussian kernel.
     """
 
     __doc__ = __doc__ + lmfit.models.COMMON_INIT_DOC
diff --git a/src/erlab/analysis/fit/spline.py b/src/erlab/analysis/fit/spline.py
index 47e449e7..741032d4 100644
--- a/src/erlab/analysis/fit/spline.py
+++ b/src/erlab/analysis/fit/spline.py
@@ -26,7 +26,6 @@ def xcsaps(arr: xr.DataArray, **kwargs) -> tuple[xr.DataArray, csaps.ISmoothingS
     spl : csaps.ISmoothingSpline
         The spline object.
     """
-
     kwargs.setdefault("normalizedsmooth", True)
     coords = [arr[d].values for d in arr.dims]
     spl = csaps.csaps(coords, arr.values, **kwargs)
diff --git a/src/erlab/analysis/gold.py b/src/erlab/analysis/gold.py
index fb2a5bb3..f987eca3 100644
--- a/src/erlab/analysis/gold.py
+++ b/src/erlab/analysis/gold.py
@@ -212,7 +212,6 @@ def edge(
         `True`.
 
     """
-
     if fast:
         params = lmfit.create_params()
         model_cls: lmfit.Model = StepEdgeModel
diff --git a/src/erlab/analysis/image.py b/src/erlab/analysis/image.py
index ae8dbe79..af096c61 100644
--- a/src/erlab/analysis/image.py
+++ b/src/erlab/analysis/image.py
@@ -7,7 +7,7 @@
 unlike the scipy default of 'reflect'.
 """
 
-from collections.abc import Collection, Mapping, Sequence, Sized, Hashable
+from collections.abc import Collection, Hashable, Mapping, Sequence, Sized
 
 import numpy as np
 import numpy.typing as npt
@@ -205,7 +205,6 @@ def gaussian_laplace(
     :func:`scipy.ndimage.gaussian_laplace` : The underlying function used to apply the
         filter.
     """
-
     if isinstance(sigma, Mapping):
         sigma_dict = dict(sigma)
     elif np.isscalar(sigma):
@@ -279,7 +278,6 @@ def gradient_magnitude(
     The gradient magnitude is then calculated as the square root of the sum of the
     squared differences.
     """
-
     dxy = np.sqrt(dx**2 + dy**2)
     dist = np.array([[dxy, dy, dxy], [dx, 0.0, dx], [dxy, dy, dxy]]).flatten()
 
@@ -378,7 +376,6 @@ def minimum_gradient(
     - The input array is assumed to be regularly spaced.
     - Any zero gradient values are replaced with NaN.
     """
-
     if darr.ndim != 2:
         raise ValueError("DataArray must be 2D")
 
@@ -505,7 +502,6 @@ def curvature(darr: xr.DataArray, a0: float = 1.0, factor: float = 1.0) -> xr.Da
     ----
     The input array is assumed to be regularly spaced.
     """
-
     if darr.ndim != 2:
         raise ValueError("DataArray must be 2D")
 
diff --git a/src/erlab/analysis/interpolate.py b/src/erlab/analysis/interpolate.py
index 991aa43f..b351bf7d 100644
--- a/src/erlab/analysis/interpolate.py
+++ b/src/erlab/analysis/interpolate.py
@@ -1,163 +1,21 @@
-"""Fast linear interpolation."""
+"""Utilities for interpolation."""
 
-__all__ = ["FastInterpolator", "interpn"]
+__all__ = ["FastInterpolator", "interpn", "slice_along_path"]
 
+import itertools
 import math
 import warnings
-from collections.abc import Sequence
+from collections.abc import Hashable, Mapping, Sequence
+from typing import cast
 
 import numba
 import numpy as np
+import numpy.typing as npt
 import scipy.interpolate
 import xarray as xr
 import xarray.core.missing
 
-
-@numba.njit(nogil=True, inline="always")
-def _do_interp1(x, v0, v1):
-    return v0 * (1 - x) + v1 * x
-
-
-@numba.njit(nogil=True, inline="always")
-def _do_interp2(x, y, v0, v1, v2, v3):
-    return _do_interp1(y, _do_interp1(x, v0, v1), _do_interp1(x, v2, v3))
-
-
-@numba.njit(nogil=True, inline="always")
-def _do_interp3(x, y, z, v0, v1, v2, v3, v4, v5, v6, v7):
-    return _do_interp1(
-        z, _do_interp2(x, y, v0, v1, v2, v3), _do_interp2(x, y, v4, v5, v6, v7)
-    )
-
-
-# @numba.njit
-# def _do_interp2(x, y, v0, v1, v2, v3):
-#     return _do_interp1(y, v0 * (1 - x) + v1 * x, v2 * (1 - x) + v3 * x)
-
-
-# @numba.njit
-# def _do_interp3(x, y, z, v0, v1, v2, v3, v4, v5, v6, v7):
-#     return (
-#         _do_interp2(x, y, v0, v1, v2, v3) * (1 - z)
-#         + _do_interp2(x, y, v4, v5, v6, v7) * z
-#     )
-
-
-@numba.njit(nogil=True, inline="always")
-def _calc_interp1(values, v0):
-    i0 = math.floor(v0)
-    n0 = values.size
-    j0 = min(i0 + 1, n0 - 1)
-    return _do_interp1(v0 - i0, values[i0], values[j0])
-
-
-@numba.njit(nogil=True, inline="always")
-def _calc_interp2(values, v0, v1):
-    i0, i1 = math.floor(v0), math.floor(v1)
-    n0, n1 = values.shape
-    j0, j1 = min(i0 + 1, n0 - 1), min(i1 + 1, n1 - 1)
-    return _do_interp2(
-        v0 - i0,
-        v1 - i1,
-        values[i0, i1],
-        values[j0, i1],
-        values[i0, j1],
-        values[j0, j1],
-    )
-
-
-@numba.njit(nogil=True, inline="always")
-def _calc_interp3(values, v0, v1, v2):
-    i0, i1, i2 = math.floor(v0), math.floor(v1), math.floor(v2)
-    n0, n1, n2 = values.shape
-    j0, j1, j2 = min(i0 + 1, n0 - 1), min(i1 + 1, n1 - 1), min(i2 + 1, n2 - 1)
-    return _do_interp3(
-        v0 - i0,
-        v1 - i1,
-        v2 - i2,
-        values[i0, i1, i2],
-        values[j0, i1, i2],
-        values[i0, j1, i2],
-        values[j0, j1, i2],
-        values[i0, i1, j2],
-        values[j0, i1, j2],
-        values[i0, j1, j2],
-        values[j0, j1, j2],
-    )
-
-
-@numba.njit(nogil=True, inline="always")
-def _val2ind(val, coord):
-    if val > coord[-1] or val < coord[0]:
-        return np.nan
-    else:
-        return np.divide(val - coord[0], coord[1] - coord[0])
-
-
-@numba.njit(nogil=True, parallel=True)
-def _interp1(x, values, xc, fill_value=np.nan):
-    n = len(xc)
-
-    arr_new = np.empty(n, values.dtype)
-
-    for m in numba.prange(n):
-        v0 = _val2ind(xc[m], x)
-
-        if np.isnan(v0):
-            arr_new[m] = fill_value
-        else:
-            arr_new[m] = _calc_interp1(values, v0)
-    return arr_new
-
-
-@numba.njit(nogil=True, parallel=True)
-def _interp2(x, y, values, xc, yc, fill_value=np.nan):
-    n = len(xc)
-
-    arr_new = np.empty(n, values.dtype)
-
-    for m in numba.prange(n):
-        v0, v1 = _val2ind(xc[m], x), _val2ind(yc[m], y)
-
-        if np.isnan(v0) or np.isnan(v1):
-            arr_new[m] = fill_value
-        else:
-            arr_new[m] = _calc_interp2(values, v0, v1)
-    return arr_new
-
-
-@numba.njit(nogil=True, parallel=True)
-def _interp3(x, y, z, values, xc, yc, zc, fill_value=np.nan):
-    n = len(xc)
-
-    arr_new = np.empty(n, values.dtype)
-
-    for m in numba.prange(n):
-        v0, v1, v2 = _val2ind(xc[m], x), _val2ind(yc[m], y), _val2ind(zc[m], z)
-
-        if np.isnan(v0) or np.isnan(v1) or np.isnan(v2):
-            arr_new[m] = fill_value
-        else:
-            arr_new[m] = _calc_interp3(values, v0, v1, v2)
-    return arr_new
-
-
-def _get_interp_func(ndim):
-    if ndim == 3:
-        return _interp3
-    elif ndim == 2:
-        return _interp2
-    elif ndim == 1:
-        return _interp1
-    else:
-        raise ValueError("Fast interpolation only supported for 2D or 3D")
-
-
-def _check_even(arr):
-    dif = np.diff(arr)
-    if dif.size == 0:
-        return False
-    return np.allclose(dif, dif[0])
+from erlab.accessors.utils import either_dict_or_kwargs
 
 
 class FastInterpolator(scipy.interpolate.RegularGridInterpolator):
@@ -333,6 +191,267 @@ def interpn(
     return interp(xi)
 
 
+@numba.njit(nogil=True, inline="always")
+def _do_interp1(x, v0, v1):
+    return v0 * (1 - x) + v1 * x
+
+
+@numba.njit(nogil=True, inline="always")
+def _do_interp2(x, y, v0, v1, v2, v3):
+    return _do_interp1(y, _do_interp1(x, v0, v1), _do_interp1(x, v2, v3))
+
+
+@numba.njit(nogil=True, inline="always")
+def _do_interp3(x, y, z, v0, v1, v2, v3, v4, v5, v6, v7):
+    return _do_interp1(
+        z, _do_interp2(x, y, v0, v1, v2, v3), _do_interp2(x, y, v4, v5, v6, v7)
+    )
+
+
+@numba.njit(nogil=True, inline="always")
+def _calc_interp1(values, v0):
+    i0 = math.floor(v0)
+    n0 = values.size
+    j0 = min(i0 + 1, n0 - 1)
+    return _do_interp1(v0 - i0, values[i0], values[j0])
+
+
+@numba.njit(nogil=True, inline="always")
+def _calc_interp2(values, v0, v1):
+    i0, i1 = math.floor(v0), math.floor(v1)
+    n0, n1 = values.shape
+    j0, j1 = min(i0 + 1, n0 - 1), min(i1 + 1, n1 - 1)
+    return _do_interp2(
+        v0 - i0,
+        v1 - i1,
+        values[i0, i1],
+        values[j0, i1],
+        values[i0, j1],
+        values[j0, j1],
+    )
+
+
+@numba.njit(nogil=True, inline="always")
+def _calc_interp3(values, v0, v1, v2):
+    i0, i1, i2 = math.floor(v0), math.floor(v1), math.floor(v2)
+    n0, n1, n2 = values.shape
+    j0, j1, j2 = min(i0 + 1, n0 - 1), min(i1 + 1, n1 - 1), min(i2 + 1, n2 - 1)
+    return _do_interp3(
+        v0 - i0,
+        v1 - i1,
+        v2 - i2,
+        values[i0, i1, i2],
+        values[j0, i1, i2],
+        values[i0, j1, i2],
+        values[j0, j1, i2],
+        values[i0, i1, j2],
+        values[j0, i1, j2],
+        values[i0, j1, j2],
+        values[j0, j1, j2],
+    )
+
+
+@numba.njit(nogil=True, inline="always")
+def _val2ind(val, coord):
+    if val > coord[-1] or val < coord[0]:
+        return np.nan
+    else:
+        return np.divide(val - coord[0], coord[1] - coord[0])
+
+
+@numba.njit(nogil=True, parallel=True)
+def _interp1(x, values, xc, fill_value=np.nan):
+    n = len(xc)
+
+    arr_new = np.empty(n, values.dtype)
+
+    for m in numba.prange(n):
+        v0 = _val2ind(xc[m], x)
+
+        if np.isnan(v0):
+            arr_new[m] = fill_value
+        else:
+            arr_new[m] = _calc_interp1(values, v0)
+    return arr_new
+
+
+@numba.njit(nogil=True, parallel=True)
+def _interp2(x, y, values, xc, yc, fill_value=np.nan):
+    n = len(xc)
+
+    arr_new = np.empty(n, values.dtype)
+
+    for m in numba.prange(n):
+        v0, v1 = _val2ind(xc[m], x), _val2ind(yc[m], y)
+
+        if np.isnan(v0) or np.isnan(v1):
+            arr_new[m] = fill_value
+        else:
+            arr_new[m] = _calc_interp2(values, v0, v1)
+    return arr_new
+
+
+@numba.njit(nogil=True, parallel=True)
+def _interp3(x, y, z, values, xc, yc, zc, fill_value=np.nan):
+    n = len(xc)
+
+    arr_new = np.empty(n, values.dtype)
+
+    for m in numba.prange(n):
+        v0, v1, v2 = _val2ind(xc[m], x), _val2ind(yc[m], y), _val2ind(zc[m], z)
+
+        if np.isnan(v0) or np.isnan(v1) or np.isnan(v2):
+            arr_new[m] = fill_value
+        else:
+            arr_new[m] = _calc_interp3(values, v0, v1, v2)
+    return arr_new
+
+
+def _get_interp_func(ndim):
+    if ndim == 3:
+        return _interp3
+    elif ndim == 2:
+        return _interp2
+    elif ndim == 1:
+        return _interp1
+    else:
+        raise ValueError("Fast interpolation only supported for 2D or 3D")
+
+
+def _check_even(arr):
+    dif = np.diff(arr)
+    if dif.size == 0:
+        return False
+    return np.allclose(dif, dif[0])
+
+
+def slice_along_path(
+    darr: xr.DataArray,
+    vertices: Mapping[Hashable, Sequence],
+    step_size: float | None = None,
+    dim_name: str = "path",
+    interp_kwargs: dict | None = None,
+    **vertices_kwargs,
+) -> xr.DataArray:
+    """Interpolate a DataArray along a path defined by a sequence of vertices.
+
+    Parameters
+    ----------
+    darr
+        The data array to interpolate.
+    vertices
+        Dictionary specifying the vertices of the path along which to interpolate the
+        DataArray. The keys of the dictionary should correspond to the dimensions of the
+        DataArray along which to interpolate.
+    step_size
+        The step size to use for the interpolation. This determines the number of points
+        along the path at which the data array will be interpolated. If None, the step
+        size is determined automatically as the smallest step size of the coordinates
+        along the dimensions of the vertices if all coordinates are evenly spaced. If
+        there exists a dimension where the coordinates are not evenly spaced,
+        `step_size` must be specified.
+    dim_name
+        The name of the new dimension that corresponds to the distance along the
+        interpolated path. Default is "path".
+    interp_kwargs
+        Additional keyword arguments passed to `xarray.DataArray.interp`.
+    **vertices_kwargs
+        The keyword arguments form of `vertices`. One of `vertices` or `vertices_kwargs`
+        must be provided.
+
+    Returns
+    -------
+    interpolated : DataArray
+        New dataarray on the new coordinate.
+
+    Examples
+    --------
+    >>> import numpy as np
+    >>> import xarray as xr
+    >>> from erlab.analysis.interpolate import slice_along_path
+    >>> x = np.linspace(0, 10, 11)
+    >>> y = np.linspace(0, 10, 11)
+    >>> z = np.linspace(0, 10, 11)
+    >>> data = np.random.rand(11, 11, 11)
+    >>> darr = xr.DataArray(data, coords={"x": x, "y": y, "z": z}, dims=["x", "y", "z"])
+    >>> vertices = {"x": [0, 5, 10], "y": [0, 5, 10], "z": [0, 5, 10]}
+    >>> interp = slice_along_path(darr, vertices)
+
+    See Also
+    --------
+    xarray.DataArray.interp
+
+    """
+    interp_kwargs = interp_kwargs or {}
+    vertices = dict(
+        either_dict_or_kwargs(vertices, vertices_kwargs, "slice_along_path")
+    )
+
+    for dim, vert in vertices.items():
+        if dim not in darr.dims:
+            raise ValueError(f"Dimension {dim} not found in data array")
+
+        # Convert to numpy array
+        vert_arr = np.asarray(vert)
+        if vert_arr.ndim != 1:
+            raise ValueError("Each vertex must be 1-dimensional")
+        vertices[dim] = vert_arr
+
+    if not all(len(p) == len(next(iter(vertices.values()))) for p in vertices.values()):
+        raise ValueError("Vertices must have the same length along all dimensions")
+
+    if step_size is None:
+        step_size = np.inf
+        for dim in vertices:
+            dif = np.diff(darr[dim].values)
+            if np.allclose(dif, dif[0], equal_nan=True):
+                step_size = min(step_size, np.abs(dif[0]))
+
+        if not np.isfinite(cast(float, step_size)):
+            raise ValueError("Could not determine step size automatically")
+    else:
+        if not np.isfinite(step_size):
+            raise ValueError("Step size must be finite")
+        if step_size <= 0:
+            raise ValueError("Step size must be positive")
+
+    points: npt.NDArray[np.floating] = np.array(list(vertices.values())).T
+
+    # Calculate number of points for each segment
+    num_points = [
+        round(np.linalg.norm(p2 - p1) / step_size) - 1  # type: ignore[call-overload,operator]
+        for p1, p2 in itertools.pairwise(points)
+    ]
+    # Generate points for each segment
+    segments = [
+        np.linspace(points[i], points[i + 1], num, axis=-1, endpoint=False)
+        for i, num in enumerate(num_points[:-1])
+    ]
+    segments.append(
+        np.linspace(points[-2], points[-1], num_points[-1] + 1, axis=-1, endpoint=True)
+    )
+
+    # Concatenate the points
+    points_arr = np.concatenate(segments, axis=-1)
+    points_arr = np.atleast_2d(points_arr)  # Ensure 2D for 1D paths
+
+    # Distance between each pair of consecutive points
+    distances = np.linalg.norm(np.diff(points_arr, axis=-1), axis=0)
+
+    # Cumulative sum of the distances
+    path_coord = np.concatenate(([0], np.cumsum(distances)))
+
+    interp_coords = [
+        xr.DataArray(p, dims=dim_name, coords={dim_name: path_coord})
+        for p in points_arr
+    ]
+    return darr.interp(
+        dict(zip(vertices.keys(), interp_coords, strict=False)), **interp_kwargs
+    )
+
+
+# Monkey patch xarray to make fast interpolator available
+
 _get_interpolator_nd_original = xarray.core.missing._get_interpolator_nd
 _get_interpolator_original = xarray.core.missing._get_interpolator
 
diff --git a/src/erlab/analysis/kspace.py b/src/erlab/analysis/kspace.py
index ef875ae6..73008357 100644
--- a/src/erlab/analysis/kspace.py
+++ b/src/erlab/analysis/kspace.py
@@ -20,9 +20,9 @@
 
 
 class AxesConfiguration(enum.IntEnum):
-    """
-    Enum class representing different types of axes configurations. See Ref.
-    :cite:p:`ishida2018kconv`.
+    """Enum class representing different types of axes configurations.
+
+    See Ref. :cite:p:`ishida2018kconv`.
 
     """
 
@@ -43,20 +43,21 @@ def _kz_func_inv(kz, inner_potential, kx, ky):
 
 
 def _kperp_func(k_tot_sq, kx, ky):
-    r""":math:`\sqrt{k^2 - k_x^2 - k_y^2}`"""
+    r""":math:`\sqrt{k^2 - k_x^2 - k_y^2}`."""
     return np.sqrt(np.clip(k_tot_sq - kx**2 - ky**2, a_min=0, a_max=None))
 
 
 def kz_func(kinetic_energy, inner_potential, kx, ky):
-    r"""
-    Calculate :math:`k_z` from the given kinetic energy, inner potential, and
-    :math:`k_x`, :math:`k_y` with
+    r"""Calculate the out-of-plane momentum.
+
+    :math:`k_z` is computed from the given kinetic energy :math:`E_k`, inner potential
+    :math:`V_0`, and in-plane momenta :math:`k_x`, and :math:`k_y` by
 
     .. math::
 
         k_z = \sqrt{k^2 - k_x^2 - k_y^2 + \frac{2 m_e V_0}{\hbar^2}}
 
-    where :math:`k = \frac{\sqrt{2 m_e E_k}}{\hbar}`.
+    where :math:`k =\sqrt{2 m_e E_k}/\hbar`.
     """
     k_tot = erlab.constants.rel_kconv * np.sqrt(kinetic_energy)
     k_perp_sq = k_tot**2 - kx**2 - ky**2
@@ -69,9 +70,9 @@ def get_kconv_func(
     configuration: AxesConfiguration,
     angle_params: dict[str, float],
 ) -> tuple[Callable, Callable]:
-    """
-    Returns the appropriate momentum conversion functions based on the given
-    configuration and kinetic energy.
+    """Return appropriate momentum conversion functions.
+
+    The appropriate function is created by the given configuration and kinetic energy.
 
     Parameters
     ----------
diff --git a/src/erlab/analysis/mask/__init__.py b/src/erlab/analysis/mask/__init__.py
index c54032c0..7be20744 100644
--- a/src/erlab/analysis/mask/__init__.py
+++ b/src/erlab/analysis/mask/__init__.py
@@ -155,8 +155,7 @@ def polygon_mask_points(
 def mask_with_hex_bz(
     kxymap: xr.DataArray, a: float = 3.54, rotate: float = 0.0, invert: bool = False
 ) -> xr.DataArray:
-    """Returns map masked with a hexagonal BZ."""
-
+    """Return map masked with a hexagonal BZ."""
     if "kx" in kxymap.dims or "qx" in kxymap.dims:
         dims = ("kx", "ky")
 
@@ -177,8 +176,7 @@ def hex_bz_mask_points(
     reciprocal: bool = False,
     invert: bool = False,
 ) -> npt.NDArray[np.bool_]:
-    """Returns a mask for given points."""
-
+    """Return a mask for given points."""
     if reciprocal:
         d = 2 * np.pi / (a * 3)
     else:
diff --git a/src/erlab/analysis/mask/polygon.py b/src/erlab/analysis/mask/polygon.py
index 42e59843..df2f3116 100644
--- a/src/erlab/analysis/mask/polygon.py
+++ b/src/erlab/analysis/mask/polygon.py
@@ -1,5 +1,6 @@
-"""
-Point-in-polygon algorithm adapted from the `CGAL C++ library
+"""Point-in-polygon algorithm.
+
+The implementation has been adapted from the `CGAL C++ library
 <https://doc.cgal.org/5.3.2/Polygon/index.html>`_.
 
 """
@@ -58,7 +59,7 @@ def _get_argmax_all(arr):
 
 @numba.njit(nogil=True, cache=True)
 def left_vertex(points: Annotated[npt.NDArray[np.float64], Literal[1, 2]]) -> int:
-    """Returns the index of the leftmost point of a polygon.
+    """Return the index of the leftmost point of a polygon.
 
     In case of a tie, the point with the smallest y-coordinate is taken.
 
@@ -78,7 +79,7 @@ def left_vertex(points: Annotated[npt.NDArray[np.float64], Literal[1, 2]]) -> in
 
 @numba.njit(nogil=True, cache=True)
 def right_vertex(points: Annotated[npt.NDArray[np.float64], Literal[1, 2]]):
-    """Returns the index of the rightmost point of a polygon.
+    """Return the index of the rightmost point of a polygon.
 
     In case of a tie, the point with the largest y-coordinate is taken.
 
@@ -136,7 +137,7 @@ def which_side_in_slab(point, low, high, points):
 def bounded_side_bool(
     points: npt.NDArray[np.float64], point: tuple[float, float], boundary: bool = True
 ) -> bool:
-    """Computes whether a point lies inside a polygon using `bounded_side`.
+    """Compute whether a point lies inside a polygon using `bounded_side`.
 
     Parameters
     ----------
@@ -166,7 +167,7 @@ def bounded_side_bool(
 
 @numba.njit(nogil=True, cache=True)
 def bounded_side(points: npt.NDArray[np.float64], point: tuple[float, float]) -> Side:
-    """Computes if a point is inside, outside, or on the boundary of a polygon.
+    """Compute if a point is inside, outside, or on the boundary of a polygon.
 
     The polygon is defined by the sequence of points [first,last). Being inside is
     defined by the odd-even rule. If the point is on a polygon edge, a special value is
@@ -194,7 +195,6 @@ def bounded_side(points: npt.NDArray[np.float64], point: tuple[float, float]) ->
     the bottom vertex is not. (Segments are half-closed).
 
     """
-
     last = len(points) - 1
     if last < 2:
         return Side.ON_UNBOUNDED_SIDE
diff --git a/src/erlab/analysis/utilities.py b/src/erlab/analysis/utilities.py
index b3bb5131..17ed118a 100644
--- a/src/erlab/analysis/utilities.py
+++ b/src/erlab/analysis/utilities.py
@@ -69,7 +69,6 @@ def shift(
     nan 1.0 2.0 4.0 5.0 6.0 nan nan 7.0
     Dimensions without coordinates: x, y
     """
-
     shift_kwargs.setdefault("order", 1)
     shift_kwargs.setdefault("mode", "constant")
     if shift_kwargs["mode"] == "constant":
diff --git a/src/erlab/characterization/__init__.py b/src/erlab/characterization/__init__.py
index cbf7ee98..238c186f 100644
--- a/src/erlab/characterization/__init__.py
+++ b/src/erlab/characterization/__init__.py
@@ -1,5 +1,6 @@
 import warnings
-from erlab.io.characterization import xrd, resistance  # noqa: F401
+
+from erlab.io.characterization import resistance, xrd  # noqa: F401
 
 warnings.warn(
     "`erlab.characterization` is deprecated. Use `erlab.io.characterization` instead",
diff --git a/src/erlab/interactive/__init__.py b/src/erlab/interactive/__init__.py
index d5e55349..0c85c2bb 100644
--- a/src/erlab/interactive/__init__.py
+++ b/src/erlab/interactive/__init__.py
@@ -1,5 +1,5 @@
 """
-Interactive plotting based on Qt and pyqtgraph
+Interactive plotting based on Qt and pyqtgraph.
 
 .. currentmodule:: erlab.interactive
 
diff --git a/src/erlab/interactive/bzplot.py b/src/erlab/interactive/bzplot.py
index b1fcfabb..3eea3f8d 100644
--- a/src/erlab/interactive/bzplot.py
+++ b/src/erlab/interactive/bzplot.py
@@ -1,5 +1,5 @@
 import sys
-from typing import Literal
+from typing import Literal, cast
 
 import mpl_toolkits
 import numpy as np
@@ -15,27 +15,31 @@
 
 
 class BZPlotter(QtWidgets.QMainWindow):
+    """
+    Interactive Brillouin zone plotter.
+
+    Parameters
+    ----------
+    params :
+        Input parameter for plotting. If `param_type` is 'lattice', it must be a
+        6-tuple containing the lengths a, b, c and angles alpha, beta, gamma.
+        Otherwise, it must be a 3 by 3 numpy array with each row vector containing
+        each real/reciprocal lattice vector. If not provided, a hexagonal lattice is
+        shown by default.
+    param_type
+        Specifies the param_type of the input parameters. Valid param_types are
+        `'lattice'`, `'avec'`, `'bvec'`. By default, `'bvec'` is assumed.
+    """
+
     def __init__(
         self,
         params: tuple[float, ...] | npt.NDArray[np.float64] | None = None,
         param_type: Literal["lattice", "avec", "bvec"] | None = None,
     ) -> None:
-        """
-        Parameters
-        ----------
-        params :
-            Input parameter for plotting. If `param_type` is 'lattice', it must be a
-            6-tuple containing the lengths a, b, c and angles alpha, beta, gamma.
-            Otherwise, it must be a 3 by 3 numpy array with each row vector containing
-            each real/reciprocal lattice vector. If not provided, a hexagonal lattice is
-            shown by default.
-        param_type
-            Specifies the param_type of the input parameters. Valid param_types are
-            `'lattice'`, `'avec'`, `'bvec'`. By default, `'bvec'` is assumed.
-        """
-        self.qapp = QtCore.QCoreApplication.instance()
+        self.qapp = cast(QtWidgets.QApplication, QtWidgets.QApplication.instance())
         if not self.qapp:
             self.qapp = QtWidgets.QApplication(sys.argv)
+
         super().__init__()
 
         if params is None:
@@ -64,7 +68,6 @@ def __init__(
             elif param_type == "bvec":
                 bvec = params
 
-        self.controls = None
         self.plot = BZPlotWidget(bvec)
         self.setCentralWidget(self.plot)
 
@@ -279,7 +282,7 @@ def __init__(self, bvec):
 
     def closeEvent(self, event):
         super().closeEvent(event)
-        QtCore.QCoreApplication.instance().quit()
+        QtWidgets.QApplication.instance().quit()
 
     def block_params_signals(self, b: bool):
         self.params_latt.blockSignals(b)
diff --git a/src/erlab/interactive/colors.py b/src/erlab/interactive/colors.py
index a506a3de..9b5e3daa 100644
--- a/src/erlab/interactive/colors.py
+++ b/src/erlab/interactive/colors.py
@@ -101,17 +101,17 @@ def setPopupMinimumWidthForItems(self):
         view = self.view()
         fm = self.fontMetrics()
         maxWidth = max(
-            [fm.boundingRect(self.itemText(i)).width() for i in range(self.count())]
+            fm.boundingRect(self.itemText(i)).width() for i in range(self.count())
         )
-        if maxWidth:
+        if maxWidth and view is not None:
             view.setMinimumWidth(maxWidth)
 
     @QtCore.Slot()
     def nextIndex(self):
         self.wheelEvent(
             QtGui.QWheelEvent(
-                QtCore.QPoint(0, 0),
-                QtCore.QPoint(0, 0),
+                QtCore.QPointF(0, 0),
+                QtCore.QPointF(0, 0),
                 QtCore.QPoint(0, 0),
                 QtCore.QPoint(0, -15),
                 QtCore.Qt.MouseButton.NoButton,
@@ -125,8 +125,8 @@ def nextIndex(self):
     def previousIndex(self):
         self.wheelEvent(
             QtGui.QWheelEvent(
-                QtCore.QPoint(0, 0),
-                QtCore.QPoint(0, 0),
+                QtCore.QPointF(0, 0),
+                QtCore.QPointF(0, 0),
                 QtCore.QPoint(0, 0),
                 QtCore.QPoint(0, 15),
                 QtCore.Qt.MouseButton.NoButton,
@@ -597,7 +597,7 @@ def pg_colormap_names(
 
 
 def pg_colormap_from_name(name: str, skipCache: bool = True) -> pg.ColorMap:
-    """Gets a :class:`pyqtgraph.ColorMap` from its name.
+    """Get a :class:`pyqtgraph.ColorMap` from its name.
 
     Parameters
     ----------
@@ -673,7 +673,7 @@ def mapping_fn(x):
 def pg_colormap_to_QPixmap(
     cmap: str | pg.ColorMap, w: int = 64, h: int = 16, skipCache: bool = True
 ) -> QtGui.QPixmap:
-    """Converts a :class:`pyqtgraph.ColorMap` to a ``w``-by-``h`` QPixmap thumbnail.
+    """Convert a :class:`pyqtgraph.ColorMap` to a ``w``-by-``h`` QPixmap thumbnail.
 
     Parameters
     ----------
@@ -690,7 +690,6 @@ def pg_colormap_to_QPixmap(
     PySide6.QtGui.QPixmap
 
     """
-
     if isinstance(cmap, str):
         cmap = pg_colormap_from_name(cmap, skipCache=skipCache)
     # cmap_arr = np.reshape(cmap.getColors()[:, None], (1, -1, 4), order='C')
diff --git a/src/erlab/interactive/imagetool/__init__.py b/src/erlab/interactive/imagetool/__init__.py
index cd24f379..ca7f942d 100644
--- a/src/erlab/interactive/imagetool/__init__.py
+++ b/src/erlab/interactive/imagetool/__init__.py
@@ -21,7 +21,7 @@
 
 import gc
 import sys
-from typing import TYPE_CHECKING, Any, cast
+from typing import TYPE_CHECKING, Any, Literal, cast
 
 import numpy as np
 import numpy.typing as npt
@@ -36,6 +36,7 @@
 )
 from erlab.interactive.imagetool.core import ImageSlicerArea, SlicerLinkProxy
 from erlab.interactive.utilities import DictMenuBar, copy_to_clipboard
+from erlab.io.plugins.merlin import BL403Loader
 
 if TYPE_CHECKING:
     from collections.abc import Callable, Sequence
@@ -87,7 +88,6 @@ def itool(
     >>> itool(data, cmap="gray", gamma=0.5)
     >>> itool(data_list, link=True)
     """
-
     qapp = QtWidgets.QApplication.instance()
     if not qapp:
         qapp = QtWidgets.QApplication(sys.argv)
@@ -407,7 +407,9 @@ def refreshMenus(self):
             self.colorAct, ["reversed", "highContrast", "zeroCentered"], strict=True
         ):
             ca.blockSignals(True)
-            ca.setChecked(cmap_props[k])
+            ca.setChecked(
+                cmap_props[cast(Literal["reversed", "highContrast", "zeroCentered"], k)]
+            )
             ca.blockSignals(False)
 
     def _set_colormap_options(self):
@@ -424,10 +426,11 @@ def _copy_cursor_idx(self):
         copy_to_clipboard(str(self.slicer_area.array_slicer._indices))
 
     def _open_file(self):
+        merlin_loader = cast(BL403Loader, erlab.io.loaders["merlin"])
         valid_loaders: dict[str, tuple[Callable, dict]] = {
             "xarray HDF5 Files (*.h5)": (erlab.io.load_hdf5, {}),
-            "ALS BL4.0.3 Raw Data (*.pxt)": (erlab.io.loaders["merlin"].load, {}),
-            "ALS BL4.0.3 Live (*.ibw)": (erlab.io.loaders["merlin"].load_live, {}),
+            "ALS BL4.0.3 Raw Data (*.pxt)": (merlin_loader.load, {}),
+            "ALS BL4.0.3 Live (*.ibw)": (merlin_loader.load_live, {}),
             "DA30 Raw Data (*.ibw *.pxt *.zip)": (erlab.io.loaders["da30"].load, {}),
             "SSRL BL5-2 Raw Data (*.h5)": (erlab.io.loaders["ssrl52"].load, {}),
             "NetCDF Files (*.nc *.nc4 *.cdf)": (xr.load_dataarray, {}),
diff --git a/src/erlab/interactive/imagetool/_deprecated/imagetool_mpl.py b/src/erlab/interactive/imagetool/_deprecated/imagetool_mpl.py
index 6c1fee6d..df16c621 100644
--- a/src/erlab/interactive/imagetool/_deprecated/imagetool_mpl.py
+++ b/src/erlab/interactive/imagetool/_deprecated/imagetool_mpl.py
@@ -1,5 +1,5 @@
 """
-.. deprecated:: 0.1
+.. deprecated:: 0.1.
 
     This module is deprecated, and is only kept for reference purposes.
     Use `erlab.interactive.imagetool` instead.
@@ -30,7 +30,7 @@
 
 
 def qt_style_names():
-    """Return a list of styles, default platform style first"""
+    """Return a list of styles, default platform style first."""
     default_style_name = QtWidgets.QApplication.style().objectName().lower()
     result = []
     for style in QtWidgets.QStyleFactory.keys():
@@ -46,7 +46,7 @@ def change_style(style_name):
 
 
 def colormap_to_QPixmap(name: str, h=64):
-    """Convert matplotlib colormap to a 256-by-`h` QPixmap"""
+    """Convert matplotlib colormap to a 256-by-`h` QPixmap."""
     cmap = plt.colormaps[name]
     cmap_arr = cmap(np.tile(np.linspace(0, 1, 256), (h, 1))) * 255
     img = QtGui.QImage(
@@ -171,8 +171,7 @@ def setColor(self, color):
         self.colorChanged.emit(color.getRgbF())
         if self._color:
             self.setStyleSheet(
-                "QWidget { background-color: %s; border: 0; }"
-                % self._color.name(QtGui.QColor.HexArgb)
+                f"QWidget {{ background-color: {self._color.name(QtGui.QColor.HexArgb)}; border: 0; }}"
             )
         else:
             self.setStyleSheet("")
@@ -1049,9 +1048,10 @@ def __init__(self, canvas, parent, coordinates=True):
         NavigationToolbar2QT.__init__(self, canvas, parent, coordinates=coordinates)
 
     def _icon(self, name):
-        """
-        Construct a `.QIcon` from an image file *name*, including the extension
-        and relative to Matplotlib's "images" data directory.
+        """Construct a `.QIcon` from an image file name.
+
+        The name must include the extension and be given relative to Matplotlib's
+        "images" data directory.
         """
         name = name.replace(".png", "")
         icons_dict = {
diff --git a/src/erlab/interactive/imagetool/_deprecated/imagetool_old.py b/src/erlab/interactive/imagetool/_deprecated/imagetool_old.py
index 17c805e0..5663e984 100644
--- a/src/erlab/interactive/imagetool/_deprecated/imagetool_old.py
+++ b/src/erlab/interactive/imagetool/_deprecated/imagetool_old.py
@@ -1,5 +1,5 @@
 """
-.. deprecated:: 0.1
+.. deprecated:: 0.1.
 
     This module is deprecated, and is only kept for reference purposes.
     Use `erlab.interactive.imagetool` instead.
@@ -8,6 +8,7 @@
 
 import colorsys
 import enum
+import importlib
 import sys
 import weakref
 from itertools import chain, compress
@@ -15,7 +16,6 @@
 
 import matplotlib.mathtext
 import numba
-import numbagg
 import numpy as np
 import pyqtgraph as pg
 import qtawesome as qta
@@ -38,6 +38,12 @@
 # pg.setConfigOption('background', 'w')
 # pg.setConfigOption('foreground', 'k')
 
+if importlib.util.find_spec("numbagg"):
+    import numbagg
+
+    _general_nanmean_func = numbagg.nanmean
+else:
+    _general_nanmean_func = np.nanmean
 
 __all__ = ["itool_", "pg_itool"]
 
@@ -243,7 +249,7 @@ def __init__(self, *args, **kwargs):
 
 
 def qt_style_names():
-    """Return a list of styles, default platform style first"""
+    """Return a list of styles, default platform style first."""
     default_style_name = QtWidgets.QApplication.style().objectName().lower()
     result = []
     for style in QtWidgets.QStyleFactory.keys():
@@ -362,8 +368,7 @@ def setColor(self, color):
         self.colorChanged.emit(color.getRgbF())
         if self._color:
             self.setStyleSheet(
-                "QWidget { background-color: %s; border: 0; }"
-                % self._color.name(QtGui.QColor.HexArgb)
+                f"QWidget {{ background-color: {self._color.name(QtGui.QColor.HexArgb)}; border: 0; }}"
             )
         else:
             self.setStyleSheet("")
@@ -504,7 +509,7 @@ def changeEvent(self, evt):
 
 
 def get_pixmap_label(s: str, prop=None, dpi=300, **text_kw):
-    """Creates a QtGui.QPixmap from a mathtext string.
+    """Create a QtGui.QPixmap from a mathtext string.
 
     Parameters
     ----------
@@ -520,11 +525,9 @@ def get_pixmap_label(s: str, prop=None, dpi=300, **text_kw):
 
     Returns
     -------
-
     A QtGui.QPixmap object.
 
     """
-
     parser = matplotlib.mathtext.MathTextParser("path")
     if prop is None:
         prop = FontProperties(size=9)
@@ -545,7 +548,7 @@ def get_pixmap_label(s: str, prop=None, dpi=300, **text_kw):
 def get_svg_label(
     s: str, outfile: QtCore.QTemporaryFile, prop=None, dpi=300, **text_kw
 ):
-    """Creates an SVG image from a mathtext string.
+    """Create an SVG image from a mathtext string.
 
     Parameters
     ----------
@@ -567,7 +570,6 @@ def get_svg_label(
         Name of the output file containing the rendered SVG.
 
     """
-
     parser = matplotlib.mathtext.MathTextParser("path")
     if prop is None:
         prop = FontProperties(size=12)
@@ -888,9 +890,6 @@ class pg_itool(pg.GraphicsLayoutWidget):
 
     """
 
-    # !TODO: ctrl + A to view all
-    # !TODO: auto adjust limits on transpose
-
     sigDataChanged = QtCore.Signal(object)  #: :meta private:
     sigIndexChanged = QtCore.Signal(list, list)  #: :meta private:
 
@@ -1424,7 +1423,7 @@ def toggle_axes(self, axis):
             self.addItem(self.axes[group[1]], *anchors[1], *ref_dims[axis][2:])
 
     def set_labels(self, labels=None):
-        """labels: list or tuple of str"""
+        """labels: list or tuple of str."""
         if labels is None:
             labels = self.data_dims
         # 0: default, 1: svg, 2: pixmap
@@ -1581,7 +1580,7 @@ def _measure_fps(self):
         self._fpsLastUpdate = now
         w = 0.8
         self._avg_fps = self._avg_fps * (1 - w) + fps * w
-        self.axes[1].setTitle("%0.2f fps" % self._avg_fps)
+        self.axes[1].setTitle(f"{self._avg_fps:0.2f} fps")
 
     def labelify(self, text):
         """Prettify some frequently used axis labels."""
@@ -1698,7 +1697,7 @@ def _get_binned_data(self, axis):
                 + (self._get_bin_slice(axis + 1),)
             ].squeeze(axis=axis)
         else:
-            return numbagg.nanmean(
+            return _general_nanmean_func(
                 self.data_vals_T[
                     (slice(None),) * (axis % self.data_ndim)
                     + (self._get_bin_slice(axis + 1),)
@@ -1712,11 +1711,11 @@ def _binned_profile(self, avg_axis):
         slices = tuple(self._get_bin_slice(ax) for ax in avg_axis)
         return self._block_slice_avg(avg_axis, slices)
         # self._slice_block = self._block_slicer(avg_axis, slices)
-        # return numbagg.nanmean(self._slice_block, axis=[(ax - 1) for ax in avg_axis])
+        # return _general_nanmean_func(self._slice_block, axis=[(ax - 1) for ax in avg_axis])
 
     def _block_slice_avg(self, axis=None, slices=None):
         axis = [(ax - 1) % self.data_ndim for ax in axis]
-        return numbagg.nanmean(
+        return _general_nanmean_func(
             self.data_vals_T[
                 tuple(
                     slices[axis.index(d)] if d in axis else slice(None)
@@ -1797,7 +1796,7 @@ def get_key_modifiers(self):
         return mods
 
     def _get_mouse_datapos(self, plot, pos):
-        """Returns mouse position in data coords"""
+        """Return mouse position in data coords."""
         mouse_point = plot.vb.mapSceneToView(pos)
         return mouse_point.x(), mouse_point.y()
 
@@ -2313,7 +2312,6 @@ def fast_isocurve(data, level, connected=False, extendToEdge=False, path=False):
                    vertex coordinates.
     ============== =========================================================
     """
-
     if path is True:
         connected = True
     np.nan_to_num(data, copy=False)
@@ -3054,7 +3052,7 @@ def showPopup(self):
     def setPopupMinimumWidthForItems(self):
         view = self.view()
         fm = self.fontMetrics()
-        maxWidth = max([fm.width(self.itemText(i)) for i in range(self.count())])
+        maxWidth = max(fm.width(self.itemText(i)) for i in range(self.count()))
         if maxWidth:
             view.setMinimumWidth(maxWidth)
 
diff --git a/src/erlab/interactive/imagetool/controls.py b/src/erlab/interactive/imagetool/controls.py
index e069c93a..aa42f3b6 100644
--- a/src/erlab/interactive/imagetool/controls.py
+++ b/src/erlab/interactive/imagetool/controls.py
@@ -1,26 +1,31 @@
+"""Widgets for controlling `ImageSlicerArea`."""
+
 from __future__ import annotations
 
 __all__ = [
-    "ColorControls",
     "ItoolBinningControls",
+    "ItoolColorControls",
     "ItoolColormapControls",
+    "ItoolControlsBase",
     "ItoolCrosshairControls",
 ]
 
-from typing import TYPE_CHECKING
+import types
+from typing import TYPE_CHECKING, cast
 
 import numpy as np
 import pyqtgraph as pg
 import qtawesome as qta
 from qtpy import QtCore, QtGui, QtWidgets
-import types
+
 from erlab.interactive.colors import ColorMapComboBox, ColorMapGammaWidget
 from erlab.interactive.utilities import BetterSpinBox
 
 if TYPE_CHECKING:
-    import xarray as xr
     from collections.abc import Mapping
 
+    import xarray as xr
+
     from erlab.interactive.imagetool.core import ImageSlicerArea
     from erlab.interactive.imagetool.slicer import ArraySlicer
 
@@ -138,10 +143,10 @@ def current_cursor(self) -> int:
         return self.slicer_area.current_cursor
 
     def initialize_layout(self):
-        self.setLayout(QtWidgets.QHBoxLayout(self))
-        self.layout().setContentsMargins(0, 0, 0, 0)
-
-        self.layout().setSpacing(3)
+        layout = QtWidgets.QHBoxLayout(self)
+        self.setLayout(layout)
+        layout.setContentsMargins(0, 0, 0, 0)
+        layout.setSpacing(3)
 
     def initialize_widgets(self):
         for ctrl in self.sub_controls:
@@ -182,19 +187,19 @@ def slicer_area(self) -> ImageSlicerArea:
 
     @slicer_area.setter
     def slicer_area(self, value: ImageSlicerArea):
-        """
-        Initially, the goal was to be able to control multiple `ImageSlicerArea`s
-        with a single control widget, so the control widgets were designed with easy
-        connection and disconnection of signals in mind. However, this has become
-        largely unnecessary since we now have `SlicerLinkProxy` as a translation layer
-        between multiple `ImageSlicerArea`s with individual control widgets. Hence, this
-        method will remain unused for the time being.
+        """Set the `ImageSlicerArea` instance for the control widget.
+
+        Initially, the goal was to be able to control multiple `ImageSlicerArea`s with a
+        single control widget, so the control widgets were designed with easy connection
+        and disconnection of signals in mind. However, this has become largely
+        unnecessary since we now have `SlicerLinkProxy` as a translation layer between
+        multiple `ImageSlicerArea`s with individual control widgets. Hence, this method
+        will remain unused for the time being.
 
         Also, in principle, most of the control widgets along with the menu bar should
         be re-written to use QActions...
 
         """
-
         # ignore until https://bugreports.qt.io/browse/PYSIDE-229 is fixed
         try:
             self.disconnect_signals()
@@ -305,12 +310,15 @@ def initialize_widgets(self):
             self.values_layouts[0].addWidget(self.btn_snap, 0, 3, 1, 1)
             self.values_layouts[0].addWidget(self.cb_cursors, 0, 0, 1, 1)
             self.values_layouts[0].addWidget(self.spin_dat, 0, 4, 1, 1)
-        self.layout().addWidget(self.values_groups[0])
+        cast(QtWidgets.QLayout, self.layout()).addWidget(self.values_groups[0])
 
         # info widgets
         self.label_dim = tuple(
-            QtWidgets.QPushButton(grp, checkable=True) for grp in self.values_groups[1:]
+            QtWidgets.QPushButton(grp) for grp in self.values_groups[1:]
         )
+        for lab in self.label_dim:
+            lab.setCheckable(True)
+
         self.spin_idx = tuple(
             BetterSpinBox(
                 grp,
@@ -364,7 +372,7 @@ def initialize_widgets(self):
                 self.values_layouts[i + 1].addWidget(self.spin_idx[i], 0, 2, 1, 1)
                 self.values_layouts[i + 1].addWidget(self.spin_val[i], 0, 3, 1, 1)
 
-            self.layout().addWidget(self.values_groups[i + 1])
+            cast(QtWidgets.QLayout, self.layout()).addWidget(self.values_groups[i + 1])
 
     def _transpose_axes(self, idx):
         if self.data.ndim == 4:
@@ -526,7 +534,7 @@ def setActiveCursor(self, value: str):
         self.slicer_area.set_current_cursor(self.cb_cursors.findText(value))
 
 
-class ColorControls(ItoolControlsBase):
+class ItoolColorControls(ItoolControlsBase):
     def __init__(self, *args, **kwargs):
         super().__init__(*args, **kwargs)
 
@@ -558,10 +566,11 @@ def initialize_widgets(self):
         self.btn_contrast.toggled.connect(self.update_colormap)
         self.btn_zero.toggled.connect(self.update_colormap)
 
-        self.layout().addWidget(self.btn_reverse)
-        self.layout().addWidget(self.btn_contrast)
-        self.layout().addWidget(self.btn_zero)
-        self.layout().addWidget(self.btn_lock)
+        layout = cast(QtWidgets.QLayout, self.layout())
+        layout.addWidget(self.btn_reverse)
+        layout.addWidget(self.btn_contrast)
+        layout.addWidget(self.btn_zero)
+        layout.addWidget(self.btn_lock)
 
     def update(self):
         self.btn_reverse.blockSignals(True)
@@ -614,9 +623,10 @@ def initialize_layout(self):
             self.setLayout(QtWidgets.QVBoxLayout(self))
         else:
             self.setLayout(QtWidgets.QHBoxLayout(self))
-        self.layout().setContentsMargins(0, 0, 0, 0)
 
-        self.layout().setSpacing(3)
+        layout = cast(QtWidgets.QBoxLayout, self.layout())
+        layout.setContentsMargins(0, 0, 0, 0)
+        layout.setSpacing(3)
 
     def initialize_widgets(self):
         super().initialize_widgets()
@@ -631,11 +641,12 @@ def initialize_widgets(self):
             QtWidgets.QSizePolicy.Policy.Fixed, QtWidgets.QSizePolicy.Policy.Fixed
         )
 
-        self.misc_controls = self.add_control(ColorControls(self))
+        self.misc_controls = self.add_control(ItoolColorControls(self))
 
-        self.layout().addWidget(self.cb_colormap)
-        self.layout().addWidget(self.gamma_widget)
-        self.layout().addWidget(self.misc_controls)
+        layout = cast(QtWidgets.QBoxLayout, self.layout())
+        layout.addWidget(self.cb_colormap)
+        layout.addWidget(self.gamma_widget)
+        layout.addWidget(self.misc_controls)
 
     def update(self):
         super().update()
@@ -666,9 +677,10 @@ def __init__(self, *args, **kwargs):
 
     def initialize_layout(self):
         self.setLayout(QtWidgets.QGridLayout(self))
-        self.layout().setContentsMargins(0, 0, 0, 0)
 
-        self.layout().setSpacing(3)
+        layout = cast(QtWidgets.QGridLayout, self.layout())
+        layout.setContentsMargins(0, 0, 0, 0)
+        layout.setSpacing(3)
 
     def initialize_widgets(self):
         super().initialize_widgets()
@@ -698,11 +710,12 @@ def initialize_widgets(self):
             toolTip="Apply for all cursors",
         )
 
+        layout = cast(QtWidgets.QGridLayout, self.layout())
         for i in range(self.data.ndim):
-            self.layout().addWidget(self.labels[i], 0, i, 1, 1)
-            self.layout().addWidget(self.spins[i], 1, i, 1, 1)
-        self.layout().addWidget(self.reset_btn, 2, 0, 1, 1)
-        self.layout().addWidget(self.all_btn, 2, 1, 1, 1)
+            layout.addWidget(self.labels[i], 0, i, 1, 1)
+            layout.addWidget(self.spins[i], 1, i, 1, 1)
+        layout.addWidget(self.reset_btn, 2, 0, 1, 1)
+        layout.addWidget(self.all_btn, 2, 1, 1, 1)
         # for spin in self.spins:
         # spin.setMinimumWidth(60)
 
diff --git a/src/erlab/interactive/imagetool/core.py b/src/erlab/interactive/imagetool/core.py
index 56998947..aea84340 100644
--- a/src/erlab/interactive/imagetool/core.py
+++ b/src/erlab/interactive/imagetool/core.py
@@ -115,9 +115,7 @@ def link_slicer(
     steps: bool = False,
     color: bool = False,
 ):
-    """
-    An internal decorator for choosing which functions to sync accross multiple
-    instances of `ImageSlicerArea`.
+    """Sync decorated methods across multiple `ImageSlicerArea` instances.
 
     Parameters
     ----------
@@ -153,9 +151,15 @@ def wrapped(*args, **kwargs):
                     )
                     all_args.apply_defaults()
                     obj: ImageSlicerArea = all_args.arguments.pop("self")
-                    obj._linking_proxy.sync(
-                        obj, func.__name__, all_args.arguments, indices, steps, color
-                    )
+                    if obj._linking_proxy is not None:
+                        obj._linking_proxy.sync(
+                            obj,
+                            func.__name__,
+                            all_args.arguments,
+                            indices,
+                            steps,
+                            color,
+                        )
             return out
 
         return wrapped
@@ -205,7 +209,7 @@ def sync(
         steps: bool,
         color: bool,
     ):
-        """The core method that propagates changes across multiple `ImageSlicerArea`s.
+        """Propagate changes across multiple `ImageSlicerArea`s.
 
         This method is invoked every time a method decorated with :func:`link_slicer` in
         a linked `ImageSlicerArea` is called.
@@ -333,22 +337,22 @@ class ImageSlicerArea(QtWidgets.QWidget):
 
     @property
     def sigCursorCountChanged(self) -> QtCore.SignalInstance:
-        """:meta private:"""
+        """:meta private:"""  # noqa: D400
         return self.array_slicer.sigCursorCountChanged
 
     @property
     def sigIndexChanged(self) -> QtCore.SignalInstance:
-        """:meta private:"""
+        """:meta private:"""  # noqa: D400
         return self.array_slicer.sigIndexChanged
 
     @property
     def sigBinChanged(self) -> QtCore.SignalInstance:
-        """:meta private:"""
+        """:meta private:"""  # noqa: D400
         return self.array_slicer.sigBinChanged
 
     @property
     def sigShapeChanged(self) -> QtCore.SignalInstance:
-        """:meta private:"""
+        """:meta private:"""  # noqa: D400
         return self.array_slicer.sigShapeChanged
 
     def __init__(
@@ -405,7 +409,7 @@ def __init__(
         self._colorbar.setVisible(False)
 
         pkw = {"image_cls": image_cls, "plotdata_cls": plotdata_cls}
-        self.manual_limits: dict[str, list[list[float]]] = {}
+        self.manual_limits: dict[str, list[float]] = {}
         self._plots: tuple[ItoolGraphicsLayoutWidget, ...] = (
             ItoolGraphicsLayoutWidget(self, image=True, display_axis=(0, 1), **pkw),
             ItoolGraphicsLayoutWidget(self, display_axis=(0,), **pkw),
@@ -458,7 +462,7 @@ def __init__(
     def on_close(self):
         self.array_slicer.clear_cache()
         self.data.close()
-        if hasattr(self, "_data"):
+        if hasattr(self, "_data") and self._data is not None:
             self._data.close()
             del self._data
 
@@ -480,7 +484,8 @@ def add_link(self, proxy: SlicerLinkProxy):
         proxy.add(self)
 
     def remove_link(self):
-        self._linking_proxy.remove(self)
+        if self.is_linked:
+            cast(SlicerLinkProxy, self._linking_proxy).remove(self)
 
     @property
     def is_linked(self) -> bool:
@@ -514,7 +519,7 @@ def profiles(self) -> tuple[ItoolPlotItem, ...]:
 
     @property
     def main_image(self) -> ItoolPlotItem:
-        """returns the main PlotItem"""
+        """Return the main PlotItem."""
         return self.get_axes(0)
 
     @property
@@ -523,7 +528,7 @@ def images(self) -> tuple[ItoolPlotItem, ...]:
 
     @property
     def axes(self) -> tuple[ItoolPlotItem, ...]:
-        """Currently valid subset of self._plots"""
+        """Currently valid subset of self._plots."""
         return self.images + self.profiles
 
     @property
@@ -618,6 +623,22 @@ def set_current_cursor(self, cursor: int, update=True):
     def set_data(
         self, data: xr.DataArray | npt.ArrayLike, rad2deg: bool | Iterable[str] = False
     ):
+        """Set the data to be displayed.
+
+        Parameters
+        ----------
+        data
+            The data to be displayed. If a `xarray.DataArray` is given, the
+            dimensions and coordinates are used to determine the axes of the plots. If a
+            :class:`xarray.Dataset` is given, the first data variable is used. If a
+            :class:`numpy.ndarray` is given, it is converted to a `xarray.DataArray`
+            with default dimensions.
+        rad2deg
+            If `True`, converts coords along dimensions that have angle-like names to
+            degrees. If an iterable of strings is given, coordinates for dimensions that
+            correspond to the given strings are converted.
+
+        """
         if hasattr(self, "_array_slicer"):
             n_cursors_old = self.n_cursors
             if isinstance(self._data, xr.DataArray):
@@ -670,6 +691,52 @@ def set_data(
         self._colorbar.cb.setImageItem()
         self.lock_levels(False)
 
+    def update_values(self, values: npt.NDArray | xr.DataArray, update: bool = True):
+        """Update only the values of the data.
+
+        The coords and shape of the data array are not changed.
+
+        Parameters
+        ----------
+        values
+            The new values to be set. If a `xarray.DataArray` is given, the dimensions
+            must match the current data array. If a `numpy.ndarray` is given, the shape
+            must match the current data array. Note that if the user has transposed the
+            current data array, passing a `numpy.ndarray` with the original shape will
+            fail.
+        update
+            If `True`, the plots are updated after setting the new values.
+
+        Note
+        ----
+        This method only checks for matching dimension name and shape, and does not
+        check for equal coordinate values.
+
+        """
+        if isinstance(values, xr.DataArray):
+            if self.array_slicer._obj.ndim != values.ndim:
+                raise ValueError("DataArray dimensions do not match")
+            if set(self.array_slicer._obj.dims) != set(values.dims):
+                raise ValueError("DataArray dimensions do not match")
+
+            if self.array_slicer._obj.dims != values.dims:
+                values = values.transpose(*self.array_slicer._obj.dims)
+            if self.array_slicer._obj.shape != values.shape:
+                raise ValueError("DataArray shape does not match")
+
+            values = values.values
+        else:
+            if self.array_slicer._obj.shape != values.shape:
+                raise ValueError(
+                    "Data shape does not match. Array is "
+                    f"{self.array_slicer._obj.shape} but {values.shape} given"
+                )
+        self.array_slicer._obj[:] = values
+
+        if update:
+            self.array_slicer.clear_val_cache(include_vals=True)
+            self.refresh_all(only_plots=True)
+
     @QtCore.Slot(int, int)
     @link_slicer
     def swap_axes(self, ax1: int, ax2: int):
@@ -848,7 +915,7 @@ def adjust_layout(
         font_size: float = 11.0,
         r: tuple[float, float, float, float] = (1.2, 1.5, 3.0, 1.0),
     ):
-        """Determines the padding and aspect ratios.
+        """Determine the padding and aspect ratios.
 
         Parameters
         ----------
@@ -879,7 +946,6 @@ def adjust_layout(
                   r[3] * r[2]
 
         """
-
         font = QtGui.QFont()
         font.setPointSizeF(float(font_size))
 
@@ -1233,7 +1299,7 @@ def refresh_manual_range(self):
         ):
             if dim is not None:
                 if auto:
-                    self.slicer_area.manual_limits.pop("dim", None)
+                    self.slicer_area.manual_limits.pop(dim, None)
                 else:
                     self.slicer_area.manual_limits[dim] = rng
 
@@ -1510,6 +1576,10 @@ def set_active_cursor(self, index: int):
                 item.setVisible(i == index)
 
     def save_current_data(self, fileName=None):
+        default_name = "data"
+        if self.slicer_area._data is not None:
+            default_name = str(self.slicer_area._data.name)
+
         if fileName is None:
             self.fileDialog = QtWidgets.QFileDialog()
             self.fileDialog.setFileMode(QtWidgets.QFileDialog.FileMode.AnyFile)
@@ -1517,13 +1587,11 @@ def save_current_data(self, fileName=None):
             self.fileDialog.setNameFilter("xarray HDF5 Files (*.h5)")
             if pg.PlotItem.lastFileDir is not None:
                 self.fileDialog.setDirectory(
-                    os.path.join(
-                        pg.PlotItem.lastFileDir, f"{self.slicer_area._data.name}.h5"
-                    )
+                    os.path.join(pg.PlotItem.lastFileDir, f"{default_name}.h5")
                 )
             else:
                 self.fileDialog.setDirectory(
-                    os.path.join(os.getcwd(), f"{self.slicer_area._data.name}.h5")
+                    os.path.join(os.getcwd(), f"{default_name}.h5")
                 )
             self.fileDialog.show()
             self.fileDialog.fileSelected.connect(self.save_current_data)
diff --git a/src/erlab/interactive/imagetool/fastbinning.py b/src/erlab/interactive/imagetool/fastbinning.py
index 15990525..d34dba80 100644
--- a/src/erlab/interactive/imagetool/fastbinning.py
+++ b/src/erlab/interactive/imagetool/fastbinning.py
@@ -1,20 +1,30 @@
-"""
-This module provides a fast, parallelized version of nanmean that supports multiple axes
-based on numba. Enables efficient real-time multidimensional binning.
+"""Fast parallelized averaging for multidimensional arrays.
+
+This module provides a numba-based fast, parallelized version of nanmean that supports
+multiple axes. This enables efficient real-time multidimensional binning.
 
 """
 
-__all__ = ["fast_nanmean"]
+__all__ = ["NANMEAN_FUNCS", "fast_nanmean", "fast_nanmean_skipcheck"]
 
-from collections.abc import Collection
+import importlib
+from collections.abc import Callable, Collection
 
 import numba
 import numba.core.registry
 import numba.typed
-import numbagg
 import numpy as np
 import numpy.typing as npt
 
+if importlib.util.find_spec("numbagg"):
+    import numbagg
+
+    _general_nanmean_func: Callable = numbagg.nanmean
+else:
+    _general_nanmean_func = np.nanmean
+
+# _SIG_N_M: List of signatures that reduces from N to M dimensions.
+
 _SIG_2_1 = [
     numba.types.Array(numba.float64, 1, "C")(numba.types.Array(numba.float32, 2, "C")),
     numba.types.Array(numba.float64, 1, "C")(numba.types.Array(numba.float32, 2, "A")),
@@ -273,7 +283,7 @@ def _nanmean_4_123(a: npt.NDArray[np.float32 | np.float64]) -> npt.NDArray[np.fl
     return output
 
 
-nanmean_funcs = {
+NANMEAN_FUNCS: dict[int, dict[int | frozenset[int], Callable]] = {
     2: {
         0: _nanmean_2_0,
         1: _nanmean_2_1,
@@ -311,13 +321,13 @@ def _nanmean_4_123(a: npt.NDArray[np.float32 | np.float64]) -> npt.NDArray[np.fl
         frozenset({0, 2, 3}): _nanmean_4_023,
         frozenset({1, 2, 3}): _nanmean_4_123,
     },
-}
+}  #: Mapping from array dimensions to axis combinations to corresponding functions.
 
 
 def fast_nanmean(
     a: npt.NDArray[np.float32 | np.float64], axis: int | Collection[int] | None = None
 ) -> npt.NDArray[np.float32 | np.float64] | np.float64:
-    """A fast, parallelized arithmetic mean for floating point arrays that ignores NaNs.
+    """Compute the mean for floating point arrays while ignoring NaNs.
 
     Parameters
     ----------
@@ -332,43 +342,56 @@ def fast_nanmean(
     numpy.ndarray or float
         The calculated mean. The output array is always C-contiguous.
 
-    Note
-    ----
-    Parallelization is only applied for :code:`N`-dimensional arrays with :code:`N <= 4`
-    and :code:`len(axis) < N`. For bigger :code:`N`, :obj:`numbagg.nanmean` is used. For
-    calculating the average of a flattened array (:code:`axis = None` or
-    :code:`len(axis) == N`), the :obj:`numba` implemenation of :obj:`numpy.nanmean` is
-    used. This function does not keep the input dimensions, i.e., the output is
-    squeezed.
+    Notes
+    -----
+    - Parallelization is only applied for ``N``-dimensional arrays with ``N <= 4`` and
+      ``len(axis) < N``.
+
+    - For calculating the average of a flattened array (``axis = None`` or ``len(axis)
+      == N``), the `numba` implemenation of `numpy.nanmean` is used.
+
+    - For bigger ``N``, ``numbagg.nanmean`` is used if `numbagg
+      <https://github.com/numbagg/numbagg>`_ is installed. Otherwise, the calculation
+      falls back to `numpy.nanmean`.
+
+    - This function does not keep the input dimensions, i.e., the output is squeezed.
+
+    - For single precision input, the calculation is performed in double precision and
+      converted back to single precision. This may lead to different results compared to
+      `numpy.nanmean`.
 
     """
     if a.ndim == 1 or axis is None:
         return _nanmean_all(a)
     elif a.ndim > 4:
-        return np.ascontiguousarray(numbagg.nanmean(a, axis))  # type: ignore[arg-type]
+        return np.ascontiguousarray(
+            _general_nanmean_func(a.astype(np.float64), axis), dtype=a.dtype
+        )
     if isinstance(axis, Collection):
         if len(axis) == a.ndim:
             return _nanmean_all(a)
         axis = frozenset(x % a.ndim for x in axis)
     else:
         axis = axis % a.ndim
-    return nanmean_funcs[a.ndim][axis](a).astype(a.dtype)
+    return NANMEAN_FUNCS[a.ndim][axis](a).astype(a.dtype)
 
 
-def _fast_nanmean_skipcheck(
+def fast_nanmean_skipcheck(
     a: npt.NDArray[np.float32 | np.float64], axis: int | Collection[int]
 ) -> npt.NDArray[np.float32 | np.float64] | np.float64:
-    """A version of `fast_nanmean` with near-zero overhead. Meant for internal use.
+    """Compute the mean for specific floating point arrays while ignoring NaNs.
 
-    Strict assumptions on the input parameters allow skipping some checks.
+    This is a version of `fast_nanmean` with near-zero overhead meant for internal use.
+    Strict assumptions on the input parameters allow skipping some checks. Failure to
+    meet these assumptions may lead to undefined behavior.
 
     Parameters
     ----------
     a
-        A numpy array of floats. :code:`a.ndim` must be one of 2, 3, and 4.
+        A numpy array of floats. ``a.ndim`` must be one of 2, 3, and 4.
     axis
         Axis or iterable of axis along which the means are computed. All elements must
-        be nonnegative integers that are less than or equal to :code:`a.ndim`, i.e.,
+        be nonnegative integers that are less than or equal to ``a.ndim``, i.e.,
         negative indexing is not allowed.
 
     Returns
@@ -381,4 +404,4 @@ def _fast_nanmean_skipcheck(
         if len(axis) == a.ndim:
             return _nanmean_all(a)
         axis = frozenset(axis)
-    return nanmean_funcs[a.ndim][axis](a).astype(a.dtype)
+    return NANMEAN_FUNCS[a.ndim][axis](a).astype(a.dtype)
diff --git a/src/erlab/interactive/imagetool/slicer.py b/src/erlab/interactive/imagetool/slicer.py
index 633d81db..a092837a 100644
--- a/src/erlab/interactive/imagetool/slicer.py
+++ b/src/erlab/interactive/imagetool/slicer.py
@@ -12,10 +12,10 @@
 import numpy.typing as npt
 from qtpy import QtCore
 
-from erlab.interactive.imagetool.fastbinning import _fast_nanmean_skipcheck
+from erlab.interactive.imagetool.fastbinning import fast_nanmean_skipcheck
 
 if TYPE_CHECKING:
-    from collections.abc import Sequence, Hashable
+    from collections.abc import Hashable, Sequence
 
     import xarray as xr
 
@@ -229,12 +229,12 @@ def absnanmin(self) -> float:
 
     @property
     def limits(self) -> tuple[float, float]:
-        """Returns the global minima and maxima of the data."""
+        """Return the global minima and maxima of the data."""
         return self.nanmin, self.nanmax
 
     @staticmethod
     def validate_array(data: xr.DataArray) -> xr.DataArray:
-        """Validates a given :class:`xarray.DataArray`.
+        """Validate a given :class:`xarray.DataArray`.
 
         If data has two momentum axes (``kx`` and ``ky``), set them (and ``eV`` if
         exists) as the first two (or three) dimensions. Then, checks the data for
@@ -282,7 +282,7 @@ def validate_array(data: xr.DataArray) -> xr.DataArray:
     def reset_property_cache(self, propname: str) -> None:
         self.__dict__.pop(propname, None)
 
-    def clear_dim_cache(self):
+    def clear_dim_cache(self, include_vals: bool = False):
         for prop in (
             "coords",
             "coords_uniform",
@@ -290,17 +290,22 @@ def clear_dim_cache(self):
             "incs_uniform",
             "lims",
             "lims_uniform",
-            "data_vals_T",
         ):
             self.reset_property_cache(prop)
 
-    def clear_val_cache(self):
+        if include_vals:
+            self.reset_property_cache("data_vals_T")
+
+    def clear_val_cache(self, include_vals: bool = False):
         for prop in ("nanmax", "nanmin", "absnanmax", "absnanmin"):
             self.reset_property_cache(prop)
 
+        if include_vals:
+            self.reset_property_cache("data_vals_T")
+
     def clear_cache(self):
         self.clear_dim_cache()
-        self.clear_val_cache()
+        self.clear_val_cache(include_vals=True)
 
     def set_array(
         self, xarray_obj: xr.DataArray, validate: bool = True, reset: bool = False
@@ -316,9 +321,9 @@ def set_array(
             i for i, d in enumerate(self._obj.dims) if str(d).endswith("_idx")
         ]
 
-        self.clear_dim_cache()
+        self.clear_dim_cache(include_vals=True)
         if validate:
-            self.clear_val_cache()
+            self.clear_val_cache(include_vals=False)
 
         if reset:
             self._bins: list[list[int]] = [[1] * self._obj.ndim]
@@ -399,7 +404,7 @@ def set_bins(
     @QtCore.Slot(int, int, int, bool, result=list)
     def set_bin(
         self, cursor: int, axis: int, value: int, update: bool = True
-    ) -> list[int | None]:
+    ) -> list[int]:
         if value is None:
             return []
         if int(value) != value:
@@ -703,7 +708,7 @@ def _bin_along_axis(
                 (slice(None),) * axis_val + (self._bin_slice(cursor, axis),)
             ].squeeze(axis=axis_val)
         else:
-            return _fast_nanmean_skipcheck(
+            return fast_nanmean_skipcheck(
                 self.data_vals_T[
                     (slice(None),) * axis_val + (self._bin_slice(cursor, axis),)
                 ],
@@ -725,6 +730,6 @@ def _bin_along_multiaxis(
             )
         ]
         if any(self.get_binned(cursor)):
-            return _fast_nanmean_skipcheck(selected, axis=axis)
+            return fast_nanmean_skipcheck(selected, axis=axis)
         else:
             return selected
diff --git a/src/erlab/interactive/utilities.py b/src/erlab/interactive/utilities.py
index 6296bcda..81890617 100644
--- a/src/erlab/interactive/utilities.py
+++ b/src/erlab/interactive/utilities.py
@@ -6,7 +6,7 @@
 import sys
 import types
 import warnings
-from typing import TYPE_CHECKING, Any, Literal, cast
+from typing import TYPE_CHECKING, Any, Literal, cast, no_type_check
 
 import numpy as np
 import numpy.typing as npt
@@ -20,6 +20,8 @@
 if TYPE_CHECKING:
     from collections.abc import Mapping
 
+    from pyqtgraph.GraphicsScene.mouseEvents import MouseDragEvent
+
 __all__ = [
     "AnalysisWidgetBase",
     "AnalysisWindow",
@@ -35,7 +37,7 @@
 ]
 
 
-def parse_data(data):
+def parse_data(data) -> xr.DataArray:
     if isinstance(data, xr.Dataset):
         raise TypeError(
             "input argument data must be a xarray.DataArray or a "
@@ -60,7 +62,7 @@ def array_rect(data):
 
 
 def copy_to_clipboard(content: str | list[str]) -> str:
-    """Convenience function for clipboard handling.
+    """Copy content to the clipboard.
 
     Parameters
     ----------
@@ -71,7 +73,6 @@ def copy_to_clipboard(content: str | list[str]) -> str:
     -------
     str
         The copied content.
-
     """
     if isinstance(content, list):
         content = "\n".join(content)
@@ -89,7 +90,7 @@ def process_arg(arg):
 
 
 def gen_single_function_code(funcname: str, *args: tuple, **kwargs):
-    """Generates the string for a Python function call.
+    """Generate the string for a Python function call.
 
     The first argument is the name of the function, and subsequent arguments are
     passed as positional arguments. Keyword arguments are also supported.
@@ -137,7 +138,7 @@ def gen_single_function_code(funcname: str, *args: tuple, **kwargs):
 
 
 def gen_function_code(copy: bool = True, **kwargs):
-    r"""Copies the Python code for function calls to the clipboard.
+    r"""Copy the Python code for function calls to the clipboard.
 
     The result can be copied to your clipboard in a form that can be pasted into an
     interactive Python session or Jupyter notebook cell.
@@ -228,7 +229,7 @@ def __init__(
         self.setDecimals(decimals)
 
         self._value = value
-        self._lastvalue = None
+        self._lastvalue: float | None = None
         self._min = -np.inf
         self._max = np.inf
         self._step = 1 if self._only_int else 0.01
@@ -399,7 +400,9 @@ def setValue(self, val):
         self._lastvalue, self._value = self._value, val
 
         self.valueChanged.emit(self.value())
-        self.lineEdit().setText(self.text())
+        line = self.lineEdit()
+        if line is not None:
+            line.setText(self.text())
         self.textChanged.emit(self.text())
 
     def fixup(self, input):
@@ -424,12 +427,18 @@ def validate(self, strn, pos):
         return (ret, strn, pos)
 
     def editingFinishedEvent(self):
-        self.setValue(self.valueFromText(self.lineEdit().text()))
+        line = self.lineEdit()
+        if line is not None:
+            self.setValue(self.valueFromText(line.text()))
 
     def keyPressEvent(self, evt):
-        if evt == QtGui.QKeySequence.StandardKey.Copy:
-            if (not evt.isAutoRepeat()) and self.lineEdit().hasSelectedText():
-                copy_to_clipboard(self.lineEdit().selectedText())
+        line = self.lineEdit()
+        if line is not None:
+            if evt == QtGui.QKeySequence.StandardKey.Copy:
+                if (not evt.isAutoRepeat()) and line.hasSelectedText():
+                    copy_to_clipboard(line.selectedText())
+            else:
+                super().keyPressEvent(evt)
         else:
             super().keyPressEvent(evt)
 
@@ -438,7 +447,8 @@ def _updateHeight(self):
             self.setFixedHeight(QtGui.QFontMetrics(self.font()).height() + 3)
 
     def _get_offset(self):
-        spin = QtWidgets.QDoubleSpinBox(self, minimum=0, maximum=0)
+        spin = QtWidgets.QDoubleSpinBox(self)
+        spin.setRange(0, 0)
         w = (
             spin.minimumSizeHint().width()
             - QtGui.QFontMetrics(spin.font())
@@ -494,7 +504,7 @@ def tickStrings(self, values, scale, spacing):
         for v in values:
             vs = v * scale
             if abs(vs) < 0.001 or abs(vs) >= 10000:
-                vstr = "%g" % vs
+                vstr = f"{vs:g}"
             else:
                 vstr = ("%%0.%df" % places) % vs
             strings.append(vstr.replace("-", "−"))
@@ -510,14 +520,13 @@ def labelString(self):
                 #     r"10<sup>\1\2</sup>",
                 #     f"(×{1.0 / self.autoSIPrefixScale:.3G})",
                 # )#.replace("-", "−")
+                search = re.search(
+                    r"1E\+?(\-?)0?(\d?\d)",
+                    f"{1.0 / self.autoSIPrefixScale:.3G}",
+                )
+                if search is not None:
+                    units = "".join(search.groups())
 
-                try:
-                    units = "".join(
-                        re.search(
-                            r"1E\+?(\-?)0?(\d?\d)",
-                            f"{1.0 / self.autoSIPrefixScale:.3G}",
-                        ).groups()
-                    )
                     for k, v in zip(
                         ("0", "1", "2", "3", "4", "5", "6", "7", "8", "9", "-"),
                         ("⁰", "¹", "²", "³", "⁴", "⁵", "⁶", "⁷", "⁸", "⁹", "⁻"),
@@ -525,7 +534,7 @@ def labelString(self):
                     ):
                         units = units.replace(k, v)
                     units = f"10{units}"
-                except AttributeError:
+                else:
                     units = f"{1.0 / self.autoSIPrefixScale:.3G}"
                 units = f"(×{units})"
 
@@ -547,7 +556,7 @@ def setLabel(self, text=None, units=None, unitPrefix=None, **args):
         self.labelUnits = units or ""
         self.labelUnitPrefix = unitPrefix or ""
         if len(args) > 0:
-            self.labelStyle = args
+            self.labelStyle: dict = args
         # Account empty string and `None` for units and text
         visible = bool(text or units)
         if text == units == "":
@@ -679,9 +688,10 @@ def param_dict(self):
 
 
 class xImageItem(BetterImageItem):
-    """
-    :class:`pyqtgraph.ImageItem` with additional functionality, including
-    :class:`xarray.DataArray` support and auto limits based on histogram analysis.
+    """:class:`pyqtgraph.ImageItem` with additional functionality.
+
+    This class provides :class:`xarray.DataArray` support and auto limits based on
+    histogram analysis.
 
     Parameters
     ----------
@@ -700,7 +710,7 @@ class xImageItem(BetterImageItem):
 
     def __init__(self, image: npt.NDArray | None = None, **kwargs):
         super().__init__(image, **kwargs)
-        self.cut_tolerance = (30, 30)
+        self.cut_tolerance = [30, 30]
         self.data_array: None | xr.DataArray = None
 
     def set_cut_tolerance(self, cut_tolerance):
@@ -749,6 +759,7 @@ def setDataArray(self, data: xr.DataArray, update_labels: bool = True, **kargs):
         self.setImage(img, rect=rect, **kargs)
         self.data_array = data
 
+    @no_type_check
     def getMenu(self):
         if self.menu is None:
             if not self.removable:
@@ -808,7 +819,6 @@ class ParameterGroup(QtWidgets.QGroupBox):
 
     Examples
     --------
-
     >>> ParameterGroup(
         **{
             "a": QtWidgets.QDoubleSpinBox(range=(0, 1), singleStep=0.01, value=0.2),
@@ -911,7 +921,7 @@ def getParameterWidget(
         **kwargs,
     ):
         """
-        Initializes the :class:`PySide6.QtWidgets.QWidget` corresponding to ``qwtype``.
+        Initialize the :class:`PySide6.QtWidgets.QWidget` corresponding to ``qwtype``.
 
         Parameters
         ----------
@@ -1001,9 +1011,12 @@ def getParameterWidget(
 
     def set_values(self, **kwargs):
         for k, v in kwargs.items():
-            self.widgets[k].blockSignals(True)
-            self.widgets[k].setValue(v)
-            self.widgets[k].blockSignals(False)
+            widget = self.widgets[k]
+            widget.blockSignals(True)
+            if hasattr(widget, "setValue"):
+                widget.setValue(v)
+            widget.blockSignals(False)
+
         self.sigParameterChanged.emit(kwargs)
 
     def widget_value(self, widget: str | QtWidgets.QWidget):
@@ -1148,6 +1161,7 @@ def roi_limits(self):
         x1, y1 = x0 + w, y0 + h
         return x0, y0, x1, y1
 
+    @no_type_check
     def update_pos(self):
         self.widgets["x0"].setMaximum(self.widgets["x1"].value())
         self.widgets["y0"].setMaximum(self.widgets["y1"].value())
@@ -1177,12 +1191,13 @@ def draw_mode(self, toggle):
             vb.setMouseMode(self._state_old)
             vb.rbScaleBox.setPen(pg.mkPen((255, 255, 100), width=1))
             vb.rbScaleBox.setBrush(pg.mkBrush(255, 255, 0, 100))
-            vb.setCursor(QtGui.QCursor(QtCore.Qt.ArrowCursor))
+
+            vb.setCursor(QtGui.QCursor(QtCore.Qt.CursorShape.ArrowCursor))
             return
 
-        vb.setCursor(QtGui.QCursor(QtCore.Qt.CrossCursor))
-        self._state_old = vb.state["mouseMode"]
-        self._drag_evt_old = vb.mouseDragEvent
+        vb.setCursor(QtGui.QCursor(QtCore.Qt.CursorShape.CrossCursor))
+        self._state_old: int = vb.state["mouseMode"]
+        self._drag_evt_old: MouseDragEvent = vb.mouseDragEvent
 
         vb.setMouseMode(pg.ViewBox.RectMode)
         vb.rbScaleBox.setPen(pg.mkPen((255, 255, 255), width=1))
@@ -1256,7 +1271,7 @@ def __init__(
         *args,
         **kwargs,
     ):
-        self.qapp = QtCore.QCoreApplication.instance()
+        self.qapp = cast(QtWidgets.QApplication, QtWidgets.QApplication.instance())
         if not self.qapp:
             self.qapp = QtWidgets.QApplication(sys.argv)
         self.qapp.setStyle("Fusion")
@@ -1273,10 +1288,12 @@ def __init__(
 
         self.controlgroup = QtWidgets.QWidget()
         if layout == "vertical":
-            self.layout = QtWidgets.QVBoxLayout(self._main)
-            self.controls = QtWidgets.QHBoxLayout(self.controlgroup)
+            layout = QtWidgets.QVBoxLayout(self._main)
+            self.controls: QtWidgets.QBoxLayout = QtWidgets.QHBoxLayout(
+                self.controlgroup
+            )
         elif layout == "horizontal":
-            self.layout = QtWidgets.QHBoxLayout(self._main)
+            layout = QtWidgets.QHBoxLayout(self._main)
             self.controls = QtWidgets.QVBoxLayout(self.controlgroup)
         else:
             raise ValueError("Layout must be 'vertical' or 'horizontal'.")
@@ -1302,14 +1319,14 @@ def __init__(
             # "refresh_all",
         ]:
             setattr(self, n, getattr(self.aw, n))
-        self.layout.addWidget(self.aw)
-        self.layout.addWidget(self.controlgroup)
+        layout.addWidget(self.aw)
+        layout.addWidget(self.controlgroup)
 
-        self.layout.setContentsMargins(QtCore.QMargins(0, 0, 0, 0))
-        self.layout.setSpacing(0)
+        layout.setContentsMargins(QtCore.QMargins(0, 0, 0, 0))
+        layout.setSpacing(0)
         self.controlgroup.setContentsMargins(QtCore.QMargins(0, 0, 0, 0))
         if data_is_input:
-            self.set_input(data)
+            self.aw.set_input(data)
 
     def __post_init__(self, execute=None):
         self.show()
@@ -1340,7 +1357,7 @@ def closeEvent(self, event: QtGui.QCloseEvent | None) -> None:
 
 
 class AnalysisWidgetBase(pg.GraphicsLayoutWidget):
-    """AnalysisWidgetBase
+    """AnalysisWidgetBase.
 
     Parameters
     ----------
@@ -1372,7 +1389,7 @@ def __init__(
             raise ValueError("Orientation must be 'vertical' or 'horizontal'.")
         self.cut_to_data = cut_to_data
 
-        self.input = None
+        self.input: None | xr.DataArray = None
 
         self.initialize_layout(num_ax)
 
@@ -1422,10 +1439,10 @@ def setStretchFactor(self, i, factor):
     def set_input(self, data=None):
         if data is not None:
             self.input = parse_data(data)
-        self.images[0].setDataArray(
-            self.input,
-            cut_to_data=self.cut_to_data in ("in", "both"),
-        )
+            self.images[0].setDataArray(
+                self.input,
+                cut_to_data=self.cut_to_data in ("in", "both"),
+            )
 
     def add_roi(self, i):
         self.roi[i] = pg.ROI(
@@ -1447,7 +1464,8 @@ def add_roi(self, i):
 
 class ComparisonWidget(AnalysisWidgetBase):
     def __init__(self, *args, **kwargs):
-        super().__init__(*args, num_ax=2, **kwargs)
+        kwargs["num_ax"] = 2
+        super().__init__(*args, **kwargs)
         self.prefunc = lambda x: x
         self.mainfunc = lambda x: x
         self.prefunc_only_values = False
@@ -1457,8 +1475,10 @@ def __init__(self, *args, **kwargs):
 
     def call_prefunc(self, x):
         if self.prefunc_only_values:
-            x = np.asarray(x)
-        return self.prefunc(x, **self.prefunc_kwargs)
+            xval = np.asarray(x)
+        else:
+            xval = x
+        return self.prefunc(xval, **self.prefunc_kwargs)
 
     def set_input(self, data=None):
         if data is not None:
@@ -1466,12 +1486,20 @@ def set_input(self, data=None):
 
         self.input = self.call_prefunc(self.input_)
         if self.prefunc_only_values:
+            if not isinstance(self.input, np.ndarray):
+                raise TypeError(
+                    "Pre-function must return a numpy array when `only_values` is True."
+                )
             self.images[0].setImage(
                 np.ascontiguousarray(self.input),
                 rect=array_rect(self.input_),
                 cut_to_data=self.cut_to_data in ("in", "both"),
             )
         else:
+            if not isinstance(self.input, xr.DataArray):
+                raise TypeError(
+                    "Pre-function must return a DataArray when `only_values` is False."
+                )
             self.images[0].setDataArray(
                 self.input,
                 cut_to_data=self.cut_to_data in ("in", "both"),
diff --git a/src/erlab/io/__init__.py b/src/erlab/io/__init__.py
index 41d7177f..c8ca5777 100644
--- a/src/erlab/io/__init__.py
+++ b/src/erlab/io/__init__.py
@@ -30,7 +30,6 @@
 
 Examples
 --------
-
 - View all registered loaders:
 
   >>> erlab.io.loaders
diff --git a/src/erlab/io/characterization/resistance.py b/src/erlab/io/characterization/resistance.py
index 69fdb543..ad83abb3 100644
--- a/src/erlab/io/characterization/resistance.py
+++ b/src/erlab/io/characterization/resistance.py
@@ -8,6 +8,7 @@
 import os
 import re
 from io import StringIO
+from pathlib import Path
 
 import numpy as np
 import pandas as pd
@@ -17,7 +18,7 @@
 
 
 def load_resistance_physlab(path: str, **kwargs) -> xr.Dataset:
-    """Loads resistance measurement acquired with physics lab III equipment.
+    """Load resistance measurement acquired with physics lab III equipment.
 
     Parameters
     ----------
@@ -50,7 +51,7 @@ def load_resistance_physlab(path: str, **kwargs) -> xr.Dataset:
 def _load_resistance_physlab_old(
     path: str, encoding: str = "windows-1252", **kwargs
 ) -> xr.Dataset:
-    """Loads resistance measurement acquired with physics lab III equipment.
+    """Load resistance measurement acquired with physics lab III equipment.
 
     Parameters
     ----------
@@ -86,14 +87,14 @@ def _load_resistance_physlab_old(
     >>> data.res.plot()
 
     """
-    with open(path, encoding=encoding) as file:
-        content = re.sub(
-            r"(e[-+]\d{3}) {2,3}(-?)",
-            "\\g<1>\\t \\g<2>",
-            file.read(),
-            count=0,
-            flags=re.MULTILINE,
-        )
+    content = re.sub(
+        r"(e[-+]\d{3}) {2,3}(-?)",
+        "\\g<1>\\t \\g<2>",
+        Path(path).read_text(encoding),
+        count=0,
+        flags=re.MULTILINE,
+    )
+
     content = content.replace("-1.#IO", "   nan")
     data = np.genfromtxt(
         StringIO(content),
diff --git a/src/erlab/io/dataloader.py b/src/erlab/io/dataloader.py
index 1000570c..60b1b524 100644
--- a/src/erlab/io/dataloader.py
+++ b/src/erlab/io/dataloader.py
@@ -1,4 +1,4 @@
-r"""Base functionality for implementing data loaders
+r"""Base functionality for implementing data loaders.
 
 This module provides a base class `LoaderBase` for implementing data loaders. Data
 loaders are plugins used to load data from various file formats. Each data loader that
@@ -48,15 +48,15 @@ def _is_monotonic(arr: npt.NDArray) -> np.bool_:
 
 
 class ValidationError(Exception):
-    """This exception is raised when the loaded data fails validation checks."""
+    """Raised when the loaded data fails validation checks."""
 
 
 class ValidationWarning(UserWarning):
-    """This warning is issued when the loaded data fails validation checks."""
+    """Issued when the loaded data fails validation checks."""
 
 
 class LoaderNotFoundError(Exception):
-    """This exception is raised when a loader is not found in the registry."""
+    """Raised when a loader is not found in the registry."""
 
     def __init__(self, key: str):
         super().__init__(f"Loader for name or alias {key} not found in the registry")
@@ -115,7 +115,11 @@ class LoaderBase:
 
     @property
     def name_map_reversed(self) -> dict[str, str]:
-        """A reversed version of the name_map dictionary."""
+        """A reversed version of the name_map dictionary.
+
+        This property is useful for mapping original names to new names.
+
+        """
         return self.reverse_mapping(self.name_map)
 
     @staticmethod
@@ -221,7 +225,6 @@ def formatter(cls, val: object):
 
         Examples
         --------
-
         >>> formatter(np.array([0.1, 0.15, 0.2]))
         '0.1→0.2 (0.05, 3)'
 
@@ -424,10 +427,11 @@ def summarize(
         display: bool = True,
         **kwargs,
     ) -> pandas.DataFrame | pandas.io.formats.style.Styler | None:
-        """
-        Takes a path to a directory and summarizes the data in the directory to a table
-        `pandas.DataFrame`, much like a log file. This is useful for quickly inspecting
-        the contents of a directory.
+        """Summarize the data in the given directory.
+
+        Takes a path to a directory and summarizes the data in the directory to a table,
+        much like a log file. This is useful for quickly inspecting the contents of a
+        directory.
 
         The dataframe is formatted using the style from :meth:`get_styler
         <erlab.io.dataloader.LoaderBase.get_styler>` and displayed in the IPython shell.
@@ -454,14 +458,20 @@ def summarize(
         Returns
         -------
         df : pandas.DataFrame or pandas.io.formats.style.Styler or None
-            Summary of the data in the directory. If `display` is `False`, the DataFrame
-            is returned. If `display` is `True` and the IPython shell is detected, the
-            summary will be displayed, and `None` will be returned. If `display` is
-            `True` but the IPython shell is not detected, the styler for the summary
-            DataFrame will be returned.
+            Summary of the data in the directory.
 
-        """
+            - If `display` is `False`, the summary DataFrame is returned.
+
+            - If `display` is `True` and the IPython shell is detected, the summary will
+              be displayed, and `None` will be returned.
+
+              * If `ipywidgets` is installed, an interactive widget will be returned
+                instead of `None`.
 
+            - If `display` is `True` but the IPython shell is not detected, the styler
+              for the summary DataFrame will be returned.
+
+        """
         if not os.path.isdir(data_dir):
             raise FileNotFoundError(f"Directory {data_dir} not found")
 
@@ -502,7 +512,7 @@ def summarize(
                     display(styled)  # type: ignore[misc]
 
                 if importlib.util.find_spec("ipywidgets"):
-                    display(self.isummarize(df=df))  # type: ignore[misc]
+                    return self._isummarize(df)
 
                 return None
 
@@ -539,6 +549,9 @@ def isummarize(self, df: pandas.DataFrame | None = None, **kwargs):
             kwargs["display"] = False
             df = cast(pandas.DataFrame, self.summarize(**kwargs))
 
+        self._isummarize(df)
+
+    def _isummarize(self, df: pandas.DataFrame):
         import matplotlib.pyplot as plt
         from ipywidgets import (
             HTML,
@@ -596,13 +609,18 @@ def _update_data(_, *, full: bool = False):
                         if n_scans > 1 and not full:
                             full_button.disabled = False
 
-                self._temp_data = self.load(path, single=not full)
+                out = self.load(path, single=not full)
+                if isinstance(out, xr.DataArray):
+                    self._temp_data = out
+                del out
 
                 data_info.value = _format_data_info(series)
+            if self._temp_data is None:
+                return
 
             if self._temp_data.ndim == 4:
                 # If the data is 4D, average over the last dimension, making it 3D
-                self._temp_data = self._temp_data.mean(self._temp_data.dims[-1])
+                self._temp_data = self._temp_data.mean(str(self._temp_data.dims[-1]))
 
             if self._temp_data.ndim == 3:
                 dim_sel.unobserve(_update_sliders, "value")
@@ -643,7 +661,7 @@ def _update_sliders(_):
             coord_sel.unobserve(_update_plot, "value")
 
             coord_sel.step = abs(scan_coords[1] - scan_coords[0])
-            coord_sel.max = 1e100  # To ensure max > min always
+            coord_sel.max = 1e100  # To ensure max > min before setting bounds
             coord_sel.min = scan_coords.min()
             coord_sel.max = scan_coords.max()
 
@@ -661,24 +679,25 @@ def _update_plot(_):
             out.clear_output(wait=True)
             with out:
                 plot_data.qplot(ax=plt.gca())
-                # Remove automatic title from xarray
-                plt.title("")
+                plt.title("")  # Remove automatically generated title
 
                 # Add line at Fermi level if the data is 2D and has an energy dimension
                 if plot_data.ndim == 2 and "eV" in plot_data.dims:
+                    # Check if binding
                     if plot_data["eV"].values[0] * plot_data["eV"].values[-1] < 0:
                         eplt.fermiline(
                             orientation="h" if plot_data.dims[0] == "eV" else "v"
                         )
-
                 show_inline_matplotlib_plots()
 
         def _next(_):
+            # Select next row
             idx = list(df.index).index(data_select.value)
             if idx + 1 < len(df.index):
                 data_select.value = list(df.index)[idx + 1]
 
         def _prev(_):
+            # Select previous row
             idx = list(df.index).index(data_select.value)
             if idx - 1 >= 0:
                 data_select.value = list(df.index)[idx - 1]
@@ -696,7 +715,7 @@ def _prev(_):
         if not self.always_single:
             buttons.append(full_button)
 
-        data_select = Select(options=list(df.index), value=next(iter(df.index)))
+        data_select = Select(options=list(df.index), value=next(iter(df.index)), rows=8)
         data_select.observe(_update_data, "value")
 
         data_info = HTML()
@@ -757,6 +776,9 @@ def identify(
         the list should contain only one file path and coordinates must be an empty
         dictionary.
 
+        The keys of the coordinates must be transformed to new names prior to returning
+        by using the mapping returned by the `name_map_reversed` property.
+
         Parameters
         ----------
         num
@@ -780,8 +802,11 @@ def identify(
     def infer_index(self, name: str) -> tuple[int | None, dict[str, Any]]:
         """Infer the index for the given file name.
 
-        This method takes a file name and tries to infer the scan index from it. If the
-        index can be inferred, it is returned; otherwise, `None` should be returned.
+        This method takes a file name with the path and extension stripped, and tries to
+        infer the scan index from it. If the index can be inferred, it is returned along
+        with additional keyword arguments that should be passed to `load`. If the index
+        is not found, `None` should be returned for the index, and an empty dictionary
+        for additional keyword arguments.
 
         Parameters
         ----------
@@ -805,9 +830,11 @@ def infer_index(self, name: str) -> tuple[int | None, dict[str, Any]]:
         raise NotImplementedError("method must be implemented in the subclass")
 
     def generate_summary(self, data_dir: str | os.PathLike) -> pandas.DataFrame:
-        """Takes a path to a directory and summarizes the data in the directory to a
-        pandas DataFrame, much like a log file. This is useful for quickly inspecting
-        the contents of a directory.
+        """Generate a dataframe summarizing the data in the given directory.
+
+        Takes a path to a directory and summarizes the data in the directory to a pandas
+        DataFrame, much like a log file. This is useful for quickly inspecting the
+        contents of a directory.
 
         Parameters
         ----------
@@ -1006,7 +1033,7 @@ def __new__(cls):
 
     @classmethod
     def instance(cls) -> Self:
-        """Returns the registry instance."""
+        """Return the registry instance."""
         return cls()
 
 
@@ -1113,7 +1140,6 @@ def loader_context(
           >>> dat_ssrl_2 = erlab.io.load(...)
 
         """
-
         if loader is None and data_dir is None:
             raise ValueError(
                 "At least one of loader or data_dir must be specified in the context"
diff --git a/src/erlab/io/exampledata.py b/src/erlab/io/exampledata.py
index c7ac81f6..85b3b480 100644
--- a/src/erlab/io/exampledata.py
+++ b/src/erlab/io/exampledata.py
@@ -2,7 +2,8 @@
 
 __all__ = ["generate_data", "generate_data_angles"]
 
-from collections.abc import Sequence
+from collections.abc import Hashable, Sequence
+from typing import cast
 
 import numpy as np
 import scipy.ndimage
@@ -224,7 +225,6 @@ def generate_data_angles(
         The generated data with coordinates for alpha, beta, and eV.
 
     """
-
     if isinstance(angrange, dict):
         alpha = np.linspace(*angrange["alpha"], shape[0])
         beta = np.linspace(*angrange["beta"], shape[1])
@@ -375,7 +375,7 @@ def generate_gold_edge(
         data[:] = rng.poisson(data).astype(float)
 
     data = erlab.analysis.image.gaussian_filter(
-        data, sigma={"eV": Eres, "alpha": angres}
+        data, sigma=cast(dict[Hashable, float], {"eV": Eres, "alpha": angres})
     )
 
     if noise:
diff --git a/src/erlab/io/igor.py b/src/erlab/io/igor.py
index 480aa7ce..acae569e 100644
--- a/src/erlab/io/igor.py
+++ b/src/erlab/io/igor.py
@@ -1,10 +1,10 @@
 import os
+from typing import Any
 
 import h5netcdf
 import igor2.binarywave
 import igor2.packed
 import igor2.record
-from typing import Any
 import numpy as np
 import xarray as xr
 
@@ -67,7 +67,7 @@ def load_experiment(
     recursive: bool = False,
     **kwargs,
 ) -> xr.Dataset:
-    """Loads waves from an igor experiment (`.pxp`) file.
+    """Load waves from an igor experiment (`.pxp`) file.
 
     Parameters
     ----------
@@ -145,6 +145,7 @@ def load_wave(
     data_dir
         The directory where the wave file is located. This parameter is only used if
         `wave` is a string or `PathLike` object. If `None`, `wave` must be a valid path.
+
     Returns
     -------
     xarray.DataArray
diff --git a/src/erlab/io/plugins/da30.py b/src/erlab/io/plugins/da30.py
index 6e734ac7..df252a34 100644
--- a/src/erlab/io/plugins/da30.py
+++ b/src/erlab/io/plugins/da30.py
@@ -1,6 +1,8 @@
-"""
-Base class for implementing loaders that can load data acquired with Scienta Omicron's
-DA30L analyzer using ``SES.exe``. Must be subclassed to implement the actual loading.
+"""Loader for Scienta Omicron DA30L analyzer with SES.
+
+Provides a base class for implementing loaders that can load data acquired with Scienta
+Omicron's DA30L analyzer using ``SES.exe``. Must be subclassed to implement the actual
+loading.
 """
 
 import configparser
@@ -8,6 +10,7 @@
 import tempfile
 import zipfile
 from typing import ClassVar
+
 import numpy as np
 import xarray as xr
 
diff --git a/src/erlab/io/utilities.py b/src/erlab/io/utilities.py
index 502f463c..f3ebbd70 100644
--- a/src/erlab/io/utilities.py
+++ b/src/erlab/io/utilities.py
@@ -13,12 +13,11 @@
 
 import numpy as np
 import numpy.typing as npt
-
 import xarray as xr
 
 
 def showfitsinfo(path: str | os.PathLike):
-    """Prints raw metadata from a ``.fits`` file.
+    """Print raw metadata from a ``.fits`` file.
 
     Parameters
     ----------
@@ -42,7 +41,7 @@ def get_files(
     contains: str | None = None,
     notcontains: str | None = None,
 ) -> list[str]:
-    """Returns a list of files in a directory with the given extensions.
+    """Return a list of files in a directory with the given extensions.
 
     Parameters
     ----------
@@ -61,7 +60,6 @@ def get_files(
         List of files in the directory.
 
     """
-
     files = []
 
     for f in os.listdir(directory):
@@ -167,7 +165,7 @@ def save_as_hdf5(
     igor_compat: bool = True,
     **kwargs,
 ):
-    """Saves data in ``HDF5`` format.
+    """Save data in ``HDF5`` format.
 
     Parameters
     ----------
@@ -220,7 +218,7 @@ def save_as_hdf5(
 
 
 def save_as_netcdf(data: xr.DataArray, filename: str | os.PathLike, **kwargs):
-    """Saves data in ``netCDF4`` format.
+    """Save data in ``netCDF4`` format.
 
     Discards invalid ``netCDF4`` attributes and produces a warning.
 
@@ -241,8 +239,3 @@ def save_as_netcdf(data: xr.DataArray, filename: str | os.PathLike, **kwargs):
         encoding={var: {"zlib": True, "complevel": 5} for var in data.coords},
         **kwargs,
     )
-
-
-def save_as_fits():
-    # TODO
-    raise NotImplementedError
diff --git a/src/erlab/lattice.py b/src/erlab/lattice.py
index 0c199f6e..5a361ec4 100644
--- a/src/erlab/lattice.py
+++ b/src/erlab/lattice.py
@@ -13,7 +13,7 @@
 
 
 def angle_between(v1: npt.NDArray[np.float64], v2: npt.NDArray[np.float64]) -> float:
-    """Returns the angle between two vectors.
+    """Return the angle between two vectors.
 
     Parameters
     ----------
@@ -31,7 +31,7 @@ def angle_between(v1: npt.NDArray[np.float64], v2: npt.NDArray[np.float64]) -> f
 def abc2avec(
     a: float, b: float, c: float, alpha: float, beta: float, gamma: float
 ) -> npt.NDArray:
-    """Constructs lattice vectors from lattice parameters."""
+    """Construct lattice vectors from lattice parameters."""
     alpha, beta, gamma = np.deg2rad(alpha), np.deg2rad(beta), np.deg2rad(gamma)
     sa, ca, sb, cb, cg = (
         np.sin(alpha),
@@ -55,7 +55,7 @@ def abc2avec(
 def avec2abc(
     avec: npt.NDArray[np.float64],
 ) -> tuple[np.floating, np.floating, np.floating, float, float, float]:
-    """Determines lattice parameters from lattice vectors."""
+    """Determine lattice parameters from lattice vectors."""
     a, b, c = (np.linalg.norm(x) for x in avec)
     alpha = angle_between(avec[1] / b, avec[2] / c)
     beta = angle_between(avec[2] / c, avec[0] / a)
@@ -64,7 +64,7 @@ def avec2abc(
 
 
 def to_reciprocal(avec: npt.NDArray[np.float64]) -> npt.NDArray[np.float64]:
-    """Constructs the reciprocal lattice vectors from real lattice vectors.
+    """Construct the reciprocal lattice vectors from real lattice vectors.
 
     Parameters
     ----------
@@ -81,7 +81,7 @@ def to_reciprocal(avec: npt.NDArray[np.float64]) -> npt.NDArray[np.float64]:
 
 
 def to_real(bvec: npt.NDArray[np.float64]) -> npt.NDArray[np.float64]:
-    """Constructs the real lattice vectors from reciprocal lattice vectors.
+    """Construct the real lattice vectors from reciprocal lattice vectors.
 
     Parameters
     ----------
diff --git a/src/erlab/parallel.py b/src/erlab/parallel.py
index 272aef85..a7ba8c6c 100644
--- a/src/erlab/parallel.py
+++ b/src/erlab/parallel.py
@@ -18,9 +18,7 @@
 
 @contextlib.contextmanager
 def joblib_progress(file=None, **kwargs):
-    """Context manager to patch joblib to report into tqdm progress bar given as
-    argument"""
-
+    """Patches joblib to report into a tqdm progress bar."""
     if file is None:
         file = sys.stdout
 
diff --git a/src/erlab/plotting/annotations.py b/src/erlab/plotting/annotations.py
index f7864318..49fd14ad 100644
--- a/src/erlab/plotting/annotations.py
+++ b/src/erlab/plotting/annotations.py
@@ -90,6 +90,56 @@
     "yocto",
 )  #: Names of the SI prefixes.
 
+PRETTY_NAMES: dict[str, tuple[str, str]] = {
+    "temperature": ("Temperature", "Temperature"),
+    "T": (r"\ensuremath{T}", r"$T$"),
+    "beta": (r"\ensuremath{\beta}", r"$\beta$"),
+    "theta": (r"\ensuremath{\theta}", r"$\theta$"),
+    "chi": (r"\ensuremath{\chi}", r"$\chi$"),
+    "alpha": (r"\ensuremath{\alpha}", r"$\alpha$"),
+    "psi": (r"\ensuremath{\psi}", r"$\psi$"),
+    "phi": (r"\ensuremath{\phi}", r"$\phi$"),
+    "xi": (r"\ensuremath{\xi}", r"$\xi$"),
+    "Eb": (r"\ensuremath{E}", r"$E$"),
+    "Ek": (r"\ensuremath{E_{\text{kin}}}", r"$E_{\text{kin}}$"),
+    "eV": (r"\ensuremath{E-E_F}", r"$E-E_F$"),
+    "kx": (r"\ensuremath{k_{x}}", r"$k_x$"),
+    "ky": (r"\ensuremath{k_{y}}", r"$k_y$"),
+    "kz": (r"\ensuremath{k_{z}}", r"$k_z$"),
+    "kp": (r"\ensuremath{k_{\parallel}}", r"$k_\parallel$"),
+    "hv": (r"\ensuremath{h\nu}", r"$h\nu$"),
+}
+"""Pretty names for labeling plots.
+
+The first element is for LaTeX, and the second is for plain text.
+
+"""
+
+PRETTY_UNITS: dict[str, tuple[str, str]] = {
+    "temperature": (r"K", r"K"),
+    "T": (r"K", r"K"),
+    "theta": (r"deg", r"deg"),
+    "beta": (r"deg", r"deg"),
+    "psi": (r"deg", r"deg"),
+    "chi": (r"deg", r"deg"),
+    "alpha": (r"deg", r"deg"),
+    "phi": (r"deg", r"deg"),
+    "xi": (r"deg", r"deg"),
+    "Eb": (r"eV", r"eV"),
+    "Ek": (r"eV", r"eV"),
+    "eV": (r"eV", r"eV"),
+    "hv": (r"eV", r"eV"),
+    "kx": (r"Å\ensuremath{{}^{-1}}", r"Å${}^{-1}$"),
+    "ky": (r"Å\ensuremath{{}^{-1}}", r"Å${}^{-1}$"),
+    "kz": (r"Å\ensuremath{{}^{-1}}", r"Å${}^{-1}$"),
+    "kp": (r"Å\ensuremath{{}^{-1}}", r"Å${}^{-1}$"),
+}
+"""Pretty units for labeling plots.
+
+The first element is for LaTeX, and the second is for plain text.
+
+"""
+
 
 def _alph_label(val, prefix, suffix, numeric, capital):
     """Generate labels from string or integer."""
@@ -123,7 +173,7 @@ def _unit_from_label(label: str) -> str | None:
 
 
 def get_si_str(si: int) -> str:
-    """Returns the SI prefix string to be plotted by :mod:`matplotlib`.
+    """Return the SI prefix string to be plotted by :mod:`matplotlib`.
 
     Parameters
     ----------
@@ -145,66 +195,33 @@ def get_si_str(si: int) -> str:
             raise ValueError("Invalid SI prefix.") from e
 
 
-def name_for_dim(dim_name, escaped=True):
-    name = {
-        "temperature": ("Temperature", "Temperature"),
-        "T": (r"\ensuremath{T}", r"$T$"),
-        "beta": (r"\ensuremath{\beta}", r"$\beta$"),
-        "theta": (r"\ensuremath{\theta}", r"$\theta$"),
-        "chi": (r"\ensuremath{\chi}", r"$\chi$"),
-        "alpha": (r"\ensuremath{\alpha}", r"$\alpha$"),
-        "psi": (r"\ensuremath{\psi}", r"$\psi$"),
-        "phi": (r"\ensuremath{\phi}", r"$\phi$"),
-        "xi": (r"\ensuremath{\xi}", r"$\xi$"),
-        "Eb": (r"\ensuremath{E}", r"$E$"),
-        "eV": (r"\ensuremath{E-E_F}", r"$E-E_F$"),
-        "kx": (r"\ensuremath{k_{x}}", r"$k_x$"),
-        "ky": (r"\ensuremath{k_{y}}", r"$k_y$"),
-        "kz": (r"\ensuremath{k_{z}}", r"$k_z$"),
-        "kp": (r"\ensuremath{k_{\parallel}}", r"$k_\parallel$"),
-        "hv": (r"\ensuremath{h\nu}", r"$h\nu$"),
-    }.get(dim_name)
+def name_for_dim(dim_name: str, escaped: bool = True) -> str:
+    names: tuple[str, str] | None = PRETTY_NAMES.get(dim_name)
 
-    if name is None:
+    if names is None:
         name = dim_name
     else:
-        name = name[0] if plt.rcParams["text.usetex"] else name[1]
+        name = names[0] if plt.rcParams["text.usetex"] else names[1]
 
     if not escaped:
         name = name.replace("$", "")
     return name
 
 
-def unit_for_dim(dim_name, deg2rad=False):
-    unit = {
-        "temperature": (r"K", r"K"),
-        "T": (r"K", r"K"),
-        "theta": (r"deg", r"deg"),
-        "beta": (r"deg", r"deg"),
-        "psi": (r"deg", r"deg"),
-        "chi": (r"deg", r"deg"),
-        "alpha": (r"deg", r"deg"),
-        "phi": (r"deg", r"deg"),
-        "xi": (r"deg", r"deg"),
-        "Eb": (r"eV", r"eV"),
-        "eV": (r"eV", r"eV"),
-        "hv": (r"eV", r"eV"),
-        "kx": (r"Å\ensuremath{{}^{-1}}", r"Å${}^{-1}$"),
-        "ky": (r"Å\ensuremath{{}^{-1}}", r"Å${}^{-1}$"),
-        "kz": (r"Å\ensuremath{{}^{-1}}", r"Å${}^{-1}$"),
-        "kp": (r"Å\ensuremath{{}^{-1}}", r"Å${}^{-1}$"),
-    }.get(dim_name)
+def unit_for_dim(dim_name: str, deg2rad: bool = False) -> str:
+    units: tuple[str, str] | None = PRETTY_UNITS.get(dim_name)
 
-    if unit is None:
+    if units is None:
         unit = ""
     else:
-        unit = unit[0] if plt.rcParams["text.usetex"] else unit[1]
+        unit = units[0] if plt.rcParams["text.usetex"] else units[1]
+
     if deg2rad:
         unit = unit.replace("deg", "rad")
     return unit
 
 
-def label_for_dim(dim_name, deg2rad=False, escaped=True):
+def label_for_dim(dim_name: str, deg2rad: bool = False, escaped: bool = True) -> str:
     name = name_for_dim(dim_name, escaped=escaped)
     unit = unit_for_dim(dim_name, deg2rad=deg2rad)
     if unit == "":
@@ -213,15 +230,16 @@ def label_for_dim(dim_name, deg2rad=False, escaped=True):
         return f"{name} ({unit})"
 
 
-def parse_special_point(name):
+def parse_special_point(name: str) -> str:
     special_points = {"G": r"\Gamma", "D": r"\Delta"}
-    try:
+
+    if name in special_points.keys():
         return special_points[name]
-    except KeyError:
+    else:
         return name
 
 
-def parse_point_labels(name: str, roman=True, bar=False):
+def parse_point_labels(name: str, roman: bool = True, bar: bool = False) -> str:
     name = parse_special_point(name)
 
     if name.endswith("*"):
@@ -463,7 +481,6 @@ def label_subplots(
         documentation for a list of all possible arguments.
 
     """
-
     kwargs["fontweight"] = fontweight
     if plt.rcParams["text.usetex"] & (fontweight == "bold"):
         prefix = "\\textbf{" + prefix
@@ -583,7 +600,6 @@ def label_subplots_nature(
         documentation for a list of all possible arguments.
 
     """
-
     kwargs["fontweight"] = fontweight
     if plt.rcParams["text.usetex"] & (fontweight == "bold"):
         prefix = "\\textbf{" + prefix
@@ -764,9 +780,10 @@ def mark_points_y(pts, labels, roman=True, bar=False, ax=None):
         label_ax.set_frame_on(False)
 
 
-# TODO: fix format using name_for_dim and unit_for_dim
 def plot_hv_text(ax, val, x=0.025, y=0.975, **kwargs):
-    s = "$h\\nu=" + str(val) + "$~eV"
+    name = name_for_dim("hv", escaped=False)
+    unit = unit_for_dim("hv")
+    s = f"${name}={val}$ {unit}"
     ax.text(
         x,
         y,
@@ -780,7 +797,9 @@ def plot_hv_text(ax, val, x=0.025, y=0.975, **kwargs):
 
 
 def plot_hv_text_right(ax, val, x=1 - 0.025, y=0.975, **kwargs):
-    s = "$h\\nu=" + str(val) + "$~eV"
+    name = name_for_dim("hv", escaped=False)
+    unit = unit_for_dim("hv")
+    s = f"${name}={val}$ {unit}"
     ax.text(
         x,
         y,
@@ -839,7 +858,12 @@ def __init__(self, si: int = 0, *args, **kwargs):
     def __call__(self, x, pos=None):
         self.orderOfMagnitude += self._si_exponent
 
-        sigfigs = int(re.match(r".*{%1.(\d+)f}", self.format).group(1))
+        match_format = re.match(r".*{%1.(\d+)f}", self.format)
+        if match_format is None:
+            # Match failed, may be due to changes in matplotlib
+            raise RuntimeError("Failed to match format string. Please report this bug")
+
+        sigfigs = int(match_format.group(1))
         self.format = self.format.replace(
             f"%1.{sigfigs}f", f"%1.{max(0, sigfigs + self._si_exponent)}f"
         )
@@ -859,8 +883,9 @@ def scale_units(
 ):
     """Rescales ticks and adds an SI prefix to the axis label.
 
-    Useful when you want to rescale the ticks. For example, when plotting a cut from a
-    low pass energy scan, you might want to convert the energy units from eV to meV.
+    Useful when you want to rescale the ticks without actually rescaling the data. For
+    example, when plotting a cut from a low pass energy scan, you might want to convert
+    the energy units from eV to meV.
 
     Using this function on an axis where the major locator is not the default formatter
     `matplotlib.ticker.ScalarFormatter` will result in undefined behavior.
@@ -948,7 +973,7 @@ def sizebar(
     frameon: bool = False,
     **kwargs,
 ):
-    """
+    """Add a size bar to an axes.
 
     Parameters
     ----------
@@ -981,7 +1006,6 @@ def sizebar(
         `mpl_toolkits.axes_grid1.anchored_artists.AnchoredSizeBar`.
 
     """
-
     size = value / resolution
     unit = get_si_str(si) + unit
     value = np.around(value / 10**si, decimals=decimals)
diff --git a/src/erlab/plotting/atoms.py b/src/erlab/plotting/atoms.py
index b187d959..50ceb073 100644
--- a/src/erlab/plotting/atoms.py
+++ b/src/erlab/plotting/atoms.py
@@ -1,4 +1,5 @@
-"""
+"""Plot atoms.
+
 Classes and functions for plotting atoms and bonds in a crystal structure using
 matplotlib's 3D plotting capabilities.
 
@@ -26,7 +27,7 @@
 
 
 def sunflower_sphere(n: int = 100) -> npt.NDArray[np.float64]:
-    """Returns n points on a sphere using the sunflower algorithm."""
+    """Return n points on a sphere using the sunflower algorithm."""
     indices = np.arange(0, n, dtype=float) + 0.5
 
     phi = np.arccos(1 - 2 * indices / n)
@@ -97,8 +98,11 @@ def projected_length_pos(ax: mpl_toolkits.mplot3d.Axes3D, length, position):
 
 
 def _zalpha(colors, zs):
-    """Modified from `mpl_toolkits.mplot3d.art3d` to modifies the brightness of the
-    color based on depth, rather than setting the transparency."""
+    """Set alpha based on zdepth.
+
+    Modified from `mpl_toolkits.mplot3d.art3d` to modifies the brightness of the
+    color based on depth, rather than setting the transparency.
+    """
     if len(colors) == 0 or len(zs) == 0:
         return np.zeros((0, 4))
     norm = plt.Normalize(min(zs), max(zs))
@@ -191,9 +195,7 @@ def draw(self, renderer):
 
 
 class Bond3DCollection(mpl_toolkits.mplot3d.art3d.Line3DCollection):
-    """
-    A subclass of `mpl_toolkits.mplot3d.art3d.Line3DCollection` representing a
-    collection of 3D bonds.
+    """A collection of 3D lines that represents bonds in a crystal structure.
 
     Parameters
     ----------
@@ -255,10 +257,11 @@ def __init__(
         mask: Callable | None = None,
         r_factor: float = 0.4,
     ):
-        """
-        A class that stores the information required to plot a three dimensional crystal
-        structure. The crystal can be repeated along the lattice vectors, and can be
-        masked. Bonds can be added between atoms.
+        """Properties of a crystal structure for plotting.
+
+        Stores the information required to plot a three dimensional crystal structure.
+        The crystal can be repeated along the lattice vectors, and can be masked. Bonds
+        can be added between atoms.
 
         Parameters
         ----------
@@ -471,7 +474,6 @@ def plot(
             plot the atoms.
 
         """
-
         if ax is None:
             ax = plt.gcf().add_subplot(projection="3d")
         ax = cast(mpl_toolkits.mplot3d.Axes3D, ax)
diff --git a/src/erlab/plotting/bz.py b/src/erlab/plotting/bz.py
index 0031bab7..852039a9 100644
--- a/src/erlab/plotting/bz.py
+++ b/src/erlab/plotting/bz.py
@@ -3,6 +3,7 @@
 __all__ = ["get_bz_edge", "plot_hex_bz"]
 
 import itertools
+from typing import Any
 
 import matplotlib.pyplot as plt
 import numpy as np
@@ -12,10 +13,10 @@
 
 from erlab.plotting.colors import axes_textcolor
 
-abbrv_kws = {
-    "facecolor": ["fc", "none"],
-    "linestyle": ["ls", "--"],
-    "linewidth": ["lw", 0.5],
+abbrv_kws: dict[str, tuple[str, Any]] = {
+    "facecolor": ("fc", "none"),
+    "linestyle": ("ls", "--"),
+    "linewidth": ("lw", 0.5),
 }
 
 
@@ -24,7 +25,7 @@ def get_bz_edge(
     reciprocal: bool = True,
     extend: tuple[int, ...] | None = None,
 ) -> tuple[npt.NDArray[np.float64], npt.NDArray[np.float64]]:
-    """Calculates the edge of the first Brillouin zone (BZ) from lattice vectors.
+    """Calculate the edge of the first Brillouin zone (BZ) from lattice vectors.
 
     Parameters
     ----------
@@ -117,9 +118,7 @@ def get_bz_edge(
 def plot_hex_bz(
     a=3.54, rotate=0.0, offset=(0.0, 0.0), reciprocal=True, ax=None, **kwargs
 ):
-    """
-    Plots a 2D hexagonal BZ overlay on the specified axes.
-    """
+    """Plot a 2D hexagonal BZ overlay on the specified axes."""
     kwargs.setdefault("zorder", 5)
     for k, v in abbrv_kws.items():
         kwargs[k] = kwargs.pop(k, kwargs.pop(*v))
diff --git a/src/erlab/plotting/colors.py b/src/erlab/plotting/colors.py
index 7035d4f7..137484db 100644
--- a/src/erlab/plotting/colors.py
+++ b/src/erlab/plotting/colors.py
@@ -60,34 +60,35 @@
 import matplotlib.transforms
 import numpy as np
 import numpy.typing as npt
-from matplotlib.typing import ColorType
+from matplotlib.typing import ColorType, RGBColorType
 
 
 class InversePowerNorm(matplotlib.colors.Normalize):
-    r"""
+    r"""Inverse power-law normalization.
+
     Linearly map a given value to the 0-1 range and then apply an inverse power-law
     normalization over that range.
 
-    For values :math:`x`, `matplotlib.colors.PowerNorm` calculates
-    :math:`x^\gamma`, whereas `InversePowerNorm` calculates :math:`1-x^{1/\gamma}`.
-    This provides higher contrast for values closer to ``vmin``.
+    For values :math:`x`, `matplotlib.colors.PowerNorm` calculates :math:`x^\gamma`,
+    whereas `InversePowerNorm` calculates :math:`1-x^{1/\gamma}`. This provides higher
+    contrast for values closer to ``vmin``.
 
     Parameters
     ----------
     gamma
         Power law normalization parameter. If equal to 1, the colormap is linear.
     vmin, vmax
-        If ``vmin`` and/or ``vmax`` is not given, they are initialized from the
-        minimum and maximum value, respectively, of the first input
-        processed; i.e., ``__call__(A)`` calls ``autoscale_None(A)``
+        If ``vmin`` and/or ``vmax`` is not given, they are initialized from the minimum
+        and maximum value, respectively, of the first input processed; i.e.,
+        ``__call__(A)`` calls ``autoscale_None(A)``
     clip
-        If ``True`` values falling outside the range ``[vmin, vmax]``,
-        are mapped to 0 or 1, whichever is closer, and masked values are
-        set to 1.  If ``False`` masked values remain masked.
+        If ``True`` values falling outside the range ``[vmin, vmax]``, are mapped to 0
+        or 1, whichever is closer, and masked values are set to 1.  If ``False`` masked
+        values remain masked.
 
-        Clipping silently defeats the purpose of setting the over, under,
-        and masked colors in a colormap, so it is likely to lead to
-        surprises; therefore the default is ``clip=False``.
+        Clipping silently defeats the purpose of setting the over, under, and masked
+        colors in a colormap, so it is likely to lead to surprises; therefore the
+        default is ``clip=False``.
 
     """
 
@@ -108,6 +109,7 @@ def __call__(self, value, clip=None):
         result, is_scalar = self.process_value(value)
 
         self.autoscale_None(result)
+
         gamma = self.gamma
         vmin, vmax = self.vmin, self.vmax
         if vmin > vmax:
@@ -271,9 +273,7 @@ def __init__(
         self._func_i = _diverging_powernorm_inv
 
     def __call__(self, value, clip=None):
-        """
-        Map value to the interval [0, 1]. The clip argument is unused.
-        """
+        """Map value to the interval [0, 1]. The clip argument is unused."""
         result, is_scalar = self.process_value(value)
 
         self.autoscale_None(result)
@@ -344,9 +344,7 @@ def __init__(
         self._func_i = _diverging_powernorm_inv
 
     def __call__(self, value, clip=None):
-        """
-        Map value to the interval [0, 1].
-        """
+        """Map value to the interval [0, 1]."""
         if clip is None:
             clip = self.clip
 
@@ -462,7 +460,7 @@ def __init__(
 def get_mappable(
     ax: matplotlib.axes.Axes, image_only: bool = False, silent: bool = False
 ) -> matplotlib.cm.ScalarMappable | None:
-    """Gets the `matplotlib.cm.ScalarMappable` from a given `matplotlib.axes.Axes`.
+    """Get the `matplotlib.cm.ScalarMappable` from a given `matplotlib.axes.Axes`.
 
     Parameters
     ----------
@@ -553,7 +551,7 @@ def proportional_colorbar(
     **kwargs,
 ) -> matplotlib.colorbar.Colorbar:
     """
-    Replaces the current colorbar or creates a new colorbar with proportional spacing.
+    Replace the current colorbar or creates a new colorbar with proportional spacing.
 
     The default behavior of colorbars in `matplotlib` does not support colors
     proportional to data in different norms. This function circumvents this behavior.
@@ -632,9 +630,6 @@ def proportional_colorbar(
     kwargs.setdefault("ticks", ticks)
     kwargs.setdefault("cmap", mappable.cmap)
     kwargs.setdefault("norm", mappable.norm)
-    kwargs.setdefault("pad", 0.05)
-    kwargs.setdefault("fraction", 0.05)
-    kwargs.setdefault("aspect", 25)
 
     cbar = plt.colorbar(
         mappable=mappable,
@@ -759,7 +754,7 @@ def __call__(self, ax, renderer):
             matplotlib.transforms.Bbox.from_bounds(*self._size_to_bounds(ax)),
             self._transAxes
             + matplotlib.transforms.ScaledTranslation(
-                *self.pads, ax.figure.dpi_scale_trans
+                self.pads[0], self.pads[1], ax.figure.dpi_scale_trans
             )
             - ax.figure.transSubfigure,
         )
@@ -827,11 +822,10 @@ def _gen_cax(ax, width=4.0, aspect=7.0, pad=3.0, horiz=False, **kwargs):
     return cax
 
 
-# TODO: fix colorbar size properly
 def nice_colorbar(
     ax: matplotlib.axes.Axes | Iterable[matplotlib.axes.Axes] | None = None,
     mappable: matplotlib.cm.ScalarMappable | None = None,
-    width: float = 5.0,
+    width: float = 8.0,
     aspect: float = 5.0,
     pad: float = 3.0,
     minmax: bool = False,
@@ -840,7 +834,7 @@ def nice_colorbar(
     ticklabels: Sequence[str] | None = None,
     **kwargs,
 ):
-    r"""Creates a colorbar with fixed width and aspect to ensure uniformity of plots.
+    r"""Create a colorbar with fixed width and aspect to ensure uniformity of plots.
 
     Parameters
     ----------
@@ -869,7 +863,6 @@ def nice_colorbar(
         The created colorbar.
 
     """
-
     is_horizontal = orientation == "horizontal"
 
     if ax is None:
@@ -900,9 +893,10 @@ def nice_colorbar(
                 ax = np.array(ax, dtype=object)
             bbox = matplotlib.transforms.Bbox.union(
                 [
-                    x.get_window_extent().transformed(
-                        x.figure.dpi_scale_trans.inverted()
-                    )
+                    x.get_position(original=True)
+                    .frozen()
+                    .transformed(x.figure.transFigure)
+                    .transformed(x.figure.dpi_scale_trans.inverted())
                     for x in ax.flat
                 ]
             )
@@ -912,7 +906,12 @@ def nice_colorbar(
             if fig is None:
                 raise RuntimeError("Axes is not attached to a figure")
 
-            bbox = ax.get_window_extent().transformed(fig.dpi_scale_trans.inverted())
+            bbox = (
+                ax.get_position(original=True)
+                .frozen()
+                .transformed(fig.transFigure)
+                .transformed(fig.dpi_scale_trans.inverted())
+            )
 
         if orientation == "horizontal":
             kwargs["anchor"] = (1, 1)
@@ -954,7 +953,7 @@ def nice_colorbar(
     return cbar
 
 
-def flatten_transparency(rgba: npt.NDArray, background: Sequence[float] | None = None):
+def flatten_transparency(rgba: npt.NDArray, background: RGBColorType | None = None):
     """
     Flatten the transparency of an RGBA image by blending it with a background color.
 
@@ -962,12 +961,15 @@ def flatten_transparency(rgba: npt.NDArray, background: Sequence[float] | None =
     ----------
     rgba
         The input RGBA image as a numpy array.
-    background
-        The background color to blend with. Defaults to white ``(1, 1, 1)``.
+    background : RGBColorType, optional
+        The background color to blend with. Defaults to white.
 
     """
     if background is None:
         background = (1, 1, 1)
+    else:
+        background = matplotlib.colors.to_rgb(background)
+
     original_shape = rgba.shape
     rgba = rgba.reshape(-1, 4)
     rgb = rgba[:, :-1]
@@ -1098,7 +1100,7 @@ def gen_2d_colormap(
         The normalization for the lightness axes.
     cnorm
         The normalization for the color axes.
-    background
+    background : ColorType, optional
         The background color. If `None`, it is set to white.
     N
         The number of levels in the colormap. Default is 256. The resulting colormap
@@ -1162,13 +1164,13 @@ def color_distance(c1: ColorType, c2: ColorType) -> float:
 
     Parameters
     ----------
-    c1, c2
+    c1, c2 : ColorType
         Color to calculate the distance between in any format that
         :func:`matplotlib.colors.to_rgb` can handle.
 
     Returns
     -------
-    float
+    distance : float
         The color distance between the two colors.
 
     Note
@@ -1196,7 +1198,7 @@ def close_to_white(c: ColorType) -> bool:
 
     Parameters
     ----------
-    c
+    c : ColorType
         Color in any format that :func:`matplotlib.colors.to_rgb` can handle.
 
     Returns
@@ -1228,10 +1230,9 @@ def prominent_color(im: matplotlib.image._ImageBase | matplotlib.collections.Qua
 def image_is_light(
     im: matplotlib.image._ImageBase | matplotlib.collections.QuadMesh,
 ) -> bool:
-    """
-    Determines if an image is *light* or *dark* by checking whether the prominent color
-    is closer to white than black.
+    """Determine if an image is *light* or *dark*.
 
+    Checks whether the prominent color is closer to white than black.
     """
     return close_to_white(prominent_color(im))
 
@@ -1243,13 +1244,13 @@ def axes_textcolor(
 
     Parameters
     ----------
-    ax
+    ax : matplotlib.axes.Axes
         The axes object for which the text color needs to be determined.
-    light
+    light : ColorType
         The *light* color, returned when :func:`image_is_light
         <erlab.plotting.colors.image_is_light>` returns `False`. Default is ``'w'``
         (white).
-    dark
+    dark : ColorType
         The *dark* color, returned when :func:`image_is_light
         <erlab.plotting.colors.image_is_light>` returns `True`. Default is ``'k'``
         (black).
diff --git a/src/erlab/plotting/general.py b/src/erlab/plotting/general.py
index d3cf85f1..0acd4b71 100644
--- a/src/erlab/plotting/general.py
+++ b/src/erlab/plotting/general.py
@@ -17,10 +17,12 @@
 
 import contextlib
 import copy
+import warnings
 from collections.abc import Iterable
 from typing import TYPE_CHECKING, Any, Literal, Union, cast
 
 import matplotlib
+import matplotlib.colorbar
 import matplotlib.colors
 import matplotlib.image
 import matplotlib.patches
@@ -43,6 +45,8 @@
 if TYPE_CHECKING:
     from collections.abc import Callable, Collection, Sequence
 
+    from matplotlib.typing import ColorType
+
 figure_width_ref = {
     "aps": [3.4, 7.0],
     "aip": [3.37, 6.69],
@@ -166,7 +170,7 @@ def __init__(
         self.needclear = False
 
     def clear(self, event):
-        """Internal event handler to clear the cursor."""
+        """Clear the cursor."""
         if self.ignore(event):
             return
         if self.useblit:
@@ -176,7 +180,7 @@ def clear(self, event):
         self.label.set_visible(False)
 
     def onmove(self, event):
-        """Internal event handler to draw the cursor when the mouse moves."""
+        """Draw the cursor when the mouse moves."""
         if self.ignore(event):
             return
         if not self.canvas.widgetlock.available(self):
@@ -256,27 +260,43 @@ def place_inset(
 
 
     """
-
     return _ez_inset(parent_axes, width, height, pad, loc, **kwargs)
 
 
-def array_extent(data: xr.DataArray) -> tuple[float, float, float, float]:
-    """
-    Gets the extent of a :class:`xarray.DataArray` to be used in
-    :func:`matplotlib.pyplot.imshow`.
+def array_extent(
+    darr: xr.DataArray, rtol: float = 1.0e-5, atol: float = 1.0e-8
+) -> tuple[float, float, float, float]:
+    """Get the extent of a :class:`xarray.DataArray`.
+
+    The extent can be used as the `extent` argument in :func:`matplotlib.pyplot.imshow`.
 
     Parameters
     ----------
-    data
+    darr
         A two-dimensional :class:`xarray.DataArray`.
+    rtol, atol
+        Tolerance for the coordinates to be considered evenly spaced. The default values
+        are consistent with `numpy.isclose`.
 
     Returns
     -------
     x0, x1, y0, y1 : float
 
     """
+    if darr.ndim != 2:
+        raise ValueError("Input array must be 2D")
+
+    data_coords = tuple(darr[dim].values for dim in darr.dims)
+    for dim, coord in zip(darr.dims, data_coords, strict=True):
+        dif = np.diff(coord)
+        if not np.allclose(dif, dif[0], rtol=rtol, atol=atol):
+            warnings.warn(
+                f"Coordinates for {dim} are not evenly spaced, and the plot may not be "
+                "accurate. Use `DataArray.plot`, `xarray.plot.pcolormesh` or "
+                "`matplotlib.pyplot.pcolormesh` for non-evenly spaced data.",
+                stacklevel=2,
+            )
 
-    data_coords = tuple(data[dim].values for dim in data.dims)
     data_incs = tuple(coord[1] - coord[0] for coord in data_coords)
     data_lims = tuple((coord[0], coord[-1]) for coord in data_coords)
     y0, x0 = data_lims[0][0] - 0.5 * data_incs[0], data_lims[1][0] - 0.5 * data_incs[1]
@@ -301,14 +321,16 @@ def plot_array(
     rad2deg: bool | Iterable[str] = False,
     func: Callable | None = None,
     func_args: dict | None = None,
+    rtol: float = 1.0e-5,
+    atol: float = 1.0e-8,
     **improps,
 ) -> matplotlib.image.AxesImage:
-    """Plots a 2D :class:`xarray.DataArray` using :func:`matplotlib.pyplot.imshow`.
+    """Plot a 2D :class:`xarray.DataArray` using :func:`matplotlib.pyplot.imshow`.
 
     Parameters
     ----------
     arr
-        A two-dimensional :class:`xarray.DataArray` to be plotted.
+        A two-dimensional :class:`xarray.DataArray` with evenly spaced coordinates.
     ax
         The target :class:`matplotlib.axes.Axes`.
     colorbar
@@ -328,6 +350,10 @@ def plot_array(
         same shape as the input.
     func_args
         Keyword arguments passed onto `func`.
+    rtol, atol
+        By default, the input array is checked for evenly spaced coordinates. ``rtol``
+        and ``atol`` are the tolerances for the coordinates to be considered evenly
+        spaced. The default values are consistent with `numpy.isclose`.
     **improps
         Keyword arguments passed onto :func:`matplotlib.axes.Axes.imshow`.
 
@@ -379,7 +405,7 @@ def plot_array(
 
     improps_default = {
         "interpolation": "none",
-        "extent": array_extent(arr),
+        "extent": array_extent(arr, rtol, atol),
         "aspect": "auto",
         "origin": "lower",
         "rasterized": True,
@@ -424,14 +450,84 @@ def plot_array_2d(
     cmap: matplotlib.colors.Colormap | str | None = None,
     lnorm: matplotlib.colors.Normalize | None = None,
     cnorm: matplotlib.colors.Normalize | None = None,
-    background: Any = None,
+    background: ColorType | None = None,
     colorbar: bool = True,
     cax: matplotlib.axes.Axes | None = None,
     colorbar_kw: dict | None = None,
     imshow_kw: dict | None = None,
     N: int = 256,
+    rtol: float = 1.0e-5,
+    atol: float = 1.0e-8,
     **indexers_kwargs,
-):
+) -> tuple[matplotlib.image.AxesImage, matplotlib.colorbar.Colorbar | None]:
+    """
+    Plot a 2D array with associated color array.
+
+    The lightness array represents the intensity values, while the color array
+    represents some other property. The arrays must have the same shape.
+
+    Parameters
+    ----------
+    larr
+        The 2D array representing the lightness values.
+    carr
+        The 2D array representing the color values.
+    ax
+        The axes on which to plot the array. If None, the current axes will be used.
+    normalize_with_larr
+        Whether to normalize the color array with the lightness array. Default is False.
+    xlim
+        The x-axis limits for the plot. If a float, it represents the symmetric limits
+        around 0. If a tuple, it represents the lower and upper limits. If None, the
+        limits are determined from the data.
+    ylim
+        The y-axis limits for the plot. If a float, it represents the symmetric limits
+        around 0. If a tuple, it represents the lower and upper limits. If None, the
+        limits are determined from the data.
+    cmap
+        The colormap to use for the color array. If None, a linear segmented colormap
+        consisting of blue, black, and red is used.
+    lnorm
+        The normalization object for the lightness array.
+    cnorm
+        The normalization object for the color array.
+    background
+        The background color to use for the plot. If None, white is used.
+    colorbar
+        Whether to create a colorbar. Default is `True`.
+    cax
+        The axes on which to create the colorbar if `colorbar` is `True`. If None, a new
+        axes will be created for the colorbar.
+    colorbar_kw
+        Additional keyword arguments to pass to `matplotlib.pyplot.colorbar`.
+    imshow_kw
+        Additional keyword arguments to pass to `matplotlib.pyplot.imshow`.
+    N
+        The number of levels in the colormap. Default is 256.
+    rtol, atol
+        By default, the input array is checked for evenly spaced coordinates. ``rtol``
+        and ``atol`` are the tolerances for the coordinates to be considered evenly
+        spaced. The default values are consistent with `numpy.isclose`.
+    **indexers_kwargs : dict
+        Additional keyword arguments to pass to `qsel` to select the data to plot. Note
+        that the resulting data after the selection must be 2D.
+
+    Returns
+    -------
+    im : matplotlib.image.AxesImage
+        The plotted image.
+    cb : matplotlib.colorbar.Colorbar or None
+        The colorbar associated with the plot. If `colorbar` is False, None is returned.
+
+    Example
+    -------
+    >>> import erlab.plotting.erplot as eplt
+    >>> import matplotlib.pyplot as plt
+    >>> import xarray as xr
+    >>> larr = xr.DataArray([[1, 2, 3], [4, 5, 6]])
+    >>> carr = xr.DataArray([[0.1, 0.2, 0.3], [0.4, 0.5, 0.6]])
+    >>> eplt.plot_array_2d(larr, carr)
+    """
     if lnorm is None:
         lnorm = plt.Normalize()
     else:
@@ -457,6 +553,10 @@ def plot_array_2d(
 
     larr = larr.qsel(**indexers_kwargs).copy(deep=True)
     carr = carr.qsel(**indexers_kwargs).copy(deep=True)
+
+    if larr.shape != carr.shape:
+        raise ValueError("Lightness and color array must have the same shape")
+
     sel_kw = {}
 
     if xlim is not None:
@@ -517,7 +617,7 @@ def plot_array_2d(
             aspect="auto",
         )
 
-    im = ax.imshow(img, extent=array_extent(larr), **imshow_kw)
+    im = ax.imshow(img, extent=array_extent(larr, rtol, atol), **imshow_kw)
     ax.set_xlabel(str(larr.dims[0]))
     ax.set_ylabel(str(larr.dims[1]))
     fancy_labels(ax)
@@ -539,7 +639,7 @@ def gradient_fill(
     ax: matplotlib.axes.Axes | None = None,
     **kwargs,
 ) -> matplotlib.image.AxesImage:
-    """Applies a gradient fill to a line plot.
+    """Apply a gradient fill to a line plot.
 
     Parameters
     ----------
@@ -1001,7 +1101,7 @@ def fermiline(
     orientation: Literal["h", "v"] = "h",
     **kwargs,
 ) -> matplotlib.lines.Line2D | MultipleLine2D:
-    """Plots a constant energy line to denote the Fermi level.
+    """Plot a constant energy line to denote the Fermi level.
 
     Parameters
     ----------


Test patch:
diff --git a/tests/accessors/test_kspace.py b/tests/accessors/test_kspace.py
index 5451f7f5..13ba792b 100644
--- a/tests/accessors/test_kspace.py
+++ b/tests/accessors/test_kspace.py
@@ -1,8 +1,8 @@
 import pytest
 import xarray
 import xarray.testing
-from erlab.io.exampledata import generate_data_angles
 from erlab.analysis.kspace import AxesConfiguration
+from erlab.io.exampledata import generate_data_angles
 
 
 @pytest.fixture()
diff --git a/tests/analysis/test_fastbinning.py b/tests/analysis/test_fastbinning.py
index da533255..047e6f4b 100644
--- a/tests/analysis/test_fastbinning.py
+++ b/tests/analysis/test_fastbinning.py
@@ -1,15 +1,25 @@
 import numpy as np
-from erlab.interactive.imagetool.fastbinning import fast_nanmean, nanmean_funcs
+from erlab.interactive.imagetool.fastbinning import NANMEAN_FUNCS, fast_nanmean
 
 
 def test_fast_nanmean():
-    for nd, funcs in nanmean_funcs.items():
+    for nd, funcs in NANMEAN_FUNCS.items():
         x = np.random.RandomState(42).randn(*((10,) * nd))
+        x64 = x.astype(np.float64)
+        x32 = x.astype(np.float32)
         for axis, func in funcs.items():
             if isinstance(axis, frozenset):
                 axis = tuple(axis)
-            if not np.allclose(np.nanmean(x, axis), fast_nanmean(x, axis)):
+            if not np.allclose(np.nanmean(x64, axis), fast_nanmean(x64, axis)):
                 raise AssertionError(
-                    f"fast_nanmean failed for {nd}D array "
+                    f"fast_nanmean failed for {nd}D float64 array "
+                    f"with axis {axis} using {func}."
+                )
+            if not np.allclose(
+                np.nanmean(x32.astype(np.float64), axis).astype(np.float32),
+                fast_nanmean(x32, axis),
+            ):
+                raise AssertionError(
+                    f"fast_nanmean failed for {nd}D float32 array "
                     f"with axis {axis} using {func}."
                 )
diff --git a/tests/analysis/test_fit_functions_general.py b/tests/analysis/test_fit_functions_general.py
index e09075fa..c8d8d596 100644
--- a/tests/analysis/test_fit_functions_general.py
+++ b/tests/analysis/test_fit_functions_general.py
@@ -55,7 +55,13 @@ def testfunc(x):
 
     resolution = 1.0
     expected = np.array(
-        [-6.61744490e-24, 5.87148334e00, 1.17429667e01, 1.76144500e01, 2.34859334e01]
+        [
+            -6.61744490e-24,
+            5.87148334e00,
+            1.17429667e01,
+            1.76144500e01,
+            2.34859334e01,
+        ]
     )
     result = do_convolve(x, testfunc, resolution)
     assert np.allclose(result, expected)
@@ -68,7 +74,13 @@ def testfunc(x):
 
     resolution = 0.5
     expected = np.array(
-        [0.00000000e00, 2.95132925e00, 3.61433592e-16, -2.95132925e00, -7.22867185e-16]
+        [
+            0.00000000e00,
+            2.95132925e00,
+            3.61433592e-16,
+            -2.95132925e00,
+            -7.22867185e-16,
+        ]
     )
     result = do_convolve(x, testfunc, resolution)
     assert np.allclose(result, expected)
diff --git a/tests/analysis/test_fit_models.py b/tests/analysis/test_fit_models.py
index 7ccca56b..bc2b1d34 100644
--- a/tests/analysis/test_fit_models.py
+++ b/tests/analysis/test_fit_models.py
@@ -1,8 +1,7 @@
+import erlab.analysis.fit.models as models
 import numpy as np
 import xarray as xr
 
-import erlab.analysis.fit.models as models
-
 
 def test_fermi_edge_model():
     # Create a test input DataArray
diff --git a/tests/analysis/test_interpolate.py b/tests/analysis/test_interpolate.py
index 6b21b8be..48c459c9 100644
--- a/tests/analysis/test_interpolate.py
+++ b/tests/analysis/test_interpolate.py
@@ -1,6 +1,7 @@
 import numpy as np
 import scipy.interpolate
-from erlab.analysis.interpolate import interpn
+import xarray as xr
+from erlab.analysis.interpolate import interpn, slice_along_path
 
 
 def value_func_1d(x):
@@ -58,3 +59,66 @@ def test_interpn_3d():
             points, values, point, method="linear", bounds_error=False
         ),
     )
+
+
+def test_slice_along_path():
+    x = np.linspace(0, 10, 11)
+    y = np.linspace(0, 10, 11)
+    z = np.linspace(0, 10, 11)
+    data = np.random.default_rng(1).random((11, 11, 11))
+    darr = xr.DataArray(data, coords={"x": x, "y": y, "z": z}, dims=["x", "y", "z"])
+    vertices = {"x": [0, 5, 10], "y": [0, 5, 10], "z": [0, 5, 10]}
+    interp = slice_along_path(darr, vertices, step_size=1.0)
+
+    np.testing.assert_allclose(
+        interp.values,
+        np.array(
+            [
+                0.51182162,
+                0.71780756,
+                0.59044693,
+                0.2016458,
+                0.40237743,
+                0.44487846,
+                0.65585751,
+                0.52695829,
+                0.77584909,
+                0.73380143,
+                0.65321432,
+                0.80464538,
+                0.68497213,
+                0.583326,
+                0.58202596,
+                0.60825019,
+                0.21265952,
+            ]
+        ),
+    )
+
+    np.testing.assert_allclose(interp.x.values, interp.y.values)
+    np.testing.assert_allclose(interp.x.values, interp.z.values)
+
+    np.testing.assert_allclose(
+        interp.x.values,
+        np.array(
+            [
+                0.0,
+                0.625,
+                1.25,
+                1.875,
+                2.5,
+                3.125,
+                3.75,
+                4.375,
+                5.0,
+                5.625,
+                6.25,
+                6.875,
+                7.5,
+                8.125,
+                8.75,
+                9.375,
+                10.0,
+            ]
+        ),
+    )
diff --git a/tests/io/test_dataloader.py b/tests/io/test_dataloader.py
index 7b438f4c..1f4314d9 100644
--- a/tests/io/test_dataloader.py
+++ b/tests/io/test_dataloader.py
@@ -64,7 +64,7 @@ def test_loader():
     for i, beta in enumerate(beta_coords):
         erlab.io.save_as_hdf5(
             make_data(beta=beta, temp=20.0, hv=50.0),
-            filename=f"{tmp_dir.name}/data_001_S{str(i+1).zfill(3)}.h5",
+            filename=f"{tmp_dir.name}/data_001_S{str(i + 1).zfill(3)}.h5",
             igor_compat=False,
         )