diff --git a/CHANGES.rst b/CHANGES.rst index 80fe9a14..847723a6 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -6,6 +6,9 @@ New Features - Added a ``disp_bounds`` argument to ``tracing.FitTrace``. The argument allows for adjusting the dispersion-axis window from which the trace peaks are estimated. +- Added a new ``specreduce.wavecal1d.WavelengthCalibration1D`` class for one-dimensional wavelength + calibration. The old ``specreduce.wavelength_calibration.WavelengthCalibration1D`` is + deprecated and will be removed in v. 2.0. 1.6.0 (2025-06-18) ------------------ diff --git a/docs/_static/specreduce.css b/docs/_static/specreduce.css index 229d8787..e69de29b 100644 --- a/docs/_static/specreduce.css +++ b/docs/_static/specreduce.css @@ -1,15 +0,0 @@ -@import url("bootstrap-astropy.css"); - -div.topbar a.brand { - background: transparent url("logo_icon.png") no-repeat 8px 3px; - background-image: url("logo_icon.png"), none; - background-size: 32px 32px; -} - -#logotext1 { - color: #519EA8; -} - -#logotext2 { - color: #FF5000; -} diff --git a/docs/api.rst b/docs/api.rst index e4b4de6b..d0f52177 100644 --- a/docs/api.rst +++ b/docs/api.rst @@ -1,7 +1,7 @@ .. _api_index: -API Index -========= +API Reference +============= .. automodapi:: specreduce :no-inheritance-diagram: @@ -21,6 +21,12 @@ API Index .. automodapi:: specreduce.extract :no-inheritance-diagram: +.. automodapi:: specreduce.wavecal1d + :no-inheritance-diagram: + +.. automodapi:: specreduce.wavesol1d + :no-inheritance-diagram: + .. automodapi:: specreduce.calibration_data :no-inheritance-diagram: :include-all-objects: diff --git a/docs/background.rst b/docs/background.rst new file mode 100644 index 00000000..28dc7c46 --- /dev/null +++ b/docs/background.rst @@ -0,0 +1,29 @@ +Background correction +===================== + +The `specreduce.background` module generates and subtracts a background image from +the input 2D spectral image. The `~specreduce.background.Background` object is +defined by one or more windows, and can be generated with: + +* `~specreduce.background.Background` +* `Background.one_sided ` +* `Background.two_sided ` + +The center of the window can either be passed as a float/integer or as a trace + +.. code-block:: python + + bg = specreduce.background.Background.one_sided(image, trace, separation=5, width=2) + +or, equivalently + +.. code-block:: python + + bg = specreduce.background.Background.one_sided(image, 15, separation=5, width=2) + +The background image can be accessed via `~specreduce.background.Background.bkg_image` +and the background-subtracted image via `~specreduce.background.Background.sub_image` +(or ``image - bg``). + +The background and trace steps can be done iteratively, to refine an automated +trace using the background-subtracted image as input. \ No newline at end of file diff --git a/docs/conf.py b/docs/conf.py index cf78023f..17afccf5 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -27,21 +27,20 @@ import sys import datetime - import sphinx from specreduce import __version__ try: - from sphinx_astropy.conf.v1 import * # noqa + from sphinx_astropy.conf.v2 import * # noqa + from sphinx_astropy.conf.v2 import extensions # noqa except ImportError: - print('ERROR: the documentation requires the sphinx-astropy package to be installed') + print("ERROR: the documentation requires the sphinx-astropy package to be installed") sys.exit(1) # xref: https://github.com/sphinx-doc/sphinx/issues/13232#issuecomment-2608708175 if sys.version_info[:2] >= (3, 13) and sphinx.version_info[:2] < (8, 2): import pathlib - from sphinx.util.typing import _INVALID_BUILTIN_CLASSES _INVALID_BUILTIN_CLASSES[pathlib.Path] = "pathlib.Path" @@ -49,10 +48,10 @@ # -- General configuration ---------------------------------------------------- # By default, highlight as Python 3. -highlight_language = 'python3' +highlight_language = "python3" # If your documentation needs a minimal Sphinx version, state it here. -#needs_sphinx = '1.2' +needs_sphinx = "3.0" # To perform a Sphinx version check that needs to be more specific than # major.minor, call `check_sphinx_version("x.y.z")` here. @@ -60,83 +59,94 @@ # List of patterns, relative to source directory, that match files and # directories to ignore when looking for source files. -exclude_patterns.append('_templates') +exclude_patterns.append("_templates") # This is added to the end of RST files - a good place to put substitutions to # be used globally. -rst_epilog += """ -""" +#rst_epilog += """ +#.. _Astropy: https://www.astropy.org/ +#""" + +extensions.extend( + [ + "sphinx_design", + "nbsphinx", + ] +) # -- Project information ------------------------------------------------------ # This does not *have* to match the package name, but typically does project = "specreduce" author = "Astropy Specreduce contributors" -copyright = '{0}, {1}'.format( - datetime.datetime.now().year, author) +copyright = "{0}, {1}".format(datetime.datetime.now().year, author) # The version info for the project you're documenting, acts as replacement for # |version| and |release|, also used in various other places throughout the # built documents. # The short X.Y version. -version = __version__.split('-', 1)[0] +version = __version__.split("-", 1)[0] # The full version, including alpha/beta/rc tags. release = __version__ # -- Options for HTML output -------------------------------------------------- -# A NOTE ON HTML THEMES -# The global astropy configuration uses a custom theme, 'bootstrap-astropy', -# which is installed along with astropy. A different theme can be used or -# the options for this theme can be modified by overriding some of the -# variables set in the global configuration. The variables set in the -# global configuration are listed below, commented out. - - -# Add any paths that contain custom themes here, relative to this directory. -# To use a different custom theme, add the directory containing the theme. -#html_theme_path = [] - -# The theme to use for HTML and HTML Help pages. See the documentation for -# a list of builtin themes. To override the custom theme, set this to the -# name of a builtin theme or the name of a custom theme in html_theme_path. -#html_theme = None -html_static_path = ['_static'] # html_theme = None -html_style = 'specreduce.css' +html_static_path = ["_static"] # html_theme = None +html_style = "specreduce.css" - -html_theme_options = { - 'logotext1': 'spec', # white, semi-bold - 'logotext2': 'reduce', # orange, light - 'logotext3': ':docs' # white, light +html_theme_options.update( + { + "github_url": "https://github.com/astropy/specreduce", + "use_edit_page_button": False, + "navigation_with_keys": False, + "logo": { + "text": f"{project}", + "image_light": "_static/logo_icon.png", + "image_dark": "_static/logo_icon.png", + }, + "secondary_sidebar_items": {"**": ["page-toc"], "index": []}, } +) + +html_context = { + "default_mode": "light", + "version_slug": os.environ.get("READTHEDOCS_VERSION") or "", + "to_be_indexed": ["stable", "latest"], + "github_user": "astropy", + "github_repo": "specreduce", + "github_version": "main", + "doc_path": "docs", + "edit_page_url_template": "{{ astropy_custom_edit_url(github_user, github_repo, github_version, doc_path, file_name, default_edit_page_url_template) }}", + "default_edit_page_url_template": "https://github.com/{github_user}/{github_repo}/edit/{github_version}/{doc_path}{file_name}", + # Tell Jinja2 templates the build is running on Read the Docs + "READTHEDOCS": os.environ.get("READTHEDOCS", "") == "True", +} # Custom sidebar templates, maps document names to template names. -#html_sidebars = {} -html_sidebars['**'] = ['localtoc.html'] -html_sidebars['index'] = ['globaltoc.html', 'localtoc.html'] +html_sidebars = {} +html_sidebars['index'] = [] +html_sidebars["contributing"] = [] -# The name of an image file (relative to this directory) to place at the top -# of the sidebar. -#html_logo = '' +# html_sidebars['**'] = ['localtoc.html'] +# html_sidebars['index'] = [] #['globaltoc.html', 'localtoc.html'] # The name of an image file (within the static path) to use as favicon of the # docs. This file should be a Windows icon file (.ico) being 16x16 or 32x32 # pixels large. -html_favicon = '_static/logo_icon.ico' +html_favicon = "_static/logo_icon.ico" # If not '', a 'Last updated on:' timestamp is inserted at every page bottom, # using the given strftime format. -#html_last_updated_fmt = '' +# html_last_updated_fmt = '' # The name for this set of Sphinx documents. If None, it defaults to # " v documentation". -html_title = '{0} v{1}'.format(project, release) +html_title = "{0} v{1}".format(project, release) # Output file base name for HTML help builder. -htmlhelp_basename = project + 'doc' +htmlhelp_basename = project + "doc" # Prefixes that are ignored for sorting the Python module index modindex_common_prefix = ["specreduce."] @@ -145,16 +155,17 @@ # Grouping the document tree into LaTeX files. List of tuples # (source start file, target name, title, author, documentclass [howto/manual]). -latex_documents = [('index', project + '.tex', project + u' Documentation', - author, 'manual')] +latex_documents = [("index", project + ".tex", project + " Documentation", author, "manual")] # -- Options for manual page output ------------------------------------------- # One entry per manual page. List of tuples # (source start file, name, description, authors, manual section). -man_pages = [('index', project.lower(), project + u' Documentation', - [author], 1)] +man_pages = [("index", project.lower(), project + " Documentation", [author], 1)] + +# -- Options for numpydoc extension ------------------------------------------- +numpydoc_xref_param_type = False # -- Options for the edit_on_github extension --------------------------------- @@ -164,10 +175,10 @@ nitpicky = True intersphinx_mapping.update( { - 'astropy': ('https://docs.astropy.org/en/stable/', None), - 'ccdproc': ('https://ccdproc.readthedocs.io/en/stable/', None), - 'specutils': ('https://specutils.readthedocs.io/en/stable/', None), - 'gwcs': ('https://gwcs.readthedocs.io/en/stable/', None) + "astropy": ("https://docs.astropy.org/en/stable/", None), + "ccdproc": ("https://ccdproc.readthedocs.io/en/stable/", None), + "specutils": ("https://specutils.readthedocs.io/en/stable/", None), + "gwcs": ("https://gwcs.readthedocs.io/en/stable/", None), } ) # diff --git a/docs/contributing.rst b/docs/contributing.rst new file mode 100644 index 00000000..44bee7fc --- /dev/null +++ b/docs/contributing.rst @@ -0,0 +1,106 @@ +.. _contributors_guide: + +Contributing +============ + +As with all `Astropy `_ coordinated packages, **Specreduce** is +developed *by* the astrophysics community, *for* the astrophysics community. We warmly welcome +contributions of all kinds, whether it’s sharing feedback, reporting bugs, suggesting new +features, or improving code and documentation. Every contribution, big or small, helps make +Specreduce more useful and reliable for everyone. + +Specreduce follows the `Astropy development and contribution guidelines +`_ and the +`Astropy Community Code of Conduct `_. +By contributing, you agree to uphold these community standards, which help ensure that the +project remains a welcoming and inclusive space for all contributors. + +Getting Started +--------------- + +If you’re new to open-source development or to the Astropy ecosystem, do not worry. The best place +to start is by visiting the +`Specreduce GitHub repository `_, +where you can find current issues, pull requests, and the latest development activity. + +Before you begin contributing code, you may want to: + +* Read through the Astropy developer documentation linked above. +* Explore the existing documentation and tutorials to get a sense of how Specreduce works. +* Comment on an issue to let others know you’d like to work on it — or open a new issue to + discuss an idea or feature you’d like to propose. + +Roadmap +------- + +.. image:: + roadmap.png + :align: center + :width: 100% + +Contribute feedback +------------------- + +The Specreduce team values feedback from all users. If something doesn’t work as expected, +if you encounter a bug, or if you have an idea for an improvement, please let us know by opening +a new issue on the +`Specreduce GitHub issue tracker `_. + +For bug reports, please include: + +* A short description of the problem. +* The version of Specreduce and Python you are using. +* A minimal example (if possible) that reproduces the issue. + +For feature requests or usability feedback, describe your idea clearly and explain how it would +improve the user experience. Even short notes are valuable and help guide development priorities. + + +Contribute Code or Documentation +-------------------------------- + +If you see an open issue you’d like to work on, or if you’ve spotted an error or missing detail +in the documentation, you can submit your own improvements through GitHub. +To get started: + +1. Fork the Specreduce repository and create a new branch for your changes. +2. Make your edits or additions, whether it’s a bug fix, a new feature, or a documentation update. +3. Commit your changes with a clear message describing what you did. +4. Open a pull request to the main repository. + +.. If you’re unsure how to start, you can look for issues labeled +.. ``good first issue`` or ``help wanted`` — these are designed to be approachable even for +.. first-time contributors. + +Contribute Tutorials and Examples +--------------------------------- + +Tutorials and worked examples are among the most valuable contributions. They help other users +learn from real data and see how different tools fit together in practice. + +While the main steps of spectroscopic reduction (for example, tracing, extraction, and wavelength +calibration) are similar across most instruments, the best workflow can still depend on the +setup and science goals. In the long term, we aim to build a library of example reduction recipes +for different instruments that users can adapt for their own observations. + +If you have a reduction example, a notebook, or a teaching resource that might help others, +please share it, either by opening a pull request or by discussing it in an issue first. +We’re happy to help with formatting and integration into the documentation. + +Staying in Touch +---------------- + +Development discussions happen mainly on GitHub, but you can also connect with the wider Astropy +community through the `Astropy Discussion Forum `_, +where you can ask questions, share ideas, and get advice from other developers and users. + +Your contributions help make Specreduce and the Astropy ecosystem stronger and more sustainable. +Whether you fix a typo, improve a function, or share a new example, you are helping build a tool +that benefits the entire astrophysics community. Thank you for being part of it! + +.. toctree:: + :maxdepth: 1 + :hidden: + + process/index + terms \ No newline at end of file diff --git a/docs/extraction_quickstart.rst b/docs/extraction.rst similarity index 50% rename from docs/extraction_quickstart.rst rename to docs/extraction.rst index febd4198..de1e6d50 100644 --- a/docs/extraction_quickstart.rst +++ b/docs/extraction.rst @@ -1,69 +1,5 @@ -.. _extraction_quickstart: - -=============================== -Spectral Extraction Quick Start -=============================== - -Specreduce provides flexible functionality for extracting a 1D spectrum from a -2D spectral image, including steps for determining the trace of a spectrum, -background subtraction, and extraction. - - -Tracing -======= - -The `specreduce.tracing` module defines the trace of a spectrum on the 2D image. -These traces can either be determined semi-automatically or manually, and are -provided as the inputs for the remaining steps of the extraction process. -Supported trace types include: - -* `~specreduce.tracing.ArrayTrace` -* `~specreduce.tracing.FlatTrace` -* `~specreduce.tracing.FitTrace` - - -Each of these trace classes takes the 2D spectral image as input, as well as -additional information needed to define or determine the trace (see the API docs -above for required parameters for each of the available trace classes):: - - trace = specreduce.tracing.FlatTrace(image, 15) - -.. note:: - The fit for `~specreduce.tracing.FitTrace` may be adversely affected by noise where the spectrum - is faint. Narrowing the window parameter or lowering the order of the fitting function may - improve the result for noisy data. - - -Background -========== - -The `specreduce.background` module generates and subtracts a background image from -the input 2D spectral image. The `~specreduce.background.Background` object is -defined by one or more windows, and can be generated with: - -* `~specreduce.background.Background` -* `Background.one_sided ` -* `Background.two_sided ` - -The center of the window can either be passed as a float/integer or as a trace:: - - bg = specreduce.background.Background.one_sided(image, trace, separation=5, width=2) - - -or, equivalently:: - - bg = specreduce.background.Background.one_sided(image, 15, separation=5, width=2) - - -The background image can be accessed via `~specreduce.background.Background.bkg_image` -and the background-subtracted image via `~specreduce.background.Background.sub_image` -(or ``image - bg``). - -The background and trace steps can be done iteratively, to refine an automated -trace using the background-subtracted image as input. - -Extraction -========== +Spectrum Extraction +=================== The `specreduce.extract` module extracts a 1D spectrum from an input 2D spectrum (likely a background-extracted spectrum from the previous step) and a defined @@ -72,24 +8,30 @@ window, using one of the following implemented methods: * `~specreduce.extract.BoxcarExtract` * `~specreduce.extract.HorneExtract` -Each of these takes the input image and trace as inputs (see the API above for -other required and optional parameters):: +Each of these takes the input image and trace as inputs (see the :ref:`api_index` for +other required and optional parameters) + +.. code-block:: python extract = specreduce.extract.BoxcarExtract(image-bg, trace, width=3) -or:: +or + +.. code-block:: python extract = specreduce.extract.HorneExtract(image-bg, trace) For the Horne algorithm, the variance array is required. If the input image is -an ``astropy.NDData`` object with ``image.uncertainty`` provided, -then this will be used. Otherwise, the ``variance`` parameter must be set.:: +an `~astropy.nddata.NDData` object with ``image.uncertainty`` provided, +then this will be used. Otherwise, the ``variance`` parameter must be set. + +.. code-block:: python extract = specreduce.extract.HorneExtract(image-bg, trace, variance=var_array) -An optional mask array for the image may be supplied to HorneExtract as well. +An optional mask array for the image may be supplied to HorneExtract as well. This follows the same convention and can either be attached to ``image`` if it -is an ``astropy.NDData`` object, or supplied as a keyword argument. +is an `~astropy.nddata.NDData` object, or supplied as a keyword argument. The extraction methods automatically detect non-finite pixels in the input image and combine them with the user-supplied mask to prevent them from biasing the @@ -102,17 +44,24 @@ differently as documented in the API. The previous examples in this section show how to initialize the BoxcarExtract or HorneExtract objects with their required parameters. To extract the 1D -spectrum:: +spectrum + +.. code-block:: python spectrum = extract.spectrum The ``extract`` object contains all the set options. The extracted 1D spectrum can be accessed via the ``spectrum`` property or by calling (e.g ``extract()``) -the ``extract`` object (which also allows temporarily overriding any values):: +the ``extract`` object (which also allows temporarily overriding any values) + +.. code-block:: python spectrum2 = extract(width=6) -or, for example to override the original ``trace_object``:: +or, for example to override the original ``trace_object`` + +.. code-block:: python + spectrum2 = extract(trace_object=new_trace) Spatial profile options @@ -126,7 +75,6 @@ supplied as well (default is a 2D Polynomial) to account for residual background in the spatial profile. This option is not supported for ``interpolated_profile``. - If the ``interpolated_profile`` option is used, the image will be sampled in various wavelength bins (set by ``n_bins_interpolated_profile``), averaged in those bins, and samples are then interpolated between (linear by default, interpolation degree can @@ -138,38 +86,25 @@ and to use the defaults for bins and interpolation degree, or to override these defaults a dictionary can be passed in. For example, to use the ``interpolated_profile`` option with default bins and -interpolation degree:: +interpolation degree - interp_profile_extraction = extract(spatial_profile='interpolated_profile') +.. code-block:: python -Or, to override the default of 10 samples and use 20 samples:: + interp_profile_extraction = extract(spatial_profile='interpolated_profile') - interp_profile_extraction = extract(spatial_profile={'name': 'interpolated_profile', - 'n_bins_interpolated_profile': 20) +Or, to override the default of 10 samples and use 20 samples -Or, to do a cubic interpolation instead of the default linear:: +.. code-block:: python interp_profile_extraction = extract(spatial_profile={'name': 'interpolated_profile', - 'interp_degree_interpolated_profile': 3) - -As usual, parameters can either be set when instantiating the HorneExtraxt object, -or supplied/overridden when calling the extraction method on that object. - -Example Workflow -================ + 'n_bins_interpolated_profile': 20) -This will produce a 1D spectrum, with flux in units of the 2D spectrum. The -wavelength units will be pixels. Wavelength and flux calibration steps are not -included here. +Or, to do a cubic interpolation instead of the default linear -Putting all these steps together, a simple extraction process might look -something like:: +.. code-block:: python - from specreduce.tracing import FlatTrace - from specreduce.background import Background - from specreduce.extract import BoxcarExtract + interp_profile_extraction = extract(spatial_profile={'name': 'interpolated_profile', + 'interp_degree_interpolated_profile': 3) - trace = FlatTrace(image, 15) - bg = Background.two_sided(image, trace, separation=5, width=2) - extract = BoxcarExtract(image-bg, trace, width=3) - spectrum = extract.spectrum +As usual, parameters can either be set when instantiating the HorneExtraxt object, +or supplied/overridden when calling the extraction method on that object. \ No newline at end of file diff --git a/docs/getting_started/index.rst b/docs/getting_started/index.rst new file mode 100644 index 00000000..bf4c86cb --- /dev/null +++ b/docs/getting_started/index.rst @@ -0,0 +1,10 @@ +.. _getting-started: + +Getting started +=============== + +.. toctree:: + :maxdepth: 2 + + installation.rst + quickstart.ipynb \ No newline at end of file diff --git a/docs/getting_started/installation.rst b/docs/getting_started/installation.rst new file mode 100644 index 00000000..dd531292 --- /dev/null +++ b/docs/getting_started/installation.rst @@ -0,0 +1,24 @@ +.. _installation: + +Installation +============ + +Specreduce can be installed from PyPI using ``pip`` + +.. code-block:: console + + pip install specreduce + +from the conda-forge repository using ``conda`` + +.. code-block:: console + + conda install conda-forge::specreduce + +or by cloning the repository from `GitHub `_ + +.. code-block:: console + + git clone https://github.com/astropy/specreduce.git + cd specreduce + pip install . \ No newline at end of file diff --git a/docs/getting_started/quickstart.ipynb b/docs/getting_started/quickstart.ipynb new file mode 100644 index 00000000..b1d10a43 --- /dev/null +++ b/docs/getting_started/quickstart.ipynb @@ -0,0 +1,754 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "d43db469-374f-494f-be9d-aa15caf30045", + "metadata": {}, + "source": [ + "# Spectrum Extraction Quickstart\n", + "\n", + "Specreduce provides a flexible toolset for extracting a 1D spectrum from a\n", + "2D spectral image, including steps for determining the trace of a spectrum,\n", + "background subtraction, extraction, and wavelength calibration.\n", + "\n", + "This quickstart covers the basic spectrum extraction and calibration steps for\n", + "a synthetic science spectrum accompanied by He and Ne arc-lamp\n", + "spectra for wavelength calibration.\n", + "\n", + "This tutorial is written as a Jupyter Notebook. You can run it\n", + "by copying it (and the `quickstart.py` module) from GitHub,\n", + "or from your own specreduce installation under the `docs/getting_started` directory.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1d462ab0-6bc7-4e8b-9bbb-702b43a39241", + "metadata": { + "ExecuteTime": { + "end_time": "2025-10-21T20:37:50.961773Z", + "start_time": "2025-10-21T20:37:49.431235Z" + } + }, + "outputs": [], + "source": [ + "import warnings\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "import astropy.units as u\n", + "\n", + "from astropy.modeling import models\n", + "\n", + "from specreduce.tracing import FitTrace\n", + "from specreduce.background import Background\n", + "from specreduce.extract import BoxcarExtract, HorneExtract\n", + "from specreduce.wavecal1d import WavelengthCalibration1D\n", + "\n", + "from quickstart import make_science_and_arcs, plot_2d_spectrum\n", + "\n", + "plt.rcParams['figure.constrained_layout.use'] = True" + ] + }, + { + "cell_type": "markdown", + "id": "1ba4ac1f-0742-4a0b-81dc-6197f82a9c17", + "metadata": {}, + "source": [ + "## Data preparation\n", + "\n", + "First, we read in our science and arc spectra. We use synthetic data\n", + "created by the `quickstart.make_science_and_arcs` utility function, which itself calls\n", + "`specreduce.utils.synth_data.make_2d_arc_image` and `specreduce.utils.synth_data.make_2d_spec_image`.\n", + "The science frame (`sci`) and the two arc frames (`arc_he` and `arc_ne`)\n", + "are returned as `astropy.nddata.CCDData` objects with `astropy.nddata.StdDevUncertainty`\n", + "uncertainties.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2e616221-6e7d-426d-8fe0-9568845d27f4", + "metadata": { + "ExecuteTime": { + "end_time": "2025-10-21T21:13:17.955893Z", + "start_time": "2025-10-21T21:13:17.049304Z" + } + }, + "outputs": [], + "source": [ + "sci2d, (arc2d_he, arc2d_ne) = make_science_and_arcs(1000, 300)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4ce01aed-67c2-4d86-9e9a-5c230e276271", + "metadata": { + "ExecuteTime": { + "end_time": "2025-10-21T21:13:19.389105Z", + "start_time": "2025-10-21T21:13:18.997903Z" + } + }, + "outputs": [], + "source": [ + "fig, axs = plt.subplots(3, 1, figsize=(6, 4), sharex='all')\n", + "plot_2d_spectrum(sci2d, ax=axs[0], label='Science')\n", + "for i, (arc, label) in enumerate(zip((arc2d_he, arc2d_ne), \"HeI NeI\".split())):\n", + " plot_2d_spectrum(arc, ax=axs[i+1], label=label)\n", + "plt.setp(axs[:-1], xlabel='')\n", + "plt.setp(axs, ylabel='C.D. axis [pix]');" + ] + }, + { + "cell_type": "markdown", + "id": "abbbccbd06a57afd", + "metadata": {}, + "source": [ + "## Spectrum Tracing\n", + "\n", + "The `specreduce.tracing` module provides classes for calculating the trace of a spectrum\n", + "on a 2D image. Traces can be determined semi-automatically or manually and are used as inputs\n", + "for the remaining extraction steps. Supported trace types include:\n", + "\n", + "- `specreduce.tracing.ArrayTrace`\n", + "- `specreduce.tracing.FlatTrace`\n", + "- `specreduce.tracing.FitTrace`\n", + "\n", + "Each trace class takes the 2D spectral image as input, along with any additional information\n", + "needed to define or determine the trace (see the API docs for required parameters).\n", + "\n", + "In this example, the spectrum exhibits significant curvature along the cross-dispersion axis,\n", + "so `FlatTrace` is not suitable. Instead, we use `FitTrace`, which fits an Astropy model to the trace.\n", + "\n", + "We choose to fit a third-degree polynomial and bin the spectrum into 40 bins along the dispersion\n", + "axis. `FitTrace` estimates the PSF centroid along the cross-dispersion axis for each bin, and then\n", + "fits the model to those centroids. Binning along the dispersion axis helps stabilize the fit against\n", + "noise if the spectrum is faint." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cab67d99-23c8-44a9-9c73-7638a804bfff", + "metadata": { + "ExecuteTime": { + "end_time": "2025-10-21T21:13:27.808661Z", + "start_time": "2025-10-21T21:13:27.796122Z" + } + }, + "outputs": [], + "source": [ + "trace = FitTrace(sci2d, bins=40, guess=200, trace_model=models.Polynomial1D(3))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3bde4d81-8df5-4742-93c0-cbd35d399679", + "metadata": { + "ExecuteTime": { + "end_time": "2025-10-21T21:13:28.172408Z", + "start_time": "2025-10-21T21:13:27.984481Z" + } + }, + "outputs": [], + "source": [ + "fig, ax = plot_2d_spectrum(sci2d)\n", + "ax.plot(trace.trace, 'k', ls='--');" + ] + }, + { + "cell_type": "markdown", + "id": "69dc8cb2-1b42-4406-9869-e354c765ec14", + "metadata": {}, + "source": [ + "## Background Subtraction\n", + "\n", + "The `specreduce.background` module contains tools to generate and subtract a \n", + "background image from the input 2D spectral image. The `specreduce.background.Background` \n", + "object is defined by one or more windows, and can be generated with:\n", + "\n", + "* `specreduce.background.Background.one_sided`\n", + "* `specreduce.background.Background.two_sided`\n", + " \n", + "The center of the window can either be passed as a float, integer, or trace\n", + "\n", + "```Python\n", + "bg = Background.one_sided(sci, trace, separation=5, width=2)\n", + "```\n", + "\n", + "or, equivalently\n", + "\n", + "```Python\n", + "bg = Background.one_sided(sci, 15, separation=5, width=2)\n", + "```\n", + "\n", + "The estimated background image can be accessed via `specreduce.background.Background.bkg_image`\n", + "and the background-subtracted image via `specreduce.background.Background.sub_image`.\n", + "\n", + "Let's calculate and remove a two-sided background estimate." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "178de389-ca6c-4468-8f86-4b056ae82d00", + "metadata": { + "ExecuteTime": { + "end_time": "2025-10-21T21:13:35.434329Z", + "start_time": "2025-10-21T21:13:35.401616Z" + } + }, + "outputs": [], + "source": [ + "background = Background.two_sided(sci2d, trace, separation=50, width=30)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "af59b997-0692-4855-8e1a-5d0b6ee9a0dd", + "metadata": { + "ExecuteTime": { + "end_time": "2025-10-21T20:37:53.279056Z", + "start_time": "2025-10-21T20:37:53.085413Z" + } + }, + "outputs": [], + "source": [ + "fig, ax = plot_2d_spectrum(sci2d)\n", + "ax.plot(trace.trace, 'k--')\n", + "for bkt in background.traces:\n", + " ax.fill_between(np.arange(sci2d.shape[1]),\n", + " bkt.trace - background.width / 2,\n", + " bkt.trace + background.width / 2,\n", + " fc='w', alpha=0.2)\n", + " ax.plot(bkt.trace+background.width/2, 'w--', lw=1)\n", + " ax.plot(bkt.trace-background.width/2, 'w--', lw=1)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "13510b0a-6f9e-4bda-aaf2-431ae42a2564", + "metadata": { + "ExecuteTime": { + "end_time": "2025-10-21T21:14:37.598262Z", + "start_time": "2025-10-21T21:14:37.589326Z" + } + }, + "outputs": [], + "source": [ + "sci2d_clean = background.sub_image()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9e9b21c9-d948-43bc-882f-7c7b0a3aa86d", + "metadata": { + "ExecuteTime": { + "end_time": "2025-10-21T21:14:37.986257Z", + "start_time": "2025-10-21T21:14:37.801907Z" + } + }, + "outputs": [], + "source": [ + "fig, ax = plot_2d_spectrum(sci2d_clean.flux.value)" + ] + }, + { + "cell_type": "markdown", + "id": "6b68f282-99d8-41ab-b6cd-5080f855bc69", + "metadata": {}, + "source": [ + "## Spectrum Extraction\n", + "\n", + "The `specreduce.extract` module extracts a 1D spectrum from an input 2D spectrum\n", + "(likely a background-extracted spectrum from the previous step) and a defined\n", + "window, using one of the following implemented methods:\n", + "\n", + "* `specreduce.extract.BoxcarExtract`\n", + "* `specreduce.extract.HorneExtract`\n", + "\n", + "The methods take the input image and trace as inputs, an optional mask treatment\n", + "agrumtn, and a set of method-specific arguments fine-tuning the behavior of each\n", + "method.\n", + "\n", + "### Boxcar Extraction\n", + "\n", + "Boxcar extraction requires a 2D spectrum, a trace, and the extraction aperture width." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b7b1f2f2-0fe7-48e2-9952-205977c76d62", + "metadata": { + "ExecuteTime": { + "end_time": "2025-10-21T21:14:38.982259Z", + "start_time": "2025-10-21T21:14:38.969168Z" + } + }, + "outputs": [], + "source": [ + "aperture_width = 20\n", + "e = BoxcarExtract(sci2d_clean, trace, width=aperture_width)\n", + "sci1d_boxcar = e.spectrum" + ] + }, + { + "cell_type": "markdown", + "id": "6ff894f4-88d2-4431-a09f-17005fcb4053", + "metadata": {}, + "source": [ + "The extracted spectrum can be accessed via the `spectrum` property or by calling\n", + "the `extract` object, which also allows you to override the values the object\n", + "was initialized with." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3f5e6c0c-621e-4852-9141-921b6b663793", + "metadata": { + "ExecuteTime": { + "end_time": "2025-10-21T21:14:42.631810Z", + "start_time": "2025-10-21T21:14:42.374783Z" + } + }, + "outputs": [], + "source": [ + "fig, axs = plt.subplots(2, 1, figsize=(6, 4), sharex='all')\n", + "plot_2d_spectrum(sci2d_clean.flux.value, ax=axs[0])\n", + "axs[0].fill_between(np.arange(sci2d_clean.shape[1]),\n", + " trace.trace - aperture_width / 2,\n", + " trace.trace + aperture_width / 2,\n", + " fc='w', alpha=0.25, ec='k', ls='--')\n", + "\n", + "axs[1].plot(sci1d_boxcar.flux)\n", + "plt.setp(axs[1],\n", + " ylim = (0, 1500),\n", + " xlabel=f'Wavelength [{sci1d_boxcar.spectral_axis.unit}]',\n", + " ylabel=f'Flux [{sci1d_boxcar.flux.unit}]');\n", + "fig.align_ylabels()" + ] + }, + { + "cell_type": "markdown", + "id": "d92ea556-6ae6-4584-825d-78ecaabce8b3", + "metadata": {}, + "source": [ + "### Horne Extraction\n", + "\n", + "\n", + "Horne extraction (a.k.a. optimal extraction) fits the source’s spatial profile across the\n", + "cross-dispersion direction using one of two approaches:\n", + "- A Gaussian profile (default). Optionally, you may include a background model (default is\n", + " a 2D polynomial) to account for residual background in the spatial profile. This\n", + " background-model option is not supported with the interpolated profile.\n", + "- An empirical interpolated profile, enabled via `interpolated_profile`.\n", + "\n", + "Using the `interpolated_profile` option:\n", + "- The image is binned along the wavelength axis (number of bins set by `n_bins_interpolated_profile`),\n", + " averaged within each bin, and then interpolated (linear by default; the interpolation degree in\n", + " x and y can be set via `interp_degree_interpolated_profile`) to form an empirical spatial profile.\n", + "- You can select this mode by passing a string to use default settings, or a dictionary to override the defaults.\n", + "\n", + "Examples:\n", + "Use the interpolated profile with default bins and interpolation degree:\n", + " \n", + " interp_profile_extraction = extract(spatial_profile='interpolated_profile')\n", + "\n", + "Use 20 bins instead of the default of 10:\n", + "\n", + " interp_profile_extraction = extract(spatial_profile={\n", + " 'name': 'interpolated_profile',\n", + " 'n_bins_interpolated_profile': 20\n", + " })\n", + "\n", + "Use cubic interpolation instead of the default linear:\n", + "\n", + " interp_profile_extraction = extract(spatial_profile={\n", + " 'name': 'interpolated_profile',\n", + " 'interp_degree_interpolated_profile': 3\n", + " })\n", + "\n", + "The Horne extraction algorithm requires a variance array. If the input image is\n", + "an `astropy.nddata.NDData` object with `image.uncertainty` provided, that uncertainty\n", + "will be used. Otherwise, you must supply the `variance` parameter.\n", + "\n", + " extract = HorneExtract(image - bg, trace, variance=var_array)\n", + "\n", + "As usual, parameters can be set when creating the `HorneExtract` object or overridden when calling the object’s extraction method.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c356be3a-7358-4f3f-888d-b75832e48548", + "metadata": { + "ExecuteTime": { + "end_time": "2025-10-21T21:14:58.294487Z", + "start_time": "2025-10-21T21:14:58.135177Z" + } + }, + "outputs": [], + "source": [ + "e = HorneExtract(sci2d_clean, trace)\n", + "sci1d_horne = e.spectrum" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c05bcd15-ac47-4e81-a2e5-d725e86099a6", + "metadata": { + "ExecuteTime": { + "end_time": "2025-10-21T21:14:58.729932Z", + "start_time": "2025-10-21T21:14:58.639898Z" + } + }, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(1, 1, figsize=(6, 2), sharex='all')\n", + "ax.plot(sci1d_horne.flux)\n", + "plt.setp(ax,\n", + " ylim = (0, 1500),\n", + " xlabel=f'Wavelength [{sci1d_horne.spectral_axis.unit}]',\n", + " ylabel=f'Flux [{sci1d_horne.flux.unit}]')\n", + "ax.autoscale(axis='x', tight=True)" + ] + }, + { + "cell_type": "markdown", + "id": "d9fa857a-b6bb-4057-8476-1ffc328a0b24", + "metadata": {}, + "source": [ + "## Wavelength Calibration\n", + "\n", + "The `specreduce.wavecal1d.WavelengthCalibration1D` class provides tools for\n", + "one-dimensional wavelength calibration given a number of arc calibration spectra and\n", + "corresponding catalog line lists.\n", + "\n", + "First, we need to extract the arc spectra from the 2D arc frames using the trace\n", + "we calculated from our science spectrum. Next, we instantiate a `WavelengthCalibration1D`\n", + "object and pass the arc spectra and the line lists to use. We also pass it the wavelength\n", + "unit we want to use, the line list bounds (we know that our instrument setup gives\n", + "a spectrum that covers approximately the range of 500-1000 nm), and the number of strongest\n", + "lines to include from the catalog line lists." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2f5a6050-d396-45bd-92e7-120a8a0d570d", + "metadata": { + "ExecuteTime": { + "end_time": "2025-10-21T21:15:32.791564Z", + "start_time": "2025-10-21T21:15:32.768331Z" + } + }, + "outputs": [], + "source": [ + "arc1d_he = BoxcarExtract(arc2d_he, trace, width=aperture_width).spectrum\n", + "arc1d_ne = BoxcarExtract(arc2d_ne, trace, width=aperture_width).spectrum" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "05e0d6e5-b493-4ca7-a25d-e22d7ecfc9d2", + "metadata": { + "ExecuteTime": { + "end_time": "2025-10-21T21:15:33.337471Z", + "start_time": "2025-10-21T21:15:33.329630Z" + } + }, + "outputs": [], + "source": [ + "wc = WavelengthCalibration1D([arc1d_he, arc1d_ne],\n", + " line_lists=[\"HeI\", \"NeI\"],\n", + " line_list_bounds=(500, 1000),\n", + " n_strogest_lines=15,\n", + " unit=u.nm)" + ] + }, + { + "cell_type": "markdown", + "id": "734f8364a6c86f1e", + "metadata": {}, + "source": [ + "### Line Finding\n", + "\n", + "Next, we wse `WavelengthCalibration1D.find_lines` to detect emission lines in the\n", + "arc spectra. Then, we visualize the catalog lines alongside the arc spectra and the\n", + "detected lines with `WavelengthCalibration1D.plot_fit`." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "844e3994-367e-4936-8f4f-245215a4ece7", + "metadata": { + "ExecuteTime": { + "end_time": "2025-10-21T21:15:35.352789Z", + "start_time": "2025-10-21T21:15:35.300727Z" + } + }, + "outputs": [], + "source": [ + "with warnings.catch_warnings():\n", + " warnings.simplefilter('ignore')\n", + " wc.find_lines(3, noise_factor=25)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "dceca5d1-1872-42ef-b6b4-5eacca755c41", + "metadata": { + "ExecuteTime": { + "end_time": "2025-10-21T21:15:36.018727Z", + "start_time": "2025-10-21T21:15:35.515189Z" + } + }, + "outputs": [], + "source": [ + "fig = wc.plot_fit(figsize=(6, 6))\n", + "plt.setp(fig.axes[::2], xlim=(490, 1000));" + ] + }, + { + "cell_type": "markdown", + "id": "8c97a1c9-0d2c-4947-96b2-1288b128ee38", + "metadata": {}, + "source": [ + "### Wavelength Solution Fitting\n", + "\n", + "With the observed and catalog lines identified, we can now compute the wavelength\n", + "solution, represented by a `specreduce.wavesol1d.WavelengthSolution1D` object. The\n", + "wavelength calibration class provides several methods; here we use the interactive\n", + "`fit_lines` method. This method:\n", + "- takes a list of matched pixel positions and their corresponding wavelengths,\n", + "- fits a polynomial to these pairs, and\n", + "- refines the solution by incorporating all observed and catalog lines within a\n", + " specified pixel distance of the initial fit." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4b77b125-5f15-4e49-b6e4-29c6607a28d7", + "metadata": { + "ExecuteTime": { + "end_time": "2025-10-21T21:15:38.870556Z", + "start_time": "2025-10-21T21:15:38.863758Z" + } + }, + "outputs": [], + "source": [ + "ws = wc.fit_lines([72, 295, 403, 772], [588, 668, 707, 838], \n", + " degree=3, refine_max_distance=5)" + ] + }, + { + "cell_type": "markdown", + "id": "e10c3e233717766d", + "metadata": {}, + "source": [ + "Now we can visualize the observed and catalog lines again. Set `obs_to_wav=True`\n", + "to display the observed arc spectra and lines in wavelength space. Matched fitted\n", + "observed and catalog lines are shown with solid lines; unmatched lines are shown\n", + "with dashed lines.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a62b4b71-789e-4df5-8ec2-61ea4d350ceb", + "metadata": { + "ExecuteTime": { + "end_time": "2025-10-21T21:15:40.398662Z", + "start_time": "2025-10-21T21:15:39.874303Z" + } + }, + "outputs": [], + "source": [ + "fig = wc.plot_fit(figsize=(6, 6), obs_to_wav=True)\n", + "plt.setp(fig.axes, xlim=ws.p2w([0, 1000]));" + ] + }, + { + "cell_type": "markdown", + "id": "ad7df493fdef8432", + "metadata": {}, + "source": [ + "We can also visualize the residuals between the observed and catalog lines\n", + "(either in wavelength or pixel space) to assess the quality of the fit.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "77063157-9cd8-4ec8-8082-80b836efaa24", + "metadata": { + "ExecuteTime": { + "end_time": "2025-10-21T21:15:43.910595Z", + "start_time": "2025-10-21T21:15:43.693771Z" + } + }, + "outputs": [], + "source": [ + "fig, axs = plt.subplots(2, 1, figsize=(6, 4), constrained_layout=True)\n", + "wc.plot_residuals(ax=axs[0], space=\"wavelength\")\n", + "wc.plot_residuals(ax=axs[1], space=\"pixel\")\n", + "fig.align_ylabels()" + ] + }, + { + "cell_type": "markdown", + "id": "ddac0171bbb9b7a7", + "metadata": {}, + "source": [ + "## Spectrum Resampling\n", + "\n", + "With a satisfactory wavelength solution in hand, we can finalize the reduction. You have two options:\n", + "- Attach the `gwcs` object from the `WavelengthSolution1D` to the science spectrum, or\n", + "- Resample the science spectrum onto a defined wavelength grid using `WavelengthSolution1D.resample`.\n", + "\n", + "Here, we choose to resample the spectrum.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f800abf7-442d-4919-8652-1f4eba3d22dd", + "metadata": { + "ExecuteTime": { + "end_time": "2025-10-21T21:10:08.038359Z", + "start_time": "2025-10-21T21:10:08.026830Z" + } + }, + "outputs": [], + "source": [ + "sci1d_resampled = ws.resample(sci1d_boxcar)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "549f06a81dfeb7b5", + "metadata": { + "ExecuteTime": { + "end_time": "2025-10-21T21:10:29.633337Z", + "start_time": "2025-10-21T21:10:29.514489Z" + } + }, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(figsize=(6, 3))\n", + "ax.plot(sci1d_resampled.spectral_axis, sci1d_resampled.flux)\n", + "plt.setp(ax,\n", + " ylim = (0, 1500),\n", + " xlabel=f\"Wavelength [{sci1d_resampled.spectral_axis.unit}]\",\n", + " ylabel=f\"Flux density [{sci1d_resampled.flux.unit}]\")\n", + "ax.autoscale(axis='x', tight=True)" + ] + }, + { + "cell_type": "markdown", + "id": "d21011efe913e431", + "metadata": {}, + "source": [ + "By default, the `resample` method creates a linear wavelength grid spanning the\n", + "spectrum’s range, with as many bins as there are input pixels. However, you can\n", + "also fully customize the grid by supplying explicit `bin_edges`." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6ecfdfd0eb9b87ed", + "metadata": { + "ExecuteTime": { + "end_time": "2025-10-21T21:16:29.849432Z", + "start_time": "2025-10-21T21:16:29.839531Z" + } + }, + "outputs": [], + "source": [ + "bin_edges = np.concatenate([np.geomspace(ws.pix_to_wav(0), 650, num=100), \n", + " np.linspace(660, 800, 50),\n", + " np.linspace(810, ws.pix_to_wav(999), 400)])\n", + "\n", + "sci1d_resampled = ws.resample(sci1d_boxcar, bin_edges=bin_edges)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "68b062f209cde4a9", + "metadata": { + "ExecuteTime": { + "end_time": "2025-10-21T21:16:31.032973Z", + "start_time": "2025-10-21T21:16:30.920578Z" + } + }, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(figsize=(6, 3))\n", + "ax.plot(sci1d_resampled.spectral_axis, sci1d_resampled.flux, '-')\n", + "plt.setp(ax,\n", + " ylim = (0, 1500),\n", + " xlabel=f\"Wavelength [{sci1d_resampled.spectral_axis.unit}]\",\n", + " ylabel=f\"Flux density [{sci1d_resampled.flux.unit}]\")\n", + "ax.autoscale(axis='x', tight=True)" + ] + }, + { + "cell_type": "markdown", + "id": "9e52e79e95c73038", + "metadata": {}, + "source": [ + "Finally, we can save the spectrum as fits." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8ab90caa-0453-4c36-b66f-5b28dea205db", + "metadata": { + "ExecuteTime": { + "end_time": "2025-10-21T21:00:44.342319Z", + "start_time": "2025-10-21T21:00:44.319720Z" + } + }, + "outputs": [], + "source": [ + "sci1d_resampled.write('science_spectrum.fits', overwrite=True)" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.10" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/docs/getting_started/quickstart.py b/docs/getting_started/quickstart.py new file mode 100644 index 00000000..67637575 --- /dev/null +++ b/docs/getting_started/quickstart.py @@ -0,0 +1,199 @@ +from typing import Sequence + +import numpy as np +import matplotlib.pyplot as plt +from astropy.modeling import models, Model + +from astropy.wcs import WCS +from astropy.nddata import StdDevUncertainty, CCDData, VarianceUncertainty +import astropy.units as u +from photutils.datasets import apply_poisson_noise + +from specreduce.utils.synth_data import make_2d_arc_image, make_2d_trace_image + + +def make_2d_spec_image( + nx: int = 1000, + ny: int = 300, + wcs: WCS | None = None, + extent: Sequence[int | float] = (6500, 9500), + wave_unit: u.Unit = u.Angstrom, + wave_air: bool = False, + background: int | float = 5, + line_fwhm: float = 5.0, + linelists: list[str] = ("OH_GMOS"), + airglow_amplitude: float = 1.0, + spectrum_amplitude: float = 1.0, + tilt_func: Model = models.Legendre1D(degree=0), + trace_center: int | float | None = None, + trace_order: int = 3, + trace_coeffs: None | dict[str, int | float] = None, + source_profile: Model = models.Moffat1D(amplitude=10, alpha=0.1), + add_noise: bool = True, +) -> CCDData: + """ + Generate a simulated 2D spectroscopic image. + + This function creates a two-dimensional synthetic spectroscopic image, combining + arc (wavelength calibration) data and trace (spatial profile) data. The image simulates + a realistic spectral profile with contributions from background, airglow, and a modeled + source profile. Noise can optionally be added to the resulting image. It employs models + for traces and tilts, and allows customization of wavelength range, amplitude scales, + and noise addition. + + Parameters + ---------- + nx : int, optional + Number of columns (spatial dimension) in the generated image, by default 3000. + ny : int, optional + Number of rows (dispersion dimension) in the generated image, by default 1000. + wcs : WCS or None, optional + World Coordinate System (WCS) object for the image. Specifies spectral coordinates. + If None, WCS is not assigned, by default None. + extent : Sequence[int or float], optional + Wavelength range for the generated image in units determined by `wave_unit`. + Defined as a tuple (min, max), by default (6500, 9500). + wave_unit : Unit, optional + Unit of the generated wavelength axis, by default astropy.units.Angstrom. + wave_air : bool, optional + If True, use air wavelengths. If False, use vacuum wavelengths, by default False. + background : int or float, optional + Constant background level added to the image, by default 5. + line_fwhm : float, optional + Full-width half maximum (in pixels) for spectral lines in the image, by default 5.0. + linelists : list of str, optional + Names of line lists (e.g., emission lines) to simulate in the arc image, + by default ["OH_GMOS"]. + airglow_amplitude : float, optional + Scaling factor for the airglow contribution to the image, by default 1.0. + spectrum_amplitude : float, optional + Scaling factor for the primary spectrum contribution (trace image), by default 1.0. + amplitude_scale : float, optional + Global amplitude scaling factor applied to both arcs and traces, by default 1.0. + tilt_func : Model, optional + Astropy model representing the spectral tilt across the spatial axis. By default, + a zero-degree Legendre polynomial (no tilt) is used. + trace_center : int, float, or None, optional + Central position of the trace in spatial pixels. If None, defaults to trace modeled + by coefficients in `trace_coeffs`, by default None. + trace_order : int, optional + Polynomial order to model the trace profile spatial variation, by default 3. + trace_coeffs : None or dict of str to int or float, optional + Coefficients for modeling the trace position (spatial axis). If None, defaults to + {"c0": 0, "c1": 50, "c2": 100}, by default None. + source_profile : Model, optional + Astropy model to simulate the source profile along the trace, by default Moffat1D with + amplitude=10 and alpha=0.1. + add_noise : bool, optional + If True, adds Poisson noise to the generated spectral image, by default True. + + Returns + ------- + CCDData + A CCDData object containing the simulated 2D spectroscopic image. The data includes + contributions from arc lines, traces, airglow, and background, with optional noise + added. The WCS information, if provided, is preserved. + """ + + if trace_coeffs is None: + trace_coeffs = {"c0": 0, "c1": 50, "c2": 100} + + arc_image = make_2d_arc_image( + nx=nx, + ny=ny, + wcs=wcs, + extent=extent, + wave_unit=wave_unit, + wave_air=wave_air, + background=0, + line_fwhm=line_fwhm, + linelists=linelists, + tilt_func=tilt_func, + add_noise=False, + ) + + trace_image = make_2d_trace_image( + nx=nx, + ny=ny, + background=0, + trace_center=trace_center, + trace_order=trace_order, + trace_coeffs=trace_coeffs, + profile=source_profile, + add_noise=False, + ) + + wl = wcs.spectral.pixel_to_world(np.arange(nx)).to(u.nm).value + signal = 0.8 + 0.2 * np.abs(np.sin((wl - 650) / 100 * 2 * np.pi)) ** 5 + + n = lambda a: a / a.max() + spec_image = ( + airglow_amplitude * n(arc_image.data) + + spectrum_amplitude * n(trace_image.data) * signal + + background + ) + + if add_noise: + from photutils.datasets import apply_poisson_noise + + spec_image = apply_poisson_noise(spec_image) + + return CCDData( + spec_image, + unit=u.count, + wcs=arc_image.wcs, + uncertainty=StdDevUncertainty(np.sqrt(spec_image)), + ) + + +def make_science_and_arcs(ndisp: int = 1000, ncross: int = 300): + refx = ndisp // 2 + wcs = WCS( + header={ + "CTYPE1": "AWAV-GRA", # Grating dispersion function with air wavelengths + "CUNIT1": "Angstrom", # Dispersion units + "CRPIX1": refx, # Reference pixel [pix] + "CRVAL1": 7410, # Reference value [Angstrom] + "CDELT1": 3 * 1.19, # Linear dispersion [Angstrom/pix] + "PV1_0": 5.0e5, # Grating density [1/m] + "PV1_1": 1, # Diffraction order + "PV1_2": 8.05, # Incident angle [deg] + "CTYPE2": "PIXEL", # Spatial detector coordinates + "CUNIT2": "pix", # Spatial units + "CRPIX2": 1, # Reference pixel + "CRVAL2": 0, # Reference value + "CDELT2": 1, # Spatial units per pixel + } + ) + + science = make_2d_spec_image( + ndisp, + ncross, + add_noise=True, + background=10, + wcs=wcs, + airglow_amplitude=40.0, + spectrum_amplitude=160.0, + trace_coeffs={"c0": 0, "c1": 30, "c2": 40}, + source_profile=models.Moffat1D(amplitude=1, alpha=0.3), + ) + + arcargs = dict(wcs=wcs, line_fwhm=3, background=0, add_noise=False) + arcs = [] + for linelist in ["HeI", "NeI"]: + arc = make_2d_arc_image(ndisp, ncross, linelists=[linelist], **arcargs) + arc.data = apply_poisson_noise(200*(arc.data / arc.data.max()) + 10) - 10 + arcs.append(arc) + return science, arcs + + +def plot_2d_spectrum(spec, ax=None, label=None): + if ax is None: + fig, ax = plt.subplots(figsize=(6, 2), constrained_layout=True) + else: + fig = ax.figure + ax.imshow(spec, origin="lower", aspect="auto") + if label is not None: + ax.text(0.98, 0.9, label, va="top", ha="right", transform=ax.transAxes, c="w") + plt.setp(ax, xlabel="Dispersion axis [pix]", ylabel="Cross-disp. axis [pix]") + return fig, ax diff --git a/docs/index.rst b/docs/index.rst index 40727b1a..f4c0b82d 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -1,82 +1,125 @@ -######################## -Specreduce Documentation -######################## - -The `specreduce `_ package -aims to provide a data reduction toolkit for optical -and infrared spectroscopy, on which applications such as pipeline processes for -specific instruments can be built. The scope of its functionality is limited to -basic spectroscopic reduction, currently encompassing the following three tasks: - -#. Determining the trace of a spectrum dispersed in a 2D image, either by setting a flat - trace, providing a custom trace array, or fitting a spline, polynomial, or other model - to the positions of the dispersed spectrum. -#. Generating a background based on a region on one or both sides of this trace, and making - available the background image, 1D spectrum of the background, and the - background-subtracted image. -#. Performing either a Horne (a.k.a. "optimal") or boxcar extraction on either the original - or background-subtracted 2D spectrum, using the trace generated by the first task to - generate a 1D spectrum. - -Further details about these capabilities are detailed in the sections linked below. -Beyond these tasks, basic *image* processing steps (such as bias subtraction) are covered by -`ccdproc `_ -and other packages, data analysis by `specutils `_, -and visualization by `matplotlib `_. A few -examples will be provided that feature ``specreduce`` in conjunction with these -complementary packages. +.. _docroot: -.. note:: - - Specreduce is currently an incomplete work-in-progress and is liable to - change. Please feel free to contribute code and suggestions through github. +########## +Specreduce +########## +| **Version**: |release| +| **Date**: |today| -.. _spectral-extraction: +:ref:`Specreduce ` is an `Astropy `_ +`coordinated package `_ that +provides a toolkit for reducing optical and infrared spectroscopic data. It +offers the building blocks for basic spectroscopic reduction, and is designed to serve as a +foundation upon which instrument-specific pipelines and analysis tools can be built. -******************* -Spectral Extraction -******************* +Specreduce includes tools for determining and modeling spectral traces, performing +background subtraction, extracting one-dimensional spectra using both optimal and boxcar methods, +and applying wavelength correction derived from calibration data. +Beyond these tasks, basic image processing steps, data analysis, and visualisation are covered by +other Astropy ecosystem packages like `ccdproc `_, +`specutils `_, and +`matplotlib `_. The documentation includes examples demonstrating how these +tools can be combined to create complete spectroscopic workflows. -.. toctree:: - :maxdepth: 2 +.. image:: + roadmap.png - extraction_quickstart.rst +.. note:: -*********** -Calibration -*********** + Specreduce is an active, community-driven project, and we warmly welcome contributions of all kinds, + from reporting bugs and suggesting new features to improving documentation or writing code. + Whether you are an experienced developer or just getting started, your input helps make Specreduce + better for everyone. If you would like to get involved, see the :ref:`contributor's guide + ` for details on how to participate and become part of the wider Astropy + community. .. toctree:: :maxdepth: 1 + :hidden: - wavelength_calibration.rst - extinction.rst - specphot_standards.rst - mask_treatment/mask_treatment.rst + getting_started/index + user_guide.rst + contributing.rst + api.rst -***** -Index -***** +.. grid:: 2 + :gutter: 2 3 4 4 -.. toctree:: - :maxdepth: 1 + .. grid-item-card:: + :text-align: center - api.rst + **Getting Started** + ^^^^^^^^^^^^^^^^^^^ -* :ref:`genindex` -* :ref:`modindex` -* :ref:`search` + New to Specreduce? Check out the getting started guides. + +++ -**************** -Development Docs -**************** + .. button-ref:: getting_started/index + :expand: + :color: primary + :click-parent: -.. toctree:: - :maxdepth: 1 + To the getting started guides + + .. grid-item-card:: + :text-align: center + + **User Guide** + ^^^^^^^^^^^^^^ + + The user guide provides in-depth information on the key concepts + of Specreduce with useful background information and explanation. + + +++ + + .. button-ref:: user_guide + :expand: + :color: primary + :click-parent: + + To the user guide + + .. grid-item-card:: + :text-align: center + + **API Reference** + ^^^^^^^^^^^^^^^^^ + + The API reference contains a detailed description of the + functions, modules, and objects included in Specreduce. It + assumes that you have an understanding of the key concepts. + + +++ + + .. button-ref:: api + :expand: + :color: primary + :click-parent: + + To the API reference + + .. grid-item-card:: + :text-align: center + + **Contributor's Guide** + ^^^^^^^^^^^^^^^^^^^^^^^ + + Want to contribute to specreduce? Found a bug? The contributing guidelines will + show you how to improve specreduce. + + +++ + + .. button-ref:: contributing + :expand: + :color: primary + :click-parent: + + To the contributor's guide - process/index - terms +.. * :ref:`genindex` +.. * :ref:`modindex` +.. * :ref:`search` diff --git a/docs/roadmap.png b/docs/roadmap.png new file mode 100644 index 00000000..4f75dc40 Binary files /dev/null and b/docs/roadmap.png differ diff --git a/docs/terms.rst b/docs/terms.rst index 815233e0..304bdfac 100644 --- a/docs/terms.rst +++ b/docs/terms.rst @@ -476,10 +476,9 @@ Processing Steps what redshift?) ----- Mentioned but not defined -------------------------- +========================= - WCS & Database archive - Cloud archiving diff --git a/docs/trace.rst b/docs/trace.rst new file mode 100644 index 00000000..7bb86063 --- /dev/null +++ b/docs/trace.rst @@ -0,0 +1,166 @@ +Tracing +======= + +The `specreduce.tracing` module provides three ``Trace`` classes that are used to define the +spatial position (trace) of a spectrum across a 2D detector image: `~specreduce.tracing.FlatTrace`, +`~specreduce.tracing.FitTrace`, and `~specreduce.tracing.ArrayTrace`. Each trace class requires +the 2D spectral image as input, along with trace-specific parameters that control how the trace +is determined. + +Flat trace +---------- + +`~specreduce.tracing.FlatTrace` assumes that the spectrum follows a straight line across the +detector, and is best for well-aligned spectrographs with minimal optical distortion. To +initialize a `~specreduce.tracing.FlatTrace`, we need to specify the fixed cross-dispersion +pixel position: + +.. code-block:: python + + trace = specreduce.tracing.FlatTrace(image, trace_pos=6) + + +.. plot:: + + import numpy as np + import matplotlib.pyplot as plt + from scipy.stats import norm + + fw = 10 + nd, ncd = 31, 13 + xd, xcd = np.arange(nd), np.arange(ncd) + + spectrum = np.zeros((ncd, nd)) + spectrum[:,:] = norm(6.0, 1.5).pdf(xcd)[:, None] + + plt.rc('font', size=13) + fig, ax = plt.subplots(figsize=(fw, fw*(ncd/nd)), constrained_layout=True) + ax.imshow(spectrum, origin='lower') + ax.plot((0, nd-1), (6, 6), c='k') + ax.set_xticks(xd+0.5, minor=True) + ax.set_yticks(xcd+0.5, minor=True) + ax.grid(alpha=0.25, lw=1, which='minor') + plt.setp(ax, xlabel='Dispersion axis', ylabel='Cross-dispersion axis') + fig.show() + +FitTrace +-------- + +`~specreduce.tracing.FitTrace` fits a polynomial function to automatically detected spectrum +positions, and is suitable for typical spectra with smooth, continuous trace profiles. The trace +model can be chosen from `~astropy.modeling.polynomial.Chebyshev1D`, +`~astropy.modeling.polynomial.Legendre1D`, `~astropy.modeling.polynomial.Polynomial1D`, +or `~astropy.modeling.spline.Spline1D`, and the fitting can be optimized by binning the spectrum +along the dispersion axis. + +.. code-block:: python + + trace = specreduce.tracing.FitTrace(image, bins=10, trace_model=Polynomial1D(3)) + +.. plot:: + + import numpy as np + import matplotlib.pyplot as plt + from scipy.stats import norm + plt.rc('font', size=13) + + fw = 10 + nd, ncd = 31, 13 + xd, xcd = np.arange(nd), np.arange(ncd) + + tr = np.poly1d([-0.01, 0.2, 7.0]) + spectrum = np.zeros((ncd, nd)) + for i,x in enumerate(xd): + spectrum[:,i] = norm(tr(x), 1.0).pdf(xcd) + + fig, ax = plt.subplots(figsize=(fw, fw*(ncd/nd)), constrained_layout=True) + ax.imshow(spectrum, origin='lower') + ax.plot(xd, tr(xd), 'k') + ax.set_xticks(xd+0.5, minor=True) + ax.set_yticks(xcd+0.5, minor=True) + ax.grid(alpha=0.25, lw=1, which='minor') + plt.setp(ax, xlabel='Dispersion axis', ylabel='Cross-dispersion axis') + fig.show() + +The method works by (optionally) binning the 2D spectrum along the dispersion axis, finding +the PSF peak position along the cross-dispersion for each bin, and then fitting a 1D polynomial to +the cross-dispersion and dispersion axis positions. Binning is optional, and the native image +resolution is used by default. Binning can significantly increase the reliability of the fitting +with low SNR spectra, and always increases the speed. + +The peak detection method can be chosen from ``max``, ``centroid``, and ``gaussian``. Of these +methods, ``max`` is the fastest but yields an integer pixel precision. Both ``centroid`` and +``gaussian`` can be used when sub-pixel precision is required, and ``gaussian``, while being the +slowest method of the three, is the best option if the data is significantly contaminated by +non-finite values. + +.. plot:: + + import numpy as np + import matplotlib.pyplot as plt + from scipy.stats import norm + from numpy.random import seed, normal + from astropy.modeling.models import Gaussian1D + from astropy.modeling.fitting import DogBoxLSQFitter + plt.rc('font', size=13) + seed(5) + + fw = 10 + nd, ncd = 31, 13 + xd, xcd = np.arange(nd), np.arange(ncd) + + psf = norm(5.4, 1.5).pdf(xcd) + normal(0, 0.01, ncd) + fitter = DogBoxLSQFitter() + m = fitter(Gaussian1D(), xcd, psf) + + fig, ax = plt.subplots(figsize=(fw, fw*(ncd/nd)), constrained_layout=True) + ax.step(xcd, psf, where='mid', c='k') + ax.axvline(xcd[np.argmax(psf)], label='max') + ax.axvline(np.average(xcd, weights=psf), ls='--', label='centroid') + ax.axvline(m.mean.value, ls=':', label='gaussian') + ax.plot(xcd, m(xcd), ls=':') + ax.legend() + plt.setp(ax, yticks=[], ylabel='Flux', xlabel='Cross-dispersion axis [pix]', xlim=(0, ncd-1)) + fig.show() + +ArrayTrace +---------- + +`~specreduce.tracing.ArrayTrace` uses a pre-defined array of positions for maximum flexibility, +and is ideal for complex or unusual trace shapes that are difficult to model mathematically. +`~specreduce.tracing.ArrayTrace` initialization requires an array of cross-dispersion pixel +positions. The size of the array must match the number of dispersion-axis pixels in the image. + +.. code-block:: python + + trace = specreduce.tracing.ArrayTrace(image, positions) + +.. plot:: + + import numpy as np + import matplotlib.pyplot as plt + from scipy.stats import norm + plt.rc('font', size=13) + + fw = 10 + nd, ncd = 31, 13 + xd, xcd = np.arange(nd), np.arange(ncd) + + tr = np.full_like(xd, 6) + tr[:6] = 4 + tr[15:23] = 8 + + spectrum = np.zeros((ncd, nd)) + + for i,x in enumerate(xd): + spectrum[:,i] = norm(tr[i], 1.0).pdf(xcd) + + plt.rc('font', size=13) + fig, ax = plt.subplots(figsize=(fw, fw*(ncd/nd)), constrained_layout=True) + ax.imshow(spectrum, origin='lower') + ax.plot(xd, tr, 'k') + ax.set_xticks(xd+0.5, minor=True) + ax.set_yticks(xcd+0.5, minor=True) + ax.grid(alpha=0.25, lw=1, which='minor') + plt.setp(ax, xlabel='Dispersion axis', ylabel='Cross-dispersion axis') + fig.show() diff --git a/docs/user_guide.rst b/docs/user_guide.rst new file mode 100644 index 00000000..2fab4e82 --- /dev/null +++ b/docs/user_guide.rst @@ -0,0 +1,20 @@ +User Guide +========== + +.. toctree:: + :maxdepth: 1 + :caption: Core Functionality + + trace + background + extraction + wavelength_calibration/wavelength_calibration + + +.. toctree:: + :maxdepth: 1 + :caption: Additional Topics + + extinction + specphot_standards + mask_treatment/mask_treatment \ No newline at end of file diff --git a/docs/wavelength_calibration.rst b/docs/wavelength_calibration.rst deleted file mode 100644 index f34d7dec..00000000 --- a/docs/wavelength_calibration.rst +++ /dev/null @@ -1,53 +0,0 @@ -.. _wavelength_calibration: - -Wavelength Calibration -====================== - -Wavelength calibration is currently supported for 1D spectra. Given a list of spectral -lines with known wavelengths and estimated pixel positions on an input calibration -spectrum, you can currently use ``specreduce`` to: - -#. Fit an ``astropy`` model to the wavelength/pixel pairs to generate a spectral WCS - solution for the dispersion. -#. Apply the generated spectral WCS to other `~specutils.Spectrum1D` objects. - -1D Wavelength Calibration -------------------------- - -The `~specreduce.wavelength_calibration.WavelengthCalibration1D` class can be used -to fit a dispersion model to a list of line positions and wavelengths. Future development -will implement catalogs of known lamp spectra for use in matching observed lines. In the -example below, the line positions (``pixel_centers``) have already been extracted from -``lamp_spectrum``:: - - import astropy.units as u - from specreduce import WavelengthCalibration1D - pixel_centers = [10, 22, 31, 43] - wavelengths = [5340, 5410, 5476, 5543]*u.AA - test_cal = WavelengthCalibration1D(lamp_spectrum, line_pixels=pixel_centers, - line_wavelengths=wavelengths) - calibrated_spectrum = test_cal.apply_to_spectrum(science_spectrum) - -The example above uses the default model (`~astropy.modeling.functional_models.Linear1D`) -to fit the input spectral lines, and then applies the calculated WCS solution to a second -spectrum (``science_spectrum``). Any other 1D ``astropy`` model can be provided as the -input ``model`` parameter to the `~specreduce.wavelength_calibration.WavelengthCalibration1D`. -In the above example, the model fit and WCS construction is all done as part of the -``apply_to_spectrum()`` call, but you could also access the `~gwcs.wcs.WCS` object itself -by calling:: - - test_cal.wcs - -The calculated WCS is a cached property that will be cleared if the ``line_list``, ``model``, -or ``input_spectrum`` properties are updated, since these will alter the calculated dispersion -fit. - -You can also provide the input pixel locations and wavelengths of the lines as an -`~astropy.table.QTable` with (at minimum) columns ``pixel_center`` and ``wavelength``, -using the ``matched_line_list`` input argument:: - - from astropy.table import QTable - pixels = [10, 20, 30, 40]*u.pix - wavelength = [5340, 5410, 5476, 5543]*u.AA - line_list = QTable([pixels, wavelength], names=["pixel_center", "wavelength"]) - test_cal = WavelengthCalibration1D(lamp_spectrum, matched_line_list=line_list) \ No newline at end of file diff --git a/docs/wavelength_calibration/wavelength_calibration.rst b/docs/wavelength_calibration/wavelength_calibration.rst new file mode 100644 index 00000000..4e76d1db --- /dev/null +++ b/docs/wavelength_calibration/wavelength_calibration.rst @@ -0,0 +1,251 @@ +.. _wavelength_calibration: + +Wavelength Calibration +====================== + +In spectroscopy, the raw data from a detector typically records the intensity of light +at different pixel positions. For scientific analysis, we need to know the physical +wavelength (or frequency, or energy) corresponding to each pixel. **Wavelength +calibration** is the process of determining this mapping, creating a model that +converts pixel coordinates to wavelength values. + +This is often achieved by observing a calibration source with well-known emission +lines (e.g., an arc lamp). By identifying the pixel positions of these lines, we can fit +a mathematical model describing the spectrograph’s dispersion. + +The ``specreduce`` library provides two complementary classes for 1D calibration: + +- :class:`~specreduce.wavecal1d.WavelengthCalibration1D`: a high-level calibration workflow + to detect and match arc lines and fit a dispersion model. +- `~specreduce.wavesol1d.WavelengthSolution1D`: a lightweight wavelength-solution + container that holds the pixel→wavelength transform with its inverse, a WCS solution, and + a method for flux-conserving resampling. + +Tools for calibrating 2D spectra will be added in later versions. + +1D Wavelength Calibration +------------------------- + +The :class:`~specreduce.wavecal1d.WavelengthCalibration1D` class encapsulates the data +and methods needed to compute a 1D wavelength solution, and the +:class:`~specreduce.wavesol1d.WavelengthSolution1D` class provides the tools to use the +computed solution. A typical workflow involves: + +1. **Initialization**: Create an instance of the class with either arc spectra or + pre-identified observed line positions, and provide one or more line lists of known + wavelengths. Multiple arcs and multiple catalogs are supported. +2. **Line Detection (if using arc spectra)**: Identify line positions in the arc spectra. +3. **Fitting the Wavelength Solution**: Use either direct fitting with known pixel– + wavelength pairs (:meth:`~specreduce.wavecal1d.WavelengthCalibration1D.fit_lines`) + or automated global optimization + (:meth:`~specreduce.wavecal1d.WavelengthCalibration1D.fit_dispersion`). +4. **Refinement and Inspection**: Optionally refine the solution and check fit quality + using RMS statistics and plots. +5. **Applying the Solution**: Access the fitted solution as a + :class:`~specreduce.wavesol1d.WavelengthSolution1D` to resample spectra, or attach the + solution to a `specutils.Spectrum` as a `~gwcs.wcs.WCS` object. + +Quickstart +---------- + +1. Initialization +***************** + +You instantiate the :class:`~specreduce.wavecal1d.WavelengthCalibration1D` by providing basic +information about your setup and data. You must provide *either* a list of arc spectra (or a +single arc spectrum) *or* a list of known observed line positions: + +* **Using an Arc Spectrum**: Provide the arc spectrum as a `specutils.Spectrum` + object via the ``arc_spectra`` argument. You also need to provide a ``line_lists`` argument, + which can be a list of known catalog wavelengths or the name(s) of standard line lists + recognized by `specreduce` (e.g., ``"HeI"``). You can query the available line lists by using + the `~specreduce.calibration_data.get_available_line_catalogs` function. + + .. code-block:: python + + wc = WavelengthCalibration1D(arc_spectra=arc_he, line_lists=["HeI"]) + +* **Using multiple Arc Spectra**: + + .. code-block:: python + + wc = WavelengthCalibration1D(arc_spectra=[arc_he, arc_hg_ar], + line_lists=[["HeI"], ["HgI", "ArI"]]) + +* **Using Observed Line Positions**: If you have already identified the pixel centroids of + lines in your calibration spectrum, you can provide them directly via the ``obs_lines`` + argument (as a list of NumPy arrays). In this case, you *must* also provide the detector's pixel + boundaries using ``pix_bounds`` (a tuple like ``(min_pixel, max_pixel)``) and a reference + pixel, ``ref_pixel``, which serves as the anchor point for the polynomial fit. You also + need to provide the ``line_lists`` containing the potential matching catalog wavelengths. + + .. code-block:: python + + obs_he = np.array([105.3, 210.8, 455.1, 512.5, 680.2]) + pixel_bounds = (0, 1024) + + wc = WavelengthCalibration1D(ref_pixel=512, + obs_lines=obs_he, + line_lists=["HeI"], + pix_bounds=pixel_bounds) + +* **Using Observed Line Positions From Multiple Arcs**: + + .. code-block:: python + + obs_he = np.array([105.3, 210.8, 455.1, 512.5, 680.2]) + obs_hg_ar = np.array([234.2, 534.1, 768.2, 879.6]) + pixel_bounds = (0, 1024) + + wc = WavelengthCalibration1D(ref_pixel=512, + obs_lines=[obs_he, obs_hg_ar], + line_lists=[["HeI"], ["HgI", "ArI"]], + pix_bounds=pixel_bounds) + +2. Finding Observed Lines +************************* + +If you initialized the class with ``arc_spectra``, you need to detect the lines in it. Use the +:meth:`~specreduce.wavecal1d.WavelengthCalibration1D.find_lines` method: + +.. code-block:: python + + wc.find_lines(fwhm=3.5, noise_factor=5) + +This populates the `~specreduce.wavecal1d.WavelengthCalibration1D.observed_line_locations` +attribute. + +.. code-block:: python + + print(wc.observed_line_locations) + +3. Matching and Fitting the Solution +************************************ + +The core of the calibration process is fitting a model that maps pixels to wavelengths. + +* **Global Fitting for Automated Pipelines**: If you have + `~specreduce.wavecal1d.WavelengthCalibration1D.observed_lines` (either found automatically or + provided initially) and + `~specreduce.wavecal1d.WavelengthCalibration1D.catalog_lines` (from ``line_lists``), but don't + know the exact pixel-wavelength pairs, you can use + :meth:`~specreduce.wavecal1d.WavelengthCalibration1D.fit_dispersion`. This method applies the + `differential evolution optimization algorithm `_ + to find the polynomial coefficients that minimize the distance between observed line + positions (transformed to wavelength space) and the nearest catalog lines. + + The fit is anchored around the reference pixel (``ref_pixel``), which defines the centre + of the polynomial model. If not explicitly set at initialization, it defaults to the middle + of the detector when arc spectra are supplied. You must provide estimated bounds for both + the wavelength and the dispersion (dλ/dx) at ``ref_pixel``. + + The polynomial degree is set with the ``degree`` argument. A higher degree can capture + more complex dispersion behaviour but risks overfitting if too few calibration lines are + available. + + The ``higher_order_limits`` argument optionally constrains the absolute values of the + higher-order polynomial coefficients, reducing the risk of unrealistic solutions. + + The ``popsize`` argument controls the population size in the differential evolution search. + Larger values generally improve robustness at the expense of longer computation times. + + + .. code-block:: python + + ws = wc.fit_dispersion(wavelength_bounds=(7450, 7550), + dispersion_bounds=(1.8, 2.2), + higher_order_limits=[0.001, 1e-5], + degree=4, + popsize=30, + refine_fit=True) + + The fitted `~specreduce.wavesol1d.WavelengthSolution1D` is returned by the method and also + stored in ``WavelengthCalibration1D.solution``. Setting ``refine_fit=True`` + automatically performs a least-squares refinement after the global fit finds an initial + solution and matches lines. + +* **Fitting Known Pairs for an Interactive Workflow**: If you have already established explicit + pairs of observed pixel centers and their corresponding known wavelengths, you can use + :meth:`~specreduce.wavecal1d.WavelengthCalibration1D.fit_lines` to perform a direct + least-squares fit. + + .. code-block:: python + + ws = wc.fit_lines(pixels=[105.3, 512.5, 780.1], + wavelengths=[6965.43, 7503.87, 7723.76], + degree=2, + refine_fit=True) + + When ``refine_fit=True`` is set, the method automatically identifies matching pairs between + observed and catalog lines, then performs a least-squares refinement using **all matching lines**. + This goes beyond the subset of lines provided to :meth:`~specreduce.wavecal1d.WavelengthCalibration1D.fit_lines`, + resulting in an improved wavelength solution. + +4. Inspecting the Fit +********************* + +Several tools help assess the quality of the wavelength solution: + +* **RMS Error**: Calculate the root-mean-square error of the fit in wavelength or pixel units + using :meth:`~specreduce.wavecal1d.WavelengthCalibration1D.rms`. + + .. code-block:: python + + rms_wave = wc.rms(space='wavelength') + rms_pix = wc.rms(space='pixel') + print(f"Fit RMS (wavelength): {rms_wave}") + print(f"Fit RMS (pixel): {rms_pix}") + +* **Plotting**: Visualize the fit and residuals: + + * :meth:`~specreduce.wavecal1d.WavelengthCalibration1D.plot_fit`: Shows the observed line + positions mapped to the wavelength axis, overlaid with the catalog lines and the fitted + solution. Also shows the fit residuals (observed - fitted wavelength) vs. pixel. + * :meth:`~specreduce.wavecal1d.WavelengthCalibration1D.plot_residuals`: Plots residuals vs. + pixel or vs. wavelength. + * :meth:`~specreduce.wavecal1d.WavelengthCalibration1D.plot_observed_lines`: Plots the + identified observed line positions (in pixels or mapped to wavelengths). Can optionally + overlay the arc spectrum. + * :meth:`~specreduce.wavecal1d.WavelengthCalibration1D.plot_catalog_lines`: Plots the catalog + line positions (in wavelengths or mapped to pixels). + +5. Using the Wavelength Solution +******************************** + +The fitted wavelength solution is stored as a :class:`~specreduce.wavesol1d.WavelengthSolution1D` +instance that you can use to transform and resample spectra. The fitting methods return the +solution, but it is also stored in ``WavelengthCalibration1D.solution``. + +* **Convert Coordinates**: Use :meth:`~specreduce.wavesol1d.WavelengthSolution1D.pix_to_wav` and + :meth:`~specreduce.wavesol1d.WavelengthSolution1D.wav_to_pix` to convert between pixel and + wavelength coordinates. + + .. code-block:: python + + pixels = np.array([100, 500, 900]) + wavelengths = ws.pix_to_wav(pixels) + print(wavelengths) + +* **Get WCS Object**: Access the `~gwcs.wcs.WCS` object representing the solution via the + :attr:`~specreduce.wavesol1d.WavelengthSolution1D.gwcs` attribute. This is particularly + useful for attaching the calibration to a :class:`~specutils.Spectrum` object. + +* **Rebin Spectrum**: Resample a spectrum onto a new wavelength grid using + :meth:`~specreduce.wavesol1d.WavelengthSolution1D.resample`. The rebinning is + flux-conserving, meaning the integrated flux in the output spectrum matches the integrated flux + in the input spectrum. + + .. code-block:: python + + resampled_arc = ws.resample(arc_spectrum, nbins=1000) + print(resampled_arc.spectral_axis) + + By default, :meth:`~specreduce.wavesol1d.WavelengthSolution1D.resample` produces a spectrum on + a linear wavelength grid. Alternatively, you can pass an arbitrary 1D array of wavelength bin + edges via ``bin_edges`` to define a custom grid. + + .. code-block:: python + + resampled_arc = ws.resample(arc_spectrum, bin_edges=np.geomspace(4500, 7500, 1000)) + print(resampled_arc.spectral_axis) \ No newline at end of file diff --git a/pyproject.toml b/pyproject.toml index 52565e49..1640fc88 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -13,6 +13,7 @@ dependencies = [ "astropy>=5.3", "scipy>=1.10", "specutils>=1.9.1", + "matplotlib>=3.10", "gwcs", ] @@ -23,13 +24,16 @@ test = [ "tox", ] docs = [ - "sphinx-astropy", + "sphinx-astropy[confv2]", + "sphinx-copybutton", + "sphinx-design", "matplotlib>=3.7", "photutils>=1.0", "synphot", + "nbsphinx", + "ipykernel" ] all = [ - "matplotlib>=3.7", "photutils>=1.0", "synphot", ] diff --git a/specreduce/__init__.py b/specreduce/__init__.py index ad1f3f06..d703b86b 100644 --- a/specreduce/__init__.py +++ b/specreduce/__init__.py @@ -3,7 +3,6 @@ from specreduce.core import * # noqa from specreduce.wavelength_calibration import * # noqa - try: from .version import version as __version__ except ImportError: diff --git a/specreduce/background.py b/specreduce/background.py index 2390f199..2288fb31 100644 --- a/specreduce/background.py +++ b/specreduce/background.py @@ -210,13 +210,13 @@ def two_sided(cls, image, trace_object, separation, **kwargs): Image with 2-D spectral image data. Assumes cross-dispersion (spatial) direction is axis 0 and dispersion (wavelength) direction is axis 1. - trace_object: `~specreduce.tracing.Trace` + trace_object : `~specreduce.tracing.Trace` estimated trace of the spectrum to center the background traces - separation: float + separation : float separation from ``trace_object`` for the background regions width : float width of each background aperture in pixels - statistic: string + statistic : string statistic to use when computing the background. 'average' will account for partial pixel weights, 'median' will include all partial pixels. @@ -258,21 +258,21 @@ def one_sided(cls, image, trace_object, separation, **kwargs): Image with 2-D spectral image data. Assumes cross-dispersion (spatial) direction is axis 0 and dispersion (wavelength) direction is axis 1. - trace_object: `~specreduce.tracing.Trace` - estimated trace of the spectrum to center the background traces - separation: float - separation from ``trace_object`` for the background, positive will be + trace_object : `~specreduce.tracing.Trace` + Estimated trace of the spectrum to center the background traces + separation : float + Separation from ``trace_object`` for the background, positive will be above the trace, negative below. width : float - width of each background aperture in pixels - statistic: string - statistic to use when computing the background. 'average' will + Width of each background aperture in pixels + statistic : string + Statistic to use when computing the background. 'average' will account for partial pixel weights, 'median' will include all partial pixels. disp_axis : int - dispersion axis + Dispersion axis crossdisp_axis : int - cross-dispersion axis + Cross-dispersion axis mask_treatment : string The method for handling masked or non-finite data. Choice of ``filter``, ``omit``, or ``zero_fill``. If `filter` is chosen, masked/non-finite data @@ -344,17 +344,17 @@ def bkg_spectrum(self, image=None, bkg_statistic=None): def sub_image(self, image=None): """ - Subtract the computed background from ``image``. + Subtract the computed background from image. Parameters ---------- image : nddata-compatible image or None - image with 2-D spectral image data. If None, will extract + image with 2D spectral image data. If None, will extract the background from ``image`` used to initialize the class. Returns ------- - spec : `~specutils.Spectrum1D` + spec : `~specutils.Spectrum` Spectrum object with same shape as ``image``. """ image = self._parse_image(image) @@ -376,14 +376,14 @@ def sub_spectrum(self, image=None): Parameters ---------- - image : nddata-compatible image or None - image with 2-D spectral image data. If None, will extract + image : `~astropy.nddata.NDData`-compatible image or None + image with 2D spectral image data. If None, will extract the background from ``image`` used to initialize the class. Returns ------- - spec : `~specutils.Spectrum1D` - The background 1-D spectrum, with flux expressed in the same + spec : `~specutils.Spectrum` + The background 1D spectrum, with flux expressed in the same units as the input image (or u.DN if none were provided) and the spectral axis expressed in pixel units. """ diff --git a/specreduce/calibration_data.py b/specreduce/calibration_data.py index 8823037d..72d831d6 100644 --- a/specreduce/calibration_data.py +++ b/specreduce/calibration_data.py @@ -211,7 +211,7 @@ def load_MAST_calspec( Returns ------- - spectrum : If the spectrum can be loaded, return it as a `~specutils.Spectrum1D`. + spectrum : If the spectrum can be loaded, return it as a `~specutils.Spectrum`. Otherwise return None. The spectral_axis units are Å and the flux units are milli-Janskys. """ filename = Path(filename) @@ -266,7 +266,7 @@ def load_onedstds( Returns ------- - spectrum : If the spectrum can be loaded, return it as a `~specutils.Spectrum1D`. + spectrum : If the spectrum can be loaded, return it as a `~specutils.Spectrum`. Otherwise return None. The spectral_axis units are Å and the flux units are milli-Janskys. """ if dataset not in SPECPHOT_DATASETS: diff --git a/specreduce/compat.py b/specreduce/compat.py index 8d9b42fd..274098d6 100644 --- a/specreduce/compat.py +++ b/specreduce/compat.py @@ -1,7 +1,7 @@ import specutils from astropy.utils import minversion -__all__ = [] +__all__ = ["Spectrum"] SPECUTILS_LT_2 = not minversion(specutils, "2.0.dev") diff --git a/specreduce/extract.py b/specreduce/extract.py index 9db444f9..ea1add7a 100644 --- a/specreduce/extract.py +++ b/specreduce/extract.py @@ -420,12 +420,6 @@ def _parse_image(self, image, variance=None, mask=None, unit=None, disp_axis=1): if image.uncertainty.uncertainty_type == "var": variance = image.uncertainty.array elif image.uncertainty.uncertainty_type == "std": - warnings.warn( - "image NDData object's uncertainty " - "interpreted as standard deviation. if " - "incorrect, use VarianceUncertainty when " - "assigning image object's uncertainty." - ) variance = image.uncertainty.array**2 elif image.uncertainty.uncertainty_type == "ivar": variance = 1 / image.uncertainty.array diff --git a/specreduce/line_matching.py b/specreduce/line_matching.py index aa9ff5c6..ed48a1c4 100644 --- a/specreduce/line_matching.py +++ b/specreduce/line_matching.py @@ -1,8 +1,10 @@ import warnings +from copy import deepcopy from typing import Sequence import astropy.units as u import numpy as np +from astropy.nddata import StdDevUncertainty from astropy.stats import gaussian_fwhm_to_sigma, gaussian_sigma_to_fwhm from astropy.modeling import models from astropy.table import QTable @@ -53,6 +55,10 @@ def find_arc_lines( if fwhm.unit != spectrum.spectral_axis.unit: raise ValueError("fwhm must have the same units as spectrum.spectral_axis.") + if spectrum.uncertainty is None: + spectrum = deepcopy(spectrum) + spectrum.uncertainty = StdDevUncertainty(np.sqrt(np.abs(spectrum.flux.value))) + detected_lines = find_lines_threshold(spectrum, noise_factor=noise_factor) detected_lines = detected_lines[detected_lines['line_type'] == 'emission'] diff --git a/specreduce/tests/test_extract.py b/specreduce/tests/test_extract.py index 4dd88a49..069996e3 100644 --- a/specreduce/tests/test_extract.py +++ b/specreduce/tests/test_extract.py @@ -2,7 +2,7 @@ import pytest from astropy import units as u from astropy.modeling import models -from astropy.nddata import NDData, VarianceUncertainty, UnknownUncertainty, StdDevUncertainty +from astropy.nddata import NDData, VarianceUncertainty, UnknownUncertainty from astropy.tests.helper import assert_quantity_allclose from specreduce.background import Background @@ -137,12 +137,6 @@ def test_horne_image_validation(mk_test_img): with pytest.raises(ValueError, match=r".*NDData object lacks uncertainty"): ext = extract(image=image) - # a warning should be raised if uncertainty is StdDevUncertainty - with pytest.warns(UserWarning, match="image NDData object's uncertainty"): - err = StdDevUncertainty(np.ones_like(image)) - image.uncertainty = err - ext = extract(image=image) - # an NDData-type image's uncertainty must be of type VarianceUncertainty # or type StdDevUncertainty with pytest.raises(ValueError, match=r".*unexpected uncertainty type.*"): diff --git a/specreduce/tests/test_wavecal1d.py b/specreduce/tests/test_wavecal1d.py new file mode 100644 index 00000000..30b2fe77 --- /dev/null +++ b/specreduce/tests/test_wavecal1d.py @@ -0,0 +1,351 @@ +# specreduce/tests/test_wavecal1d.py +import astropy.units as u +import numpy as np +import pytest +from astropy.modeling import models +from astropy.nddata import StdDevUncertainty +from matplotlib import pyplot as plt +from matplotlib.figure import Figure +from numpy import array + +from specreduce.compat import Spectrum +from specreduce.wavecal1d import WavelengthCalibration1D + +ref_pixel = 250 +pix_bounds = (0, 500) +p2w = models.Shift(-ref_pixel) | models.Polynomial1D(degree=3, c0=1, c1=0.2, c2=0.001) + + +@pytest.fixture +def mk_lines(): + obs_lines = ref_pixel + array([1, 2, 5, 8, 10]) + cat_lines = p2w(ref_pixel + array([1, 3, 5, 7, 8, 10])) + return obs_lines, cat_lines + + +@pytest.fixture +def mk_matched_lines(): + obs_lines = ref_pixel + array([1, 2, 5, 8, 10]) + cat_lines = p2w(obs_lines) + return obs_lines, cat_lines + + +@pytest.fixture +def mk_wc(mk_lines): + obs_lines, cat_lines = mk_lines + return WavelengthCalibration1D( + obs_lines=obs_lines, + line_lists=cat_lines, + pix_bounds=(0, 10), + ref_pixel=ref_pixel, + ) + + +@pytest.fixture +def mk_good_wc_with_transform(mk_lines): + obs_lines, cat_lines = mk_lines + wc = WavelengthCalibration1D( + obs_lines=obs_lines, + line_lists=cat_lines, + pix_bounds=pix_bounds, + ref_pixel=ref_pixel, + ) + wc.solution.p2w = p2w + return wc + + +@pytest.fixture +def mk_arc(): + return Spectrum( + flux=np.ones(pix_bounds[1]) * u.DN, + spectral_axis=np.arange(pix_bounds[1]) * u.pix, + uncertainty=StdDevUncertainty(np.ones(pix_bounds[1])), + ) + + +def test_init(mk_arc, mk_lines): + arc = mk_arc + obs_lines, cat_lines = mk_lines + WavelengthCalibration1D(arc_spectra=arc, line_lists=cat_lines, ref_pixel=ref_pixel) + WavelengthCalibration1D( + obs_lines=obs_lines, + line_lists=cat_lines, + pix_bounds=(0, 10), + ref_pixel=ref_pixel, + ) + + with pytest.raises(ValueError, match="Only one of arc_spectra or obs_lines can be provided."): + WavelengthCalibration1D(arc_spectra=arc, obs_lines=obs_lines, ref_pixel=ref_pixel) + + arc = Spectrum( + flux=np.array([[1, 2, 3, 4, 5]]) * u.DN, + spectral_axis=np.array([1, 2, 3, 4, 5]) * u.angstrom, + ) + with pytest.raises(ValueError, match="The arc spectrum must be one dimensional."): + WavelengthCalibration1D(arc_spectra=arc, ref_pixel=ref_pixel) + + with pytest.raises(ValueError, match="Must give pixel bounds when"): + WavelengthCalibration1D(obs_lines=obs_lines, ref_pixel=ref_pixel) + + +def test_init_line_list(mk_arc): + arc = mk_arc + WavelengthCalibration1D(arc_spectra=arc, line_lists=["ArI"]) + WavelengthCalibration1D(arc_spectra=arc, line_lists="ArI") + WavelengthCalibration1D(arc_spectra=arc, line_lists=[array([0.1])]) + WavelengthCalibration1D(arc_spectra=arc, line_lists=array([0.1])) + WavelengthCalibration1D(arc_spectra=[arc, arc], line_lists=[["ArI"], ["ArI"]]) + WavelengthCalibration1D(arc_spectra=[arc, arc], line_lists=["ArI", ["ArI"]]) + WavelengthCalibration1D(arc_spectra=[arc, arc], line_lists=["ArI", "ArI"]) + WavelengthCalibration1D(arc_spectra=[arc, arc], line_lists=[array([0.1]), array([0.1])]) + WavelengthCalibration1D(arc_spectra=[arc, arc], line_lists=[array([0.1, 0.3]), ["ArI"]]) + with pytest.raises(ValueError, match="The number of line lists"): + WavelengthCalibration1D(arc_spectra=[arc, arc], line_lists=[["ArI"]]) + + +def test_find_lines(mocker, mk_arc): + wc = WavelengthCalibration1D(arc_spectra=mk_arc) + mock_find_arc_lines = mocker.patch("specreduce.wavecal1d.find_arc_lines") + mock_find_arc_lines.return_value = {"centroid": np.array([5.0]) * u.angstrom} + wc.find_lines(fwhm=2.0, noise_factor=1.5) + assert wc._obs_lines is not None + assert len(wc._obs_lines) == 1 + + wc = WavelengthCalibration1D() + with pytest.raises(ValueError, match="Must provide arc spectra to find lines."): + wc.find_lines(fwhm=2.0, noise_factor=1.5) + + +def test_fit_lines(mk_matched_lines): + lo, lc = mk_matched_lines + wc = WavelengthCalibration1D(pix_bounds=pix_bounds, ref_pixel=ref_pixel) + wc.fit_lines(pixels=lo, wavelengths=lc) + assert wc.solution.p2w is not None + assert wc.solution.p2w[1].degree == wc.degree + + wc = WavelengthCalibration1D( + obs_lines=lo, line_lists=lc, pix_bounds=pix_bounds, ref_pixel=ref_pixel + ) + wc.fit_lines(pixels=lo, wavelengths=lc, match_cat=True, match_obs=True) + assert wc.solution.p2w is not None + assert wc.solution.p2w[1].degree == wc.degree + + wc = WavelengthCalibration1D(pix_bounds=pix_bounds, ref_pixel=ref_pixel) + with pytest.warns(UserWarning, match="The degree of the polynomial"): + wc.fit_lines(degree=5, pixels=lo[:3], wavelengths=lc[:3]) + with pytest.raises(ValueError, match="Cannot fit without catalog"): + wc.fit_lines(pixels=lo, wavelengths=lc, match_cat=True, match_obs=True) + with pytest.raises(ValueError, match="Cannot fit without observed"): + wc.fit_lines(pixels=lo, wavelengths=lc, match_cat=False, match_obs=True) + + +def test_observed_lines(mk_lines): + wc = WavelengthCalibration1D() + assert wc.observed_lines is None + obs_lines, cat_lines = mk_lines + wc = WavelengthCalibration1D(obs_lines=obs_lines, pix_bounds=pix_bounds, ref_pixel=ref_pixel) + assert len(wc.observed_lines) == 1 + assert wc.observed_lines[0].shape == (len(obs_lines), 2) + assert wc.observed_lines[0].mask.shape == (len(obs_lines), 2) + np.testing.assert_allclose(wc.observed_lines[0].data[:, 0], obs_lines) + assert np.all(wc.observed_lines[0].mask[:, 0] == 0) + + wc.observed_lines = obs_lines + assert len(wc.observed_lines) == 1 + assert wc.observed_lines[0].mask.shape == (len(obs_lines), 2) + np.testing.assert_allclose(wc.observed_lines[0].data[:, 0], obs_lines) + assert np.all(wc.observed_lines[0].mask[:, 0] == 0) + + wc.observed_lines = wc.observed_lines + assert len(wc.observed_lines) == 1 + assert wc.observed_lines[0].mask.shape == (len(obs_lines), 2) + np.testing.assert_allclose(wc.observed_lines[0].data[:, 0], obs_lines) + assert np.all(wc.observed_lines[0].mask[:, 0] == 0) + + # Line locations and amplitudes + assert wc.observed_line_locations[0].shape == (len(obs_lines),) + np.testing.assert_allclose(wc.observed_line_locations[0], obs_lines) + assert wc.observed_line_amplitudes[0].shape == (len(obs_lines),) + np.testing.assert_allclose(wc.observed_line_amplitudes[0], 0.0) + + +def test_catalog_lines(mk_lines): + wc = WavelengthCalibration1D(ref_pixel=ref_pixel) + assert wc.catalog_lines is None + obs_lines, cat_lines = mk_lines + wc = WavelengthCalibration1D( + obs_lines=obs_lines, line_lists=cat_lines, pix_bounds=pix_bounds, ref_pixel=ref_pixel + ) + assert len(wc.catalog_lines) == 1 + assert wc.catalog_lines[0].shape == (len(cat_lines), 2) + np.testing.assert_allclose(wc.catalog_lines[0].data[:, 0], cat_lines) + assert np.all(wc.catalog_lines[0].data[:, 1] == 0) + + wc.catalog_lines = cat_lines + assert len(wc.catalog_lines) == 1 + np.testing.assert_allclose(wc.catalog_lines[0].data[:, 0], cat_lines) + assert np.all(wc.catalog_lines[0].data[:, 1] == 0) + + wc.catalog_lines = wc.catalog_lines + assert len(wc.catalog_lines) == 1 + np.testing.assert_allclose(wc.catalog_lines[0].data[:, 0], cat_lines) + assert np.all(wc.catalog_lines[0].data[:, 1] == 0) + + +def test_fit_lines_raises_error_for_mismatched_sizes(): + pixels = array([2, 4, 6]) + wavelengths = p2w(array([2, 4, 5, 6])) + wc = WavelengthCalibration1D(pix_bounds=pix_bounds, ref_pixel=ref_pixel) + with pytest.raises(ValueError, match="The sizes of pixel and wavelength arrays must match."): + wc.fit_lines(pixels=pixels, wavelengths=wavelengths) + + +def test_fit_lines_raises_error_for_insufficient_lines(): + pixels = [5] + wavelengths = p2w(pixels) + wc = WavelengthCalibration1D(pix_bounds=pix_bounds, ref_pixel=ref_pixel) + with pytest.raises(ValueError, match="Need at least two lines for a fit"): + wc.fit_lines(pixels=pixels, wavelengths=wavelengths) + + +def test_fit_lines_raises_error_for_missing_pixel_bounds(): + pixels = [2, 4, 6, 8] + wavelengths = p2w(pixels) + wc = WavelengthCalibration1D(ref_pixel=ref_pixel) + with pytest.raises(ValueError, match="Cannot fit without pixel bounds set."): + wc.fit_lines(pixels=pixels, wavelengths=wavelengths) + + +def test_fit_global(): + p2w = models.Shift(-ref_pixel) | models.Polynomial1D(degree=3, c0=650, c1=50.0, c2=-0.001) + + lines_obs = ref_pixel + array([2, 4, 4.5, 5, 6, 6.3, 8]) + lines_cat = p2w(ref_pixel + array([1, 2, 2.3, 4, 5, 6, 6.3, 7, 8, 9, 10])) + wavelength_bounds = (649, 651) + dispersion_bounds = (49, 51) + wc = WavelengthCalibration1D( + obs_lines=lines_obs, line_lists=lines_cat, pix_bounds=pix_bounds, ref_pixel=ref_pixel + ) + ws = wc.fit_dispersion(wavelength_bounds, dispersion_bounds, popsize=10, refine_fit=True) + np.testing.assert_allclose(ws.p2w[1].parameters, [650.0, 50.0, -0.001, 0.0], atol=1e-4) + assert wc._fit is not None + assert wc._fit.success + assert wc.solution.p2w is not None + + wc = WavelengthCalibration1D( + obs_lines=lines_obs, line_lists=lines_cat, pix_bounds=pix_bounds, ref_pixel=ref_pixel + ) + wc.fit_dispersion(wavelength_bounds, dispersion_bounds, popsize=10, refine_fit=False) + + +def test_rms(mk_good_wc_with_transform): + wc = mk_good_wc_with_transform + assert np.isclose(wc.rms(space="wavelength"), 0) # Perfect match, so RMS should be zero + assert np.isclose(wc.rms(space="pixel"), 0) # Perfect match, so RMS should be zero + with pytest.raises(ValueError, match="Space must be either 'pixel' or 'wavelength'"): + wc.rms(space="wavelenght") + + +def test_remove_unmatched_lines(mk_good_wc_with_transform): + wc = mk_good_wc_with_transform + wc.match_lines() + wc.remove_unmatched_lines() + assert wc.catalog_lines[0].size == wc.observed_lines[0].size + + +def test_plot_lines_with_valid_input(): + wc = WavelengthCalibration1D(ref_pixel=ref_pixel) + wc.observed_lines = [np.ma.masked_array([100, 200, 300], mask=[False, True, False])] + wc._cat_lines = wc.observed_lines + fig = wc._plot_lines(kind="observed", frames=0, figsize=(8, 4), plot_labels=True) + assert isinstance(fig, Figure) + assert fig.axes[0].has_data() + + fig = wc._plot_lines(kind="catalog", frames=0, figsize=(8, 4), plot_labels=True) + assert isinstance(fig, Figure) + assert fig.axes[0].has_data() + + fig, ax = plt.subplots(1, 1) + fig = wc._plot_lines(kind="catalog", frames=0, axes=ax, plot_labels=True) + assert isinstance(fig, Figure) + assert fig.axes[0].has_data() + + fig, axs = plt.subplots(1, 2) + fig = wc._plot_lines(kind="catalog", frames=0, axes=axs, plot_labels=True) + assert isinstance(fig, Figure) + assert fig.axes[0].has_data() + + fig = wc._plot_lines(kind="observed", frames=0, axes=axs, plot_labels=True) + assert isinstance(fig, Figure) + assert fig.axes[0].has_data() + + +def test_plot_lines_raises_for_missing_transform(mk_wc): + wc = mk_wc + with pytest.raises(ValueError, match="Cannot map between pixels and"): + wc._plot_lines(kind="observed", map_x=True) + + +def test_plot_lines_calls_transform_correctly(mk_good_wc_with_transform): + wc = mk_good_wc_with_transform + wc._plot_lines(kind="observed", map_x=True) + wc._plot_lines(kind="catalog", map_x=True) + + +def test_plot_catalog_lines(mk_wc): + wc = mk_wc + wc.catalog_lines = [np.ma.masked_array([400, 500, 600], mask=[False, True, False])] + fig = wc.plot_catalog_lines(frames=0, figsize=(10, 6), plot_labels=True, map_to_pix=False) + assert isinstance(fig, Figure) + assert fig.axes[0].has_data() + + fig, ax = plt.subplots(1, 1) + fig = wc.plot_catalog_lines(frames=0, axes=ax, plot_labels=True) + assert isinstance(fig, Figure) + assert fig.axes[0].has_data() + + fig, axs = plt.subplots(1, 2) + fig = wc.plot_catalog_lines(frames=[0], axes=axs, plot_labels=False) + assert isinstance(fig, Figure) + assert fig.axes[0].has_data() + + +def test_plot_observed_lines(mk_good_wc_with_transform, mk_arc): + wc = mk_good_wc_with_transform + wc.observed_lines = [np.ma.masked_array([100, 200, 300], mask=[False, True, False])] + wc.arc_spectra = [mk_arc] + for frames in [None, 0]: + fig = wc.plot_observed_lines(frames=frames, figsize=(10, 5), plot_labels=True) + assert isinstance(fig, Figure) + assert fig.axes[0].has_data() + assert len(fig.axes) == 1 + + +def test_plot_fit(mk_arc, mk_good_wc_with_transform): + wc = mk_good_wc_with_transform + wc.arc_spectra = [mk_arc] + for frames in [None, 0]: + fig = wc.plot_fit(frames=frames, figsize=(12, 6), plot_labels=True) + assert isinstance(fig, Figure) + assert len(fig.axes) == 2 + assert fig.axes[0].has_data() + assert fig.axes[1].has_data() + wc.plot_fit(frames=frames, figsize=(12, 6), plot_labels=True, obs_to_wav=True) + + +def test_plot_residuals(mk_good_wc_with_transform): + wc = mk_good_wc_with_transform + + fig = wc.plot_residuals(space="pixel", figsize=(8, 4)) + assert isinstance(fig, Figure) + assert fig.axes[0].has_data() + + fig = wc.plot_residuals(space="wavelength", figsize=(8, 4)) + assert isinstance(fig, Figure) + assert fig.axes[0].has_data() + + fig, ax = plt.subplots(1, 1) + wc.plot_residuals(ax=ax, space="wavelength", figsize=(8, 4)) + + with pytest.raises(ValueError, match="Invalid space specified"): + wc.plot_residuals(space="wavelenght", figsize=(8, 4)) diff --git a/specreduce/tests/test_wavesol1d.py b/specreduce/tests/test_wavesol1d.py new file mode 100644 index 00000000..d2cd42a2 --- /dev/null +++ b/specreduce/tests/test_wavesol1d.py @@ -0,0 +1,132 @@ +import astropy.units as u +import numpy as np +import pytest +from astropy.modeling import models +from astropy.modeling.polynomial import Polynomial1D +from astropy.nddata import StdDevUncertainty +from gwcs import wcs +from specreduce.wavesol1d import _diff_poly1d, WavelengthSolution1D +from specreduce.compat import Spectrum + + +ref_pixel = 250.0 +p2w = models.Shift(ref_pixel) | models.Polynomial1D(degree=3, c0=1, c1=0.2, c2=0.001) +pix_bounds = (0, 500) +wav_bounds = p2w(pix_bounds) + + +@pytest.fixture +def mk_ws_without_transform(): + return WavelengthSolution1D(None, pix_bounds, u.angstrom) + + +@pytest.fixture +def mk_ws_with_transform(): + return WavelengthSolution1D(p2w, pix_bounds, u.angstrom) + + +@pytest.fixture +def mk_spectrum(): + return Spectrum( + flux=np.ones(pix_bounds[1]) * u.DN, + spectral_axis=np.arange(pix_bounds[1]) * u.pix, + uncertainty=StdDevUncertainty(np.ones(pix_bounds[1])), + ) + + +def test_diff_poly1d(): + p = _diff_poly1d(Polynomial1D(3, c0=1.0, c1=2.0, c2=3.0, c3=4.0)) + np.testing.assert_array_equal(p.parameters, [2.0, 6.0, 12.0]) + + +def test_init(): + ws = WavelengthSolution1D(p2w, pix_bounds, u.angstrom) + assert ws._p2w is p2w + assert ws.bounds_pix == pix_bounds + assert ws.unit == u.angstrom + assert "w2p" not in ws.__dict__ + assert "p2d_dldx" not in ws.__dict__ + assert "gwcs" not in ws.__dict__ + + # Test that the cached properties are created correctly + ws.w2p(0.5) + assert "w2p" in ws.__dict__ + ws.p2w_dldx(pix_bounds[0]) + assert "p2w_dldx" in ws.__dict__ + wcs = ws.gwcs # noqa: F841 + assert "gwcs" in ws.__dict__ + + # Test that the cached properties are deleted correctly + ws.p2w = p2w + assert "w2p" not in ws.__dict__ + assert "p2d_dldx" not in ws.__dict__ + assert "gwcs" not in ws.__dict__ + + ws = WavelengthSolution1D(p2w, pix_bounds, u.micron) + assert ws.unit == u.micron + + ws = WavelengthSolution1D(None, pix_bounds, u.angstrom) + assert ws._p2w is None + + +def test_resample(mk_spectrum, mk_ws_with_transform, mk_ws_without_transform): + ws = mk_ws_with_transform + spectrum = mk_spectrum + + # Resample a spectrum with uncertainty + resampled = ws.resample(spectrum, nbins=50) + assert resampled is not None + assert len(resampled.flux) == 50 + assert resampled.flux.unit == u.DN / u.angstrom + + pix_edges = np.arange(spectrum.spectral_axis.size + 1) - 0.5 + f0 = (spectrum.flux.value * np.diff(ws._p2w(pix_edges))).sum() + f1 = (resampled.flux.value * np.diff(resampled.spectral_axis.value)[0]).sum() + np.testing.assert_approx_equal(f0, f1, 5) + + resampled = ws.resample(spectrum, wlbounds=wav_bounds) + resampled = ws.resample(spectrum, bin_edges=np.linspace(*wav_bounds, num=50)) + + # Resample a spectrum without uncertainty + spectrum.uncertainty = None + resampled = ws.resample(spectrum, nbins=50) + assert resampled.uncertainty is not None + + ws = mk_ws_without_transform + with pytest.raises(ValueError, match="Wavelength solution not set."): + ws.resample(mk_spectrum) + + ws = mk_ws_with_transform + with pytest.raises(ValueError, match="Number of bins must be positive"): + ws.resample(mk_spectrum, nbins=-5) + + +def test_pix_to_wav(mk_ws_with_transform): + ws = mk_ws_with_transform + pix = np.array([1, 2, 3, 4, 5]) + np.testing.assert_array_equal(ws.pix_to_wav(pix), p2w(pix)) + + pix = np.ma.masked_array([1, 2, 3], mask=[0, 1, 0]) + wav = ws.pix_to_wav(pix) + np.testing.assert_array_equal(wav.data, p2w(pix.data)) + np.testing.assert_array_equal(wav.mask, np.array([0, 1, 0])) + + +def test_wav_to_pix(mk_ws_with_transform): + ws = mk_ws_with_transform + pix_values_orig = np.array([1, 2, 3, 4, 5]) + pix_values_tran = ws.wav_to_pix(ws.pix_to_wav(pix_values_orig)) + np.testing.assert_array_almost_equal(pix_values_orig, pix_values_tran) + + pix_values_orig = np.ma.masked_array([1, 2, 3, 4, 5], mask=[0, 1, 0, 1, 0]) + pix_values_tran = ws.wav_to_pix(ws.pix_to_wav(pix_values_orig)) + np.testing.assert_array_almost_equal(pix_values_orig.data, pix_values_tran.data) + np.testing.assert_array_almost_equal(pix_values_orig.mask, pix_values_tran.mask) + + +def test_wcs_creates_valid_gwcs_object(mk_ws_with_transform): + wc = mk_ws_with_transform + wcs_obj = wc.gwcs + assert wcs_obj is not None + assert isinstance(wcs_obj, wcs.WCS) + assert wcs_obj.output_frame.unit[0] == u.angstrom diff --git a/specreduce/utils/synth_data.py b/specreduce/utils/synth_data.py index 64b3dacb..a270db93 100644 --- a/specreduce/utils/synth_data.py +++ b/specreduce/utils/synth_data.py @@ -24,7 +24,7 @@ def make_2d_trace_image( background: int | float = 5, trace_center: int | float | None = None, trace_order: int = 3, - trace_coeffs: dict[str, int | float] = {'c0': 0, 'c1': 50, 'c2': 100}, + trace_coeffs: None | dict[str, int | float] = None, profile: Model = models.Moffat1D(amplitude=10, alpha=0.1), add_noise: bool = True ) -> CCDData: @@ -35,27 +35,37 @@ def make_2d_trace_image( Parameters ---------- - nx : Size of image in X axis which is assumed to be the dispersion axis + nx + Size of image in X axis which is assumed to be the dispersion axis - ny : Size of image in Y axis which is assumed to be the spatial axis + ny + Size of image in Y axis which is assumed to be the spatial axis - background : Level of constant background in counts + background + Level of constant background in counts - trace_center : Zeropoint of the trace. If None, then use center of Y (spatial) axis. + trace_center + Zeropoint of the trace. If None, then use center of Y (spatial) axis. - trace_order : Order of the Chebyshev polynomial used to model the source's trace + trace_order + Order of the Chebyshev polynomial used to model the source's trace - trace_coeffs : Dict containing the Chebyshev polynomial coefficients to use in the trace model + trace_coeffs + Dict containing the Chebyshev polynomial coefficients to use in the trace model - profile : Model to use for the source's spatial profile - - add_noise : If True, add Poisson noise to the image + profile + Model to use for the source's spatial profile + add_noise + If True, add Poisson noise to the image Returns ------- - ccd_im : CCDData instance containing synthetic 2D spectroscopic image + ccd_im + `~astropy.nddata.CCDData` instance containing synthetic 2D spectroscopic image """ + if trace_coeffs is None: + trace_coeffs = {'c0': 0, 'c1': 50, 'c2': 100} x = np.arange(nx) y = np.arange(ny) xx, yy = np.meshgrid(x, y) @@ -87,7 +97,7 @@ def make_2d_arc_image( wave_air: bool = False, background: int | float = 5, line_fwhm: float = 5., - linelists: list[str] = ['HeI'], + linelists: list[str] = ('HeI',), amplitude_scale: float = 1., tilt_func: Model = models.Legendre1D(degree=0), add_noise: bool = True @@ -101,36 +111,50 @@ def make_2d_arc_image( Parameters ---------- - nx : Size of image in X axis which is assumed to be the dispersion axis + nx + Size of image in X axis which is assumed to be the dispersion axis - ny : Size of image in Y axis which is assumed to be the spatial axis + ny + Size of image in Y axis which is assumed to be the spatial axis - wcs : 2D WCS to apply to the image. Must have a spectral axis defined along with + wcs + 2D WCS to apply to the image. Must have a spectral axis defined along with appropriate spectral wavelength units. - extent : If ``wcs`` is not provided, this defines the beginning and end wavelengths + extent + If ``wcs`` is not provided, this defines the beginning and end wavelengths of the dispersion axis. - wave_unit : If ``wcs`` is not provided, this defines the wavelength units of the + wave_unit + If ``wcs`` is not provided, this defines the wavelength units of the dispersion axis. - wave_air : If True, convert the vacuum wavelengths used by ``pypeit`` to air wavelengths. + wave_air + If True, convert the vacuum wavelengths used by ``pypeit`` to air wavelengths. - background : Level of constant background in counts + background + Level of constant background in counts - line_fwhm : Gaussian FWHM of the emission lines in pixels + line_fwhm + Gaussian FWHM of the emission lines in pixels - linelists : Specification for linelists to load from ``pypeit`` + linelists + Specification for linelists to load from ``pypeit`` - amplitude_scale : Scale factor to apply to amplitudes provided in the linelists + amplitude_scale + Scale factor to apply to amplitudes provided in the linelists - tilt_func : The tilt function to apply along the cross-dispersion axis to simulate + tilt_func + The tilt function to apply along the cross-dispersion axis to simulate tilted or curved emission lines. - add_noise : If True, add Poisson noise to the image; requires ``photutils`` to be installed. + + add_noise + If True, add Poisson noise to the image; requires ``photutils`` to be installed. Returns ------- - ccd_im : CCDData instance containing synthetic 2D spectroscopic image + ccd_im + `~astropy.nddata.CCDData` instance containing synthetic 2D spectroscopic image Examples -------- @@ -326,12 +350,12 @@ def make_2d_spec_image( wave_air: bool = False, background: int | float = 5, line_fwhm: float = 5., - linelists: list[str] = ['OH_GMOS'], + linelists: list[str] = ('OH_GMOS',), amplitude_scale: float = 1., tilt_func: Model = models.Legendre1D(degree=0), trace_center: int | float | None = None, trace_order: int = 3, - trace_coeffs: dict[str, int | float] = {'c0': 0, 'c1': 50, 'c2': 100}, + trace_coeffs: None | dict[str, int | float] = None, source_profile: Model = models.Moffat1D(amplitude=10, alpha=0.1), add_noise: bool = True ) -> CCDData: @@ -341,47 +365,67 @@ def make_2d_spec_image( Parameters ---------- - nx : Number of pixels in the dispersion direction. + nx + Number of pixels in the dispersion direction. - ny : Number of pixels in the spatial direction. + ny + Number of pixels in the spatial direction. - wcs : 2D WCS to apply to the image. Must have a spectral axis defined along with + wcs + 2D WCS to apply to the image. Must have a spectral axis defined along with appropriate spectral wavelength units. - extent : If ``wcs`` is not provided, this defines the beginning and end wavelengths + extent + If ``wcs`` is not provided, this defines the beginning and end wavelengths of the dispersion axis. - wave_unit : If ``wcs`` is not provided, this defines the wavelength units of the + wave_unit + If ``wcs`` is not provided, this defines the wavelength units of the dispersion axis. - wave_air : If True, convert the vacuum wavelengths used by ``pypeit`` to air wavelengths. + wave_air + If True, convert the vacuum wavelengths used by ``pypeit`` to air wavelengths. - background : Constant background level in counts. + background + Constant background level in counts. - line_fwhm : Gaussian FWHM of the emission lines in pixels + line_fwhm + Gaussian FWHM of the emission lines in pixels - linelists : Specification for linelists to load from ``pypeit`` + linelists + Specification for linelists to load from ``pypeit`` - amplitude_scale : Scale factor to apply to amplitudes provided in the linelists + amplitude_scale + Scale factor to apply to amplitudes provided in the linelists - tilt_func : The tilt function to apply along the cross-dispersion axis to simulate + tilt_func + The tilt function to apply along the cross-dispersion axis to simulate tilted or curved emission lines. - trace_center : Zeropoint of the trace. If None, then use center of Y (spatial) axis. + trace_center + Zeropoint of the trace. If None, then use center of Y (spatial) axis. - trace_order : Order of the Chebyshev polynomial used to model the source's trace + trace_order + Order of the Chebyshev polynomial used to model the source's trace - trace_coeffs : Dict containing the Chebyshev polynomial coefficients to use in the trace model + trace_coeffs + Dict containing the Chebyshev polynomial coefficients to use in the trace model - source_profile : Model to use for the source's spatial profile + source_profile + Model to use for the source's spatial profile - add_noise : If True, add Poisson noise to the image; requires ``photutils`` to be installed. + add_noise + If True, add Poisson noise to the image; requires ``photutils`` to be installed. Returns ------- - ccd_im : CCDData instance containing synthetic 2D spectroscopic image + ccd_im + `~astropy.nddata.CCDData` instance containing synthetic 2D spectroscopic image """ + if trace_coeffs is None: + trace_coeffs = {'c0': 0, 'c1': 50, 'c2': 100} + arc_image = make_2d_arc_image( nx=nx, ny=ny, diff --git a/specreduce/wavecal1d.py b/specreduce/wavecal1d.py new file mode 100644 index 00000000..3c7cc382 --- /dev/null +++ b/specreduce/wavecal1d.py @@ -0,0 +1,1142 @@ +import warnings +from functools import cached_property +from math import isclose +from typing import Sequence, Literal + +import astropy.units as u +import numpy as np +from astropy.modeling import models, Model, fitting, CompoundModel +from matplotlib.pyplot import Axes, Figure, setp, subplots +from numpy.ma import MaskedArray +from numpy.typing import ArrayLike +from scipy import optimize, ndimage +from scipy.spatial import KDTree + +from specreduce.calibration_data import load_pypeit_calibration_lines +from specreduce.compat import Spectrum +from specreduce.line_matching import find_arc_lines + +__all__ = ["WavelengthCalibration1D"] + +from specreduce.wavesol1d import WavelengthSolution1D + + +def _format_linelist(lst: ArrayLike) -> MaskedArray: + """Force a line list into a MaskedArray with a shape of (n, 2) where n is the number of lines. + + Parameters + ---------- + lst + Input array of centroids or centroids with amplitudes. Must be either: + - A 1D array with a shape [n] for centroids. + - A 2D array with a shape [n, 2] for centroids and amplitudes. + + Returns + ------- + numpy.ma.MaskedArray + Formatted and standardized line list array with shape [n, 2], where each row + contains a line centroid and amplitude. + + Raises + ------ + ValueError + If the input line list does not meet the specified dimensional or shape + requirements. + """ + lst: MaskedArray = MaskedArray(lst, copy=True) + lst.mask = np.ma.getmaskarray(lst) + + if lst.ndim > 2 or lst.ndim == 2 and lst.shape[1] > 2: + raise ValueError( + "Line lists must be 1D with a shape [n] (centroids) or " + "2D with a shape [n, 2] (centroids and amplitudes)." + ) + + if lst.ndim == 1: + lst = MaskedArray(np.tile(lst[:, None], [1, 2])) + lst[:, 1] = 0.0 + lst.mask[:, :] = lst.mask.any(axis=1)[:, None] + + return lst[np.argsort(lst.data[:, 0])] + + +def _unclutter_text_boxes(labels: Sequence) -> None: + """Remove overlapping labels from the plot. + + Removes overlapping text labels from a set of matplotlib label objects. The function iterates + over all combinations of labels, checks for overlaps among their bounding boxes, and removes + the label with the lower z-order in case of an overlap. + + Parameters + ---------- + labels + A list of matplotlib.text.Text objects. + """ + to_remove = set() + for i in range(len(labels)): + for j in range(i + 1, len(labels)): + l1 = labels[i] + l2 = labels[j] + bbox1 = l1.get_window_extent() + bbox2 = l2.get_window_extent() + if bbox1.overlaps(bbox2): + if l1.zorder < l2.zorder: + to_remove.add(l1) + else: + to_remove.add(l2) + + for label in to_remove: + label.remove() + + +class WavelengthCalibration1D: + def __init__( + self, + arc_spectra: Spectrum | Sequence[Spectrum] | None = None, + obs_lines: ArrayLike | Sequence[ArrayLike] | None = None, + line_lists: ArrayLike | None = None, + unit: u.Unit = u.angstrom, + ref_pixel: float | None = None, + pix_bounds: tuple[int, int] | None = None, + line_list_bounds: None | tuple[float, float] = None, + n_strogest_lines: None | int = None, + wave_air: bool = False, + ) -> None: + """A class for wavelength calibration of one-dimensional spectra. + + This class is designed to facilitate wavelength calibration of one-dimensional spectra, + with support for both direct input of line lists and observed spectra. It uses a polynomial + model for fitting the wavelength solution and offers features to incorporate catalog lines + and observed line positions. + + Parameters + ---------- + arc_spectra + Arc spectra provided as ``Spectrum`` objects for wavelength fitting, by default + None. This parameter and ``obs_lines`` cannot be provided simultaneously. + + obs_lines + Pixel positions of observed spectral lines for wavelength fitting, by default None. This + parameter and ``arc_spectra`` cannot be provided simultaneously. + + line_lists + Catalogs of spectral line wavelengths for wavelength calibration. Provide either an + array of line wavelengths or a list of `PypeIt `_ + catalog names. If `None`, no line lists are used. You can query the list of available + catalog names via `~specreduce.calibration_data.get_available_line_catalogs`. + + unit + The unit of the wavelength calibration, by default ``astropy.units.Angstrom``. + + ref_pixel + The reference pixel in which the wavelength solution will be centered. + + pix_bounds + Lower and upper pixel bounds for fitting, defined as a tuple (min, max). If + ``obs_lines`` is provided, this parameter is mandatory. + + line_list_bounds + Wavelength bounds as a tuple (min, max) for filtering usable spectral + lines from the provided line lists. + + n_strogest_lines + The number of strongest lines to be included from the line lists. If `None`, all + are included. + + wave_air + Boolean indicating whether the input wavelengths correspond to air rather than vacuum; + by default `False`, meaning vacuum wavelengths. + """ + self.unit = unit + self._unit_str = unit.to_string("latex") + self.degree = None + self.ref_pixel = ref_pixel + self.nframes = 0 + + if ref_pixel is not None and ref_pixel < 0: + raise ValueError("Reference pixel must be positive.") + + self.arc_spectra: list[Spectrum] | None = None + self.bounds_pix: tuple[int, int] | None = pix_bounds + self.bounds_wav: tuple[float, float] | None = None + self._cat_lines: list[MaskedArray] | None = None + self._obs_lines: list[MaskedArray] | None = None + self._trees: list[KDTree] | None = None + + self._fit: optimize.OptimizeResult | None = None + self.solution = WavelengthSolution1D(None, pix_bounds, unit) + + # Read and store the observational data if given. The user can provide either a list of arc + # spectra as Spectrum objects or a list of line pixel position arrays. An attempt to give + # both raises an error. + if arc_spectra is not None and obs_lines is not None: + raise ValueError("Only one of arc_spectra or obs_lines can be provided.") + + if arc_spectra is not None: + self.arc_spectra = [arc_spectra] if isinstance(arc_spectra, Spectrum) else arc_spectra + self.nframes = len(self.arc_spectra) + for s in self.arc_spectra: + if s.data.ndim > 1: + raise ValueError("The arc spectrum must be one dimensional.") + if len(set([s.data.size for s in self.arc_spectra])) != 1: + raise ValueError("All arc spectra must have the same length.") + self.bounds_pix = (0, self.arc_spectra[0].shape[0]) + self.solution.bounds_pix = self.bounds_pix + if self.ref_pixel is None: + self.ref_pixel = self.arc_spectra[0].data.size / 2 + + elif obs_lines is not None: + self.observed_lines = obs_lines + self.nframes = len(self._obs_lines) + if self.bounds_pix is None: + raise ValueError("Must give pixel bounds when providing observed line positions.") + if self.ref_pixel is None: + raise ValueError("Must give reference pixel when providing observed lines.") + + # Read the line lists if given. The user can provide an array of line wavelength positions + # or a list of line list names (used by `load_pypeit_calibration_lines`) for each arc + # spectrum. + if line_lists is not None: + if not isinstance(line_lists, (tuple, list)): + line_lists = [line_lists] + + if len(line_lists) != self.nframes: + raise ValueError("The number of line lists must match the number of arc spectra.") + self._read_linelists( + line_lists, + line_list_bounds=line_list_bounds, + wave_air=wave_air, + n_strongest=n_strogest_lines, + ) + + def _line_match_distance(self, x: ArrayLike, model: Model, max_distance: float = 100) -> float: + """Compute the sum of distances between catalog lines and transformed observed lines. + + This function evaluates the pixel-to-wavelength model at the observed line positions, + queries the nearest catalog line (via KDTree), and sums the distances after clipping + them at `max_distance`. The result is suitable as a scalar objective for global + optimization of the wavelength solution. + + Parameters + ---------- + x + Pixel-to-wavelength model parameters (e.g., Polynomial1D coefficients c0..cN). + model + The pixel-to-wavelength model to be evaluated. + max_distance + Upper bound used to clip individual distances before summation. + + Returns + ------- + float + Sum of nearest-neighbor distances between transformed observed lines and catalog lines. + + """ + total_distance = 0.0 + for t, l in zip(self._trees, self.observed_line_locations): + transformed_lines = model.evaluate(l, -self.ref_pixel, *x)[:, None] + total_distance += np.clip(t.query(transformed_lines)[0], 0, max_distance).sum() + return total_distance + + def _read_linelists( + self, + line_lists: Sequence, + line_list_bounds: None | tuple[float, float] = None, + wave_air: bool = False, + n_strongest: None | int = None, + ) -> None: + """Read and processes line lists. + + Parameters + ---------- + line_lists + A collection of line lists that can either be arrays of wavelengths or `PypeIt + `_ + lamp names. You can query the list of available catalog names via + `~specreduce.calibration_data.get_available_line_catalogs`. + + line_list_bounds + A tuple specifying the minimum and maximum wavelength bounds. Only wavelengths + within this range are retained. + + wave_air + If True, convert the vacuum wavelengths used by `PypeIt + `_ to air wavelengths. + + n_strongest + The number of strongest lines to be used. If `None`, all lines are used. + """ + + lines = [] + for lst in line_lists: + if isinstance(lst, np.ndarray): + lines.append(lst) + else: + if isinstance(lst, str): + lst = [lst] + lines.append([]) + for ll in lst: + line_table = load_pypeit_calibration_lines(ll, wave_air=wave_air) + if n_strongest is not None: + ix = np.argsort(line_table["amplitude"].value)[::-1] + lines[-1].append(line_table[ix][:n_strongest]["wavelength"].to( + self.unit).value) + else: + lines[-1].append(line_table["wavelength"].to(self.unit).value) + lines[-1] = np.ma.masked_array(np.sort(np.concatenate(lines[-1]))) + + if line_list_bounds is not None: + for i, lst in enumerate(lines): + lines[i] = lst[(lst >= line_list_bounds[0]) & (lst <= line_list_bounds[1])] + + self.catalog_lines = lines + self._create_trees() + + def _create_trees(self) -> None: + """Initialize the KDTree instances for the current set of catalog line locations.""" + self._trees = [KDTree(lst.compressed()[:, None]) for lst in self.catalog_line_locations] + + def find_lines(self, fwhm: float, noise_factor: float = 1.0) -> None: + """Find lines in the provided arc spectra. + + Determines the spectral lines within each spectrum of the arc spectra based on the + provided initial guess for the line Full Width at Half Maximum (FWHM). + + Parameters + ---------- + fwhm + Initial guess for the FWHM for the spectral lines, used as a parameter in + the ``find_arc_lines`` function to locate and identify spectral arc lines. + + noise_factor + The factor to multiply the uncertainty by to determine the noise threshold + in the `~specutils.fitting.find_lines_threshold` routine. + """ + if self.arc_spectra is None: + raise ValueError("Must provide arc spectra to find lines.") + + line_lists = [] + for i, arc in enumerate(self.arc_spectra): + lines = find_arc_lines(arc, fwhm, noise_factor=noise_factor) + ix = np.round(lines["centroid"].value).astype(int) + if np.any((ix < 0) | (ix >= arc.shape[0])): + raise ValueError( + "Error in arc line identification. Try increasing ``noise_factor``." + ) + amplitudes = ndimage.maximum_filter1d(arc.flux.value, 5)[ix] + line_lists.append( + np.ma.masked_array(np.transpose([lines["centroid"].value, amplitudes])) + ) + self.observed_lines = line_lists + + def _create_model(self, degree: int, coeffs: None | ArrayLike = None) -> CompoundModel: + """Initialize the polynomial model with the given degree and an optional base model. + + This method sets up a polynomial transformation based on the reference pixel and degree. + If coefficients are provided, they are copied to the initialized model up to the degree + specified. + + Parameters + ---------- + degree + Degree of the polynomial model to be initialized. + + coeffs + Optional initial polynomial coefficients. + """ + self.degree = degree + + pars = {} + if coeffs is not None: + nc = min(degree + 1, len(coeffs)) + pars = {f"c{i}": c for i, c in enumerate(coeffs[:nc])} + return models.Shift(-self.ref_pixel) | models.Polynomial1D(self.degree, **pars) + + def fit_lines( + self, + pixels: ArrayLike, + wavelengths: ArrayLike, + degree: int = 3, + match_obs: bool = False, + match_cat: bool = False, + refine_fit: bool = True, + refine_max_distance: float = 5.0, + refined_fit_degree: int | None = None, + ) -> WavelengthSolution1D: + """Fit the pixel-to-wavelength model using provided line pairs. + + This method fits the pixel-to-wavelength transformation using explicitly provided pairs + of pixel coordinates and their corresponding wavelengths via a linear least-squares fit + + Optionally, the provided pixel and wavelength values can be "snapped" to the nearest + values present in the internally stored observed line list and catalog line list, + respectively. This allows the inputs to be approximate, as the snapping step selects + the nearest precise centroids and catalog values when available. + + Parameters + ---------- + pixels + An array of pixel positions corresponding to known spectral lines. + + wavelengths + An array of the same size as ``pixels``, containing the known + wavelengths corresponding to the given pixel positions. + + degree + The polynomial degree for the wavelength solution. + + match_obs + If True, snap the input ``pixels`` values to the nearest + pixel values found in ``self.observed_line_locations`` (if available). This helps + ensure the fit uses the precise centroids detected by `find_lines` + or provided initially. + + match_cat + If True, snap the input ``wavelengths`` values to the + nearest wavelength values found in ``self.catalog_line_locations`` (if available). + This ensures the fit uses the precise catalog wavelengths. + + refine_fit + If True (default), automatically call the ``refine_fit`` method + immediately after the global optimization to improve the solution + using a least-squares fit on matched lines. + + refine_max_distance + Maximum allowed separation between catalog and observed lines for them to + be considered a match during ``refine_fit``. Ignored if ``refine_fit`` is False. + + refined_fit_degree + The polynomial degree for the refined fit. Can be higher than ``degree``. If ``None``, + equals to ``degree``. + """ + pixels = np.asarray(pixels) + wavelengths = np.asarray(wavelengths) + if pixels.size != wavelengths.size: + raise ValueError("The sizes of pixel and wavelength arrays must match.") + nlines = pixels.size + if nlines < 2: + raise ValueError("Need at least two lines for a fit") + + if self.bounds_pix is None: + raise ValueError("Cannot fit without pixel bounds set.") + + # Match the input wavelengths to catalog lines. + if match_cat: + if self._cat_lines is None: + raise ValueError("Cannot fit without catalog lines set.") + tree = KDTree( + np.concatenate([c.compressed() for c in self.catalog_line_locations])[:, None] + ) + ix = tree.query(wavelengths[:, None])[1] + wavelengths = tree.data[ix][:, 0] + + # Match the input pixel values to observed pixel values. + if match_obs: + if self._obs_lines is None: + raise ValueError("Cannot fit without observed lines set.") + tree = KDTree( + np.concatenate([c.compressed() for c in self.observed_line_locations])[:, None] + ) + ix = tree.query(pixels[:, None])[1] + pixels = tree.data[ix][:, 0] + + fitter = fitting.LinearLSQFitter() + shift, model = self._create_model(degree) + if model.degree > nlines: + warnings.warn( + "The degree of the polynomial model is higher than the number of lines. " + "Fixing the higher-order coefficients to zero." + ) + for i in range(nlines, model.degree + 1): + model.fixed[f"c{i}"] = True + model = fitter(model, pixels - self.ref_pixel, wavelengths) + for i in range(model.degree + 1): + model.fixed[f"c{i}"] = False + self.solution.p2w = shift | model + + can_match = self._cat_lines is not None and self._obs_lines is not None + if refine_fit and can_match: + self.refine_fit(refined_fit_degree, max_match_distance=refine_max_distance) + else: + if can_match: + self.match_lines() + return self.solution + + def fit_dispersion( + self, + wavelength_bounds: tuple[float, float], + dispersion_bounds: tuple[float, float], + higher_order_limits: Sequence[float] | None = None, + degree: int = 3, + popsize: int = 30, + max_distance: float = 100, + refine_fit: bool = True, + refine_max_distance: float = 5.0, + refined_fit_degree: int | None = None, + ) -> WavelengthSolution1D: + """Calculate a wavelength solution using all the catalog and observed lines. + + This method estimates a wavelength solution without pre-matched pixel–wavelength + pairs, making it suitable for automated pipelines on stable, well-characterized + spectrographs. It uses differential evolution to optimize the polynomial parameters + that minimize the distance between the predicted wavelengths of the observed lines + and their nearest catalog lines. The resulting solution can optionally be refined + with a least-squares fit to automatically matched lines. + + Parameters + ---------- + wavelength_bounds + (min, max) bounds for the wavelength at ``ref_pixel``; used as an optimization + constraint. + + dispersion_bounds + (min, max) bounds for the dispersion d(wavelength)/d(pixel) at ``ref_pixel``; used + as an optimization constraint. + + higher_order_limits + Absolute limits for the higher-order polynomial coefficients. Each coefficient is + constrained to [-limit, limit]. If provided, the number of limits must equal + (polynomial degree - 1). + + degree + The polynomial degree for the wavelength solution. + + popsize + Population size for ``scipy.optimize.differential_evolution``. Larger values can + improve the chance of finding the global minimum at the cost of additional time. + + max_distance + Maximum wavelength separation used when associating observed and catalog lines in + the optimization. Distances larger than this threshold are clipped to this value + in the cost function to limit the impact of outliers. + + refine_fit + If True (default), call ``refine_fit`` after global optimization to improve the + solution using a least-squares fit on matched lines. + + refine_max_distance + Maximum allowed separation between catalog and observed lines for them to + be considered a match during ``refine_fit``. Ignored if ``refine_fit`` is False. + + refined_fit_degree + The polynomial degree for the refined fit. Can be higher than ``degree``. If ``None``, + equals to ``degree``. + """ + + # Define bounds for differential_evolution. + bounds = [np.asarray(wavelength_bounds), np.asarray(dispersion_bounds)] + model = self._create_model(degree) + + if higher_order_limits is not None: + if len(higher_order_limits) != model[1].degree - 1: + raise ValueError( + "The number of higher-order limits must match the degree of the polynomial " + "model minus one." + ) + for v in higher_order_limits: + bounds.append(np.asarray([-v, v])) + else: + for i in range(2, model[1].degree + 1): + bounds.append( + np.array([-1, 1]) * 10 ** (np.log10(np.mean(dispersion_bounds)) - 2 * i) + ) + bounds = np.array(bounds) + + self._fit = optimize.differential_evolution( + lambda x: self._line_match_distance(x, model, max_distance), + bounds=bounds, + popsize=popsize, + init="sobol", + ) + self.solution.p2w = self._create_model(degree, coeffs=self._fit.x) + + can_match = self._cat_lines is not None and self._obs_lines is not None + if refine_fit: + self.refine_fit(refined_fit_degree, max_match_distance=refine_max_distance) + else: + if can_match: + self.match_lines() + return self.solution + + def refine_fit( + self, degree: None | int = None, max_match_distance: float = 5.0, max_iter: int = 5 + ) -> WavelengthSolution1D: + """Refine the pixel-to-wavelength transformation fit. + + Fits (or re-fits) the polynomial wavelength solution using the currently + matched pixel–wavelength pairs. Optionally adjusts the polynomial degree, + filters matches by a maximum pixel-space separation, and iterates the fit. + + Parameters + ---------- + degree + The polynomial degree for the wavelength solution. If ``None``, the degree + previously set by the `~WavelengthCalibration1D.fit_lines` or + `~WavelengthCalibration1D.fit_dispersion` method will be used. + + max_match_distance + Maximum allowable distance used to identify matched pixel and wavelength + data points. Points exceeding the bound will not be considered in the fit. + + max_iter + Maximum number of fitting iterations. + """ + + # Create a new model with the current parameters if degree is specified. + if degree is not None and degree != self.degree: + model = self._create_model(degree, coeffs=self.solution.p2w[1].parameters) + else: + model = self.solution.p2w + + shift, poly = model + fitter = fitting.LinearLSQFitter() + rms = np.nan + for i in range(max_iter): + self.match_lines(max_match_distance) + matched_pix = np.ma.concatenate(self.observed_line_locations).compressed() + matched_wav = np.ma.concatenate(self.catalog_line_locations).compressed() + rms_new = np.sqrt(((matched_wav - model(matched_pix)) ** 2).mean()) + if isclose(rms_new, rms): + break + model = shift | fitter(poly, matched_pix - self.ref_pixel, matched_wav) + rms = rms_new + self.solution.p2w = model + return self.solution + + @property + def degree(self) -> None | int: + return self._degree + + @degree.setter + def degree(self, degree: int | None): + if degree is not None and degree < 1: + raise ValueError("Degree must be at least 1.") + self._degree = degree + + @property + def observed_lines(self) -> None | list[MaskedArray]: + """Pixel positions and amplitudes of the observed lines as a list of masked arrays.""" + return self._obs_lines + + @cached_property + def observed_line_locations(self) -> None | list[MaskedArray]: + """Pixel positions of the observed lines as a list of masked arrays.""" + if self._obs_lines is None: + return None + else: + return [line[:, 0] for line in self._obs_lines] + + @cached_property + def observed_line_amplitudes(self) -> None | list[MaskedArray]: + """Amplitudes of the observed lines as a list of masked arrays.""" + if self._obs_lines is None: + return None + else: + return [line[:, 1] for line in self._obs_lines] + + @observed_lines.setter + def observed_lines(self, line_lists: ArrayLike | list[ArrayLike]): + if not isinstance(line_lists, Sequence): + line_lists = [line_lists] + self._obs_lines = [] + for lst in line_lists: + self._obs_lines.append(_format_linelist(lst)) + if hasattr(self, "observed_line_locations"): + del self.observed_line_locations + if hasattr(self, "observed_line_amplitudes"): + del self.observed_line_amplitudes + + @property + def catalog_lines(self) -> None | list[MaskedArray]: + """Catalog line wavelengths as a list of masked arrays.""" + return self._cat_lines + + @cached_property + def catalog_line_locations(self) -> None | list[MaskedArray]: + """Pixel positions of the catalog lines as a list of masked arrays.""" + if self._cat_lines is None: + return None + else: + return [line[:, 0] for line in self._cat_lines] + + @cached_property + def catalog_line_amplitudes(self) -> None | list[MaskedArray]: + """Amplitudes of the catalog lines as a list of masked arrays.""" + if self._obs_lines is None: + return None + else: + return [line[:, 1] for line in self._cat_lines] + + @catalog_lines.setter + def catalog_lines(self, line_lists: ArrayLike | list[ArrayLike]): + if not isinstance(line_lists, Sequence): + line_lists = [line_lists] + self._cat_lines = [] + for lst in line_lists: + self._cat_lines.append(_format_linelist(lst)) + if hasattr(self, "catalog_line_locations"): + del self.catalog_line_locations + if hasattr(self, "catalog_line_amplitudes"): + del self.catalog_line_amplitudes + + def match_lines(self, max_distance: float = 5) -> None: + """Match the observed lines to theoretical lines. + + Parameters + ---------- + max_distance + The maximum allowed distance between the catalog and observed lines for them to be + considered a match. + """ + + for iframe, tree in enumerate(self._trees): + l, ix = tree.query( + self.solution.p2w(self.observed_line_locations[iframe].data)[:, None], + distance_upper_bound=max_distance, + ) + m = np.isfinite(l) + + # Check for observed lines that match a catalog line. + # Remove all but the nearest match. This isn't an optimal solution, + # we could also iterate the match by removing the currently matched + # lines, but this works for now. + uix, cnt = np.unique(ix[m], return_counts=True) + if any(n := cnt > 1): + for i, c in zip(uix[n], cnt[n]): + s = ix == i + r = np.zeros(c, dtype=bool) + r[np.argmin(l[s])] = True + m[s] = r + + self._cat_lines[iframe].mask[:, :] = True + self._cat_lines[iframe].mask[ix[m], :] = False + self._obs_lines[iframe].mask[:, :] = ~m[:, None] + + def remove_unmatched_lines(self) -> None: + """Remove unmatched lines from observation and catalog line data.""" + self.observed_lines = [lst.compressed().reshape([-1, 2]) for lst in self._obs_lines] + self.catalog_lines = [lst.compressed().reshape([-1, 2]) for lst in self._cat_lines] + self._create_trees() + + def rms(self, space: Literal["pixel", "wavelength"] = "wavelength") -> float: + """Compute the RMS of the residuals between matched lines in the pixel or wavelength space. + + Parameters + ---------- + space + The space in which to calculate the RMS residual. If 'wavelength', + the calculation is performed in the wavelength space. If 'pixel', + it is performed in the pixel space. Default is 'wavelength'. + + Returns + ------- + float + """ + self.match_lines() + mpix = np.ma.concatenate(self.observed_line_locations).compressed() + mwav = np.ma.concatenate(self.catalog_line_locations).compressed() + if space == "wavelength": + return np.sqrt(((mwav - self.solution.p2w(mpix)) ** 2).mean()) + elif space == "pixel": + return np.sqrt(((mpix - self.solution.w2p(mwav)) ** 2).mean()) + else: + raise ValueError("Space must be either 'pixel' or 'wavelength'") + + def _plot_lines( + self, + kind: Literal["observed", "catalog"], + frames: int | Sequence[int] | None = None, + axes: Axes | Sequence[Axes] | None = None, + figsize: tuple[float, float] | None = None, + plot_labels: bool | Sequence[bool] = True, + map_x: bool = False, + label_kwargs: dict | None = None, + ) -> Figure: + """ + Plot lines with optional features such as wavelength mapping and label customization. + + Parameters + ---------- + kind + Specifies the line list to plot. + + frames + Frame indices to plot. If None, all frames are plotted. + + axes + Axes object(s) where the lines should be plotted. If None, new Axes are generated. + + figsize + Size of the figure to use if creating new Axes. Ignored if axes are provided. + + plot_labels + Flag(s) indicating whether to display labels for the lines. If a single value is + provided, it is applied to all frames. + + map_x + If True, maps the x-axis values between pixel and wavelength space. + + label_kwargs + Additional keyword arguments to customize the label style. + + Returns + ------- + Figure + The Figure object containing the plotted spectral lines. + """ + largs = dict(backgroundcolor="w", rotation=90, size="small") + if label_kwargs is not None: + largs.update(label_kwargs) + + if frames is None: + frames = np.arange(self.nframes) + else: + frames = np.atleast_1d(frames) + + if axes is None: + fig, axes = subplots( + frames.size, 1, figsize=figsize, constrained_layout=True, sharex="all" + ) + elif isinstance(axes, Axes): + fig = axes.figure + axes = [axes] + else: + fig = axes[0].figure + axes = np.atleast_1d(axes) + + if isinstance(plot_labels, bool): + plot_labels = np.full(frames.size, plot_labels, dtype=bool) + + if map_x and self.solution.p2w is None: + raise ValueError("Cannot map between pixels and wavelengths without a fitted model.") + + if kind == "observed": + transform = self.solution.pix_to_wav if map_x else lambda x: x + linelists = self.observed_lines + spectra = self.arc_spectra + lc = "C0" + else: + transform = self.solution.wav_to_pix if map_x else lambda x: x + linelists = self.catalog_lines + spectra = None + lc = "C1" + + ypad = 1.3 + labels = [] + for iframe, (ax, frame) in enumerate(zip(axes, frames)): + if spectra is not None: + spc = self.arc_spectra[iframe] + vmax = np.nanmax(spc.flux.value) + ax.plot(transform(spc.spectral_axis.value), spc.flux.value / vmax, "k") + else: + vmax = 1.0 + + if linelists is not None: + labels.append([]) + + # Loop over individual lines in the line list. + for i in range(linelists[iframe].shape[0]): + c, a = linelists[iframe].data[i] + ls = "-" if linelists[iframe].mask[i, 0] == 0 else ":" + + ax.plot(transform([c, c]), [a / vmax + 0.1, 1.27], c=lc, ls=ls, zorder=-100) + if plot_labels[iframe]: + lloc = transform(c) + labels[-1].append( + ax.text( + lloc, + ypad, + np.round(lloc, 4 - 1 - int(np.floor(np.log10(lloc)))), + ha="center", + va="top", + **largs, + ) + ) + labels[-1][-1].set_clip_on(True) + labels[-1][-1].zorder = a + + if (kind == "observed" and not map_x) or (kind == "catalog" and map_x): + xlabel = "Pixel" + else: + xlabel = f"Wavelength {self._unit_str}" + + if kind == "catalog": + axes[0].xaxis.set_label_position("top") + axes[0].xaxis.tick_top() + setp(axes[0], xlabel=xlabel) + for ax in axes[1:]: + ax.set_xticklabels([]) + else: + setp(axes[-1], xlabel=xlabel) + for ax in axes[:-1]: + ax.set_xticklabels([]) + + xlims = np.array([ax.get_xlim() for ax in axes]) + setp(axes, xlim=(xlims[:, 0].min(), xlims[:, 1].max()), yticks=[]) + + if linelists is not None: + fig.canvas.draw() + for i in range(len(frames)): + if plot_labels[i]: + + # Calculate the label bounding box upper limits and adjust the y-axis limits. + tr_to_data = axes[i].transData.inverted() + ymax = -np.inf + for lb in labels[i]: + ymax = max(ymax, tr_to_data.transform(lb.get_window_extent().p1)[1]) + setp(axes[i], ylim=(-0.04, ymax * 1.06)) + + # Remove the overlapping labels prioritizing the high-amplitude lines. + _unclutter_text_boxes(labels[i]) + + return fig + + def plot_catalog_lines( + self, + frames: int | Sequence[int] | None = None, + axes: Axes | Sequence[Axes] | None = None, + figsize: tuple[float, float] | None = None, + plot_labels: bool | Sequence[bool] = True, + map_to_pix: bool = False, + label_kwargs: dict | None = None, + ) -> Figure: + """Plot the catalog lines. + + Parameters + ---------- + frames + Specifies the frames to be plotted. If an integer, only one frame is plotted. + If a sequence, the specified frames are plotted. If None, default selection + or all frames are plotted. + + axes + The matplotlib axes where catalog data will be plotted. If provided, the function + will plot on these axes. If None, new axes will be created. + + figsize + Specifies the dimensions of the figure as (width, height). If None, the default + dimensions are used. + + plot_labels + If True, the numerical values associated with the catalog data will be displayed + in the plot. If False, only the graphical representation of the lines will be shown. + + map_to_pix + Indicates whether the catalog data should be mapped to pixel coordinates + before plotting. If True, the data is converted to pixel coordinates. + + label_kwargs + Specifies the keyword arguments for the line label text objects. + + Returns + ------- + Figure + The matplotlib figure containing the plotted catalog lines. + + """ + return self._plot_lines( + "catalog", + frames=frames, + axes=axes, + figsize=figsize, + plot_labels=plot_labels, + map_x=map_to_pix, + label_kwargs=label_kwargs, + ) + + def plot_observed_lines( + self, + frames: int | Sequence[int] | None = None, + axes: Axes | Sequence[Axes] | None = None, + figsize: tuple[float, float] | None = None, + plot_labels: bool | Sequence[bool] = True, + map_to_wav: bool = False, + label_kwargs: dict | None = None, + ) -> Figure: + """Plot observed spectral lines for the given arc spectra. + + Parameters + ---------- + frames + Specifies the frame(s) for which the plot is to be generated. If None, all frames + are plotted. When an integer is provided, a single frame is used. For a sequence + of integers, multiple frames are plotted. + + axes + Axes object(s) to plot the spectral lines on. If None, new axes are created. + + figsize + Dimensions of the figure to be created, specified as a tuple (width, height). Ignored + if ``axes`` is provided. + + plot_labels + If True, plots the numerical values of the observed lines at their respective + locations on the graph. + + map_to_wav + Determines whether to map the x-axis values to wavelengths. + + label_kwargs + Specifies the keyword arguments for the line label text objects. + + Returns + ------- + Figure + The matplotlib figure containing the observed lines plot. + """ + + fig = self._plot_lines( + "observed", + frames=frames, + axes=axes, + figsize=figsize, + plot_labels=plot_labels, + map_x=map_to_wav, + label_kwargs=label_kwargs, + ) + + for ax in fig.axes: + ax.autoscale(True, "x", tight=True) + return fig + + def plot_fit( + self, + frames: Sequence[int] | int | None = None, + figsize: tuple[float, float] | None = None, + plot_labels: bool = True, + obs_to_wav: bool = False, + cat_to_pix: bool = False, + label_kwargs: dict | None = None, + ) -> Figure: + """Plot the fitted catalog and observed lines for the specified arc spectra. + + Parameters + ---------- + frames + The indices of the frames to plot. If `None`, all frames from 0 to + ``self.nframes - 1`` are plotted. + + figsize + Defines the width and height of the figure in inches. If `None`, the + default size is used. + + plot_labels + If `True`, print line locations over the plotted lines. Can also be a list with + the same length as ``frames``. + + obs_to_wav + If `True`, transform the x-axis of observed lines to the wavelength domain + using `self._p2w`, if available. + + cat_to_pix + If `True`, transforms catalog data points to pixel values before plotting. + + label_kwargs + Specifies the keyword arguments for the line label text objects. + + Returns + ------- + matplotlib.figure.Figure + The figure object containing the generated subplots. + """ + if frames is None: + frames = np.arange(self.nframes) + else: + frames = np.atleast_1d(frames) + + fig, axs = subplots(2 * frames.size, 1, constrained_layout=True, figsize=figsize) + self.plot_catalog_lines( + frames, + axs[0::2], + plot_labels=plot_labels, + map_to_pix=cat_to_pix, + label_kwargs=label_kwargs, + ) + self.plot_observed_lines( + frames, + axs[1::2], + plot_labels=plot_labels, + map_to_wav=obs_to_wav, + label_kwargs=label_kwargs, + ) + + xlims = np.array([ax.get_xlim() for ax in axs[::2]]) + if obs_to_wav: + setp(axs, xlim=(xlims[:, 0].min(), xlims[:, 1].max())) + else: + setp(axs[::2], xlim=(xlims[:, 0].min(), xlims[:, 1].max())) + + setp(axs[0], yticks=[], xlabel=f"Wavelength [{self._unit_str}]") + for ax in axs[1:-1]: + ax.set_xlabel("") + ax.set_xticklabels("") + + axs[0].xaxis.set_label_position("top") + axs[0].xaxis.tick_top() + return fig + + def plot_residuals( + self, + ax: Axes | None = None, + space: Literal["pixel", "wavelength"] = "wavelength", + figsize: tuple[float, float] | None = None, + ) -> Figure: + """Plot the residuals of pixel-to-wavelength or wavelength-to-pixel transformation. + + Parameters + ---------- + ax + Matplotlib Axes object to plot on. If None, a new figure and axes are created. + + space + The reference space used for plotting residuals. Options are 'pixel' for residuals + in pixel space or 'wavelength' for residuals in wavelength space. + + figsize + The size of the figure in inches, if a new figure is created. + + Returns + ------- + matplotlib.figure.Figure + """ + if ax is None: + fig, ax = subplots(figsize=figsize, constrained_layout=True) + else: + fig = ax.figure + + self.match_lines() + mpix = np.ma.concatenate(self.observed_line_locations).compressed() + mwav = np.ma.concatenate(self.catalog_line_locations).compressed() + + if space == "wavelength": + twav = self.solution.pix_to_wav(mpix) + ax.plot(mwav, mwav - twav, ".") + ax.text( + 0.98, + 0.95, + f"RMS = {np.sqrt(((mwav - twav) ** 2).mean()):4.2f} {self._unit_str}", + transform=ax.transAxes, + ha="right", + va="top", + ) + setp( + ax, + xlabel=f"Wavelength [{self._unit_str}]", + ylabel=f"Residuals [{self._unit_str}]", + ) + elif space == "pixel": + tpix = self.solution.wav_to_pix(mwav) + ax.plot(mpix, mpix - tpix, ".") + ax.text( + 0.98, + 0.95, + f"RMS = {np.sqrt(((mpix - tpix) ** 2).mean()):4.2f} pix", + transform=ax.transAxes, + ha="right", + va="top", + ) + setp(ax, xlabel="Pixel", ylabel="Residuals [pix]") + else: + raise ValueError("Invalid space specified for plotting residuals.") + ax.axhline(0, c="k", lw=1, ls="--") + return fig diff --git a/specreduce/wavesol1d.py b/specreduce/wavesol1d.py new file mode 100644 index 00000000..6bd7eda0 --- /dev/null +++ b/specreduce/wavesol1d.py @@ -0,0 +1,246 @@ +from functools import cached_property +from typing import Callable + +import astropy.units as u +import gwcs +import numpy as np +from astropy.modeling import models, CompoundModel +from astropy.nddata import VarianceUncertainty +from gwcs import coordinate_frames +from numpy.ma import MaskedArray +from numpy.typing import ArrayLike, NDArray +from scipy.interpolate import interp1d + +from specreduce.compat import Spectrum + + +__all__ = ["WavelengthSolution1D"] + + +def _diff_poly1d(m: models.Polynomial1D) -> models.Polynomial1D: + """Compute the derivative of a Polynomial1D model. + + Computes the derivative of a Polynomial1D model and returns a new Polynomial1D + model representing the derivative. The coefficients of the input model are + used to calculate the coefficients of the derivative model. For a Polynomial1D + of degree n, the derivative is a Polynomial1D of degree n-1. + + Parameters + ---------- + m + A Polynomial1D model for which the derivative is to be computed. + + Returns + ------- + A new Polynomial1D model representing the derivative of the input Polynomial1D model. + """ + coeffs = {f"c{i-1}": i * getattr(m, f"c{i}").value for i in range(1, m.degree + 1)} + return models.Polynomial1D(m.degree - 1, **coeffs) + + +class WavelengthSolution1D: + def __init__( + self, + p2w: None | CompoundModel, + bounds_pix: tuple[int, int], + unit: u.Unit, + ) -> None: + """Class defining a one-dimensional wavelength solution. + + This class manages the mapping between pixel positions and wavelength values in a 1D + spectrum, supporting both forward and reverse transformations. It provides methods for + resampling spectra in the pixel-to-wavelength space while conserving flux, and integrates + with GWCS for coordinate transformations. + + Initializes an object with pixel-to-wavelength transformation, pixel bounds, and + measurement unit. Also, converts the unit to its LaTeX string representation. + + Parameters + ---------- + p2w + The pixel-to-wavelength transformation model. If None, no transformation + will be set. + bounds_pix + The lower and upper pixel bounds defining the range of the spectrum. + unit + The wavelength unit. + """ + self.unit = unit + self._unit_str = unit.to_string("latex") + self.bounds_pix: tuple[int, int] = bounds_pix + self.bounds_wav: tuple[float, float] | None = None + self._p2w: None | CompoundModel = None + self.p2w = p2w + + @property + def p2w(self) -> None | CompoundModel: + """Pixel-to-wavelength transformation.""" + return self._p2w + + @p2w.setter + def p2w(self, m: CompoundModel) -> None: + self._p2w = m + self.ref_pixel = m[0].offset.value if m is not None else None + + if "p2w_dldx" in self.__dict__: + del self.p2w_dldx + if "w2p" in self.__dict__: + del self.w2p + if "gwcs" in self.__dict__: + del self.gwcs + + @cached_property + def p2w_dldx(self) -> CompoundModel: + """Partial derivative of the pixel-to-wavelength transformation, (d lambda) / (d pix).""" + return models.Shift(self._p2w.offset_0) | _diff_poly1d(self._p2w[1]) + + @cached_property + def w2p(self) -> Callable: + """Wavelength-to-pixel transformation.""" + p = np.arange(self.bounds_pix[0] - 2, self.bounds_pix[1] + 2) + self.bounds_wav = self.p2w(self.bounds_pix) + return interp1d(self.p2w(p), p, bounds_error=False, fill_value=np.nan) + + def pix_to_wav(self, pix: float | ArrayLike) -> float | NDArray | MaskedArray: + """Map pixel values into wavelength values. + + Parameters + ---------- + pix + The pixel value(s) to be transformed into wavelength value(s). + + Returns + ------- + Transformed wavelength value(s) corresponding to the input pixel value(s). + """ + if isinstance(pix, MaskedArray): + wav = self.p2w(pix.data) + return np.ma.masked_array(wav, mask=pix.mask) + else: + return self.p2w(pix) + + def wav_to_pix(self, wav: float | ArrayLike) -> float | NDArray | MaskedArray: + """Map wavelength values into pixel values. + + Parameters + ---------- + wav + The wavelength value(s) to be converted into pixel value(s). + + Returns + ------- + The corresponding pixel value(s) for the input wavelength(s). + """ + if isinstance(wav, MaskedArray): + pix = self.w2p(wav.data) + return np.ma.masked_array(pix, mask=wav.mask) + else: + return self.w2p(wav) + + @cached_property + def gwcs(self) -> gwcs.wcs.WCS: + """GWCS object defining the mapping between pixel and spectral coordinate frames.""" + pixel_frame = coordinate_frames.CoordinateFrame( + 1, "SPECTRAL", (0,), axes_names=["x"], unit=[u.pix] + ) + spectral_frame = coordinate_frames.SpectralFrame( + axes_names=("wavelength",), unit=[self.unit] + ) + pipeline = [(pixel_frame, self._p2w), (spectral_frame, None)] + return gwcs.wcs.WCS(pipeline) + + def resample( + self, + spectrum: "Spectrum", + nbins: int | None = None, + wlbounds: tuple[float, float] | None = None, + bin_edges: ArrayLike | None = None, + ) -> Spectrum: + """Bin the given pixel-space 1D spectrum to a wavelength space conserving the flux. + + This method bins a pixel-space spectrum to a wavelength space using the computed + pixel-to-wavelength and wavelength-to-pixel transformations and their derivatives with + respect to the spectral axis. The binning is exact and conserves the integrated flux. + + Parameters + ---------- + spectrum + A Spectrum instance containing the flux to be resampled over the wavelength + space. + + nbins + The number of bins for resampling. If not provided, it defaults to the size of the + input spectrum. + + wlbounds + A tuple specifying the starting and ending wavelengths for resampling. If not + provided, the wavelength bounds are inferred from the object's methods and the + entire flux array is used. + + bin_edges + Explicit bin edges in the wavelength space. Should be an 1D array-like [e_0, e_1, + ..., e_n] with n = nbins + 1. The bins are created as [[e_0, e_1], [e_1, e_2], ..., + [e_n-1, n]]. If provided, ``nbins`` and ``wlbounds`` are ignored. + + Returns + ------- + 1D spectrum binned to the specified wavelength bins. + """ + if nbins is not None and nbins < 0: + raise ValueError("Number of bins must be positive.") + + if self._p2w is None: + raise ValueError("Wavelength solution not set.") + + flux = spectrum.flux.value + pixels = spectrum.spectral_axis.value + + if spectrum.uncertainty is not None: + ucty = spectrum.uncertainty.represent_as(VarianceUncertainty).array + ucty_type = type(spectrum.uncertainty) + else: + ucty = np.zeros_like(flux) + ucty_type = VarianceUncertainty + npix = flux.size + nbins = npix if nbins is None else nbins + if wlbounds is None: + l1, l2 = self.p2w(pixels[[0, -1]] + np.array([-0.5, 0.5])) + else: + l1, l2 = wlbounds + + if bin_edges is not None: + bin_edges_wav = np.asarray(bin_edges) + nbins = bin_edges_wav.size - 1 + else: + bin_edges_wav = np.linspace(l1, l2, num=nbins + 1) + + bin_edges_pix = np.clip(self.w2p(bin_edges_wav) + 0.5, 0, npix - 1e-12) + bin_edge_ix = np.floor(bin_edges_pix).astype(int) + bin_edge_w = bin_edges_pix - bin_edge_ix + bin_centers_wav = 0.5 * (bin_edges_wav[:-1] + bin_edges_wav[1:]) + flux_wl = np.zeros(nbins) + ucty_wl = np.zeros(nbins) + weights = np.zeros(npix) + + dldx = np.diff(self.p2w(np.arange(pixels[0], pixels[-1] + 2) - 0.5)) + + for i in range(nbins): + i1, i2 = bin_edge_ix[i : i + 2] + weights[:] = 0.0 + if i1 != i2: + weights[i1 + 1 : i2] = 1.0 + weights[i1] = 1 - bin_edge_w[i] + weights[i2] = bin_edge_w[i + 1] + sl = slice(i1, i2 + 1) + w = weights[sl] + flux_wl[i] = (w * flux[sl] * dldx[sl]).sum() + ucty_wl[i] = (w**2 * ucty[sl] * dldx[sl]).sum() + else: + fracw = bin_edges_pix[i + 1] - bin_edges_pix[i] + flux_wl[i] = fracw * flux[i1] * dldx[i1] + ucty_wl[i] = fracw**2 * ucty[i1] * dldx[i1] + + bin_widths_wav = np.diff(bin_edges_wav) + flux_wl = flux_wl / bin_widths_wav * spectrum.flux.unit / self.unit + ucty_wl = VarianceUncertainty(ucty_wl / bin_widths_wav**2).represent_as(ucty_type) + return Spectrum(flux_wl, bin_centers_wav * self.unit, uncertainty=ucty_wl)