Using the general rotation matrix to convert the shifts in arcseconds to pixel

coalignment.py

def mapsequence_coalign_by_rotation(mc, layer_index=0, clip=True, shift=None, **kwargs):
    """
    Move the layers in a mapsequence according to the input shifts.

    If an input shift is not given, the shifts due to
    solar rotation relative to an index layer is calculated and
    applied. 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 mapsequence, animate
    the derotated mapsequence, and compare the differences you see to the
    calculated shifts.

    Parameters
    ----------
    mc : `sunpy.map.MapSequence`
        A mapsequence of shape (ny, nx, nt), where nt is the number of layers in
        the mapsequence.
    layer_index : int
        Solar derotation shifts of all maps in the mapsequence are assumed
        to be relative to the layer in the mapsequence indexed by ``layer_index``.
    clip : bool
        If True, then clip off x, y edges in the datasequence that are potentially
        affected by edges effects.
    ``**kwargs``
        These keywords are passed to
        `~sunkit_image.coalignment.calculate_solar_rotate_shift`.

    Returns
    -------
    `sunpy.map.MapSequence`
        The results of the shifts applied to the input mapsequence.

    Examples
    --------
    >>> import sunpy.data.sample  # doctest: +REMOTE_DATA
    >>> from sunkit_image.coalignment import mapsequence_coalign_by_rotation
    >>> map1 = sunpy.map.Map(sunpy.data.sample.AIA_171_IMAGE)  # doctest: +REMOTE_DATA
    >>> map2 = sunpy.map.Map(sunpy.data.sample.EIT_195_IMAGE)  # doctest: +REMOTE_DATA
    >>> mc = sunpy.map.Map([map1, map2], sequence=True)  # doctest: +REMOTE_DATA
    >>> derotated_mc = mapsequence_coalign_by_rotation(mc)  # doctest: +SKIP
    >>> derotated_mc = mapsequence_coalign_by_rotation(mc, layer_index=-1)  # doctest: +SKIP
    >>> derotated_mc = mapsequence_coalign_by_rotation(mc, clip=False)  # doctest: +SKIP
    """
    nt = len(mc.maps)
    xshift_keep = np.zeros(nt) * u.pix
    yshift_keep = np.zeros_like(xshift_keep)

    if shift is None:
        shift = calculate_solar_rotate_shift(mc, layer_index=layer_index, **kwargs)
    xshift_arcseconds = shift["x"]
    yshift_arcseconds = shift["y"]

    for i, m in enumerate(mc):
        # Convert the WCS to a header
        header = m.wcs.to_header()

        # Extract the PC matrix from the header
        pc_matrix = np.array([[header['PC1_1'], header['PC1_2']],
                              [header['PC2_1'], header['PC2_2']]])
        # Calculate the inverse of the PC matrix
        pc_inv = np.linalg.inv(pc_matrix)
        
        # Extracting the CDELT values for the scale in x and y directions.
        scale_x = header['CDELT1'] * (u.arcsec / u.pix)
        scale_y = header['CDELT2'] * (u.arcsec / u.pix)
        
        # Calculate the shifts from arcseconds to pixels
        x_shift_pixels = xshift_arcseconds[i] / scale_x
        y_shift_pixels = yshift_arcseconds[i] / scale_y
        
        # Create a 2D array
        shifts_wcs = np.array([[x_shift_pixels] , [y_shift_wcs]])
        # Calculate the actual shift pixels
        xshift_keep[i], yshift_keep[i] = np.dot(shifts_wcs, pc_inv)
        
    return apply_shifts(mc, yshift_keep, xshift_keep, clip=clip)