diff --git a/.circleci/config.yml b/.circleci/config.yml index e763e1a4..42ee92ed 100644 --- a/.circleci/config.yml +++ b/.circleci/config.yml @@ -32,7 +32,7 @@ jobs: jobname: type: string docker: - - image: cimg/python:3.12 + - image: cimg/python:3.13 environment: TOXENV=<< parameters.jobname >> steps: - run: *no-backports @@ -56,7 +56,7 @@ jobs: jobname: type: string docker: - - image: cimg/python:3.12 + - image: cimg/python:3.13 environment: TOXENV: << parameters.jobname >> GIT_SSH_COMMAND: ssh -i ~/.ssh/id_rsa_6464b6a8248237ca368fd4690777d921 @@ -97,13 +97,13 @@ workflows: matrix: parameters: jobname: - - "py312-figure" + - "py313-figure" - deploy-reference-images: name: baseline-<< matrix.jobname >> matrix: parameters: jobname: - - "py312-figure" + - "py313-figure" requires: - << matrix.jobname >> filters: 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/290.breaking.rst b/changelog/290.breaking.rst index 078720c6..eed6e479 100644 --- a/changelog/290.breaking.rst +++ b/changelog/290.breaking.rst @@ -5,4 +5,4 @@ Increased the minimum required versions: - sunpy 7.0.0 (from 6.0.0) - NumPy 1.26.0 (from 1.23.5) - Astropy 6.1.0 (from 5.3.0) -- SciPy 1.11.0 (from 1.10.1) +- SciPy 1.12.0 (from 1.10.1) diff --git a/changelog/293.breaking.rst b/changelog/293.breaking.rst new file mode 100644 index 00000000..aab64a7c --- /dev/null +++ b/changelog/293.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/asda.rst b/docs/code_ref/asda.rst index 51f5f445..8ade76c6 100644 --- a/docs/code_ref/asda.rst +++ b/docs/code_ref/asda.rst @@ -1 +1,7 @@ +ASDA (`sunkit_image.asda`) +************************** + +This package contains an implementation of the `Automated Swirl Detection +Algorithm (ASDA) `__. + .. automodapi:: sunkit_image.asda diff --git a/docs/code_ref/coalignment.rst b/docs/code_ref/coalignment.rst index 86d17ad6..1e5e47a2 100644 --- a/docs/code_ref/coalignment.rst +++ b/docs/code_ref/coalignment.rst @@ -1 +1,16 @@ +Coalignment (`sunkit_image.coalignment`) +**************************************** + +Routines to perform coalignment of solar images. + .. automodapi:: sunkit_image.coalignment + +.. automodapi:: sunkit_image.coalignment.interface + :no-inheritance-diagram: + :skip: coalign_map + +.. automodapi:: sunkit_image.coalignment.register + +.. automodapi:: sunkit_image.coalignment.match_template + +.. automodapi:: sunkit_image.coalignment.phase_cross_correlation diff --git a/docs/code_ref/enhance.rst b/docs/code_ref/enhance.rst index e5843e35..591d65a0 100644 --- a/docs/code_ref/enhance.rst +++ b/docs/code_ref/enhance.rst @@ -1 +1,6 @@ +Enhancement (`sunkit_image.enhance`) +************************************ + +This package contains enhancement routines for solar physics data. + .. automodapi:: sunkit_image.enhance diff --git a/docs/code_ref/granule.rst b/docs/code_ref/granule.rst index 71b890e2..f4c74ad1 100644 --- a/docs/code_ref/granule.rst +++ b/docs/code_ref/granule.rst @@ -1 +1,6 @@ +Granule Detection (`sunkit_image.granule`) +****************************************** + +This package contains functions that will segment images for granule detection. + .. automodapi:: sunkit_image.granule diff --git a/docs/code_ref/radial.rst b/docs/code_ref/radial.rst index 053824b9..70fc16d5 100644 --- a/docs/code_ref/radial.rst +++ b/docs/code_ref/radial.rst @@ -1 +1,7 @@ +Radial Enhancement (`sunkit_image.radial`) +****************************************** + +This package contains functions that can be used to enhance the regions above a +given radius. + .. automodapi:: sunkit_image.radial diff --git a/docs/code_ref/stara.rst b/docs/code_ref/stara.rst index 240e7241..4140f81c 100644 --- a/docs/code_ref/stara.rst +++ b/docs/code_ref/stara.rst @@ -1 +1,7 @@ +Sunspot Tracking And Recognition (`sunkit_image.stara`) +******************************************************* + +This package contains an implementation of the `Sunspot Tracking And Recognition +Algorithm (STARA) `__. + .. automodapi:: sunkit_image.stara diff --git a/docs/code_ref/time_lag.rst b/docs/code_ref/time_lag.rst index a6d858ba..a6782492 100644 --- a/docs/code_ref/time_lag.rst +++ b/docs/code_ref/time_lag.rst @@ -1 +1,9 @@ +Time lag (`sunkit_image.time_lag`) +********************************** + +This package contains functions for calculating the cross-correlation and time +lag between intensity cubes. + +Useful for understanding time variability in EUV light curves. + .. automodapi:: sunkit_image.time_lag diff --git a/docs/code_ref/trace.rst b/docs/code_ref/trace.rst index 06d9758e..e986ee94 100644 --- a/docs/code_ref/trace.rst +++ b/docs/code_ref/trace.rst @@ -1 +1,6 @@ +Loop tracing (`sunkit_image.trace`) +*********************************** + +This package contains functions that will the trace out coronal loop-like structures in an image. + .. automodapi:: sunkit_image.trace diff --git a/docs/code_ref/utils.rst b/docs/code_ref/utils.rst index ac372006..e7f28489 100644 --- a/docs/code_ref/utils.rst +++ b/docs/code_ref/utils.rst @@ -1 +1,4 @@ +Utilities (`sunkit_image.utils`) +******************************** + .. automodapi:: sunkit_image.utils diff --git a/docs/conf.py b/docs/conf.py index 2ec64e32..920b7f4b 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -39,7 +39,7 @@ elif _version.is_devrelease: version = release = f"{_version.base_version}.dev{_version.dev}" is_development = _version.is_devrelease -is_release = not(_version.is_prerelease or _version.is_devrelease) +is_release = not (_version.is_prerelease or _version.is_devrelease) project = "sunkit_image" author = "The SunPy Community" @@ -105,7 +105,7 @@ ), "matplotlib": ("https://matplotlib.org/stable", None), "sunpy": ("https://docs.sunpy.org/en/stable/", None), - "ndcube": ('https://docs.sunpy.org/projects/ndcube/en/stable/', None), + "ndcube": ("https://docs.sunpy.org/projects/ndcube/en/stable/", None), "astropy": ("https://docs.astropy.org/en/stable/", None), "dask": ("https://docs.dask.org/en/latest", None), "skimage": ("https://scikit-image.org/docs/stable/", None), 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..2c365ae8 --- /dev/null +++ b/docs/how_to_guide/adding_a_coalignment_method.rst @@ -0,0 +1,57 @@ +.. _sunkit-image-how-to-guide-add-a-new-coalignment-method: + +**************************** +Add a New Coalignment Method +**************************** + +In addition to the coalignment methods provided in sunkit-image, users can also write their own +coalignment methods and "register" these methods with sunkit-image such that they can be used through the same interface as the builtin methods without having to alter the underlying +sunkit-image package. + +At a minimum, your new coalignment function should do the following: + +1. Take the following inputs: + + - ``target_array`` : The 2D array to be coaligned. + - ``reference_array`` : The 2D array to align to. + - ``**kwargs``: Optional keyword arguments used by your method + +2. Decide the values of the affine transformation - translation, scale and rotation. In most cases, this means calculating the shifts in the x- and y-directions needed to align ``input_array`` with ``target_array``. + +3. Return an instance of `~sunkit_image.coalignment.interface.AffineParams` with the results of your coalignment procedure. + +Additionally, registered methods are expected to handled NaNs and Infs should they arise as a result of your coalignment procedure. +The :func:`~sunkit_image.coalignment.coalign_map` function does not make any attempt to filter out +these non-finite values. + +To register your new coalignment method, you can use the :func:`~sunkit_image.coalignment.register.register_coalignment_method` decorator to register your new method with a custom name. An example of how to use this decorator is shown below: + +.. code-block:: python + + from sunkit_image.coalignment.interface import AffineParams, register_coalignment_method + + @register_coalignment_method("my_custom_coalignment_method") + def my_coalignment_method(target_array, reference_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) + + +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) + +If your coalignment method has been successfully registered, you should now be able to call it +through the `~sunkit_image.coalignment.coalign_map` interface: + +.. code-block:: python + + from sunkit_image.coalignment import coalign_map + + coaligned_map = coalign_map(target_map, reference_map, method='my_custom_coalignment_method') 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 d0c4d312..97216a5f 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/advanced_wow.py b/examples/advanced_wow.py index b690a943..61b2d1b7 100644 --- a/examples/advanced_wow.py +++ b/examples/advanced_wow.py @@ -63,7 +63,9 @@ # gamma-stretched input, with weight "h". gamma = 4 -bilateral_denoised_merged_wow = enhance.wow(submap_193, bilateral=1, denoise_coefficients=denoise_coefficients, gamma=gamma, h=0.99) +bilateral_denoised_merged_wow = enhance.wow( + submap_193, bilateral=1, denoise_coefficients=denoise_coefficients, gamma=gamma, h=0.99 +) ########################################################################### # Finally, we will plot the full set of outputs created and @@ -71,19 +73,19 @@ fig = plt.figure(figsize=(8, 12)) variations = { - f'Input | Gamma = {gamma} stretch': {'map': submap_193, 'stretch': PowerStretch(1 / gamma)}, - 'WOW | Linear stretch': {'map': default_wow, 'stretch': LinearStretch()}, - 'Denoised WOW': {'map': denoised_wow, 'stretch': LinearStretch()}, - 'Edge-aware WOW': {'map': bilateral_wow, 'stretch': LinearStretch()}, - 'Edge-aware & denoised WOW': {'map': bilateral_denoised_wow, 'stretch': LinearStretch()}, - 'Merged with input': {'map': bilateral_denoised_merged_wow, 'stretch': LinearStretch()} + f"Input | Gamma = {gamma} stretch": {"map": submap_193, "stretch": PowerStretch(1 / gamma)}, + "WOW | Linear stretch": {"map": default_wow, "stretch": LinearStretch()}, + "Denoised WOW": {"map": denoised_wow, "stretch": LinearStretch()}, + "Edge-aware WOW": {"map": bilateral_wow, "stretch": LinearStretch()}, + "Edge-aware & denoised WOW": {"map": bilateral_denoised_wow, "stretch": LinearStretch()}, + "Merged with input": {"map": bilateral_denoised_merged_wow, "stretch": LinearStretch()}, } interval = AsymmetricPercentileInterval(1, 99.9) for i, (title, image) in enumerate(variations.items()): - ax = fig.add_subplot(3, 2, i + 1, projection=image['map']) - image['map'].plot(norm=ImageNormalize(image['map'].data, interval=interval, stretch=image['stretch'])) + ax = fig.add_subplot(3, 2, i + 1, projection=image["map"]) + image["map"].plot(norm=ImageNormalize(image["map"].data, interval=interval, stretch=image["stretch"])) ax.set_title(title) - ax.axis('off') + ax.axis("off") fig.tight_layout() diff --git a/examples/aligning_aia_with_eis_maps.py b/examples/aligning_aia_with_eis_maps.py new file mode 100644 index 00000000..0cb778c1 --- /dev/null +++ b/examples/aligning_aia_with_eis_maps.py @@ -0,0 +1,92 @@ +""" +===================== +Coaligning EIS to AIA +===================== + +This example shows how to coalign EIS rasters to AIA images in order correct for the pointing +uncertainty in EIS. +""" + +import matplotlib.pyplot as plt + +import astropy.units as u +from astropy.visualization import AsinhStretch, ImageNormalize + +import sunpy.map + +from sunkit_image.coalignment import coalign_map + +################################################################################### +# For this example, we will use a preprocessed EIS raster image of the Fe XII +# 195.119 Å line. +# This raster image was prepared using the `eispac `__ package. + +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" +) + +################################################################################### +# Next, let's download the AIA data we will use as a reference image. +# We want AIA data near the beginning of the EIS raster and we'll use the +# 193 Å channel of AIA as it sees plasma at approximately the same temperature as +# the 195.119 Å line in our EIS raster. +# We have stored this file on Github so we can download it directly. + +aia_map = sunpy.map.Map( + "https://github.com/sunpy/data/raw/refs/heads/main/sunkit-image/aia.lev1.193A_2014_01_08T09_57_30.84Z.image_lev1.fits" +) + +#################################################################################### +# Before coaligning the images, we first resample the AIA image to the same plate +# scale as the EIS image. This will ensure better results from our coalignment. + +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 and AIA using cross-correlation. By default, this function +# uses the "match_template" method. For details of the implementation refer to the +# documentation of `skimage.feature.match_template`. + +coaligned_eis_map = coalign_map(eis_map, aia_downsampled) + +#################################################################################### +# 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. + +fig = plt.figure(figsize=(10, 7.5), layout="constrained") +ax = fig.add_subplot(121, projection=eis_map) +eis_map.plot( + axes=ax, + title="Before coalignment", + aspect=eis_map.scale.axis2 / eis_map.scale.axis1, + cmap="Blues_r", + norm=ImageNormalize(stretch=AsinhStretch()), +) +aia_map.draw_contours([800] * aia_map.unit, axes=ax) +ax = fig.add_subplot(122, projection=coaligned_eis_map, sharex=ax, sharey=ax) +coaligned_eis_map.plot( + axes=ax, + title="After coalignment", + aspect=coaligned_eis_map.scale.axis2 / coaligned_eis_map.scale.axis1, + cmap="Blues_r", + norm=ImageNormalize(stretch=AsinhStretch()), +) +aia_map.draw_contours([800] * aia_map.unit, axes=ax) + +#################################################################################### +# We can also visualize this shift by overlaying the extent of the two EIS +# maps on the AIA image we used to perform the coalignment. + +fig = plt.figure(figsize=(8, 8)) +ax = fig.add_subplot(projection=aia_map) +aia_map.plot(axes=ax) +eis_map.draw_extent(axes=ax, color="C0", label="Original") +coaligned_eis_map.draw_extent(axes=ax, color="C1", label="Coaligned") +ax.set_xlim(1700, 2500) +ax.set_ylim(1200, 2200) +ax.legend() + +plt.show() diff --git a/examples/radial_gradient_filters.py b/examples/radial_gradient_filters.py index 2eead64c..e4cadb99 100644 --- a/examples/radial_gradient_filters.py +++ b/examples/radial_gradient_filters.py @@ -48,12 +48,7 @@ # choose them according to your requirements. These can be changed by tweaking the following keywords: ``mean_attenuation_range`` and ``std_attenuation_range``. order = 20 -base_fnrgf = radial.fnrgf( - aia_map, - radial_bin_edges=radial_bin_edges, - order=order, - application_radius=1 * u.R_sun -) +base_fnrgf = radial.fnrgf(aia_map, radial_bin_edges=radial_bin_edges, order=order, application_radius=1 * u.R_sun) ########################################################################### # Now we will also use the final filter, RHEF. diff --git a/pyproject.toml b/pyproject.toml index 344acd77..d340e62c 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -45,8 +45,6 @@ watroo = ["watroo"] tests = [ "sunkit_image[all]", "dask", -# NumPy>=2.4 not compatible with some helpers in astropy<7.2 - "numpy<2.4", "pytest-astropy", "pytest-mpl", "pytest-xdist", diff --git a/sunkit_image/_dev/scm_version.py b/sunkit_image/_dev/scm_version.py index 988debf5..d0c4c1f1 100644 --- a/sunkit_image/_dev/scm_version.py +++ b/sunkit_image/_dev/scm_version.py @@ -5,8 +5,8 @@ try: from setuptools_scm import get_version - version = get_version(root=Path('../..'), relative_to=__file__) + version = get_version(root=Path("../.."), relative_to=__file__) except ImportError: raise except Exception as e: - raise ValueError('setuptools_scm can not determine version.') from e + raise ValueError("setuptools_scm can not determine version.") from e diff --git a/sunkit_image/asda.py b/sunkit_image/asda.py index 57f8e6e4..bb16210c 100644 --- a/sunkit_image/asda.py +++ b/sunkit_image/asda.py @@ -1,8 +1,3 @@ -""" -This module contains an implementation of the Automated Swirl Detection -Algorithm (ASDA). -""" - import warnings from itertools import product diff --git a/sunkit_image/coalignment.py b/sunkit_image/coalignment.py deleted file mode 100644 index 576e2720..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__ = [ - "apply_shifts", - "calculate_match_template_shift", - "mapsequence_coalign_by_match_template", - "match_template_to_layer", -] - - -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 `__. - 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..d56dcfc8 --- /dev/null +++ b/sunkit_image/coalignment/__init__.py @@ -0,0 +1,7 @@ +# Need to register the the coalignment functions +from sunkit_image.coalignment.match_template import match_template_coalign # isort:skip +from sunkit_image.coalignment.phase_cross_correlation import phase_cross_correlation_coalign # isort:skip + +from sunkit_image.coalignment.interface import coalign_map + +__all__ = ["coalign_map"] diff --git a/sunkit_image/coalignment/interface.py b/sunkit_image/coalignment/interface.py new file mode 100644 index 00000000..39b25e77 --- /dev/null +++ b/sunkit_image/coalignment/interface.py @@ -0,0 +1,230 @@ +import warnings +from dataclasses import dataclass + +import numpy as np + +import astropy +import astropy.units as u + +from sunpy import log +from sunpy.sun.models import differential_rotation +from sunpy.util.exceptions import SunpyUserWarning + +__all__ = [ + "AffineParams", + "coalign_map", + "update_map_metadata", +] + + +@dataclass +class AffineParams: + """ + A 2-element tuple containing scale values defining the image scaling along + the x and y axes. + """ + + scale: tuple[float, float] | list[float] | np.ndarray + """ + A 2x2 matrix defining the rotation transformation of the image. + """ + rotation_matrix: tuple[tuple[float, float], tuple[float, float]] | np.ndarray + """ + A 2-element tuple stores the translation of the image along the x and y + axes. + """ + translation: tuple[float, float] | list[float] | np.ndarray + + +def update_map_metadata(target_map, reference_map, affine_params): + """ + Update the metadata of a sunpy.map.Map` based on affine transformation + parameters. + + .. warning:: + + This function is currently only designed to update the reference coordinate + metadata based on the translation component of the affine transformation. + Changes to the rotation and scale metadata are not currently supported. + If you have a use case that requires changes to the rotation or scale metadata + `please open an issue on the issue tracker `__. + + Parameters + ---------- + target_map : `sunpy.map.Map` + The original map object whose metadata is to be updated. + reference_map : `sunpy.map.Map` + The reference map object to which the target map is to be coaligned. + affine_params : `NamedTuple` + A `NamedTuple` containing the affine transformation parameters. + + Returns + ------- + `sunpy.map.Map` + A new sunpy map object with updated metadata reflecting the affine transformation. + """ + # NOTE: Currently, the only metadata updates that are supported are shifts in + # the reference coordinate. Once other updates are supported, this check can be removed. + full_msg = ( + "\nIf you have a use case that requires changes to the rotation or scale metadata, " + "please open an issue at https://github.com/sunpy/sunkit-image/issues with details " + "about your use case and the type of metadata changes you require." + ) + if not (affine_params.rotation_matrix == np.eye(2)).all(): + raise NotImplementedError("Changes to the rotation metadata are currently not supported." + full_msg) + if not (affine_params.scale == np.array([1, 1])).all(): + raise NotImplementedError("Changes to the pixel scale metadata are currently not supported." + full_msg) + # Calculate the new reference pixel. + old_reference_pixel = u.Quantity(target_map.reference_pixel).to_value("pixel") + new_reference_pixel = ( + affine_params.scale * affine_params.rotation_matrix @ old_reference_pixel + affine_params.translation + ) + new_reference_coordinate = reference_map.wcs.pixel_to_world(*new_reference_pixel) + # Create a new map with the updated metadata + log.debug( + f"Shifting reference coordinate from {target_map.reference_coordinate} to {new_reference_coordinate} by {new_reference_coordinate.Tx - target_map.reference_coordinate.Tx}, {new_reference_coordinate.Ty - target_map.reference_coordinate.Ty}" + ) + return target_map.shift_reference_coord( + new_reference_coordinate.Tx - target_map.reference_coordinate.Tx, + new_reference_coordinate.Ty - target_map.reference_coordinate.Ty, + ) + + +def _warn_user_of_separation(target_map, reference_map): + """ + Issues a warning if the separation between the ``reference_map`` and + ``target_map`` is "large". + + Parameters + ---------- + target_map : `sunpy.map.Map` + The target map to be coaligned to the reference map. + reference_map : `sunpy.map.Map` + The reference map to which the target map is to be coaligned. + """ + # Maximum angular separation allowed between the reference and target maps + tolerable_angular_separation = 1 * u.deg + if astropy.__version__ >= "6.1.0": + angular_separation = reference_map.observer_coordinate.separation( + target_map.observer_coordinate, origin_mismatch="ignore" + ) + else: + angular_separation = reference_map.observer_coordinate.separation(target_map.observer_coordinate) + if angular_separation > tolerable_angular_separation: + warnings.warn( + "The angular separation between the observer coordinates of " + "the reference and target maps is large. This can lead to spurious " + "results when calculating shift between two images. Consider choosing " + "a reference map with a similar observer position to your target map.", + SunpyUserWarning, + stacklevel=3, + ) + # Calculate time difference and convert to separation angle + time_diff = np.abs((reference_map.date - target_map.date).to("s")) + time_separation_angle = differential_rotation( + time_diff, reference_map.center.transform_to("heliographic_stonyhurst").lat, model="howard" + ) + if time_separation_angle > tolerable_angular_separation: + warnings.warn( + "The difference in observation times of the reference and target maps is large. " + "This can lead to spurious results when calculating shifts between two images." + "Consider choosing a reference map with a similar observation time to your target map.", + SunpyUserWarning, + stacklevel=3, + ) + + +def _warn_user_of_plate_scale_difference(target_map, reference_map): + """ + Issues a warning if there is a plate scale difference between the + ``reference_map`` and ``target_map``. + + Parameters + ---------- + target_map : `sunpy.map.Map` + The target map to be coaligned to the reference map. + reference_map : `sunpy.map.Map` + The reference map to which the target map is to be coaligned. + """ + if not u.allclose(target_map.scale, reference_map.scale): + warnings.warn( + "The reference and target maps have different plate scales. " + "This can lead to spurious results when calculating the shift between two arrays. " + "Consider resampling the reference map to have the same plate scale as the target map.", + SunpyUserWarning, + stacklevel=3, + ) + + +def coalign_map(target_map, reference_map, method="match_template", **kwargs): + """ + Performs image coalignment using the specified method by updating the metadata of the + target map to align it with the reference map. + + .. note:: + + This function is intended to correct maps with inaccurate metadata. + It is not designed to correct for differential rotation or changes in observer + location. + + .. warning:: + + This function is currently only designed to update the reference coordinate + metadata based on the translation component of the affine transformation. + Changes to the rotation and scale metadata are not currently supported. + If you have a use case that requires changes to the rotation or scale metadata + `please open an issue on the issue tracker `__. + + .. note:: + + This 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 + ---------- + target_map : `sunpy.map.GenericMap` + The map to be coaligned. + reference_map : `sunpy.map.GenericMap` + The map to which the target map is to be coaligned. It is expected that the pointing data of this map + is accurate. For best results, ``reference_map`` and ``target_map`` should have approximately the + same observer location and observation time. For coalignment methods which do + not account for different pixel scales or rotations, it is recommended that ``reference_map`` + and ``target_map`` are resampled and/or rotated such that they have the same orientation and + plate scale. + method : {{{coalignment_function_names}}}, optional + The name of the registered coalignment method to use. + Defaults to `~sunkit_image.coalignment.match_template.match_template_coalign`. + kwargs : `dict` + Additional keyword arguments to pass to the registered method. + + Returns + ------- + `sunpy.map.GenericMap` + 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: {','.join(REGISTERED_METHODS.keys())}. " + "If you expect this method to be present, ensure the method has been registered with the correct import." + ) + raise ValueError(msg) + _warn_user_of_separation(target_map, reference_map) + _warn_user_of_plate_scale_difference(target_map, reference_map) + + affine_params = REGISTERED_METHODS[method](target_map.data, reference_map.data, **kwargs) + return update_map_metadata(target_map, reference_map, affine_params) + + +from sunkit_image.coalignment.register import REGISTERED_METHODS # isort:skip # NOQA: E402 + +# 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_map. We cannot use the add_common_docstring decorator +# due to what would be a circular loop in definitions. +coalign_map.__doc__ = coalign_map.__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..eca5530b --- /dev/null +++ b/sunkit_image/coalignment/match_template.py @@ -0,0 +1,131 @@ +import numpy as np +from skimage.feature import match_template + +from sunpy import log + +from sunkit_image.coalignment.register import 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(target_array, reference_array, **kwargs): + """ + Perform coalignment by matching the input array to the target array. + + Coalign ``target_array`` to ``reference_array`` using normalized correlation + via `skimage.feature.match_template`. For more details on this approach, please check + the documentation of that function including the available keyword arguments and the + details of the algorithm. + + .. warning:: + + This method only returns a translation. The scale and rotation + parameters are unity. + + Parameters + ---------- + target_array : `numpy.ndarray` + The input 2D array to be coaligned. + reference_array : `numpy.ndarray` + The template 2D array to align to. + kwargs: dict + Passed to `skimage.feature.match_template`. + + Returns + ------- + `sunkit_image.coalignment.interface.AffineParams` + """ + from sunkit_image.coalignment.interface import AffineParams # NOQA: PLC0415 + + corr = match_template(np.float64(reference_array), np.float64(target_array), **kwargs) + if len(corr) < 2: + raise ValueError( + "match_template returned an array with fewer than two values, and cannot be used for coalignment. " + "This implies the cross-correlation failed to find a match." + ) + # Find the best match location + y_shift, x_shift = _find_best_match_location(corr) + if kwargs.get("pad_input"): + # For pad_input=True matches correspond to the center + y_shift -= target_array.shape[0] // 2 + x_shift -= target_array.shape[1] // 2 + log.debug(f"Match template shift: x: {x_shift}, y: {y_shift}") + # 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/phase_cross_correlation.py b/sunkit_image/coalignment/phase_cross_correlation.py new file mode 100644 index 00000000..a16f14cb --- /dev/null +++ b/sunkit_image/coalignment/phase_cross_correlation.py @@ -0,0 +1,57 @@ +import numpy as np +from skimage.registration import phase_cross_correlation + +from sunpy import log + +from sunkit_image.coalignment.register import register_coalignment_method + +__all__ = ["phase_cross_correlation_coalign"] + + +@register_coalignment_method("phase_cross_correlation") +def phase_cross_correlation_coalign(target_array, reference_array, **kwargs): + """ + Perform coalignment by phase cross correlation input array to the target array. + + .. note:: This requires both the input and target arrays to be the same size. + + Coalign ``target_array`` to ``reference_array`` using phase cross-correlation + via `skimage.registration.phase_cross_correlation`. For more details on this approach, + please check the documentation of that function including the available keyword + arguments and the details of the algorithm. + + .. warning:: + + This method only returns a translation. The scale and rotation + parameters are unity. + + Parameters + ---------- + target_array : `numpy.ndarray` + The input 2D array to be coaligned. + reference_array : `numpy.ndarray` + The template 2D array to align to. + kwargs : dict + Passed to `skimage.registration.phase_cross_correlation`. + + Returns + ------- + `sunkit_image.coalignment.interface.AffineParams` + This method only returns a translation. The scale and rotation + parameters are unity. + """ + from sunkit_image.coalignment.interface import AffineParams # NOQA: PLC0415 + + if target_array.shape != reference_array.shape: + raise ValueError("Input and target arrays must be the same shape.") + shift, _, _ = phase_cross_correlation(reference_array, target_array, **kwargs) + # Shift has axis ordering which is consistent with the axis order of the input arrays, so y, x + # See the example here: + # https://scikit-image.org/docs/stable/auto_examples/registration/plot_register_translation.html + x_shift, y_shift = shift[1], shift[0] + log.debug(f"Phase cross correlation shift: x: {x_shift}, y: {y_shift}") + # 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/register.py b/sunkit_image/coalignment/register.py new file mode 100644 index 00000000..482d4a81 --- /dev/null +++ b/sunkit_image/coalignment/register.py @@ -0,0 +1,24 @@ +__all__ = [ + "REGISTERED_METHODS", + "register_coalignment_method", +] + +# Global Dictionary to store the registered methods and their names +REGISTERED_METHODS = {} + + +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 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..065fba13 --- /dev/null +++ b/sunkit_image/coalignment/tests/test_coalignment.py @@ -0,0 +1,186 @@ +import matplotlib.pyplot as plt +import numpy as np +import pytest + +import astropy.units as u +from astropy.coordinates import SkyCoord +from astropy.io import fits +from astropy.tests.helper import assert_quantity_allclose + +import sunpy.map +from sunpy.map.maputils import all_coordinates_from_map, coordinate_is_on_solar_disk +from sunpy.util.exceptions import SunpyUserWarning + +from sunkit_image.coalignment.interface import AffineParams, coalign_map, update_map_metadata +from sunkit_image.tests.helpers import figure_test + +POSITIVE_POINTING_SHIFT = [25, 50] * u.arcsec +NEGATIVE_POINTING_SHIFT = [-25, -50] * u.arcsec +MIXED_POINTING_SHIFT = [25, -50] * u.arcsec + + +@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.fixture() +def aia193_test_map(): + # This is matched to the EIS observation time + url = "https://github.com/sunpy/data/raw/refs/heads/main/sunkit-image/aia.lev1.193A_2014_01_08T09_57_30.84Z.image_lev1.fits" + with fits.open(url) as hdul: + hdul.verify("silentfix") + return sunpy.map.Map(hdul[1].data, hdul[1].header) + + +@pytest.fixture() +def incorrect_shifted_once_map(aia171_test_map): + return aia171_test_map.shift_reference_coord(30 * u.arcsec, -40 * u.arcsec) + + +@pytest.fixture() +def cutout_map(aia171_test_map): + bottom_left = SkyCoord(0 * u.arcsec, 0 * u.arcsec, frame=aia171_test_map.coordinate_frame) + top_right = SkyCoord(900 * u.arcsec, 900 * u.arcsec, frame=aia171_test_map.coordinate_frame) + return aia171_test_map.submap(bottom_left, top_right=top_right) + + +@pytest.fixture( + params=[ + POSITIVE_POINTING_SHIFT, + NEGATIVE_POINTING_SHIFT, + MIXED_POINTING_SHIFT, + ] +) +def incorrect_pointing_map_and_shift(request, aia171_test_map): + return aia171_test_map.shift_reference_coord(*request.param), request.param + + +@pytest.fixture( + params=[ + POSITIVE_POINTING_SHIFT, + NEGATIVE_POINTING_SHIFT, + MIXED_POINTING_SHIFT, + ] +) +def incorrect_pointing_cutout_map_and_shift(request, cutout_map): + return cutout_map.shift_reference_coord(*request.param), request.param + + +@pytest.mark.remote_data() +def test_coalignment_eis_aia(eis_test_map, aia193_test_map): + 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_map(eis_test_map, aia193_test_downsampled_map, method="match_template") + # Check that correction is as expected based on known pointing offset + assert_quantity_allclose( + eis_test_map.reference_coordinate.separation(coaligned_eis_map.reference_coordinate), + 5.95935177 * u.arcsec, + ) + + +def test_coalignment_match_template_full_map(incorrect_pointing_map_and_shift, aia171_test_map): + incorrect_pointing_map, pointing_shift = incorrect_pointing_map_and_shift + + # Crop out the array that is outside the solar disk to have something to align to + hpc_coords = all_coordinates_from_map(incorrect_pointing_map) + mask = coordinate_is_on_solar_disk(hpc_coords) + masked_map = sunpy.map.Map(np.abs(incorrect_pointing_map.data * mask), incorrect_pointing_map.meta) + + # We have to pad the input map because otherwise match_template cannot find a good match + fixed_map = coalign_map(masked_map, aia171_test_map, pad_input=True) + + assert_quantity_allclose( + (incorrect_pointing_map.reference_coordinate.Tx - fixed_map.reference_coordinate.Tx), + pointing_shift[0], + atol=1 * u.arcsec, + ) + assert_quantity_allclose( + (incorrect_pointing_map.reference_coordinate.Ty - fixed_map.reference_coordinate.Ty), + pointing_shift[1], + atol=1 * u.arcsec, + ) + + +def test_coalignment_match_template_cutout(incorrect_pointing_cutout_map_and_shift, aia171_test_map): + incorrect_pointing_cutout_map, pointing_shift = incorrect_pointing_cutout_map_and_shift + fixed_cutout_map = coalign_map(incorrect_pointing_cutout_map, aia171_test_map) + assert_quantity_allclose( + (incorrect_pointing_cutout_map.reference_coordinate.Tx - fixed_cutout_map.reference_coordinate.Tx), + pointing_shift[0], + atol=0.05 * u.arcsec, + ) + assert_quantity_allclose( + (incorrect_pointing_cutout_map.reference_coordinate.Ty - fixed_cutout_map.reference_coordinate.Ty), + pointing_shift[1], + atol=0.05 * u.arcsec, + ) + + +def test_coalign_phase_cross_correlation(incorrect_pointing_map_and_shift, aia171_test_map): + incorrect_pointing_map, pointing_shift = incorrect_pointing_map_and_shift + fixed_map = coalign_map(incorrect_pointing_map, aia171_test_map, method="phase_cross_correlation") + assert_quantity_allclose( + (incorrect_pointing_map.reference_coordinate.Tx - fixed_map.reference_coordinate.Tx), + pointing_shift[0], + atol=1 * u.arcsec, + ) + assert_quantity_allclose( + (incorrect_pointing_map.reference_coordinate.Ty - fixed_map.reference_coordinate.Ty), + pointing_shift[1], + atol=1 * u.arcsec, + ) + + +@figure_test +def test_coalignment_figure(incorrect_pointing_cutout_map_and_shift, cutout_map, aia171_test_map): + # This is three separate figure tests + levels = [200, 800] * cutout_map.unit + incorrect_pointing_cutout_map, _ = incorrect_pointing_cutout_map_and_shift + fixed_cutout_map = coalign_map(incorrect_pointing_cutout_map, aia171_test_map) + fig = plt.figure(figsize=(10, 7.5)) + # Messed up map + ax = fig.add_subplot(121, projection=incorrect_pointing_cutout_map) + incorrect_pointing_cutout_map.plot(axes=ax, title="Incorrect Pointing") + cutout_map.draw_contours(levels, axes=ax) + # After coalignment + ax = fig.add_subplot(122, projection=fixed_cutout_map) + fixed_cutout_map.plot(axes=ax, title="Fixed Pointing") + cutout_map.draw_contours(levels, axes=ax) + fig.tight_layout() + return fig + + +def test_unsupported_affine_parameters(incorrect_shifted_once_map, aia171_test_map): + affine_rot = AffineParams( + scale=[1, 1], + rotation_matrix=2 * np.eye(2), + translation=[0, 0], + ) + with pytest.raises(NotImplementedError, match=r"Changes to the rotation metadata are currently not supported."): + update_map_metadata(incorrect_shifted_once_map, aia171_test_map, affine_rot) + affine_scale = AffineParams( + scale=[2, 3], + rotation_matrix=np.eye(2), + translation=[0, 0], + ) + with pytest.raises(NotImplementedError, match=r"Changes to the pixel scale metadata are currently not supported."): + update_map_metadata(incorrect_shifted_once_map, aia171_test_map, affine_scale) + + +def test_warnings_coalign_map(incorrect_shifted_once_map, aia171_test_map): + time_shift_meta = incorrect_shifted_once_map.meta.copy() + time_shift_meta["DATE-OBS"] = "2014-01-08T09:57:30.84" + time_shift = sunpy.map.Map(incorrect_shifted_once_map.data, time_shift_meta) + with pytest.warns( + SunpyUserWarning, match=r"The difference in observation times of the reference and target maps is large." + ): + with pytest.raises(ValueError, match=r"match_template returned an array with fewer than two values,"): + coalign_map(time_shift, aia171_test_map) + + resampled_map = incorrect_shifted_once_map.resample([100, 100] * u.pix) + with pytest.warns(SunpyUserWarning, match=r"The reference and target maps have different plate scales."): + coalign_map(resampled_map, aia171_test_map) 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..77631cfe --- /dev/null +++ b/sunkit_image/coalignment/tests/test_match_template.py @@ -0,0 +1,19 @@ +import numpy as np +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 + + +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/coalignment/tests/test_register.py b/sunkit_image/coalignment/tests/test_register.py new file mode 100644 index 00000000..00589591 --- /dev/null +++ b/sunkit_image/coalignment/tests/test_register.py @@ -0,0 +1,16 @@ +from sunkit_image.coalignment import coalign_map +from sunkit_image.coalignment.register import REGISTERED_METHODS, register_coalignment_method + + +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" + + +def test_coalignment_method_in_docstring(): + assert "phase_cross_correlation" in coalign_map.__doc__ diff --git a/sunkit_image/conftest.py b/sunkit_image/conftest.py index e7af65a9..e22ee77b 100644 --- a/sunkit_image/conftest.py +++ b/sunkit_image/conftest.py @@ -15,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 @@ -107,7 +108,7 @@ def pytest_runtest_teardown(item): # Clear the pyplot figure stack if it is not empty after the test # You can see these log messages by passing "-o log_cli=true" to pytest on the command line if HAVE_MATPLOTLIB and plt.get_fignums(): - msg = f"Removing {len(plt.get_fignums())} pyplot figure(s) " f"left open by {item.name}" + msg = f"Removing {len(plt.get_fignums())} pyplot figure(s) left open by {item.name}" console_logger.info(msg) plt.close("all") @@ -173,6 +174,7 @@ def granule_minimap3(): ) return sunpy.map.GenericMap(arr, header) + @pytest.fixture(params=["array", "map"]) def aia_171(request): # VerifyWarning: Invalid 'BLANK' keyword in header. @@ -199,5 +201,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/enhance.py b/sunkit_image/enhance.py index b5198bf6..d9bc2d35 100644 --- a/sunkit_image/enhance.py +++ b/sunkit_image/enhance.py @@ -1,7 +1,3 @@ -""" -This module contains enhancement routines for solar physics data. -""" - import warnings import numpy as np @@ -92,7 +88,7 @@ def mgn( Sol Phys 289, 2945-2955, 2014 `doi:10.1007/s11207-014-0523-9 `__ """ - olderr = np.seterr(all='ignore') + olderr = np.seterr(all="ignore") if sigma is None: sigma = [1.25, 2.5, 5, 10, 20, 40] if np.isnan(data).any(): diff --git a/sunkit_image/granule.py b/sunkit_image/granule.py index 405a13cf..7cfffbc5 100644 --- a/sunkit_image/granule.py +++ b/sunkit_image/granule.py @@ -1,7 +1,3 @@ -""" -This module contains functions that will segment images for granule detection. -""" - import logging import matplotlib as mpl diff --git a/sunkit_image/radial.py b/sunkit_image/radial.py index 4dfa3271..b7abaad4 100644 --- a/sunkit_image/radial.py +++ b/sunkit_image/radial.py @@ -1,8 +1,3 @@ -""" -This module contains functions that can be used to enhance the regions above a -radius. -""" - import numpy as np from scipy import stats from tqdm import tqdm @@ -101,6 +96,7 @@ def _normalize_fit_radial_intensity(radii, polynomial, normalization_radius): polynomial, ) + def _select_rank_method(method): def _percentile_ranks_scipy(arr): mask = ~np.isnan(arr) @@ -543,7 +539,9 @@ def fnrgf( data = np.ones_like(smap.data) * fill # Set attenuation coefficients - attenuation_coefficients = _set_attenuation_coefficients(order, mean_attenuation_range, std_attenuation_range, cutoff) + attenuation_coefficients = _set_attenuation_coefficients( + order, mean_attenuation_range, std_attenuation_range, cutoff + ) # Iterate over each circular ring for i in tqdm(range(nbins), desc="FNRGF: ", disable=not progress): @@ -720,9 +718,7 @@ def rhef( # Loop over each radial bin to apply the filter for i in tqdm(range(radial_bin_edges.shape[1]), desc="RHEF: ", disable=not progress): # Identify pixels within the current radial bin - here = np.logical_and( - map_r >= radial_bin_edges[0, i].to(u.R_sun), map_r < radial_bin_edges[1, i].to(u.R_sun) - ) + here = np.logical_and(map_r >= radial_bin_edges[0, i].to(u.R_sun), map_r < radial_bin_edges[1, i].to(u.R_sun)) if application_radius is not None and application_radius > 0: here = np.logical_and(here, map_r >= application_radius) # Apply ranking function diff --git a/sunkit_image/stara.py b/sunkit_image/stara.py index d41e69f2..80a0dcce 100644 --- a/sunkit_image/stara.py +++ b/sunkit_image/stara.py @@ -1,8 +1,3 @@ -""" -This module contains an implementation of the Sunspot Tracking And Recognition -Algorithm (STARA). -""" - import numpy as np import skimage from skimage.filters import median diff --git a/sunkit_image/tests/figure_hashes_mpl_3100_ft_261_sunpy_700_astropy_710.json b/sunkit_image/tests/figure_hashes_mpl_3108_ft_261_sunpy_710_astropy_720.json similarity index 77% rename from sunkit_image/tests/figure_hashes_mpl_3100_ft_261_sunpy_700_astropy_710.json rename to sunkit_image/tests/figure_hashes_mpl_3108_ft_261_sunpy_710_astropy_720.json index 72148618..15f5b35f 100644 --- a/sunkit_image/tests/figure_hashes_mpl_3100_ft_261_sunpy_700_astropy_710.json +++ b/sunkit_image/tests/figure_hashes_mpl_3108_ft_261_sunpy_710_astropy_720.json @@ -1,4 +1,7 @@ { + "sunkit_image.coalignment.tests.test_coalignment.test_coalignment_figure[incorrect_pointing_cutout_map_and_shift0]": "ce40f5667e62443a8a48f7fac2ea320d159f6026572e7c879d7746da800ae17e", + "sunkit_image.coalignment.tests.test_coalignment.test_coalignment_figure[incorrect_pointing_cutout_map_and_shift1]": "2b8134d3d32cc38f0855b4aca96e8a04219bd4205def9af0e87877232029f634", + "sunkit_image.coalignment.tests.test_coalignment.test_coalignment_figure[incorrect_pointing_cutout_map_and_shift2]": "f2c3574c80d36e0e0cd8c365aae037c5e5ef5fe74a4232663d1b8bc8db705f81", "sunkit_image.tests.test_enhance.test_mgn[array]": "d3ed034e34d03288aad08c06447fde42a8eb5c7ce4c0543a2090bb04c75fad1a", "sunkit_image.tests.test_enhance.test_mgn[map]": "50bb490a6cc2408befe13c7f2a54f7433df80c2473dd21b619ace35de7e8f250", "sunkit_image.tests.test_enhance.test_mgn_submap": "cf3948d3ffa8ed4eacc500bee170db68bb1a96d2cec1c92e48f74c01827dc397", @@ -15,4 +18,4 @@ "sunkit_image.tests.test_trace.test_occult2_fig[array]": "e59f4c476b6b27788ab8dc397cb0d2f34f8d6032eced289fd2eabeb324b39558", "sunkit_image.tests.test_trace.test_occult2_fig[map]": "e59f4c476b6b27788ab8dc397cb0d2f34f8d6032eced289fd2eabeb324b39558", "sunkit_image.tests.test_trace.test_occult2_cutout": "1b18459111342c3bab4a2c8fb08b012eef4a69767e08753d8918fe4987fa165d" -} +} \ No newline at end of file diff --git a/sunkit_image/tests/test_coalignment.py b/sunkit_image/tests/test_coalignment.py deleted file mode 100644 index 4c0e1c68..00000000 --- a/sunkit_image/tests/test_coalignment.py +++ /dev/null @@ -1,412 +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=r"The layer image has nonfinite entries."): - _check_for_nonfinite_entries(b, np.ones((3, 3))) - - with pytest.warns(SunpyUserWarning, match=r"The template image has nonfinite entries."): - _check_for_nonfinite_entries(np.ones((3, 3)), b) - - with pytest.warns(SunpyUserWarning, match=r"The layer image has nonfinite entries."): - with pytest.warns(SunpyUserWarning, match=r"The template image has nonfinite entries."): - _check_for_nonfinite_entries(b, b) - - -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=r"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=r"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=r"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])} diff --git a/sunkit_image/tests/test_enhance.py b/sunkit_image/tests/test_enhance.py index fedcfce8..2566954a 100644 --- a/sunkit_image/tests/test_enhance.py +++ b/sunkit_image/tests/test_enhance.py @@ -10,7 +10,10 @@ import sunkit_image.enhance as enhance from sunkit_image.tests.helpers import figure_test -pytestmark = [pytest.mark.filterwarnings("ignore:Missing metadata for observer"), pytest.mark.filterwarnings("ignore:Missing metadata for observation time")] +pytestmark = [ + pytest.mark.filterwarnings("ignore:Missing metadata for observer"), + pytest.mark.filterwarnings("ignore:Missing metadata for observation time"), +] @figure_test diff --git a/sunkit_image/tests/test_radial.py b/sunkit_image/tests/test_radial.py index c1e90ccd..98a1fb41 100644 --- a/sunkit_image/tests/test_radial.py +++ b/sunkit_image/tests/test_radial.py @@ -14,7 +14,10 @@ import sunkit_image.utils as utils from sunkit_image.tests.helpers import figure_test, skip_windows -pytestmark = [pytest.mark.filterwarnings("ignore:Missing metadata for observer"), pytest.mark.filterwarnings("ignore:Missing metadata for observation time")] +pytestmark = [ + pytest.mark.filterwarnings("ignore:Missing metadata for observer"), + pytest.mark.filterwarnings("ignore:Missing metadata for observation time"), +] @pytest.fixture() @@ -174,6 +177,7 @@ def test_fnrgf_errors(map_test1): cutoff=0, ) + @figure_test @pytest.mark.remote_data() def test_fig_nrgf(aia_171_map): @@ -189,7 +193,14 @@ def test_fig_fnrgf(aia_171_map): radial_bin_edges = utils.equally_spaced_bins() radial_bin_edges *= u.R_sun order = 20 - out = rad.fnrgf(aia_171_map, radial_bin_edges=radial_bin_edges, order=order, mean_attenuation_range=[1.0, 0.0], std_attenuation_range=[1.0, 0.0], cutoff=0) + out = rad.fnrgf( + aia_171_map, + radial_bin_edges=radial_bin_edges, + order=order, + mean_attenuation_range=[1.0, 0.0], + std_attenuation_range=[1.0, 0.0], + cutoff=0, + ) out.plot() @@ -240,6 +251,7 @@ def test_multifig_rhef(aia_171_map): return fig + def test_set_attenuation_coefficients(): order = 1 # Hand calculated diff --git a/sunkit_image/tests/test_stara.py b/sunkit_image/tests/test_stara.py index 5f068da5..c31dec5f 100644 --- a/sunkit_image/tests/test_stara.py +++ b/sunkit_image/tests/test_stara.py @@ -26,9 +26,9 @@ def test_stara_threshold_adjustment(hmi_map): higher_limb_filtered_result = stara(hmi_upscaled, limb_filter=30 * u.percent) # Assert that the lower threshold results in more features detected assert lower_threshold_result.sum() > higher_threshold_result.sum(), "Lower threshold should detect more features" - assert ( - lower_limb_filtered_result.sum() > higher_limb_filtered_result.sum() - ), "Lower Limb filter should detect more features" + assert lower_limb_filtered_result.sum() > higher_limb_filtered_result.sum(), ( + "Lower Limb filter should detect more features" + ) @figure_test diff --git a/sunkit_image/time_lag.py b/sunkit_image/time_lag.py index c4c7404d..e81a213b 100644 --- a/sunkit_image/time_lag.py +++ b/sunkit_image/time_lag.py @@ -1,10 +1,3 @@ -""" -This module contains functions for calculating the cross-correlation and time -lag between intensity cubes. - -Useful for understanding time variability in EUV light curves. -""" - import numpy as np import astropy.units as u diff --git a/sunkit_image/trace.py b/sunkit_image/trace.py index c15be25c..7e019688 100644 --- a/sunkit_image/trace.py +++ b/sunkit_image/trace.py @@ -1,8 +1,3 @@ -""" -This module contains functions that will the trace out coronal loop-like -structures in an image. -""" - import numpy as np from scipy import interpolate diff --git a/sunkit_image/utils/utils.py b/sunkit_image/utils/utils.py index b3126b0d..2efb8c38 100644 --- a/sunkit_image/utils/utils.py +++ b/sunkit_image/utils/utils.py @@ -24,7 +24,7 @@ "get_radial_intensity_summary", "points_in_poly", "reform2d", - "remove_duplicate" + "remove_duplicate", ] @@ -195,9 +195,7 @@ def reform2d(array, factor=1): msg = "Input array must be 2d!" raise ValueError(msg) if factor > 1: - congridx = RectBivariateSpline( - np.arange(0, array.shape[0]), np.arange(0, array.shape[1]), array, kx=1, ky=1 - ) + congridx = RectBivariateSpline(np.arange(0, array.shape[0]), np.arange(0, array.shape[1]), array, kx=1, ky=1) return congridx(np.arange(0, array.shape[0], 1 / factor), np.arange(0, array.shape[1], 1 / factor)) return array diff --git a/sunkit_image/version.py b/sunkit_image/version.py index 98de46e0..b6f8a9ab 100644 --- a/sunkit_image/version.py +++ b/sunkit_image/version.py @@ -9,11 +9,9 @@ except Exception: # NOQA: BLE001 import warnings - warnings.warn( - f'could not determine {__name__.split(".")[0]} package version; this indicates a broken installation' - ) + warnings.warn(f"could not determine {__name__.split('.')[0]} package version; this indicates a broken installation") del warnings - version = '0.0.0' + version = "0.0.0" __all__ = ["version"] diff --git a/tox.ini b/tox.ini index c82f5256..e479870d 100644 --- a/tox.ini +++ b/tox.ini @@ -49,9 +49,9 @@ deps = # Handle minimum dependencies via minimum_dependencies oldestdeps: minimum_dependencies # Figure tests need a tightly controlled environment - figure-!devdeps: astropy==7.1.0 - figure-!devdeps: matplotlib==3.10.0 - figure-!devdeps: sunpy==7.0.0 + figure-!devdeps: astropy==7.2.0 + figure-!devdeps: matplotlib==3.10.8 + figure-!devdeps: sunpy==7.1.0 # The following indicates which extras_require will be installed extras = all