diff --git a/_typos.toml b/_typos.toml
index b4283011..dc94d4ef 100644
--- a/_typos.toml
+++ b/_typos.toml
@@ -7,3 +7,6 @@ default.extend-ignore-identifiers-re = [
     "iy",
     "BA",
 ]
+
+[default.extend-words]
+eis = "eis"
diff --git a/changelog/207.breaking.rst b/changelog/207.breaking.rst
new file mode 100644
index 00000000..aab64a7c
--- /dev/null
+++ b/changelog/207.breaking.rst
@@ -0,0 +1,2 @@
+The previous coalignment API has been deleted and replaced with a new set of imports and functions.
+Please see this example: :ref:`sphx_glr_generated_gallery_aligning_aia_with_eis_maps.py`.
diff --git a/docs/code_ref/coalignment.rst b/docs/code_ref/coalignment.rst
index 86d17ad6..b66cd802 100644
--- a/docs/code_ref/coalignment.rst
+++ b/docs/code_ref/coalignment.rst
@@ -1 +1,15 @@
+Coalignment Package
+*******************
+
+`sunkit_image.coalignment` provides routines to perform coalignment of solar images.
+
+The main entry point is the `~sunkit_image.coalignment.coalign` function, which accepts a reference map, a target map, and a specified method for coalignment.
+This method returns a new map with updated metadata reflecting the applied affine transformations, such as scaling, rotation, and translation.
+
+The module supports various transformation methods registered via the `~sunkit_image.coalignment.interface.register_coalignment_method` decorator, allowing for flexible coalignment strategies based on the specific needs of the data.
+
 .. automodapi:: sunkit_image.coalignment
+
+.. automodapi:: sunkit_image.coalignment.interface
+
+.. automodapi:: sunkit_image.coalignment.match_template
diff --git a/docs/how_to_guide/adding_a_coalignment_method.rst b/docs/how_to_guide/adding_a_coalignment_method.rst
new file mode 100644
index 00000000..0ba7d6f2
--- /dev/null
+++ b/docs/how_to_guide/adding_a_coalignment_method.rst
@@ -0,0 +1,60 @@
+.. _sunkit-image-how-to-guide-add-a-new-coalignment-method:
+
+****************************
+Add a New Coalignment Method
+****************************
+
+If you want to register a new coalignment method that can be used by :func:`~sunkit_image.coalignment.coalign`, you can use :func:`~sunkit_image.coalignment.interface.register_coalignment_method`:
+
+.. code-block:: python
+
+    from sunkit_image.coalignment.interface import AffineParams, register_coalignment_method
+
+    @register_coalignment_method("my_coalign")
+    def my_coalignment_method(input_array, target_array, **kwargs):
+        # Your coalignment code goes here
+        # This should encompass calculating the shifts,
+        # handling NaN values appropriately.
+        # Return the shifts in an affine style, such as the scale, rotation and translation.
+        return AffineParams(scale, rotation, translation)
+
+Decorator Parameters
+====================
+
+Currently the decorator takes one parameter:
+
+- ``name`` : The name of your custom coalignment method, which in the above example is  "my_coalign".
+
+Function Requirements
+=====================
+
+Your coalignment function should:
+
+1. **Take Input Parameters**:
+
+   - ``input_array`` : The 2D array to be coaligned.
+   - ``target_array`` : The 2D array to align to.
+   - ``**kwargs``: So extra arguments can be passed down the stack.
+
+2. **Compute Shifts** : Calculate the shifts in the x and y directions needed to align ``input_array`` with ``target_array``.
+
+3. **Determine Affine Parameters** : Decide the values of the affine parameters - translation, scale and rotation.
+
+4. **Return** : Use the ``AffineParams`` named tuple included or provide your own that exposes the three parameters as attributes.
+
+Handling NaNs and Infs
+======================
+
+Proper handling of these values is expected to be included in the registered methods.
+The :func:`~sunkit_image.coalignment.coalign` function does not change any problematic values.
+
+Checking if the Method is Registered
+====================================
+
+To check if your method is registered, you can check if it is present in the registered methods dictionary ``REGISTERED_METHODS`` using the following code:
+
+.. code-block:: python
+
+    from sunkit_image.coalignment.interface import REGISTERED_METHODS
+
+    print(REGISTERED_METHODS)
diff --git a/docs/how_to_guide/index.rst b/docs/how_to_guide/index.rst
new file mode 100644
index 00000000..a9ff69b6
--- /dev/null
+++ b/docs/how_to_guide/index.rst
@@ -0,0 +1,10 @@
+.. _sunkit-image-how-to-reference:
+
+*************
+How To Guide
+*************
+
+.. toctree::
+   :maxdepth: 1
+
+   adding_a_coalignment_method
diff --git a/docs/index.rst b/docs/index.rst
index 79ec6ad7..9a66a828 100644
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -14,6 +14,7 @@ sunkit-image
           :maxdepth: 1
 
           code_ref/index
+          how_to_guide/index
 
     .. grid-item-card:: Other info
 
diff --git a/examples/aligning_aia_with_eis_maps.py b/examples/aligning_aia_with_eis_maps.py
new file mode 100644
index 00000000..71b63dd7
--- /dev/null
+++ b/examples/aligning_aia_with_eis_maps.py
@@ -0,0 +1,84 @@
+"""
+=====================
+Coaligning EIS to AIA
+=====================
+
+This example shows how to EIS data to AIA using cross-correlation which is implemented as the "match_template" method.
+"""
+# sphinx_gallery_thumbnail_number = 2
+
+import matplotlib.pyplot as plt
+
+import astropy.units as u
+from astropy.visualization import AsinhStretch, ImageNormalize
+
+import sunpy.map
+from sunpy.net import Fido
+from sunpy.net import attrs as a
+
+from sunkit_image.coalignment import coalign
+
+###################################################################################
+# Firstly, let us acquire the EIS and AIA data we need for this example.
+#
+# For this example, we will use the EIS data from the sunpy data repository.
+# This is a preprocessed EIS raster data.
+
+
+eis_map = sunpy.map.Map("https://github.com/sunpy/data/raw/main/sunkit-image/eis_20140108_095727.fe_12_195_119.2c-0.int.fits")
+
+fig = plt.figure()
+
+ax = fig.add_subplot(111, projection=eis_map)
+eis_map.plot(axes=ax, aspect=eis_map.meta['cdelt2'] / eis_map.meta['cdelt1'],
+             cmap='Blues_r', norm=ImageNormalize(stretch=AsinhStretch()))
+
+###################################################################################
+# Lets find the AIA image that we want to use as a reference.
+# We want to be using an image near the "date_average" of the EIS raster.
+
+query = Fido.search(a.Time(start=eis_map.meta["date_beg"], near=eis_map.meta["date_avg"], end=eis_map.meta["date_end"]), a.Instrument('aia'), a.Wavelength(193*u.angstrom))
+aia_file = Fido.fetch(query)
+aia_map = sunpy.map.Map(aia_file)
+
+####################################################################################
+# Before coaligning the images, we first downsample the AIA image to the same plate
+# scale as the EIS image. This is not done automatically.
+
+nx = (aia_map.scale.axis1 * aia_map.dimensions.x) / eis_map.scale.axis1
+ny = (aia_map.scale.axis2 * aia_map.dimensions.y) / eis_map.scale.axis2
+
+aia_downsampled = aia_map.resample(u.Quantity([nx, ny]))
+
+####################################################################################
+# Now we can coalign EIS to AIA using cross-correlation. For this we would be using the
+# "match_template" method. For details of the implementation refer to the
+# documentation of `~sunkit_image.coalignment.match_template.match_template_coalign`.
+
+coaligned_eis_map = coalign(aia_downsampled, eis_map)
+
+####################################################################################
+# To check now effective this has been, we will plot the EIS data and
+# overlap the bright regions from AIA before and after the coalignment.
+
+levels = [800] * aia_map.unit
+
+fig = plt.figure(figsize=(15, 7.5))
+
+# Before coalignment
+ax = fig.add_subplot(121, projection=eis_map)
+eis_map.plot(axes=ax, title='Before coalignment',
+             aspect=coaligned_eis_map.meta['cdelt2'] / coaligned_eis_map.meta['cdelt1'],
+             cmap='Blues_r', norm=ImageNormalize(stretch=AsinhStretch()))
+aia_map.draw_contours(levels, axes=ax)
+
+# After coalignment
+ax = fig.add_subplot(122, projection=coaligned_eis_map)
+coaligned_eis_map.plot(axes=ax, title='After coalignment',
+                       aspect=coaligned_eis_map.meta['cdelt2'] / coaligned_eis_map.meta['cdelt1'],
+                       cmap='Blues_r', norm=ImageNormalize(stretch=AsinhStretch()))
+aia_map.draw_contours(levels, axes=ax)
+
+fig.tight_layout()
+
+plt.show()
diff --git a/sunkit_image/coalignment.py b/sunkit_image/coalignment.py
deleted file mode 100644
index 62b2611f..00000000
--- a/sunkit_image/coalignment.py
+++ /dev/null
@@ -1,564 +0,0 @@
-"""
-This module provides routines for the co-alignment of images and
-`~sunpy.map.mapsequence.MapSequence` objects through both template matching and
-corrections due to solar rotation.
-"""
-
-import warnings
-from copy import deepcopy
-
-import numpy as np
-from scipy.ndimage import shift
-from skimage.feature import match_template
-
-import astropy.units as u
-
-import sunpy.map
-from sunpy.map.mapbase import GenericMap
-from sunpy.util.exceptions import SunpyUserWarning
-
-__all__ = [
-    "match_template_to_layer",
-    "apply_shifts",
-    "mapsequence_coalign_by_match_template",
-    "calculate_match_template_shift",
-]
-
-
-def _default_fmap_function(data):
-    """
-    This function ensures that the data are floats.
-
-    It is the default data manipulation function for the coalignment
-    method.
-    """
-    return np.float64(data)
-
-
-def _calculate_shift(this_layer, template):
-    """
-    Calculates the pixel shift required to put the template in the "best"
-    position on a layer.
-
-    Parameters
-    ----------
-    this_layer : `numpy.ndarray`
-        A numpy array of size ``(ny, nx)``, where the first two dimensions are
-        spatial dimensions.
-    template : `numpy.ndarray`
-        A numpy array of size ``(N, M)`` where ``N < ny`` and ``M < nx``.
-
-    Returns
-    -------
-    `tuple`
-        Pixel shifts ``(yshift, xshift)`` relative to the offset of the template
-        to the input array.
-    """
-    # Warn user if any NANs, Infs, etc are present in the layer or the template
-    _check_for_nonfinite_entries(this_layer, template)
-    # Calculate the correlation array matching the template to this layer
-    corr = match_template_to_layer(this_layer, template)
-    # Calculate the y and x shifts in pixels
-    return _find_best_match_location(corr)
-
-
-@u.quantity_input
-def _clip_edges(data, yclips: u.pix, xclips: u.pix):
-    """
-    Clips off the "y" and "x" edges of a 2D array according to a list of pixel
-    values. This function is useful for removing data at the edge of 2d images
-    that may be affected by shifts from solar de- rotation and layer co-
-    registration, leaving an image unaffected by edge effects.
-
-    Parameters
-    ----------
-    data : `numpy.ndarray`
-        A numpy array of shape ``(ny, nx)``.
-    yclips : `astropy.units.Quantity`
-        The amount to clip in the y-direction of the data. Has units of
-        pixels, and values should be whole non-negative numbers.
-    xclips : `astropy.units.Quantity`
-        The amount to clip in the x-direction of the data. Has units of
-        pixels, and values should be whole non-negative numbers.
-
-    Returns
-    -------
-    `numpy.ndarray`
-        A 2D image with edges clipped off according to ``yclips`` and ``xclips``
-        arrays.
-    """
-    ny = data.shape[0]
-    nx = data.shape[1]
-    # The purpose of the int below is to ensure integer type since by default
-    # astropy quantities are converted to floats.
-    return data[int(yclips[0].value) : ny - int(yclips[1].value), int(xclips[0].value) : nx - int(xclips[1].value)]
-
-
-@u.quantity_input
-def _calculate_clipping(y: u.pix, x: u.pix):
-    """
-    Return the upper and lower clipping values for the "y" and "x" directions.
-
-    Parameters
-    ----------
-    y : `astropy.units.Quantity`
-        An array of pixel shifts in the y-direction for an image.
-    x : `astropy.units.Quantity`
-        An array of pixel shifts in the x-direction for an image.
-
-    Returns
-    -------
-    `tuple`
-        The tuple is of the form ``([y0, y1], [x0, x1])``.
-        The number of (integer) pixels that need to be clipped off at each
-        edge in an image. The first element in the tuple is a list that gives
-        the number of pixels to clip in the y-direction. The first element in
-        that list is the number of rows to clip at the lower edge of the image
-        in y. The clipped image has "clipping[0][0]" rows removed from its
-        lower edge when compared to the original image. The second element in
-        that list is the number of rows to clip at the upper edge of the image
-        in y. The clipped image has "clipping[0][1]" rows removed from its
-        upper edge when compared to the original image. The second element in
-        the "clipping" tuple applies similarly to the x-direction (image
-        columns). The parameters ``y0, y1, x0, x1`` have the type
-        `~astropy.units.Quantity`.
-    """
-    return (
-        [_lower_clip(y.value), _upper_clip(y.value)] * u.pix,
-        [_lower_clip(x.value), _upper_clip(x.value)] * u.pix,
-    )
-
-
-def _upper_clip(z):
-    """
-    Find smallest integer bigger than all the positive entries in the input
-    array.
-    """
-    zcond = z >= 0
-    return int(np.max(np.ceil(z[zcond]))) if np.any(zcond) else 0
-
-
-def _lower_clip(z):
-    """
-    Find smallest positive integer bigger than the absolute values of the
-    negative entries in the input array.
-    """
-    zcond = z <= 0
-    return int(np.max(np.ceil(-z[zcond]))) if np.any(zcond) else 0
-
-
-def match_template_to_layer(layer, template):
-    """
-    Calculate the correlation array that describes how well the template
-    matches the layer. All inputs are assumed to be numpy arrays.
-
-    Parameters
-    ----------
-    layer : `numpy.ndarray`
-        A numpy array of size ``(ny, nx)``.
-    template : `numpy.ndarray`
-        A numpy array of size ``(N, M)`` where ``N < ny`` and ``M < nx``.
-
-    Returns
-    -------
-    `numpy.ndarray`
-        A correlation array between the layer and the template.
-        The values in the array range between 0 and 1.
-    """
-    return match_template(layer, template)
-
-
-def _find_best_match_location(corr):
-    """
-    Calculate an estimate of the location of the peak of the correlation result
-    in image pixels.
-
-    Parameters
-    ----------
-    corr : `numpy.ndarray`
-        A 2D correlation array.
-
-    Returns
-    -------
-    `~astropy.units.Quantity`
-        The shift amounts ``(y, x)`` in image pixels. Subpixel values are
-        possible.
-    """
-    # Get the index of the maximum in the correlation function
-    ij = np.unravel_index(np.argmax(corr), corr.shape)
-    cor_max_x, cor_max_y = ij[::-1]
-
-    # Get the correlation function around the maximum
-    array_maximum = corr[
-        np.max([0, cor_max_y - 1]) : np.min([cor_max_y + 2, corr.shape[0] - 1]),
-        np.max([0, cor_max_x - 1]) : np.min([cor_max_x + 2, corr.shape[1] - 1]),
-    ]
-    y_shift_maximum, x_shift_maximum = _get_correlation_shifts(array_maximum)
-
-    # Get shift relative to correlation array
-    y_shift_correlation_array = y_shift_maximum + cor_max_y * u.pix
-    x_shift_correlation_array = x_shift_maximum + cor_max_x * u.pix
-
-    return y_shift_correlation_array, x_shift_correlation_array
-
-
-def _get_correlation_shifts(array):
-    """
-    Estimate the location of the maximum of a fit to the input array. The
-    estimation in the "x" and "y" directions are done separately. The location
-    estimates can be used to implement subpixel shifts between two different
-    images.
-
-    Parameters
-    ----------
-    array : `numpy.ndarray`
-        An array with at least one dimension that has three elements. The
-        input array is at most a 3x3 array of correlation values calculated
-        by matching a template to an image.
-
-    Returns
-    -------
-    `~astropy.units.Quantity`
-        The ``(y, x)`` location of the peak of a parabolic fit, in image pixels.
-    """
-    # Check input shape
-    ny = array.shape[0]
-    nx = array.shape[1]
-    if nx > 3 or ny > 3:
-        msg = "Input array dimension should not be greater than 3 in any dimension."
-        raise ValueError(msg)
-
-    # Find where the maximum of the input array is
-    ij = np.unravel_index(np.argmax(array), array.shape)
-    x_max_location, y_max_location = ij[::-1]
-
-    # Estimate the location of the parabolic peak if there is enough data.
-    # Otherwise, just return the location of the maximum in a particular
-    # direction.
-    y_location = _parabolic_turning_point(array[:, x_max_location]) if ny == 3 else 1.0 * y_max_location
-
-    x_location = _parabolic_turning_point(array[y_max_location, :]) if nx == 3 else 1.0 * x_max_location
-
-    return y_location * u.pix, x_location * u.pix
-
-
-def _parabolic_turning_point(y):
-    """
-    Find the location of the turning point for a parabola ``y(x) = ax^2 + bx +
-    c``, given input values ``y(-1), y(0), y(1)``. The maximum is located at
-    ``x0 = -b / 2a``. Assumes that the input array represents an equally spaced
-    sampling at the locations ``y(-1), y(0) and y(1)``.
-
-    Parameters
-    ----------
-    y : `numpy.ndarray`
-        A one dimensional numpy array of shape "3" with entries that sample the
-        parabola at "-1", "0", and "1".
-
-    Returns
-    -------
-    `float`
-        A float, the location of the parabola maximum.
-    """
-    numerator = -0.5 * y.dot([-1, 0, 1])
-    denominator = y.dot([1, -2, 1])
-    return numerator / denominator
-
-
-def _check_for_nonfinite_entries(layer_image, template_image):
-    """
-    Issue a warning if there is any nonfinite entry in the layer or template
-    images.
-
-    Parameters
-    ----------
-    layer_image : `numpy.ndarray`
-        A two-dimensional `numpy.ndarray`.
-    template_image : `numpy.ndarray`
-        A two-dimensional `numpy.ndarray`.
-    """
-    if not np.all(np.isfinite(layer_image)):
-        warnings.warn(
-            "The layer image has nonfinite entries. "
-            "This could cause errors when calculating shift between two "
-            "images. Please make sure there are no infinity or "
-            "Not a Number values. For instance, replacing them with a "
-            "local mean.",
-            SunpyUserWarning,
-            stacklevel=3,
-        )
-
-    if not np.all(np.isfinite(template_image)):
-        warnings.warn(
-            "The template image has nonfinite entries. "
-            "This could cause errors when calculating shift between two "
-            "images. Please make sure there are no infinity or "
-            "Not a Number values. For instance, replacing them with a "
-            "local mean.",
-            SunpyUserWarning,
-            stacklevel=3,
-        )
-
-
-@u.quantity_input
-def apply_shifts(mc, yshift: u.pix, xshift: u.pix, *, clip=True, **kwargs):
-    """
-    Apply a set of pixel shifts to a `~sunpy.map.MapSequence`, and return a new
-    `~sunpy.map.MapSequence`.
-
-    Parameters
-    ----------
-    mc : `sunpy.map.MapSequence`
-        A `~sunpy.map.MapSequence` of shape ``(ny, nx, nt)``, where ``nt`` is the number of
-        layers in the `~sunpy.map.MapSequence`. ``ny`` is the number of pixels in the
-        "y" direction, ``nx`` is the number of pixels in the "x" direction.
-    yshift : `~astropy.units.Quantity`
-        An array of pixel shifts in the y-direction for an image.
-    xshift : `~astropy.units.Quantity`
-        An array of pixel shifts in the x-direction for an image.
-    clip : `bool`, optional
-        If `True` (default), then clip off "x", "y" edges of the maps in the sequence that are
-        potentially affected by edges effects.
-
-    Notes
-    -----
-    All other keywords are passed to `scipy.ndimage.shift`.
-
-    Returns
-    -------
-    `sunpy.map.MapSequence`
-        A `~sunpy.map.MapSequence` of the same shape as the input. All layers in
-        the `~sunpy.map.MapSequence` have been shifted according the input shifts.
-    """
-    # New mapsequence will be constructed from this list
-    new_mc = []
-
-    # Calculate the clipping
-    if clip:
-        yclips, xclips = _calculate_clipping(-yshift, -xshift)
-
-    # Shift the data and construct the mapsequence
-    for i, m in enumerate(mc):
-        shifted_data = shift(deepcopy(m.data), [yshift[i].value, xshift[i].value], **kwargs)
-        new_meta = deepcopy(m.meta)
-        # Clip if required.  Use the submap function to return the appropriate
-        # portion of the data.
-        if clip:
-            shifted_data = _clip_edges(shifted_data, yclips, xclips)
-            new_meta["naxis1"] = shifted_data.shape[1]
-            new_meta["naxis2"] = shifted_data.shape[0]
-            # Add one to go from zero-based to one-based indexing
-            new_meta["crpix1"] = m.reference_pixel.x.value + 1 + xshift[i].value - xshift[0].value
-            new_meta["crpix2"] = m.reference_pixel.y.value + 1 + yshift[i].value - yshift[0].value
-
-        new_map = sunpy.map.Map(shifted_data, new_meta)
-
-        # Append to the list
-        new_mc.append(new_map)
-
-    return sunpy.map.Map(new_mc, sequence=True)
-
-
-def calculate_match_template_shift(mc, template=None, layer_index=0, func=_default_fmap_function):
-    """
-    Calculate the arcsecond shifts necessary to co-register the layers in a
-    `~sunpy.map.MapSequence` according to a template taken from that
-    `~sunpy.map.MapSequence`.
-
-    When using this functionality, it is a good idea to check that the shifts
-    that were applied to were reasonable and expected. One way of checking this
-    is to animate the original `~sunpy.map.MapSequence`, animate the coaligned
-    `~sunpy.map.MapSequence`, and compare the differences you see to the
-    calculated shifts.
-
-    Parameters
-    ----------
-    mc : `sunpy.map.MapSequence`
-        A `~sunpy.map.MapSequence` of shape ``(ny, nx, nt)``, where ``nt`` is the number of
-        layers in the `~sunpy.map.MapSequence`.
-    template : {`None` , `~sunpy.map.Map` , `numpy.ndarray`}, optional
-        The template used in the matching. If an ~numpy.ndarray` is passed,
-        the `numpy.ndarray` has to have two dimensions.
-    layer_index : `int`,  optional
-        The template is assumed to refer to the map in the `~sunpy.map.MapSequence`
-        indexed by the value of "layer_index". Displacements of all maps in the
-        `~sunpy.map.MapSequence` are assumed to be relative to this layer. The
-        displacements of the template relative to this layer are therefore
-        ``(0, 0)``.
-    func : function, optional
-        A function which is applied to the data values before the coalignment
-        method is applied. This can be useful in coalignment, because it is
-        sometimes better to co-align on a function of the data rather than the
-        data itself. The calculated shifts are applied to the original data.
-        Examples of useful functions to consider for EUV images are the
-        logarithm or the square root. The function is of the form
-        ``func = F(data)``. The default function ensures that the data are
-        floats.
-    """
-    nt = len(mc.maps)
-
-    # Calculate a template.  If no template is passed then define one
-    # from the index layer.
-    if template is None:
-        # Size of the data
-        ny = mc.maps[layer_index].data.shape[0]
-        nx = mc.maps[layer_index].data.shape[1]
-        tplate = mc.maps[layer_index].data[int(ny / 4) : int(3 * ny / 4), int(nx / 4) : int(3 * nx / 4)]
-    elif isinstance(template, GenericMap):
-        tplate = template.data
-    elif isinstance(template, np.ndarray):
-        tplate = template
-    else:
-        msg = "Invalid template."
-        raise ValueError(msg)
-
-    # Apply the function to the template
-    tplate = func(tplate)
-
-    # Storage for the pixel shift
-    xshift_keep = np.zeros(nt) * u.pix
-    yshift_keep = np.zeros_like(xshift_keep)
-
-    # Storage for the arcsecond shift
-    xshift_arcseconds = np.zeros(nt) * u.arcsec
-    yshift_arcseconds = np.zeros_like(xshift_arcseconds)
-
-    # Match the template and calculate shifts
-    for i, m in enumerate(mc.maps):
-        # Get the next 2-d data array
-        this_layer = func(m.data)
-
-        # Calculate the y and x shifts in pixels
-        yshift, xshift = _calculate_shift(this_layer, tplate)
-
-        # Keep shifts in pixels
-        yshift_keep[i] = yshift
-        xshift_keep[i] = xshift
-
-    # Calculate shifts relative to the template layer
-    yshift_keep = yshift_keep - yshift_keep[layer_index]
-    xshift_keep = xshift_keep - xshift_keep[layer_index]
-
-    for i, m in enumerate(mc.maps):
-        # Calculate the shifts required in physical units, which are
-        # presumed to be arcseconds.
-        xshift_arcseconds[i] = xshift_keep[i] * m.scale[0]
-        yshift_arcseconds[i] = yshift_keep[i] * m.scale[1]
-
-    return {"x": xshift_arcseconds, "y": yshift_arcseconds}
-
-
-def mapsequence_coalign_by_match_template(
-    mc,
-    *,
-    template=None,
-    layer_index=0,
-    func=_default_fmap_function,
-    clip=True,
-    shift=None,
-    **kwargs,
-):
-    """
-    Co-register the layers in a `~sunpy.map.MapSequence` according to a
-    template taken from that `~sunpy.map.MapSequence`.
-
-    When using this functionality, it is a good
-    idea to check that the shifts that were applied were reasonable and
-    expected. One way of checking this is to animate the original
-    `~sunpy.map.MapSequence`, animate the coaligned `~sunpy.map.MapSequence`,
-    and compare the differences you see to the calculated shifts.
-
-    Currently this module provides image co-alignment by template matching.
-    and is partially inspired by the SSWIDL routine
-    `tr_get_disp.pro <http://www.heliodocs.com/php/xdoc_print.php?file=$SSW/trace/idl/util/tr_get_disp.pro>`__.
-    In this implementation, the template matching is handled via
-    `skimage.feature.match_template`.
-
-    Parameters
-    ----------
-    mc : `sunpy.map.MapSequence`
-        A `~sunpy.map.MapSequence` of shape ``(ny, nx, nt)``, where ``nt`` is the number of
-        layers in the `~sunpy.map.MapSequence`.
-    template : { `None` , `sunpy.map.Map` , `numpy.ndarray` }, optional
-        The template used in the matching. If an `numpy.ndarray` is passed,
-        the `numpy.ndarray` has to have two dimensions.
-    layer_index : `int`, optional
-        The template is assumed to refer to the map in the `~sunpy.map.MapSequence`
-        indexed by the value of ``layer_index``. Displacements of all maps in the
-        `~sunpy.map.MapSequence` are assumed to be relative to this layer. The
-        displacements of the template relative to this layer are therefore
-        ``(0, 0)``.
-    func : ``function``, optional
-        A function which is applied to the data values before the coalignment
-        method is applied. This can be useful in coalignment, because it is
-        sometimes better to co-align on a function of the data rather than the
-        data itself. The calculated shifts are applied to the original data.
-        Examples of useful functions to consider for EUV images are the
-        logarithm or the square root. The function is of the form
-        ``func = F(data)``. The default function ensures that the data are
-        floats.
-    clip : `bool`, optional
-        If True, then clip off x, y edges of the maps in the sequence that are
-        potentially affected by edges effects.
-    shift : `dict`, optional
-        A dictionary with two keys, 'x' and 'y'.  Key 'x' is an astropy
-        quantities array of corresponding to the amount of shift in the
-        x-direction (in arcseconds, assuming the helio-projective
-        Cartesian coordinate system) that is applied to the input
-        `~sunpy.map.MapSequence`.  Key 'y' is an `~astropy.units.Quantity` array
-        corresponding to the amount of shift in the y-direction (in arcseconds,
-        assuming the helio-projective Cartesian coordinate system) that is
-        applied to the input `~sunpy.map.MapSequence`. The number of elements in
-        each array must be the same as the number of maps in the
-        `~sunpy.map.MapSequence`.  If a shift is passed in to the function, that
-        shift is applied to the input `~sunpy.map.MapSequence` and the template
-        matching algorithm is not used.
-
-    Notes
-    -----
-    The remaining keyword arguments are sent to `~sunkit_image.coalignment.apply_shifts`.
-
-    Returns
-    -------
-    `sunpy.map.MapSequence`
-        A `~sunpy.map.MapSequence` that has co-aligned by matching the template.
-
-    Examples
-    --------
-    >>> from sunkit_image.coalignment import mapsequence_coalign_by_match_template as mc_coalign
-    >>> coaligned_mc = mc_coalign(mc)  # doctest: +SKIP
-    >>> coaligned_mc = mc_coalign(mc, layer_index=-1)  # doctest: +SKIP
-    >>> coaligned_mc = mc_coalign(mc, clip=False)  # doctest: +SKIP
-    >>> coaligned_mc = mc_coalign(mc, template=sunpy_map)  # doctest: +SKIP
-    >>> coaligned_mc = mc_coalign(mc, template=two_dimensional_ndarray)  # doctest: +SKIP
-    >>> coaligned_mc = mc_coalign(mc, func=np.log)  # doctest: +SKIP
-
-    References
-    ----------
-    * http://scribblethink.org/Work/nvisionInterface/nip.html
-    * J.P. Lewis, Fast Template Matching, Vision Interface 95, Canadian Image
-      Processing and Pattern Recognition Society, Quebec City, Canada, May 15-19,
-      1995, p. 120-123 http://www.scribblethink.org/Work/nvisionInterface/vi95_lewis.pdf.
-    """
-    # Number of maps
-    nt = len(mc.maps)
-
-    # Storage for the pixel shifts and the shifts in arcseconds
-    xshift_keep = np.zeros(nt) * u.pix
-    yshift_keep = np.zeros_like(xshift_keep)
-
-    if shift is None:
-        shifts = calculate_match_template_shift(mc, template=template, layer_index=layer_index, func=func)
-        xshift_arcseconds = shifts["x"]
-        yshift_arcseconds = shifts["y"]
-    else:
-        xshift_arcseconds = shift["x"]
-        yshift_arcseconds = shift["y"]
-
-    # Calculate the pixel shifts
-    for i, m in enumerate(mc):
-        xshift_keep[i] = xshift_arcseconds[i] / m.scale[0]
-        yshift_keep[i] = yshift_arcseconds[i] / m.scale[1]
-
-    # Apply the shifts and return the coaligned mapsequence
-    return apply_shifts(mc, -yshift_keep, -xshift_keep, clip=clip, **kwargs)
diff --git a/sunkit_image/coalignment/__init__.py b/sunkit_image/coalignment/__init__.py
new file mode 100644
index 00000000..977b6899
--- /dev/null
+++ b/sunkit_image/coalignment/__init__.py
@@ -0,0 +1,5 @@
+# This will register the function
+import sunkit_image.coalignment.match_template
+from sunkit_image.coalignment.interface import coalign
+
+__all__ = ["coalign"]
diff --git a/sunkit_image/coalignment/interface.py b/sunkit_image/coalignment/interface.py
new file mode 100644
index 00000000..f846ea2c
--- /dev/null
+++ b/sunkit_image/coalignment/interface.py
@@ -0,0 +1,186 @@
+import warnings
+from typing import NamedTuple
+
+import numpy as np
+
+import astropy
+import astropy.units as u
+from astropy.coordinates import SkyCoord
+
+from sunpy.sun.models import differential_rotation
+from sunpy.util.exceptions import SunpyUserWarning
+
+__all__ = ["AffineParams", "register_coalignment_method", "REGISTERED_METHODS"]
+
+# Global Dictionary to store the registered methods and their names
+REGISTERED_METHODS = {}
+
+
+class AffineParams(NamedTuple):
+    """
+    A 2-element tuple containing scale values defining the image scaling along
+    the x and y axes.
+    """
+    scale: tuple[float, float]
+    """
+    A 2x2 matrix defining the rotation transformation of the image.
+    """
+    rotation_matrix: tuple[tuple[float, float], tuple[float, float]]
+    """
+    A 2-element tuple stores the translation of the image along the x and y
+    axes.
+    """
+    translation: tuple[float, float]
+
+
+def register_coalignment_method(name):
+    """
+    Registers a coalignment method to be used by the coalignment interface.
+
+    Parameters
+    ----------
+    name : str
+        The name of the coalignment method.
+    """
+
+    def decorator(func):
+        REGISTERED_METHODS[name] = func
+        return func
+
+    return decorator
+
+
+def _update_fits_wcs_metadata(reference_map, target_map, affine_params):
+    """
+    Update the metadata of a sunpy.map.Map` based on affine transformation
+    parameters.
+
+    Parameters
+    ----------
+    reference_map : `sunpy.map.Map`
+        The reference map object to which the target map is to be coaligned.
+    target_map : `sunpy.map.Map`
+        The original map object whose metadata is to be updated.
+    affine_params : `NamedTuple`
+        A `NamedTuple` containing the affine transformation parameters.
+        If you want to use a custom object, it must have attributes for "translation" (dx, dy), "scale", and "rotation_matrix".
+
+    Returns
+    -------
+    `sunpy.map.Map`
+        A new sunpy map object with updated metadata reflecting the affine transformation.
+    """
+    # Updating the PC_ij matrix
+    new_pc_matrix = target_map.rotation_matrix @ affine_params.rotation_matrix
+    # Calculate the new reference pixel.
+    old_reference_pixel = np.asarray([target_map.reference_pixel.x.value, target_map.reference_pixel.y.value])
+    new_reference_pixel = affine_params.scale * affine_params.rotation_matrix @ old_reference_pixel + affine_params.translation
+    reference_coord = reference_map.wcs.pixel_to_world(new_reference_pixel[0],new_reference_pixel[1])
+    Txshift = reference_coord.Tx - target_map.reference_coordinate.Tx
+    Tyshift = reference_coord.Ty - target_map.reference_coordinate.Ty
+
+    # Create a new map with the updated metadata
+    fixed_map = target_map.shift_reference_coord(Txshift, Tyshift)
+    fixed_map.meta["PC1_1"] = new_pc_matrix[0, 0]
+    fixed_map.meta["PC1_2"] = new_pc_matrix[0, 1]
+    fixed_map.meta["PC2_1"] = new_pc_matrix[1, 0]
+    fixed_map.meta["PC2_2"] = new_pc_matrix[1, 1]
+    fixed_map.meta['cdelt1'] = (target_map.scale[0] / affine_params.scale[0]).value
+    fixed_map.meta['cdelt2'] = (target_map.scale[1] / affine_params.scale[1]).value
+    return fixed_map
+
+
+def _warn_user_of_separation(reference_map, target_map):
+    """
+    Issues a warning if the separation between the ``reference_map`` and
+    ``target_map`` is large.
+
+    Parameters
+    ----------
+    reference_map : `sunpy.map.Map`
+        The reference map to which the target map is to be coaligned.
+    target_map : `sunpy.map.Map`
+        The target map to be coaligned to the reference map.
+    """
+    # Maximum angular separation allowed between the reference and target maps
+    tolerable_angular_separation = 1*u.deg
+    ref_coord = SkyCoord(reference_map.observer_coordinate)
+    target_coord = SkyCoord(target_map.observer_coordinate)
+    if astropy.__version__ >= "6.1.0":
+        angular_separation = ref_coord.separation(target_coord, origin_mismatch="ignore")
+    else:
+        angular_separation = ref_coord.separation(target_coord)
+    if angular_separation > tolerable_angular_separation:
+        warnings.warn(
+            "The angular separation between the reference and target maps is large. "
+            "This could cause errors when calculating shift between two "
+            "images. Please make sure the maps are close in space.",
+            SunpyUserWarning,
+            stacklevel=3,
+        )
+    # Calculate time difference and convert to separation angle
+    ref_time = reference_map.date
+    target_time = target_map.date
+    time_diff = abs(ref_time - target_time)
+    time_separation_angle = differential_rotation(time_diff.to(u.day), reference_map.center.Tx, model='howard')
+    if time_separation_angle > tolerable_angular_separation:
+        warnings.warn(
+            "The time difference between the reference and target maps in time is large. "
+            "This could cause errors when calculating shift between two "
+            "images. Please make sure the maps are close in time.",
+            SunpyUserWarning,
+            stacklevel=3,
+        )
+
+
+def coalign(reference_map, target_map, method='match_template', **kwargs):
+    """
+    Performs image coalignment using the specified method.
+
+    This function updates the metadata of the target map to align it with the reference map.
+
+    .. note::
+
+        * This function is intended to correct maps with known incorrect metadata.
+          It is not designed to address issues like differential rotation or changes in observer location, which are encoded in the coordinate metadata.
+        * The function modifies the metadata of the map, not the underlying array data.
+          For adjustments that involve coordinate transformations, consider using `~sunpy.map.GenericMap.reproject_to` instead.
+
+    Parameters
+    ----------
+    reference_map : `sunpy.map.Map`
+        The reference map to which the target map is to be coaligned.
+    target_map : `sunpy.map.Map`
+        The target map to be coaligned to the reference map.
+    method : {{{coalignment_function_names}}}, optional
+        The name of the registered coalignment method to use.
+        Defaults to 'match_template'.
+    kwargs : `dict`
+        Additional keyword arguments to pass to the registered method.
+
+    Returns
+    -------
+    `sunpy.map.Map`
+        The coaligned target map with the updated metadata.
+
+    Raises
+    ------
+    ValueError
+        If the specified method is not registered.
+    """
+    if method not in REGISTERED_METHODS:
+        msg = (f"Method {method} is not a registered method: {list(REGISTERED_METHODS.keys())}."
+        "The method needs to be registered, with the correct import.")
+        raise ValueError(msg)
+    target_array = target_map.data
+    reference_array = reference_map.data
+    _warn_user_of_separation(reference_map, target_map)
+    affine_params = REGISTERED_METHODS[method](reference_array, target_array, **kwargs)
+    return _update_fits_wcs_metadata(reference_map, target_map, affine_params)
+
+
+# Generate the string with allowable coalignment-function names for use in docstrings
+_coalignment_function_names = ", ".join([f"``'{name}'``" for name in REGISTERED_METHODS])
+# Insert into the docstring for coalign. We cannot use the add_common_docstring decorator
+# due to what would be a circular loop in definitions.
+coalign.__doc__ = coalign.__doc__.format(coalignment_function_names=_coalignment_function_names)  # type: ignore
diff --git a/sunkit_image/coalignment/match_template.py b/sunkit_image/coalignment/match_template.py
new file mode 100644
index 00000000..88cb2969
--- /dev/null
+++ b/sunkit_image/coalignment/match_template.py
@@ -0,0 +1,120 @@
+import numpy as np
+from skimage.feature import match_template
+
+from sunkit_image.coalignment.interface import AffineParams, register_coalignment_method
+
+__all__ = ["match_template_coalign"]
+
+
+def _parabolic_turning_point(y):
+    """
+    Calculate the turning point of a parabola given three points.
+
+    Parameters
+    ----------
+    y : `numpy.ndarray`
+        An array of three points defining the parabola.
+
+    Returns
+    -------
+    float
+        The x-coordinate of the turning point.
+    """
+    return (-0.5 * y.dot([-1, 0, 1]))/ y.dot([1, -2, 1])
+
+
+def _get_correlation_shifts(array):
+    """
+    Calculate the shifts in x and y directions based on the correlation array.
+
+    Parameters
+    ----------
+    array : `numpy.ndarray`
+        A 2D array representing the correlation values.
+
+    Returns
+    -------
+    tuple
+        The shifts in y and x directions.
+
+    Raises
+    ------
+    ValueError
+        If the input array dimensions are greater than 3 in any direction.
+    """
+    ny, nx = array.shape
+    if nx > 3 or ny > 3:
+        msg = "Input array dimension should not be greater than 3 in any dimension."
+        raise ValueError(msg)
+
+    ij = np.unravel_index(np.argmax(array), array.shape)
+    x_max_location, y_max_location = ij[::-1]
+
+    y_location = _parabolic_turning_point(array[:, x_max_location]) if ny == 3 else 1.0 * y_max_location
+    x_location = _parabolic_turning_point(array[y_max_location, :]) if nx == 3 else 1.0 * x_max_location
+
+    return y_location, x_location
+
+
+def _find_best_match_location(corr):
+    """
+    Find the best match location in the correlation array.
+
+    Parameters
+    ----------
+    corr : `numpy.ndarray`
+        The correlation array.
+
+    Returns
+    -------
+    tuple
+        The best match location in the y and x directions.
+    """
+    ij = np.unravel_index(np.argmax(corr), corr.shape)
+    cor_max_x, cor_max_y = ij[::-1]
+
+    array_maximum = corr[
+        max(0, cor_max_y - 1) : min(cor_max_y + 2, corr.shape[0] - 1),
+        max(0, cor_max_x - 1) : min(cor_max_x + 2, corr.shape[1] - 1),
+    ]
+
+    y_shift_maximum, x_shift_maximum = _get_correlation_shifts(array_maximum)
+
+    y_shift_correlation_array = y_shift_maximum + cor_max_y
+    x_shift_correlation_array = x_shift_maximum + cor_max_x
+
+    return y_shift_correlation_array, x_shift_correlation_array
+
+
+@register_coalignment_method("match_template")
+def match_template_coalign(input_array, target_array, **kwargs):
+    """
+    Perform coalignment by matching the template array to the input array.
+
+    Parameters
+    ----------
+    input_array : `numpy.ndarray`
+        The input 2D array to be coaligned.
+    template_array : `numpy.ndarray`
+        The template 2D array to align to.
+
+    Returns
+    -------
+    `sunkit_image.coalignment.interface.AffineParams`
+        A `NamedTuple` containing the following affine transformation parameters:
+
+        - scale : `list`
+            A list of tuples representing the scale transformation as an identity matrix.
+        - rotation : `float`
+            The rotation angle in radians, which is fixed at 0.0 in this function.
+        - translation : `tuple`
+            A tuple containing the x and y translation values.
+    """
+    corr = match_template(np.float64(input_array), np.float64(target_array))
+    # Find the best match location
+    y_shift, x_shift = _find_best_match_location(corr)
+    # Particularly for this method, there is no change in the rotation or scaling,
+    # hence the hardcoded values of scale to 1.0 & rotation to identity matrix
+    scale = np.array([1.0, 1.0])
+    rotation_matrix = np.eye(2)
+    return AffineParams(scale=scale, rotation_matrix=rotation_matrix, translation=(x_shift, y_shift))
diff --git a/sunkit_image/coalignment/tests/__init__.py b/sunkit_image/coalignment/tests/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/sunkit_image/coalignment/tests/test_coalignment.py b/sunkit_image/coalignment/tests/test_coalignment.py
new file mode 100644
index 00000000..46498f73
--- /dev/null
+++ b/sunkit_image/coalignment/tests/test_coalignment.py
@@ -0,0 +1,118 @@
+import matplotlib.pyplot as plt
+import numpy as np
+import pytest
+from numpy.testing import assert_allclose
+
+import astropy.units as u
+from astropy.coordinates import SkyCoord
+from astropy.io import fits
+
+import sunpy.map
+from sunpy.net import Fido
+from sunpy.net import attrs as a
+
+from sunkit_image.coalignment import coalign
+from sunkit_image.coalignment.interface import REGISTERED_METHODS, AffineParams, register_coalignment_method
+from sunkit_image.tests.helpers import figure_test
+
+
+@pytest.fixture()
+def eis_test_map():
+    url = "https://github.com/sunpy/data/raw/main/sunkit-image/eis_20140108_095727.fe_12_195_119.2c-0.int.fits"
+    with fits.open(url) as hdul:
+        return sunpy.map.Map(hdul[0].data, hdul[0].header)
+
+
+@pytest.mark.remote_data()
+def test_coalignment(eis_test_map):
+    aia193_test_map = sunpy.map.Map(Fido.fetch(Fido.search(a.Time(start=eis_test_map.meta["date_beg"], near=eis_test_map.meta["date_avg"], end=eis_test_map.meta["date_end"]), a.Instrument('aia'), a.Wavelength(193*u.angstrom))))
+    nx = (aia193_test_map.scale.axis1 * aia193_test_map.dimensions.x) / eis_test_map.scale[0]
+    ny = (aia193_test_map.scale.axis2 * aia193_test_map.dimensions.y) / eis_test_map.scale[1]
+    aia193_test_downsampled_map = aia193_test_map.resample(u.Quantity([nx, ny]))
+    coaligned_eis_map = coalign(aia193_test_downsampled_map, eis_test_map, "match_template")
+    assert coaligned_eis_map.data.shape == eis_test_map.data.shape
+    assert_allclose(coaligned_eis_map.wcs.wcs.crval[0], aia193_test_downsampled_map.wcs.wcs.crval[0],rtol = 1e-2, atol = 0.13)
+    assert_allclose(coaligned_eis_map.wcs.wcs.crval[1], aia193_test_downsampled_map.wcs.wcs.crval[1],rtol = 1e-2, atol = 0.13)
+
+
+@pytest.fixture()
+def cutout_map(aia171_test_map):
+    aia_map = sunpy.map.Map(aia171_test_map)
+    bottom_left = SkyCoord(-300 * u.arcsec, -300 * u.arcsec, frame = aia_map.coordinate_frame)
+    top_right = SkyCoord(800 * u.arcsec, 600 * u.arcsec, frame = aia_map.coordinate_frame)
+    return aia_map.submap(bottom_left, top_right=top_right)
+
+
+def test_coalignment_reflects_pixel_shifts(cutout_map, aia171_test_map):
+    """
+    Check if coalignment adjusts world coordinates as expected based on
+    reference coordinate shifts.
+    """
+    messed_map = cutout_map.shift_reference_coord(25 * u.arcsec, 50 * u.arcsec)
+    original_world_coords = cutout_map.reference_coordinate
+    fixed_cutout_map = coalign(aia171_test_map, messed_map)
+    fixed_world_coords = fixed_cutout_map.reference_coordinate
+    # The actual shifts applied by coalignment should be equal to the expected shifts
+    assert_allclose(original_world_coords.Tx, fixed_world_coords.Tx, rtol=1e-2, atol=0.4)
+    assert_allclose(original_world_coords.Ty, fixed_world_coords.Ty, rtol=1e-2, atol=0.4)
+
+
+@figure_test
+def test_coalignment_figure(cutout_map, aia171_test_map):
+    levels = [200, 400, 500, 700, 800] * cutout_map.unit
+    messed_map = cutout_map.shift_reference_coord(25*u.arcsec, 50*u.arcsec)
+    fixed_cutout_map = coalign(aia171_test_map, messed_map)
+    fig = plt.figure(figsize=(15, 7.5))
+    # Before coalignment
+    ax1 = fig.add_subplot(131, projection=cutout_map)
+    cutout_map.plot(axes=ax1, title='Original Cutout')
+    cutout_map.draw_contours(levels, axes=ax1, alpha=0.3)
+    # Messed up map
+    ax2 = fig.add_subplot(132, projection=messed_map)
+    messed_map.plot(axes=ax2, title='Messed Cutout')
+    cutout_map.draw_contours(levels, axes=ax2, alpha=0.3)
+    # After coalignment
+    ax3 = fig.add_subplot(133, projection=fixed_cutout_map)
+    fixed_cutout_map.plot(axes=ax3, title='Fixed Cutout')
+    cutout_map.draw_contours(levels, axes=ax3, alpha=0.3)
+    fig.tight_layout()
+    return fig
+
+
+@pytest.fixture()
+@register_coalignment_method("scaling")
+def coalign_scaling(reference_map, target_map):
+    return AffineParams(scale=np.array([0.25, 0.25]), rotation_matrix=np.eye(2), translation=(0 , 0))
+
+
+def test_coalignment_reflects_scaling(cutout_map, aia171_test_map):
+    scale_factor = 4
+    rescaled_map = cutout_map.resample(u.Quantity([cutout_map.dimensions.x * scale_factor, cutout_map.dimensions.y * scale_factor]))
+    fixed_cutout_map = coalign(aia171_test_map, rescaled_map, method="scaling")
+    assert_allclose(fixed_cutout_map.scale[0].value, cutout_map.scale[0].value, rtol=1e-5, atol=0)
+    assert_allclose(fixed_cutout_map.scale[1].value, cutout_map.scale[1].value, rtol=1e-5, atol=0)
+
+
+@pytest.fixture()
+@register_coalignment_method("rotation")
+def coalign_rotation(reference_map, target_map):
+    rotation_matrix = np.array([[0.866, -0.5], [0.5, 0.866]])
+    return AffineParams(scale=np.array([1.0, 1.0]), rotation_matrix= rotation_matrix, translation=(0 , 0))
+
+
+def test_coalignment_reflects_rotation(cutout_map, aia171_test_map):
+    rotation_angle = 30 * u.deg
+    rotated_map = cutout_map.rotate(angle=rotation_angle)
+    fixed_cutout_map = coalign(aia171_test_map, rotated_map, method="rotation")
+    assert_allclose(fixed_cutout_map.rotation_matrix[0, 0], cutout_map.rotation_matrix[0, 0], rtol=1e-4, atol=0)
+    assert_allclose(fixed_cutout_map.rotation_matrix[1, 1], cutout_map.rotation_matrix[1, 1], rtol=1e-4, atol=0)
+
+
+def test_register_coalignment_method():
+    @register_coalignment_method("test_method")
+    def test_func():
+        return "Test function"
+
+    assert "test_method" in REGISTERED_METHODS
+    assert REGISTERED_METHODS["test_method"] == test_func
+    assert test_func() == "Test function"
diff --git a/sunkit_image/coalignment/tests/test_match_template.py b/sunkit_image/coalignment/tests/test_match_template.py
new file mode 100644
index 00000000..7baef8a1
--- /dev/null
+++ b/sunkit_image/coalignment/tests/test_match_template.py
@@ -0,0 +1,58 @@
+import numpy as np
+import pytest
+from numpy.testing import assert_allclose
+from skimage.feature import match_template
+
+import astropy.units as u
+
+from sunkit_image.coalignment.interface import REGISTERED_METHODS
+from sunkit_image.coalignment.match_template import (
+    _find_best_match_location,
+    _get_correlation_shifts,
+    _parabolic_turning_point,
+)
+
+
+def test_parabolic_turning_point():
+    assert _parabolic_turning_point(np.asarray([6.0, 2.0, 0.0])) == 1.5
+
+
+def test_get_correlation_shifts():
+    # Input array is 3 by 3, the most common case
+    test_array = np.zeros((3, 3))
+    test_array[1, 1] = 1
+    test_array[2, 1] = 0.6
+    test_array[1, 2] = 0.2
+    y_test, x_test = _get_correlation_shifts(test_array)
+    assert_allclose(y_test, 0.214285714286, rtol=1e-2, atol=0)
+    assert_allclose(x_test, 0.0555555555556, rtol=1e-2, atol=0)
+
+    # Input array is smaller in one direction than the other.
+    test_array = np.zeros((2, 2))
+    test_array[0, 0] = 0.1
+    test_array[0, 1] = 0.2
+    test_array[1, 0] = 0.4
+    test_array[1, 1] = 0.3
+    y_test, x_test = _get_correlation_shifts(test_array)
+    assert_allclose(y_test, 1.0, rtol=1e-2, atol=0)
+    assert_allclose(x_test, 0.0, rtol=1e-2, atol=0)
+
+    # Input array is too big in either direction
+    test_array = np.zeros((4, 3))
+    with pytest.raises(ValueError, match="Input array dimension should not be greater than 3 in any dimension."):
+        _get_correlation_shifts(test_array)
+
+    test_array = np.zeros((3, 4))
+    with pytest.raises(ValueError, match="Input array dimension should not be greater than 3 in any dimension."):
+        _get_correlation_shifts(test_array)
+
+
+def test_find_best_match_location(aia171_test_map, aia171_test_template, aia171_test_shift):
+    with np.errstate(all='ignore'):
+        result = match_template(aia171_test_map.data, aia171_test_template)
+    match_location = u.Quantity(_find_best_match_location(result))
+    assert_allclose(match_location.value, np.array(result.shape) / 2.0 - 0.5 + aia171_test_shift, rtol=1e-3, atol=0)
+
+
+def test_registered_method():
+    assert REGISTERED_METHODS["match_template"] is not None
diff --git a/sunkit_image/conftest.py b/sunkit_image/conftest.py
index 298d134d..0609495f 100644
--- a/sunkit_image/conftest.py
+++ b/sunkit_image/conftest.py
@@ -3,6 +3,7 @@
 import tempfile
 import warnings
 import importlib.util
+from pathlib import Path
 
 import numpy as np
 import pytest
@@ -14,6 +15,7 @@
 from astropy.utils.data import get_pkg_data_filename
 
 import sunpy.data.sample
+import sunpy.data.test
 import sunpy.map
 from sunpy.coordinates import Helioprojective, get_earth
 from sunpy.map.header_helper import make_fitswcs_header
@@ -198,5 +200,35 @@ def aia_171_map():
 
 @pytest.fixture()
 def hmi_map():
-    hmi_file = get_test_filepath("hmi_continuum_test_lowres_data.fits")
-    return sunpy.map.Map(hmi_file)
+    return sunpy.map.Map(get_test_filepath("hmi_continuum_test_lowres_data.fits"))
+
+
+@pytest.fixture()
+def aia171_test_clipping():
+    return np.asarray([0.2, -0.3, -1.0001])
+
+
+@pytest.fixture()
+def aia171_test_map():
+    testpath = sunpy.data.test.rootdir
+    return sunpy.map.Map(Path(testpath) / "aia_171_level1.fits")
+
+
+@pytest.fixture()
+def aia171_test_map_layer_shape(aia171_test_map_layer):
+    return aia171_test_map_layer.shape
+
+
+@pytest.fixture()
+def aia171_test_shift():
+    return np.array([3, 5])
+
+
+@pytest.fixture()
+def aia171_test_template(aia171_test_shift, aia171_test_map):
+    a1 = aia171_test_shift[0] + aia171_test_map.data.shape[0] // 4
+    a2 = aia171_test_shift[0] + 3 * aia171_test_map.data.shape[0] // 4
+    b1 = aia171_test_shift[1] + aia171_test_map.data.shape[1] // 4
+    b2 = aia171_test_shift[1] + 3 * aia171_test_map.data.shape[1] // 4
+    # Requires at least 32-bit floats to pass.
+    return aia171_test_map.data[a1:a2, b1:b2].astype("float32")
diff --git a/sunkit_image/tests/figure_hashes_mpl_390_ft_261_sunpy_600_astropy_610.json b/sunkit_image/tests/figure_hashes_mpl_390_ft_261_sunpy_600_astropy_610.json
index e159e946..40b76a36 100644
--- a/sunkit_image/tests/figure_hashes_mpl_390_ft_261_sunpy_600_astropy_610.json
+++ b/sunkit_image/tests/figure_hashes_mpl_390_ft_261_sunpy_600_astropy_610.json
@@ -1,4 +1,5 @@
 {
+  "sunkit_image.coalignment.tests.test_coalignment.test_coalignment_figure": "2339ce38ff576888472a5f527b80c24ff34f0f7e9d7631af0eeefc96accb7ac7",
   "sunkit_image.tests.test_enhance.test_mgn[array]": "f7a8e5a7377422c652a47217981d1c3bd1a66b5c211eae6ef18e4cef2660e438",
   "sunkit_image.tests.test_enhance.test_mgn[map]": "50bb490a6cc2408befe13c7f2a54f7433df80c2473dd21b619ace35de7e8f250",
   "sunkit_image.tests.test_enhance.test_mgn_submap": "cf3948d3ffa8ed4eacc500bee170db68bb1a96d2cec1c92e48f74c01827dc397",
diff --git a/sunkit_image/tests/test_coalignment.py b/sunkit_image/tests/test_coalignment.py
deleted file mode 100644
index d3a604f0..00000000
--- a/sunkit_image/tests/test_coalignment.py
+++ /dev/null
@@ -1,417 +0,0 @@
-import warnings
-from copy import deepcopy
-from pathlib import Path
-
-import numpy as np
-import pytest
-from numpy.testing import assert_allclose, assert_array_almost_equal
-from scipy.ndimage import shift as sp_shift
-
-import astropy.units as u
-from astropy.coordinates import SkyCoord
-
-import sunpy.data.test
-from sunpy.map import Map, MapSequence
-from sunpy.util import SunpyUserWarning
-
-from sunkit_image.coalignment import (
-    _calculate_clipping,
-    _check_for_nonfinite_entries,
-    _clip_edges,
-    _default_fmap_function,
-    _find_best_match_location,
-    _get_correlation_shifts,
-    _lower_clip,
-    _parabolic_turning_point,
-    _upper_clip,
-    apply_shifts,
-    calculate_match_template_shift,
-    mapsequence_coalign_by_match_template,
-    match_template_to_layer,
-)
-
-
-@pytest.fixture()
-def aia171_test_clipping():
-    return np.asarray([0.2, -0.3, -1.0001])
-
-
-@pytest.fixture()
-def aia171_test_map():
-    testpath = sunpy.data.test.rootdir
-    return sunpy.map.Map(Path(testpath) / "aia_171_level1.fits")
-
-
-@pytest.fixture()
-def aia171_test_shift():
-    return np.array([3, 5])
-
-
-@pytest.fixture()
-def aia171_test_map_layer(aia171_test_map):
-    return aia171_test_map.data.astype("float32")  # SciPy 1.4 requires at least 16-bit floats
-
-
-@pytest.fixture()
-def aia171_test_map_layer_shape(aia171_test_map_layer):
-    return aia171_test_map_layer.shape
-
-
-@pytest.fixture()
-def aia171_test_template(aia171_test_shift, aia171_test_map_layer, aia171_test_map_layer_shape):
-    # Test template
-    a1 = aia171_test_shift[0] + aia171_test_map_layer_shape[0] // 4
-    a2 = aia171_test_shift[0] + 3 * aia171_test_map_layer_shape[0] // 4
-    b1 = aia171_test_shift[1] + aia171_test_map_layer_shape[1] // 4
-    b2 = aia171_test_shift[1] + 3 * aia171_test_map_layer_shape[1] // 4
-    return aia171_test_map_layer[a1:a2, b1:b2]
-
-
-@pytest.fixture()
-def aia171_test_template_shape(aia171_test_template):
-    return aia171_test_template.shape
-
-
-def test_parabolic_turning_point():
-    assert _parabolic_turning_point(np.asarray([6.0, 2.0, 0.0])) == 1.5
-
-
-def test_check_for_nonfinite_entries():
-    with warnings.catch_warnings(record=True) as warning_list:
-        a = np.zeros((3, 3))
-        b = np.ones((3, 3))
-        _check_for_nonfinite_entries(a, b)
-
-    assert len(warning_list) == 0
-
-    for i in range(9):
-        for non_number in [np.nan, np.inf]:
-            a = np.ones(9)
-            a[i] = non_number
-            b = a.reshape(3, 3)
-
-            with pytest.warns(SunpyUserWarning, match="The layer image has nonfinite entries.") as warning_list:
-                _check_for_nonfinite_entries(b, np.ones((3, 3)))
-
-            assert len(warning_list) == 1
-
-            with pytest.warns(SunpyUserWarning, match="The template image has nonfinite entries.") as warning_list:
-                _check_for_nonfinite_entries(np.ones((3, 3)), b)
-
-            assert len(warning_list) == 1
-
-            with pytest.warns(Warning) as warning_list:
-                _check_for_nonfinite_entries(b, b)
-
-            assert len(warning_list) == 2
-
-
-def test_match_template_to_layer(
-    aia171_test_map_layer,
-    aia171_test_template,
-    aia171_test_map_layer_shape,
-    aia171_test_template_shape,
-):
-    result = match_template_to_layer(aia171_test_map_layer, aia171_test_template)
-    assert_allclose(
-        result.shape[0],
-        aia171_test_map_layer_shape[0] - aia171_test_template_shape[0] + 1,
-    )
-    assert_allclose(
-        result.shape[1],
-        aia171_test_map_layer_shape[1] - aia171_test_template_shape[1] + 1,
-    )
-    assert_allclose(np.max(result), 1.00, rtol=1e-2, atol=0)
-
-
-def test_get_correlation_shifts():
-    # Input array is 3 by 3, the most common case
-    test_array = np.zeros((3, 3))
-    test_array[1, 1] = 1
-    test_array[2, 1] = 0.6
-    test_array[1, 2] = 0.2
-    y_test, x_test = _get_correlation_shifts(test_array)
-    assert_allclose(y_test.value, 0.214285714286, rtol=1e-2, atol=0)
-    assert_allclose(x_test.value, 0.0555555555556, rtol=1e-2, atol=0)
-
-    # Input array is smaller in one direction than the other.
-    test_array = np.zeros((2, 2))
-    test_array[0, 0] = 0.1
-    test_array[0, 1] = 0.2
-    test_array[1, 0] = 0.4
-    test_array[1, 1] = 0.3
-    y_test, x_test = _get_correlation_shifts(test_array)
-    assert_allclose(y_test.value, 1.0, rtol=1e-2, atol=0)
-    assert_allclose(x_test.value, 0.0, rtol=1e-2, atol=0)
-
-    # Input array is too big in either direction
-    test_array = np.zeros((4, 3))
-    with pytest.raises(ValueError, match="Input array dimension should not be greater than 3 in any dimension."):
-        _get_correlation_shifts(test_array)
-
-    test_array = np.zeros((3, 4))
-    with pytest.raises(ValueError, match="Input array dimension should not be greater than 3 in any dimension."):
-        _get_correlation_shifts(test_array)
-
-
-def test_find_best_match_location(aia171_test_map_layer, aia171_test_template, aia171_test_shift):
-    result = match_template_to_layer(aia171_test_map_layer, aia171_test_template)
-    match_location = u.Quantity(_find_best_match_location(result))
-    assert_allclose(match_location.value, np.array(result.shape) / 2.0 - 0.5 + aia171_test_shift, rtol=1e-3, atol=0)
-
-
-def test_lower_clip(aia171_test_clipping):
-    # No element is less than zero
-    test_array = np.asarray([1.1, 0.1, 3.0])
-    assert _lower_clip(test_array) == 0
-    assert _lower_clip(aia171_test_clipping) == 2.0
-
-
-def test_upper_clip(aia171_test_clipping):
-    assert _upper_clip(aia171_test_clipping) == 1.0
-    # No element is greater than zero
-    test_array = np.asarray([-1.1, -0.1, -3.0])
-    assert _upper_clip(test_array) == 0
-
-
-def test_calculate_clipping(aia171_test_clipping):
-    answer = _calculate_clipping(aia171_test_clipping * u.pix, aia171_test_clipping * u.pix)
-    assert_array_almost_equal(answer, ([2.0, 1.0] * u.pix, [2.0, 1.0] * u.pix))
-
-
-def test_clip_edges():
-    a = np.zeros(shape=(341, 156))
-    yclip = [4, 0] * u.pix
-    xclip = [1, 2] * u.pix
-    new_a = _clip_edges(a, yclip, xclip)
-    assert a.shape[0] - (yclip[0].value + yclip[1].value) == new_a.shape[0]
-    assert a.shape[1] - (xclip[0].value + xclip[1].value) == new_a.shape[1]
-
-
-def test__default_fmap_function():
-    assert _default_fmap_function([1, 2, 3]).dtype == np.float64(1).dtype
-
-
-#
-# The following tests test functions that have mapsequences as inputs
-#
-# Setup the test mapsequences that have displacements
-# Pixel displacements have the y-displacement as the first entry
-
-
-@pytest.fixture()
-def aia171_test_mc_pixel_displacements():
-    return np.asarray([1.6, 10.1])
-
-
-@pytest.fixture()
-def aia171_mc_arcsec_displacements(aia171_test_mc_pixel_displacements, aia171_test_map):
-    return {
-        "x": np.asarray([0.0, aia171_test_mc_pixel_displacements[1] * aia171_test_map.scale[0].value]) * u.arcsec,
-        "y": np.asarray([0.0, aia171_test_mc_pixel_displacements[0] * aia171_test_map.scale[1].value]) * u.arcsec,
-    }
-
-
-@pytest.fixture()
-def aia171_test_mc(aia171_test_map, aia171_test_map_layer, aia171_test_mc_pixel_displacements):
-    # Create a map that has been shifted a known amount.
-    d1 = sp_shift(aia171_test_map_layer, aia171_test_mc_pixel_displacements)
-    m1 = Map((d1, aia171_test_map.meta))
-    # Create the mapsequence
-    return Map([aia171_test_map, m1], sequence=True)
-
-
-def test_calculate_match_template_shift(
-    aia171_test_mc,
-    aia171_mc_arcsec_displacements,
-    aia171_test_map,
-    aia171_test_map_layer,
-    aia171_test_map_layer_shape,
-):
-    # Define these local variables to make the code more readable
-    ny = aia171_test_map_layer_shape[0]
-    nx = aia171_test_map_layer_shape[1]
-
-    # Test to see if the code can recover the displacements.
-    test_displacements = calculate_match_template_shift(aia171_test_mc)
-    assert_allclose(test_displacements["x"], aia171_mc_arcsec_displacements["x"], rtol=5e-2, atol=0)
-    assert_allclose(test_displacements["y"], aia171_mc_arcsec_displacements["y"], rtol=5e-2, atol=0)
-
-    # Test setting the template as a ndarray
-    template_ndarray = aia171_test_map_layer[ny // 4 : 3 * ny // 4, nx // 4 : 3 * nx // 4]
-    test_displacements = calculate_match_template_shift(aia171_test_mc, template=template_ndarray)
-    assert_allclose(test_displacements["x"], aia171_mc_arcsec_displacements["x"], rtol=5e-2, atol=0)
-    assert_allclose(test_displacements["y"], aia171_mc_arcsec_displacements["y"], rtol=5e-2, atol=0)
-
-    # Test setting the template as GenericMap
-    submap = aia171_test_map.submap([nx / 4, ny / 4] * u.pix, top_right=[3 * nx / 4, 3 * ny / 4] * u.pix)
-    test_displacements = calculate_match_template_shift(aia171_test_mc, template=submap)
-    assert_allclose(test_displacements["x"], aia171_mc_arcsec_displacements["x"], rtol=5e-2, atol=0)
-    assert_allclose(test_displacements["y"], aia171_mc_arcsec_displacements["y"], rtol=5e-2, atol=0)
-
-    # Test setting the template as something other than a ndarray and a
-    # GenericMap. This should throw a ValueError.
-    with pytest.raises(ValueError, match="Invalid template."):
-        calculate_match_template_shift(aia171_test_mc, template="broken")
-
-
-def test_mapsequence_coalign_by_match_template(aia171_test_mc, aia171_test_map_layer_shape):
-    # Define these local variables to make the code more readable
-    ny = aia171_test_map_layer_shape[0]
-    nx = aia171_test_map_layer_shape[1]
-
-    # Get the calculated test displacements
-    test_displacements = calculate_match_template_shift(aia171_test_mc)
-
-    # Test passing in displacements
-    test_mc = mapsequence_coalign_by_match_template(aia171_test_mc, shift=test_displacements)
-
-    # Make sure the output is a mapsequence
-    assert isinstance(test_mc, MapSequence)
-
-    # Test returning with no clipping.  Output layers should have the same size
-    # as the original input layer.
-    test_mc = mapsequence_coalign_by_match_template(aia171_test_mc, clip=False)
-    assert test_mc[0].data.shape == aia171_test_map_layer_shape
-    assert test_mc[1].data.shape == aia171_test_map_layer_shape
-
-    # Test the returned mapsequence using the default - clipping on.
-    # All output layers should have the same size
-    # which is smaller than the input by a known amount
-    test_mc = mapsequence_coalign_by_match_template(aia171_test_mc)
-    x_displacement_pixels = test_displacements["x"] / test_mc[0].scale[0]
-    y_displacement_pixels = test_displacements["y"] / test_mc[0].scale[1]
-    expected_clipping = _calculate_clipping(y_displacement_pixels, x_displacement_pixels)
-    number_of_pixels_clipped = [np.sum(np.abs(expected_clipping[0])), np.sum(np.abs(expected_clipping[1]))]
-
-    assert test_mc[0].data.shape == (
-        ny - number_of_pixels_clipped[0].value,
-        nx - number_of_pixels_clipped[1].value,
-    )
-    assert test_mc[1].data.shape == (
-        ny - number_of_pixels_clipped[0].value,
-        nx - number_of_pixels_clipped[1].value,
-    )
-
-    # Test the returned mapsequence explicitly using clip=True.
-    # All output layers should have the same size
-    # which is smaller than the input by a known amount
-    test_mc = mapsequence_coalign_by_match_template(aia171_test_mc, clip=True)
-    x_displacement_pixels = test_displacements["x"] / test_mc[0].scale[0]
-    y_displacement_pixels = test_displacements["y"] / test_mc[0].scale[1]
-    expected_clipping = _calculate_clipping(y_displacement_pixels, x_displacement_pixels)
-    number_of_pixels_clipped = [np.sum(np.abs(expected_clipping[0])), np.sum(np.abs(expected_clipping[1]))]
-
-    assert test_mc[0].data.shape == (
-        ny - number_of_pixels_clipped[0].value,
-        nx - number_of_pixels_clipped[1].value,
-    )
-    assert test_mc[1].data.shape == (
-        ny - number_of_pixels_clipped[0].value,
-        nx - number_of_pixels_clipped[1].value,
-    )
-
-    # Test that the reference pixel of each map in the coaligned mapsequence is
-    # correct.
-    for im, m in enumerate(aia171_test_mc):
-        for i_s, s in enumerate(["x", "y"]):
-            assert_allclose(
-                aia171_test_mc[im].reference_pixel[i_s] - test_mc[im].reference_pixel[i_s],
-                test_displacements[s][im] / m.scale[i_s],
-                rtol=5e-2,
-                atol=0,
-            )
-
-
-def test_apply_shifts(aia171_test_map):
-    # take two copies of the AIA image and create a test mapsequence.
-    mc = Map([aia171_test_map, aia171_test_map], sequence=True)
-
-    # Pixel displacements have the y-displacement as the first entry
-    numerical_displacements = {"x": np.asarray([0.0, -2.7]), "y": np.asarray([0.0, -10.4])}
-    astropy_displacements = {
-        "x": numerical_displacements["x"] * u.pix,
-        "y": numerical_displacements["y"] * u.pix,
-    }
-
-    # Test to see if the code can detect the fact that the input shifts are not
-    # astropy quantities
-    with pytest.raises(TypeError):
-        apply_shifts(mc, numerical_displacements["y"], astropy_displacements["x"])
-    with pytest.raises(TypeError):
-        apply_shifts(mc, astropy_displacements["y"], numerical_displacements["x"])
-    with pytest.raises(TypeError):
-        apply_shifts(mc, numerical_displacements["y"], numerical_displacements["x"])
-
-    # Test returning with no extra options - the code returns a mapsequence only
-    test_output = apply_shifts(mc, astropy_displacements["y"], astropy_displacements["x"])
-    assert isinstance(test_output, MapSequence)
-
-    # Test returning with no clipping.  Output layers should have the same size
-    # as the original input layer.
-    test_mc = apply_shifts(mc, astropy_displacements["y"], astropy_displacements["x"], clip=False)
-    assert test_mc[0].data.shape == aia171_test_map.data.shape
-    assert test_mc[1].data.shape == aia171_test_map.data.shape
-
-    # Test returning with clipping.  Output layers should be smaller than the
-    # original layer by a known amount.
-    test_mc = apply_shifts(mc, astropy_displacements["y"], astropy_displacements["x"], clip=True)
-    for i in range(len(test_mc.maps)):
-        clipped = _calculate_clipping(astropy_displacements["y"], astropy_displacements["x"])
-        assert test_mc[i].data.shape[0] == mc[i].data.shape[0] - np.max(clipped[0].value)
-        assert test_mc[i].data.shape[1] == mc[i].data.shape[1] - np.max(clipped[1].value)
-
-    # Test returning with default clipping.  The default clipping is set to
-    # true, that is the mapsequence is clipped.  Output layers should be smaller
-    # than the original layer by a known amount.
-    test_mc = apply_shifts(mc, astropy_displacements["y"], astropy_displacements["x"])
-    for i in range(len(test_mc.maps)):
-        clipped = _calculate_clipping(astropy_displacements["y"], astropy_displacements["x"])
-        assert test_mc[i].data.shape[0] == mc[i].data.shape[0] - np.max(clipped[0].value)
-        assert test_mc[i].data.shape[1] == mc[i].data.shape[1] - np.max(clipped[1].value)
-
-    # Test that keywords are correctly passed
-    # Test for an individual keyword
-    test_mc = apply_shifts(mc, astropy_displacements["y"], astropy_displacements["x"], clip=False, cval=np.nan)
-    assert np.all(np.logical_not(np.isfinite(test_mc[1].data[:, -1])))
-
-    # Test for a combination of keywords, and that changing the interpolation
-    # order and how the edges are treated changes the results.
-    test_mc1 = apply_shifts(
-        mc,
-        astropy_displacements["y"],
-        astropy_displacements["x"],
-        clip=False,
-        order=2,
-        mode="reflect",
-    )
-    test_mc2 = apply_shifts(mc, astropy_displacements["y"], astropy_displacements["x"], clip=False)
-    assert np.all(test_mc1[1].data[:, -1] != test_mc2[1].data[:, -1])
-
-
-@pytest.fixture()
-def aia171_test_submap(aia171_test_map):
-    return aia171_test_map.submap(SkyCoord(((0, 0), (400, 500)) * u.arcsec, frame=aia171_test_map.coordinate_frame))
-
-
-@pytest.fixture()
-def aia171_test_mapsequence(aia171_test_submap):
-    m2header = deepcopy(aia171_test_submap.meta)
-    m2header["date-obs"] = "2011-02-15T01:00:00.34"
-    m2 = sunpy.map.Map((aia171_test_submap.data, m2header))
-    m3header = deepcopy(aia171_test_submap.meta)
-    m3header["date-obs"] = "2011-02-15T02:00:00.34"
-    m3 = sunpy.map.Map((aia171_test_submap.data, m3header))
-    return sunpy.map.Map([aia171_test_submap, m2, m3], sequence=True)
-
-
-@pytest.fixture()
-def known_displacements_layer_index0():
-    # Known displacements for these mapsequence layers when the layer index is set to 0
-    return {"x": np.asarray([0.0, -9.827465, -19.676442]), "y": np.asarray([0.0, 0.251137, 0.490014])}
-
-
-@pytest.fixture()
-def known_displacements_layer_index1():
-    # Known displacements for these mapsequence layers when the layer index is set to 1
-    return {"x": np.asarray([9.804878, 0.0, -9.827465]), "y": np.asarray([-0.263369, 0.0, 0.251137])}