diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml
new file mode 100644
index 0000000000..c158ebf179
--- /dev/null
+++ b/.github/workflows/build.yml
@@ -0,0 +1,75 @@
+name: Build
+
+on:
+  workflow_dispatch:
+  schedule:
+    - cron:  '0 0 * * 0'
+  pull_request:
+    types: [labeled, ready_for_review, reopened]
+
+jobs:
+  build:
+    if: ${{ contains(github.event.pull_request.labels.*.name, 'run-packaging') || github.event_name == 'workflow_dispatch'  || github.event_name == 'schedule' }}
+    name: Build
+    runs-on: ubuntu-latest
+    env:
+      PYTHON_VERSION: '3.11'
+    steps:
+      - uses: actions/checkout@v4
+
+      - uses: actions/setup-python@v4
+        name: Install Python
+        with:
+          python-version: ${{ env.PYTHON_VERSION }}
+
+      - name: Display version
+        run: |
+          python --version
+          pip --version
+
+      - name: Install pypa/build
+        run: |
+          pip install build
+
+      - name: Build a binary wheel and a source tarball
+        run: |
+          python -m build
+
+      - name: Display content dist folder
+        run: |
+          ls -shR dist/
+
+      - uses: actions/upload-artifact@v3
+        with:
+          path: ./dist/*
+          name: dist
+
+  test:
+    name: Test Packaging
+    needs: build
+    runs-on: ubuntu-latest
+    env:
+      PYTHON_VERSION: '3.11'
+    steps:
+      - uses: actions/setup-python@v4
+        name: Install Python
+        with:
+          python-version: ${{ env.PYTHON_VERSION }}
+
+      - uses: actions/download-artifact@v3
+
+      - name: Display content working folder
+        run: |
+          ls -R
+
+      - name: Install distribution
+        run: |
+          pip install --pre --find-links dist hyperspy[all,tests]
+
+      - name: Pip list
+        run: |
+          pip list
+
+      - name: Test distribution
+        run: |
+          pytest --pyargs hyperspy --reruns 3 -n 2
diff --git a/doc/user_guide/install.rst b/doc/user_guide/install.rst
index 2ec25c8a3b..3ca35c4318 100644
--- a/doc/user_guide/install.rst
+++ b/doc/user_guide/install.rst
@@ -172,15 +172,14 @@ use:
 
 See the following list of selectors to select the installation of optional
 dependencies required by specific functionalities:
-
+* ``ipython`` for integration with the `ipython` terminal and parallel processing using `ipyparallel`,
 * ``learning`` for some machine learning features,
 * ``gui-jupyter`` to use the `Jupyter widgets <https://ipywidgets.readthedocs.io/en/stable/>`_
   GUI elements,
 * ``gui-traitsui`` to use the GUI elements based on `traitsui <https://docs.enthought.com/traitsui/>`_,
-* ``mrcz`` to read mrcz file,
-* ``speed`` to speed up some functionalities,
-* ``usid`` to read usid file,
+* ``speed`` install numba and numexpr to speed up some functionalities,
 * ``tests`` to install required libraries to run HyperSpy's unit tests,
+* ``coverage`` to coverage statistics when running the tests,
 * ``build-doc`` to install required libraries to build HyperSpy's documentation,
 * ``dev`` to install all the above,
 * ``all`` to install all the above except the development requirements
diff --git a/hyperspy/_components/bleasdale.py b/hyperspy/_components/bleasdale.py
index a7a6276ac0..4739782a2f 100644
--- a/hyperspy/_components/bleasdale.py
+++ b/hyperspy/_components/bleasdale.py
@@ -24,30 +24,34 @@ class Bleasdale(Expression):
 
     r"""Bleasdale function component.
 
-    Also called the Bleasdale-Nelder function. Originates from the description of the yield-density relationship in crop growth.
+    Also called the Bleasdale-Nelder function. Originates from
+    the description of the yield-density relationship in crop growth.
 
     .. math::
 
         f(x) = \left(a+b\cdot x\right)^{-1/c}
 
     Parameters
-    -----------
-        a : Float
-
-        b : Float
-
-        c : Float
-
-        **kwargs
-            Extra keyword arguments are passed to the
-            :py:class:`~._components.expression.Expression` component.
-
+    ----------
+    a : float, default=1.0
+        The value of Parameter a.
+    b : float, default=1.0
+        The value of Parameter b.
+    c : float, default=1.0
+        The value of Parameter c.        
+    **kwargs
+        Extra keyword arguments are passed to
+        :py:class:`~.api.model.components1D.Expression`.
+
+    Note
+    ----
     For :math:`(a+b\cdot x)\leq0`, the component will be set to 0.
+
     """
 
-    def __init__(self, a=1., b=1., c=1., module="numexpr", **kwargs):
+    def __init__(self, a=1., b=1., c=1., module=None, **kwargs):
         super().__init__(
-            expression="where((a + b * x) > 0, (a + b * x) ** (-1 / c), 0)",
+            expression="where((a + b * x) > 0, pow(a + b * x, -1 / c), 0)",
             name="Bleasdale",
             a=a,
             b=b,
@@ -58,6 +62,12 @@ def __init__(self, a=1., b=1., c=1., module="numexpr", **kwargs):
             linear_parameter_list=['b'],
             check_parameter_linearity=False,
             **kwargs)
+        module = self._whitelist['module'][1]
+        if module in ("numpy", "scipy"):
+            # Issue with calculating component at 0...
+            raise ValueError(
+                f"{module} is not supported for this component, use numexpr instead."
+                )
 
     def grad_a(self, x):
         """
diff --git a/hyperspy/_components/exponential.py b/hyperspy/_components/exponential.py
index 3d3094f6ac..dc85245b7f 100644
--- a/hyperspy/_components/exponential.py
+++ b/hyperspy/_components/exponential.py
@@ -52,7 +52,7 @@ class Exponential(Expression):
         :py:class:`~._components.expression.Expression` component.
     """
 
-    def __init__(self, A=1., tau=1., module="numexpr", **kwargs):
+    def __init__(self, A=1., tau=1., module=None, **kwargs):
         super().__init__(
             expression="A * exp(-x / tau)",
             name="Exponential",
diff --git a/hyperspy/_components/expression.py b/hyperspy/_components/expression.py
index 9018711e49..d8e6741e82 100644
--- a/hyperspy/_components/expression.py
+++ b/hyperspy/_components/expression.py
@@ -17,6 +17,7 @@
 # along with HyperSpy. If not, see <https://www.gnu.org/licenses/#GPL>.
 
 from functools import wraps
+import logging
 import numpy as np
 import sympy
 
@@ -26,6 +27,9 @@
 from hyperspy.docstrings.parameters import FUNCTION_ND_DOCSTRING
 
 
+_logger = logging.getLogger(__name__)
+
+
 _CLASS_DOC = \
     """%s component (created with Expression).
 
@@ -87,9 +91,10 @@ class docstring.
         applicable. It enables interative adjustment of the position of the
         component in the model. For 2D components, a tuple must be passed
         with the name of the two parameters e.g. `("x0", "y0")`.
-    module : {"numpy", "numexpr", "scipy"}, default "numpy"
+    module : {"numpy", "numexpr", "scipy", None}, default "numpy"
         Module used to evaluate the function. numexpr is often faster but
         it supports fewer functions and requires installing numexpr.
+        If None, the "numexpr" will be used if installed.
     add_rotation : bool, default False
         This is only relevant for 2D components. If `True` it automatically
         adds `rotation_angle` parameter.
@@ -114,10 +119,10 @@ class docstring.
         parameters.
     check_parameter_linearity : bool
         If `True`, automatically check if each parameter is linear and set
-        its corresponding attribute accordingly. If `False`, the default is to 
+        its corresponding attribute accordingly. If `False`, the default is to
         set all parameters, except for those who are specified in
         ``linear_parameter_list``.
-        
+
     **kwargs
         Keyword arguments can be used to initialise the value of the
         parameters.
@@ -162,6 +167,17 @@ def __init__(self, expression, name, position=None, module="numpy",
                  linear_parameter_list=None, check_parameter_linearity=True,
                  **kwargs):
 
+        if module is None:
+            try:
+                import numexpr
+                module = "numexpr"
+            except ImportError:
+                module = "numpy"
+                _logger.warning(
+                    "Numexpr is not installed, falling back to numpy, "
+                    "which is slower to calculate model."
+                )
+
         if linear_parameter_list is None:
             linear_parameter_list = []
         self._add_rotation = add_rotation
@@ -224,13 +240,13 @@ def __init__(self, expression, name, position=None, module="numpy",
                 if p.name not in linear_parameter_list:
                     # _parsed_expr used "non public" parameter name and we
                     # need to use the correct parameter name by using
-                    # _rename_pars_inv 
+                    # _rename_pars_inv
                     p._linear = _check_parameter_linearity(
                         self._parsed_expr,
                         self._rename_pars_inv.get(p.name, p.name)
                         )
 
-    def compile_function(self, module="numpy", position=False):
+    def compile_function(self, module, position=False):
         """
         Compile the function and calculate the gradient automatically when
         possible.
@@ -379,7 +395,7 @@ def _separate_pseudocomponents(self):
         Separate an expression into a group of lambdified functions
         that can compute the free parts of the expression, and a single
         lambdified function that computes the fixed parts of the expression
-        
+
         Used by the _compute_expression_part method.
         """
         expr = self._str_expression
@@ -455,4 +471,3 @@ def _check_parameter_linearity(expr, name):
                       "determined automatically.", UserWarning)
         return False
     return True
-
diff --git a/hyperspy/_components/gaussian.py b/hyperspy/_components/gaussian.py
index badb96c921..6198794847 100644
--- a/hyperspy/_components/gaussian.py
+++ b/hyperspy/_components/gaussian.py
@@ -105,7 +105,7 @@ class Gaussian(Expression):
     ~._components.gaussianhf.GaussianHF
     """
 
-    def __init__(self, A=1., sigma=1., centre=0., module="numexpr", **kwargs):
+    def __init__(self, A=1., sigma=1., centre=0., module=None, **kwargs):
         super().__init__(
             expression="A * (1 / (sigma * sqrt(2*pi))) * exp(-(x - centre)**2 \
                         / (2 * sigma**2))",
diff --git a/hyperspy/_components/gaussian2d.py b/hyperspy/_components/gaussian2d.py
index cd7e33e16c..461450a1fb 100644
--- a/hyperspy/_components/gaussian2d.py
+++ b/hyperspy/_components/gaussian2d.py
@@ -70,7 +70,7 @@ class Gaussian2D(Expression):
     """
 
     def __init__(self, A=1., sigma_x=1., sigma_y=1., centre_x=0.,
-                 centre_y=0, module="numexpr", **kwargs):
+                 centre_y=0, module=None, **kwargs):
         super().__init__(
             expression="A * (1 / (sigma_x * sigma_y * 2 * pi)) * \
                        exp(-((x - centre_x) ** 2 / (2 * sigma_x ** 2) \
diff --git a/hyperspy/_components/gaussianhf.py b/hyperspy/_components/gaussianhf.py
index 4e7514e965..0827d0c900 100644
--- a/hyperspy/_components/gaussianhf.py
+++ b/hyperspy/_components/gaussianhf.py
@@ -83,7 +83,7 @@ class GaussianHF(Expression):
 
     """
 
-    def __init__(self, height=1., fwhm=1., centre=0., module="numexpr",
+    def __init__(self, height=1., fwhm=1., centre=0., module=None,
                  **kwargs):
         super().__init__(
             expression="height * exp(-(x - centre)**2 * 4 * log(2)/fwhm**2)",
diff --git a/hyperspy/_components/logistic.py b/hyperspy/_components/logistic.py
index beb556853d..b17a507dfa 100644
--- a/hyperspy/_components/logistic.py
+++ b/hyperspy/_components/logistic.py
@@ -57,7 +57,7 @@ class Logistic(Expression):
         :py:class:`~._components.expression.Expression` component.
     """
 
-    def __init__(self, a=1., b=1., c=1., origin=0., module="numexpr", **kwargs):
+    def __init__(self, a=1., b=1., c=1., origin=0., module=None, **kwargs):
         super().__init__(
             expression="a / (1 + b * exp(-c * (x - origin)))",
             name="Logistic",
diff --git a/hyperspy/_components/lorentzian.py b/hyperspy/_components/lorentzian.py
index 0e2740e1db..90bdb07005 100644
--- a/hyperspy/_components/lorentzian.py
+++ b/hyperspy/_components/lorentzian.py
@@ -94,7 +94,7 @@ class Lorentzian(Expression):
     the full-with-half-maximum and height of the distribution, respectively.
     """
 
-    def __init__(self, A=1., gamma=1., centre=0., module="numexpr", **kwargs):
+    def __init__(self, A=1., gamma=1., centre=0., module=None, **kwargs):
         # We use `_gamma` internally to workaround the use of the `gamma`
         # function in sympy
         super().__init__(
diff --git a/hyperspy/_components/polynomial.py b/hyperspy/_components/polynomial.py
index 017389aadc..7ee2d3b2b4 100644
--- a/hyperspy/_components/polynomial.py
+++ b/hyperspy/_components/polynomial.py
@@ -51,7 +51,7 @@ class Polynomial(Expression):
 
     """
 
-    def __init__(self, order=2, module="numexpr", **kwargs):
+    def __init__(self, order=2, module=None, **kwargs):
         if order == 0:
             raise ValueError("Polynomial of order 0 is not supported.")
         coeff_list = ['{}'.format(o).zfill(len(list(str(order)))) for o in
diff --git a/hyperspy/_components/power_law.py b/hyperspy/_components/power_law.py
index e994925266..1fa2b45b4f 100644
--- a/hyperspy/_components/power_law.py
+++ b/hyperspy/_components/power_law.py
@@ -61,7 +61,7 @@ class PowerLaw(Expression):
     """
 
     def __init__(self, A=10e5, r=3., origin=0., left_cutoff=0.0,
-                 module="numexpr", compute_gradients=False, **kwargs):
+                 module=None, compute_gradients=False, **kwargs):
         super().__init__(
             expression="where(left_cutoff<x, A*(-origin + x)**-r, 0)",
             name="PowerLaw",
diff --git a/hyperspy/_components/rc.py b/hyperspy/_components/rc.py
index 99044cccc2..6ca854e99f 100644
--- a/hyperspy/_components/rc.py
+++ b/hyperspy/_components/rc.py
@@ -54,7 +54,7 @@ class RC(Expression):
 
     """
 
-    def __init__(self, Vmax=1., V0=0., tau=1., module="numexpr", **kwargs):
+    def __init__(self, Vmax=1., V0=0., tau=1., module=None, **kwargs):
         super().__init__(
             expression="V0 + Vmax * (1 - exp(-x / tau))",
             name="RC",
diff --git a/hyperspy/_signals/lazy.py b/hyperspy/_signals/lazy.py
index 01bd00a4db..ea087c1db1 100644
--- a/hyperspy/_signals/lazy.py
+++ b/hyperspy/_signals/lazy.py
@@ -60,6 +60,19 @@
     _logger.info("Dask widgets not loaded (dask >=2021.11.1 is required)")
 
 
+def _get():
+    try:
+        get = dask.threaded.get
+    except AttributeError:  # pragma: no cover
+        # For pyodide
+        get = dask.get
+        _logger.warning(
+            "Dask scheduler with threads is not available in this environment. "
+            "Falling back to synchronous scheduler (single-threaded)."
+        )
+    return get
+
+
 def to_array(thing, chunks=None):
     """Accepts BaseSignal, dask or numpy arrays and always produces either
     numpy or dask array.
@@ -769,7 +782,7 @@ def _calculate_summary_statistics(self, rechunk=False):
 
     def _block_iterator(self,
                         flat_signal=True,
-                        get=dask.threaded.get,
+                        get=None,
                         navigation_mask=None,
                         signal_mask=None):
         """A function that allows iterating lazy signal data by blocks,
@@ -783,9 +796,10 @@ def _block_iterator(self,
             optionally masked elements missing. If false, returns
             the equivalent of s.inav[{blocks}].data, where masked elements are
             set to np.nan or 0.
-        get : dask scheduler
-            the dask scheduler to use for computations;
-            default `dask.threaded.get`
+        get : dask scheduler or None
+            The dask scheduler to use for computations. If ``None``, 
+            ``dask.threaded.get` will be used if possible, otherwise
+            ``dask.get`` will be used, for example in pyodide interpreter.
         navigation_mask : {BaseSignal, numpy array, dask array}
             The navigation locations marked as True are not returned (flat) or
             set to NaN or 0.
@@ -794,6 +808,8 @@ def _block_iterator(self,
             to NaN or 0.
 
         """
+        if get is None:
+            get = _get()
         self._make_lazy()
         data = self._data_aligned_with_axes
         nav_chunks = data.chunks[:self.axes_manager.navigation_dimension]
@@ -853,7 +869,7 @@ def decomposition(
         output_dimension=None,
         signal_mask=None,
         navigation_mask=None,
-        get=dask.threaded.get,
+        get=None,
         num_chunks=None,
         reproject=True,
         print_info=True,
@@ -875,9 +891,10 @@ def decomposition(
         output_dimension : int or None, default None
             Number of components to keep/calculate. If None, keep all
             (only valid for 'SVD' algorithm)
-        get : dask scheduler
-            the dask scheduler to use for computations;
-            default `dask.threaded.get`
+        get : dask scheduler or None
+            The dask scheduler to use for computations. If ``None``, 
+            ``dask.threaded.get` will be used if possible, otherwise
+            ``dask.get`` will be used, for example in pyodide interpreter.
         num_chunks : int or None, default None
             the number of dask chunks to pass to the decomposition model.
             More chunks require more memory, but should run faster. Will be
@@ -912,6 +929,8 @@ def decomposition(
         * :py:class:`~.learn.ornmf.ORNMF`
 
         """
+        if get is None:
+            get = _get()
         # Check algorithms requiring output_dimension
         algorithms_require_dimension = ["PCA", "ORPCA", "ORNMF"]
         if algorithm in algorithms_require_dimension and output_dimension is None:
diff --git a/hyperspy/component.py b/hyperspy/component.py
index 4fe37210fb..d00d16ff24 100644
--- a/hyperspy/component.py
+++ b/hyperspy/component.py
@@ -1213,8 +1213,8 @@ def as_dictionary(self, fullcopy=True):
         export_to_dictionary(self, self._whitelist, dic, fullcopy)
         from hyperspy.model import _COMPONENTS
         if self._id_name not in _COMPONENTS:
-            import dill
-            dic['_class_dump'] = dill.dumps(self.__class__)
+            import cloudpickle
+            dic['_class_dump'] = cloudpickle.dumps(self.__class__)
         return dic
 
     def _load_dictionary(self, dic):
diff --git a/hyperspy/conftest.py b/hyperspy/conftest.py
index 43868f40d3..c8e79ed797 100644
--- a/hyperspy/conftest.py
+++ b/hyperspy/conftest.py
@@ -16,6 +16,8 @@
 # You should have received a copy of the GNU General Public License
 # along with HyperSpy. If not, see <https://www.gnu.org/licenses/#GPL>.
 
+from pathlib import Path
+
 try:
     # Set traits toolkit to work in a headless system
     # Capture error when toolkit is already previously set which typically
@@ -73,3 +75,14 @@ def pytest_configure(config):
             "mpl_image_compare: dummy marker registration to allow running "
             "without the pytest-mpl plugin."
         )
+
+
+def pytest_configure(config):
+    # raise an error if the baseline images are not present
+    # which is the case when installing from a wheel
+    baseline_images_path = Path(__file__).parent / "tests" / "drawing" / "plot_signal"
+    if config.getoption("--mpl") and not baseline_images_path.exists():
+        raise ValueError(
+            "`--mpl` flag can't not be used because the "
+            "baseline images are not packaged."
+            )
diff --git a/hyperspy/decorators.py b/hyperspy/decorators.py
index 06060f2d36..87d66cc4d3 100644
--- a/hyperspy/decorators.py
+++ b/hyperspy/decorators.py
@@ -18,11 +18,16 @@
 
 from functools import wraps
 import inspect
+import logging
 from typing import Callable, Optional, Union
 import warnings
 
 import numpy as np
 
+
+_logger = logging.getLogger(__name__)
+
+
 def lazify(func, **kwargs):
     from hyperspy.signal import BaseSignal
     from hyperspy.model import BaseModel
@@ -232,4 +237,27 @@ def wrapped(*args, **kwargs):
                 )
             return func(*args, **kwargs)
 
-        return wrapped
\ No newline at end of file
+        return wrapped
+
+
+def jit_ifnumba(*args, **kwargs):
+    try:
+        import numba
+
+        if "nopython" not in kwargs:
+            kwargs["nopython"] = True
+        return numba.jit(*args, **kwargs)
+    except ImportError:
+
+        _logger.warning(
+            "Numba is not installed, falling back to "
+            "non-accelerated implementation."
+        )
+
+        def wrap1(func):
+            def wrap2(*args2, **kwargs2):
+                return func(*args2, **kwargs2)
+
+            return wrap2
+
+        return wrap1
diff --git a/hyperspy/misc/array_tools.py b/hyperspy/misc/array_tools.py
index 70185e875a..a86fec8266 100644
--- a/hyperspy/misc/array_tools.py
+++ b/hyperspy/misc/array_tools.py
@@ -21,9 +21,8 @@
 
 import dask.array as da
 import numpy as np
-from numba import njit
-
 
+from hyperspy.decorators import jit_ifnumba
 from hyperspy.misc.math_tools import anyfloatin
 from hyperspy.docstrings.utils import REBIN_ARGS
 
@@ -193,7 +192,7 @@ def rebin(a, new_shape=None, scale=None, crop=True, dtype=None):
 rebin.__doc__ %= REBIN_ARGS.replace("        ", "    ")
 
 
-@njit(cache=True)
+@jit_ifnumba(cache=True)
 def _linear_bin_loop(result, data, scale):  # pragma: no cover
     for j in range(result.shape[0]):
         # Begin by determining the upper and lower limits of a given new pixel.
@@ -355,7 +354,7 @@ def numba_histogram(data, bins, ranges):
     return _numba_histogram(data, bins, ranges)
 
 
-@njit(cache=True)
+@jit_ifnumba(cache=True)
 def _numba_histogram(data, bins, ranges):
     """
     Numba histogram computation requiring native endian datatype.
@@ -408,7 +407,7 @@ def get_signal_chunk_slice(index, chunks):
     raise ValueError("Index out of signal range.")
 
 
-@njit(cache=True)
+@jit_ifnumba(cache=True)
 def numba_closest_index_round(axis_array, value_array):
     """For each value in value_array, find the closest value in axis_array and
     return the result as a numpy array of the same shape as value_array.
@@ -434,7 +433,7 @@ def numba_closest_index_round(axis_array, value_array):
     return index_array
 
 
-@njit(cache=True)
+@jit_ifnumba(cache=True)
 def numba_closest_index_floor(axis_array, value_array):  # pragma: no cover
     """For each value in value_array, find the closest smaller value in
     axis_array and return the result as a numpy array of the same shape
@@ -460,7 +459,7 @@ def numba_closest_index_floor(axis_array, value_array):  # pragma: no cover
     return index_array
 
 
-@njit(cache=True)
+@jit_ifnumba(cache=True)
 def numba_closest_index_ceil(axis_array, value_array):  # pragma: no cover
     """For each value in value_array, find the closest larger value in
     axis_array and return the result as a numpy array of the same shape
@@ -485,7 +484,7 @@ def numba_closest_index_ceil(axis_array, value_array):  # pragma: no cover
     return index_array
 
 
-@njit(cache=True)
+@jit_ifnumba(cache=True)
 def round_half_towards_zero(array, decimals=0):  # pragma: no cover
     """
     Round input array using "half towards zero" strategy.
@@ -561,7 +560,7 @@ def get_value_at_index(array,
         return stop_
 
 
-@njit(cache=True)
+@jit_ifnumba(cache=True)
 def round_half_away_from_zero(array, decimals=0):  # pragma: no cover
     """
     Round input array using "half away from zero" strategy.
diff --git a/hyperspy/misc/eds/__init__.py b/hyperspy/misc/eds/__init__.py
deleted file mode 100644
index e69de29bb2..0000000000
diff --git a/hyperspy/misc/export_dictionary.py b/hyperspy/misc/export_dictionary.py
index 98bd003acd..c5d2d9b651 100644
--- a/hyperspy/misc/export_dictionary.py
+++ b/hyperspy/misc/export_dictionary.py
@@ -19,7 +19,7 @@
 from operator import attrgetter
 from hyperspy.misc.utils import attrsetter
 from copy import deepcopy
-import dill
+import cloudpickle
 from dask.array import Array
 
 
@@ -73,7 +73,7 @@ def export_to_dictionary(target, whitelist, dic, fullcopy=True):
         * 'fn': the targeted attribute is a function, and may be pickled. A
           tuple of (thing, value) will be exported to the dictionary,
           where thing is None if function is passed as-is, and True if
-          dill package is used to pickle the function, with the value as
+          cloudpickle package is used to pickle the function, with the value as
           the result of the pickle.
         * 'id': the id of the targeted attribute is exported (e.g. id(target.name))
         * 'sig': The targeted attribute is a signal, and will be converted to a
@@ -122,7 +122,7 @@ def export_to_dictionary(target, whitelist, dic, fullcopy=True):
                     value['data'] = deepcopy(value['data'])
         elif 'fn' in flags:
             if fullcopy:
-                value = (True, dill.dumps(value))
+                value = (True, cloudpickle.dumps(value))
             else:
                 value = (None, value)
         elif fullcopy:
@@ -156,7 +156,7 @@ def load_from_dictionary(target, dic):
         * 'init': object used for initialization of the target. Will be
           copied to the _whitelist after loading
         * 'fn': the targeted attribute is a function, and may have been
-          pickled (preferably with dill package).
+          pickled (preferably with cloudpickle package).
         * 'id': the id of the original object was exported and the
           attribute will not be set. The key has to be '_id_'
         * 'sig': The targeted attribute was a signal, and may have been
@@ -197,11 +197,11 @@ def reconstruct_object(flags, value):
             value._assign_subclass()
         return value
     if 'fn' in flags:
-        ifdill, thing = value
-        if ifdill is None:
+        ifcloudpickle, thing = value
+        if ifcloudpickle is None:
             return thing
-        if ifdill in [True, 'True', b'True']:
-            return dill.loads(thing)
+        if ifcloudpickle in [True, 'True', b'True']:
+            return cloudpickle.loads(thing)
         # should not be reached
         raise ValueError("The object format is not recognized")
     if isinstance(value, Array):
diff --git a/hyperspy/misc/holography/__init__.py b/hyperspy/misc/holography/__init__.py
deleted file mode 100644
index e69de29bb2..0000000000
diff --git a/hyperspy/misc/holography/example_signals/00_ref_Vbp_130V_0V_bin2_crop.hdf5 b/hyperspy/misc/holography/example_signals/00_ref_Vbp_130V_0V_bin2_crop.hdf5
deleted file mode 100644
index 41ce464ce7..0000000000
Binary files a/hyperspy/misc/holography/example_signals/00_ref_Vbp_130V_0V_bin2_crop.hdf5 and /dev/null differ
diff --git a/hyperspy/misc/holography/example_signals/01_holo_Vbp_130V_0V_bin2_crop.hdf5 b/hyperspy/misc/holography/example_signals/01_holo_Vbp_130V_0V_bin2_crop.hdf5
deleted file mode 100644
index 3e3ab04800..0000000000
Binary files a/hyperspy/misc/holography/example_signals/01_holo_Vbp_130V_0V_bin2_crop.hdf5 and /dev/null differ
diff --git a/hyperspy/misc/lowess_smooth.py b/hyperspy/misc/lowess_smooth.py
index 00f47107b4..fc25559980 100755
--- a/hyperspy/misc/lowess_smooth.py
+++ b/hyperspy/misc/lowess_smooth.py
@@ -18,7 +18,8 @@
 # https://gist.github.com/agramfort/850437
 
 import numpy as np
-from numba import njit
+
+from hyperspy.decorators import jit_ifnumba
 
 
 def lowess(y, x, f=2.0 / 3.0, n_iter=3):
@@ -51,7 +52,7 @@ def lowess(y, x, f=2.0 / 3.0, n_iter=3):
     return _lowess(y, x, f, n_iter)
 
 
-@njit(cache=True, nogil=True)
+@jit_ifnumba(cache=True, nogil=True)
 def _lowess(y, x, f=2.0 / 3.0, n_iter=3):  # pragma: no cover
     """Lowess smoother requiring native endian datatype (for numba).
 
diff --git a/hyperspy/model.py b/hyperspy/model.py
index de1980cbb9..fc6da82b19 100644
--- a/hyperspy/model.py
+++ b/hyperspy/model.py
@@ -25,7 +25,7 @@
 from contextlib import contextmanager
 from functools import partial
 
-import dill
+import cloudpickle
 import numpy as np
 import dask.array as da
 from dask.diagnostics import ProgressBar
@@ -124,8 +124,8 @@ def reconstruct_component(comp_dictionary, **init_args):
                 _COMPONENTS[_id]["module"]), _COMPONENTS[_id]["class"])
     elif "_class_dump" in comp_dictionary:
         # When a component is not registered using the extension mechanism,
-        # it is serialized using dill.
-        _class = dill.loads(comp_dictionary['_class_dump'])
+        # it is serialized using cloudpickle.
+        _class = cloudpickle.loads(comp_dictionary['_class_dump'])
     else:
         # For component saved with hyperspy <2.0 and moved to exspy
         if comp_dictionary["_id_name"] in EXSPY_HSPY_COMPONENTS:
diff --git a/hyperspy/samfire.py b/hyperspy/samfire.py
index 2c557ed37d..83e71f1760 100644
--- a/hyperspy/samfire.py
+++ b/hyperspy/samfire.py
@@ -20,7 +20,7 @@
 from multiprocessing import cpu_count
 import warnings
 
-import dill
+import cloudpickle
 import numpy as np
 
 from hyperspy.misc.utils import DictionaryTreeBrowser
@@ -184,7 +184,7 @@ def _setup(self, **kwargs):
         """Set up SAMFire - configure models, set up pool if necessary"""
         from hyperspy.samfire_utils.samfire_pool import SamfirePool
         self._figure = None
-        self.metadata._gt_dump = dill.dumps(self.metadata.goodness_test)
+        self.metadata._gt_dump = cloudpickle.dumps(self.metadata.goodness_test)
         self._enable_optional_components()
 
         if hasattr(self.model, '_suspend_auto_fine_structure_width'):
@@ -214,9 +214,9 @@ def start(self, **kwargs):
         if self._workers and self.pool is not None:
             self.pool.update_parameters()
         if 'min_function' in kwargs:
-            kwargs['min_function'] = dill.dumps(kwargs['min_function'])
+            kwargs['min_function'] = cloudpickle.dumps(kwargs['min_function'])
         if 'min_function_grad' in kwargs:
-            kwargs['min_function_grad'] = dill.dumps(
+            kwargs['min_function_grad'] = cloudpickle.dumps(
                 kwargs['min_function_grad'])
         self._args = kwargs
         num_of_strat = len(self.strategies)
diff --git a/hyperspy/samfire_utils/samfire_kernel.py b/hyperspy/samfire_utils/samfire_kernel.py
index a9cd51691d..3393341d63 100644
--- a/hyperspy/samfire_utils/samfire_kernel.py
+++ b/hyperspy/samfire_utils/samfire_kernel.py
@@ -124,7 +124,7 @@ def multi_kernel(
     import numpy as np
     import copy
     from hyperspy.utils.model_selection import AICc
-    import dill
+    import cloudpickle
 
     def generate_values_iterator(compnames, vals, turned_on_component_inds):
         turned_on_names = [compnames[i] for i in turned_on_component_inds]
@@ -148,7 +148,7 @@ def send_good_results(model, previous_switching, cur_p, result_q, ind):
         result['current'] = cur_p._identity
         result_q.put((ind, result, True))
 
-    test = dill.loads(test_dict)
+    test = cloudpickle.loads(test_dict)
     cur_p = current_process()
     previous_switching = []
     comb = []
diff --git a/hyperspy/samfire_utils/samfire_worker.py b/hyperspy/samfire_utils/samfire_worker.py
index 2c789e799a..7feb657759 100644
--- a/hyperspy/samfire_utils/samfire_worker.py
+++ b/hyperspy/samfire_utils/samfire_worker.py
@@ -22,7 +22,7 @@
 import sys
 from itertools import combinations, product
 from queue import Empty
-import dill
+import cloudpickle
 import numpy as np
 import matplotlib
 
@@ -167,11 +167,11 @@ def run_pixel(self, ind, value_dict):
 
         self.fitting_kwargs = self.value_dict.pop('fitting_kwargs', {})
         if 'min_function' in self.fitting_kwargs:
-            self.fitting_kwargs['min_function'] = dill.loads(
+            self.fitting_kwargs['min_function'] = cloudpickle.loads(
                 self.fitting_kwargs['min_function'])
         if 'min_function_grad' in self.fitting_kwargs and isinstance(
                 self.fitting_kwargs['min_function_grad'], bytes):
-            self.fitting_kwargs['min_function_grad'] = dill.loads(
+            self.fitting_kwargs['min_function_grad'] = cloudpickle.loads(
                 self.fitting_kwargs['min_function_grad'])
         self.model.signal.data[:] = self.value_dict.pop('signal.data')
 
@@ -239,7 +239,7 @@ def send_results(self, current=False):
         self.result_queue.put(to_send)
 
     def setup_test(self, test_string):
-        self.fit_test = dill.loads(test_string)
+        self.fit_test = cloudpickle.loads(test_string)
 
     def start_listening(self):
         self._listening = True
diff --git a/hyperspy/tests/component/test_backcompatibility.py b/hyperspy/tests/component/test_backcompatibility.py
index 22061f6f57..a3d69228f8 100644
--- a/hyperspy/tests/component/test_backcompatibility.py
+++ b/hyperspy/tests/component/test_backcompatibility.py
@@ -127,5 +127,7 @@ def test_loading_components_exspy_not_installed():
             m = s.models.restore('a')
             assert "exspy is not installed" in str(err.value)
     else:
+        # The model contains components using numexpr
+        pytest.importorskip("numexpr")
         # This should work fine
         m = s.models.restore('a')
diff --git a/hyperspy/tests/component/test_bleasdale.py b/hyperspy/tests/component/test_bleasdale.py
index 3a61a719ae..c90f7f957d 100644
--- a/hyperspy/tests/component/test_bleasdale.py
+++ b/hyperspy/tests/component/test_bleasdale.py
@@ -16,6 +16,7 @@
 # You should have received a copy of the GNU General Public License
 # along with HyperSpy. If not, see <https://www.gnu.org/licenses/#GPL>.
 
+import pytest
 
 import numpy as np
 
@@ -23,6 +24,7 @@
 
 
 def test_function():
+    pytest.importorskip("numexpr")
     g = Bleasdale()
     g.a.value = 1
     g.b.value = 2
@@ -30,9 +32,16 @@ def test_function():
     assert g.function(-0.5) == 0
     assert g.function(0) == 1
     assert g.function(12) == 0.2
-    np.testing.assert_allclose(g.function(-.48),5)
+    np.testing.assert_allclose(g.function(-.48), 5)
     assert g.grad_a(0) == -0.5
     assert g.grad_b(0) == 0
     assert g.grad_c(0) == 0
     assert g.grad_a(-1) == 0
     assert g.grad_b(-1) == 0
+
+
+def test_module_error():
+    with pytest.raises(ValueError):
+        Bleasdale(module="numpy")
+    with pytest.raises(ValueError):
+        Bleasdale(module="scipy")
diff --git a/hyperspy/tests/component/test_components.py b/hyperspy/tests/component/test_components.py
index 2974d9f575..5dfccc9d33 100644
--- a/hyperspy/tests/component/test_components.py
+++ b/hyperspy/tests/component/test_components.py
@@ -54,6 +54,9 @@ def test_creation_components1d(component_name):
         kwargs['signal1D'] = s
     elif component_name == 'Expression':
         kwargs.update({'expression': "a*x+b", "name": "linear"})
+    elif component_name == 'Bleasdale':
+        # This component only works with numexpr.
+        pytest.importorskip("numexpr")
 
     component = getattr(components1d, component_name)(**kwargs)
     component.function(np.arange(1, 100))
diff --git a/hyperspy/tests/drawing/signal_markers_hs1_7_5.hspy b/hyperspy/tests/drawing/data/signal_markers_hs1_7_5.hspy
similarity index 100%
rename from hyperspy/tests/drawing/signal_markers_hs1_7_5.hspy
rename to hyperspy/tests/drawing/data/signal_markers_hs1_7_5.hspy
diff --git a/hyperspy/tests/drawing/test_markers.py b/hyperspy/tests/drawing/test_markers.py
index c91c2d3254..4e9ce5115c 100644
--- a/hyperspy/tests/drawing/test_markers.py
+++ b/hyperspy/tests/drawing/test_markers.py
@@ -1100,7 +1100,7 @@ def test_load_old_markers():
     plt.savefig('test.png', dpi=300)
     s.save("signal_markers_hs1_7_5.hspy")
     """
-    s = hs.load(FILE_PATH / "signal_markers_hs1_7_5.hspy")
+    s = hs.load(FILE_PATH / "data" / "signal_markers_hs1_7_5.hspy")
     s.metadata.General.original_filename = ""
     s.tmp_parameters.filename = ""
     s.plot(axes_ticks=True)
diff --git a/hyperspy/tests/drawing/test_plot_signal1d.py b/hyperspy/tests/drawing/test_plot_signal1d.py
index 11a13a4a03..dabf90cff6 100644
--- a/hyperspy/tests/drawing/test_plot_signal1d.py
+++ b/hyperspy/tests/drawing/test_plot_signal1d.py
@@ -73,6 +73,7 @@ def _matplotlib_pick_event(figure, click, artist):
 
 @pytest.fixture
 def setup_teardown(request, scope="class"):
+    plot_testing = request.config.getoption("--mpl")
     try:
         import pytest_mpl
         # This option is available only when pytest-mpl is installed
@@ -84,7 +85,7 @@ def setup_teardown(request, scope="class"):
     # duplicate baseline images to match the test_name when the
     # parametrized 'test_plot_spectra' are run. For a same 'style', the
     # expected images are the same.
-    if mpl_generate_path_cmdopt is None:
+    if mpl_generate_path_cmdopt is None and plot_testing:
         for filename in _generate_filename_list(style):
             copyfile(f"{str(filename)[:-5]}.png", filename)
     yield
@@ -94,9 +95,10 @@ def setup_teardown(request, scope="class"):
     if mpl_generate_path_cmdopt:
         for filename in _generate_filename_list(style):
             copyfile(filename, f"{str(filename)[:-5]}.png")
-    # Delete the images that have been created in 'setup_class'
-    for filename in _generate_filename_list(style):
-        os.remove(filename)
+    if plot_testing:
+        # Delete the images that have been created in 'setup_class'
+        for filename in _generate_filename_list(style):
+            os.remove(filename)
 
 
 @pytest.mark.usefixtures("setup_teardown")
diff --git a/hyperspy/tests/samfire/test_samfire.py b/hyperspy/tests/samfire/test_samfire.py
index 3288beef06..c2933d0a8a 100644
--- a/hyperspy/tests/samfire/test_samfire.py
+++ b/hyperspy/tests/samfire/test_samfire.py
@@ -19,7 +19,7 @@
 import copy
 import gc
 
-import dill
+import cloudpickle
 import numpy as np
 import pytest
 
@@ -469,7 +469,7 @@ def setup_method(self, method):
         self.args = {}
         self.model_letter = 'sldkfjg'
         from hyperspy.samfire_utils.fit_tests import red_chisq_test as rct
-        self._gt_dump = dill.dumps(rct(tolerance=1.0))
+        self._gt_dump = cloudpickle.dumps(rct(tolerance=1.0))
         m_slice = m.inav[self.ind[::-1]]
         m_slice.store(self.model_letter)
         m_dict = m_slice.signal._to_dictionary(False)
diff --git a/hyperspy/tests/signals/test_find_peaks1D_ohaver.py b/hyperspy/tests/signals/test_find_peaks1D_ohaver.py
index 7160cb201c..d1f25a4241 100644
--- a/hyperspy/tests/signals/test_find_peaks1D_ohaver.py
+++ b/hyperspy/tests/signals/test_find_peaks1D_ohaver.py
@@ -28,12 +28,10 @@
 class TestFindPeaks1DOhaver:
     def setup_method(self, method):
         with pytest.warns(VisibleDeprecationWarning):
-            filepath = (
-                Path(__file__)
-                .resolve()
-                .parent.joinpath("data/test_find_peaks1D_ohaver.hdf5")
-            )
-            self.signal = load(filepath, reader='hspy')
+            self.signal = load(
+                Path(__file__).parent / "data" / "test_find_peaks1D_ohaver.hdf5",
+                reader='hspy'
+                )
 
     def test_find_peaks1D_ohaver_high_amp_thres(self):
         signal1D = self.signal
diff --git a/hyperspy/utils/peakfinders2D.py b/hyperspy/utils/peakfinders2D.py
index 0a8cd967a3..99852d8eeb 100644
--- a/hyperspy/utils/peakfinders2D.py
+++ b/hyperspy/utils/peakfinders2D.py
@@ -18,13 +18,14 @@
 
 import copy
 
-from numba import njit
 import numpy as np
 import scipy.ndimage as ndi
 from skimage.feature import blob_dog, blob_log, match_template, peak_local_max
 
+from hyperspy.decorators import jit_ifnumba
 from hyperspy.misc.machine_learning import import_sklearn
 
+
 NO_PEAKS = np.array([[np.nan, np.nan]])
 
 def _get_peak_position_and_intensity(X, f, **kwargs):
@@ -59,7 +60,7 @@ def _get_peak_position_and_intensity(X, f, **kwargs):
         return np.concatenate([peaks, intensity[:, np.newaxis]], axis=1)
 
 
-@njit(cache=True)
+@jit_ifnumba(cache=True)
 def _fast_mean(X):  # pragma: no cover
     """JIT-compiled mean of array.
 
@@ -82,7 +83,7 @@ def _fast_mean(X):  # pragma: no cover
     return np.mean(X)
 
 
-@njit(cache=True)
+@jit_ifnumba(cache=True)
 def _fast_std(X):  # pragma: no cover
     """JIT-compiled standard deviation of array.
 
diff --git a/setup.py b/setup.py
index aa6fc074e5..16d880f276 100644
--- a/setup.py
+++ b/setup.py
@@ -34,17 +34,14 @@
 
 
 install_req = [
+    'cloudpickle',
     'dask[array]>=2021.3.1',
-    'dill',
     # included in stdlib since v3.8, but this required version requires Python 3.10
     # We can remove this requirement when the minimum supported version becomes Python 3.10
     'importlib-metadata>=3.6',
     'jinja2',
     'matplotlib>=3.1.3',
     'natsort',
-    # non-uniform axis requirement
-    'numba>=0.52',
-    'numexpr',
     'numpy>=1.20.0',
     'packaging',
     'pint>=0.10',
@@ -65,6 +62,7 @@
 extras_require = {
 	"ipython": ["IPython>7.0, !=8.0", "ipyparallel"],
     "learning": ["scikit-learn>=1.0.1"],
+    "speed":["numba", "numexpr"],
     # UPDATE BEFORE RELEASE
     "gui-jupyter": ["hyperspy_gui_ipywidgets @ git+https://github.com/ericpre/hyperspy_gui_ipywidgets.git@hyperspy2.0",
                     "ipympl"],
@@ -178,7 +176,6 @@ def __exit__(self, type, value, traceback):
                   'hyperspy.tests.misc',
                   'hyperspy.models',
                   'hyperspy.misc',
-                  'hyperspy.misc.holography',
                   'hyperspy.misc.machine_learning',
                   'hyperspy.external',
                   'hyperspy.external.astropy',
@@ -197,21 +194,9 @@ def __exit__(self, type, value, traceback):
             'hyperspy':
             [
                 'tests/component/data/*.hspy',
-                'tests/drawing/*.png',
                 'tests/drawing/data/*.hspy',
-                'tests/drawing/plot_signal/*.png',
-                'tests/drawing/plot_signal1d/*.png',
-                'tests/drawing/plot_signal2d/*.png',
-                'tests/drawing/plot_markers/*.png',
-                'tests/drawing/plot_model1d/*.png',
-                'tests/drawing/plot_model/*.png',
-                'tests/drawing/plot_roi/*.png',
                 'misc/dask_widgets/*.html.j2',
-                'tests/drawing/plot_mva/*.png',
-                'tests/drawing/plot_widgets/*.png',
-                'tests/drawing/plot_signal_tools/*.png',
-                'tests/signals/data/test_find_peaks1D_ohaver.hdf5',
-                'tests/signals/data/*.hspy',
+                'tests/signals/data/*.hdf5',
                 'hyperspy_extension.yaml',
             ],
         },
diff --git a/upcoming_changes/3255.enhancements.rst b/upcoming_changes/3255.enhancements.rst
new file mode 100644
index 0000000000..631704c572
--- /dev/null
+++ b/upcoming_changes/3255.enhancements.rst
@@ -0,0 +1,6 @@
+Make HyperSpy to work on pyodide:
+- ``numba`` and ``numexpr`` optional dependencies.
+- Replace ``dill`` by ``cloudpickle``.
+- Fallback to dask synchronous scheduler on pyodide.
+- Reduce packaging size to less than 1MB.
+- Add packaging test on GitHub CI.
\ No newline at end of file