diff --git a/.gitignore b/.gitignore index 5219e90..67e14cf 100644 --- a/.gitignore +++ b/.gitignore @@ -143,6 +143,7 @@ media static # Specifc Template +docs/sg_execution_times.rst docs/_build docs/generated docs/api diff --git a/changelog/318.doc.rst b/changelog/318.doc.rst new file mode 100644 index 0000000..5a1740f --- /dev/null +++ b/changelog/318.doc.rst @@ -0,0 +1 @@ +Transformed the documentation layout. diff --git a/docs/code_ref/index.rst b/docs/code_ref/index.rst index e41a898..f1fb0af 100644 --- a/docs/code_ref/index.rst +++ b/docs/code_ref/index.rst @@ -1,3 +1,5 @@ +.. _aia_api_reference: + ============= API reference ============= diff --git a/docs/conf.py b/docs/conf.py index a0484cc..a224127 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -51,7 +51,7 @@ source_suffix = ".rst" master_doc = "index" default_role = "obj" -ogp_image = "https://gitlab.com/LMSAL_HUB/aia_hub/aiapy/-/raw/main/docs/_static/sdo.png" +ogp_image = "https://github.com/LM-SAL/aiapy/blob/main/docs/_static/sdo.png" ogp_use_first_image = True ogp_description_length = 160 ogp_custom_meta_tags = [ diff --git a/docs/develop.rst b/docs/develop.rst index 32e8fac..02feb4e 100644 --- a/docs/develop.rst +++ b/docs/develop.rst @@ -1,10 +1,10 @@ -.. _dev-guide: +.. _aiapy_dev-guide: ============ Contributing ============ -Contributing to open source projects is a fantastic way to advance one's coding skills; it's trying to create something, making mistakes, and learning from those mistakes. +Contributing to open source projects is a good way to advance one's coding skills; it's trying to create something, making mistakes, and learning from those mistakes. That's how we all improve, and we are happy to help others learn. Being an open source contributor doesn't just mean writing code, either. @@ -14,126 +14,26 @@ Some of these contributions may be the most valuable to the project as a whole, Issue Tracking -------------- -All bugs, feature requests, and other issues related to ``aiapy`` should be recorded using the GitLab issue tracker. -You can find instructions for how to create an issue `here `__. +All bugs, feature requests, and other issues related to ``aiapy`` should be recorded using the GitHub issue tracker. +You can find instructions for how to create an issue `here `__. All conversation regarding these issues should take place on the issue tracker. When a merge request resolves an issue, the issue will be closed and the appropriate merge request will be referenced. Issues will not be closed without a reason given. -Creating a fork ---------------- +Code +---- If you would like to contribute to ``aiapy``, you will first need to setup your development environment. -First create a fork of the main ``aiapy`` repository under your GitLab username. -You can find the instructions for how to do this `here `__. -If you don't already have an account on GitLab, you'll need to create one. -You can also sign into GitLab using your GitHub username. -Next, clone your fork of ``aiapy`` to your local machine, +We suggest reading through the `SunPy developer's guide`_ for a more detailed description of the development process, as we follow there process for ``aiapy``. -.. code:: shell - - git clone https://gitlab.com//``aiapy``.git - -Now add the main ``aiapy`` repository as an upstream repository, - -.. code:: shell - - git remote add upstream https://gitlab.com/LMSAL_HUB/aia_hub/``aiapy``.git - -You can now keep your fork up to date with main repository by running, - -.. code:: shell - - git pull upstream main - -Installation -------------- - -If you're using the `Miniconda Python distribution `__ (recommended), -create a new environment for ``aiapy`` development, - -.. code-block:: shell - - conda create --name ``aiapy``-dev pip - conda activate ``aiapy``-dev - -If you're using an alternate python installation, you can also use `virtual environments `__. -Next, install the needed dependencies, - -.. code-block:: shell - - cd aiapy - pip install -e .[test,docs] - -This includes all of the dependencies for the package as well as ``pytest`` for running tests and ``sphinx`` for building the docs. - -To make sure everything is working alright, you can run the tests, - -.. code-block:: shell - - pytest --remote-data=any - -See :ref:`tests` for more details regarding running the tests. - -Making a contribution ---------------------- - -If you want to add a feature or bugfix to ``aiapy``, start by first making sure the main branch of your fork is up to date with the main branch of the main repository (see above, this will help to prevent potential file conflicts). -Next, create a new branch and switch to it, - -.. code:: shell - - git checkout -b my-new-feature - -After you've made your changes, commit and push them up to GitLab, - -.. code:: shell - - git add changed_file_1.py changed_file_2.py - git commit -m "short description of my change" - git push origin my-new-feature - -Once you see the changes in GitLab, create a merge request against the main ``aiapy`` repository. -You can find instructions for how to do this `here `__. -Others will likely have comments and suggestions regarding your proposed changes. -You can make these changes using the instructions listed above. - -At least one other ``aiapy`` developer must approve your changes before the code can be merged. -Additionally, all automated tests should pass and all conversations should be resolved. -Once these steps are complete, the code can be merged and you can delete your branch ``my-new-feature``. - -.. _tests: +This will hopefully make it easier for you to contribute to ``aiapy`` and other SunPy affiliated packages in the future. +If you encounter any problems, please don't hesitate to ask for help on any of our communication channels (found on the landing page of the documentation). Testing ------- -Before committing any changes, you should ensure that the all of the tests pass locally. -To run the tests, - -.. code:: shell - - pytest --remote-data=any - -This will generate report showing which tests passed and which failed (if any). -Dropping the ``--remote-data`` flag will skip tests that require a network connection. -``aiapy`` uses the `pytest `__ framework for discovering and running all of the tests. - -Additions to the codebase should be accompanied by appropriate tests such that the test coverage of the entire package does not decrease. -You can check the test coverage by running, - -.. code:: shell - - pytest --remote-data=any --cov aiapy - -Additionally, the test suite, including the documentation build and code style checks can be run with `tox `__. -See the `SunPy developer's guide`_ for more information on running the test suite with ``tox``. - -Tests should be added to the directory in the appropriate subpackage, e.g. for ``calibrate``, the tests should be placed in ``calibrate/tests``. -Your tests can be added to an existing file or placed in a new file following the naming convention ``test_*.py``. -This organization allows the tests to be automatically discovered by pytest. - There are several tests that require a working installation of `sswidl `__ in order to compare results from IDL and Python. This is managed via the `hissw `__ package. If you'd like to run these tests, you must first tell ``hissw`` where to find your IDL and SSW installations by placing the following lines in the file: ``$HOME/.hissw/hisswrc``, @@ -148,26 +48,11 @@ where ``ssw_home`` is the path to the top of the sswidl tree and ``idl_home`` is For more details, see the `hissw documentation `__. If a working installation is not available, these tests are automatically skipped. -Documentation -------------- - -All documentation is written in `reStructuredText `__ and rendered using `Sphinx `__. -Documentation strings are automatically pulled from all modules, functions and classes to create the API documentation. -You can build and test the documentation locally by running, - -.. code:: shell - - cd docs - make html - -This will run Sphinx on the restructured text files in order to create the HTML version of the documentation. -The built documentation, in HTML format, is in ``docs/_build/html``. - -Best practices --------------- +Standards +--------- All contributors to the ``aiapy`` codebase should follow the `SunPy developer's guide`_. This guide lays out a set of best practices for contributing, reviewing, testing, and documenting code. All contributions to ``aiapy`` must adhere to the `Python in Heliophysics Community Standards `__. -.. _`SunPy developer's guide`: https://docs.sunpy.org/en/latest/dev_guide/index.html +.. _`SunPy developer's guide`: https://docs.sunpy.org/en/latest/dev_guide/contents/newcomers.html diff --git a/docs/getting_started.rst b/docs/getting_started.rst deleted file mode 100644 index 12dcec1..0000000 --- a/docs/getting_started.rst +++ /dev/null @@ -1,24 +0,0 @@ -Getting Started -=============== - -The stable version of ``aiapy`` is available through both PyPI via ``pip``, - -.. code-block:: shell - - pip install aiapy - -as well as through the `Miniconda distribution `__ via ``conda-forge``, - -.. code-block:: shell - - conda install -c conda-forge aiapy - -You can also install the development version from GitLab, - -.. code-block:: shell - - git clone https://gitlab.com/LMSAL_HUB/aia_hub/aiapy.git - cd aiapy - pip install -e ".[dev]" - -If you will be developing ``aiapy``, please see the :ref:`dev-guide`. diff --git a/docs/index.rst b/docs/index.rst index 86e33a7..fac3e37 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -1,36 +1,76 @@ +.. _aia_docs_index: + +=================== aiapy documentation =================== -``aiapy`` is a Python package for analyzing data from the Atmospheric Imaging Assembly (AIA) instrument onboard the Solar Dynamics Observatory (SDO) spacecraft. +``aiapy`` is a Python library for analyzing data from the Atmospheric Imaging Assembly (AIA) instrument onboard the Solar Dynamics Observatory (SDO) spacecraft. -``aiapy`` includes software for converting AIA images from level 1 to level 1.5, point spread function deconvolution, and computing the wavelength and temperature response functions for the EUV channels. +It allows for converting AIA images from level 1 to level 1.5, point spread function deconvolution, and computing the wavelength and temperature response functions for the EUV channels. .. grid:: 1 2 2 2 - :gutter: 3 + :gutter: 2 + + .. grid-item-card:: Topic Guide + :link: aiapy-topic-guide-index + :link-type: ref + :text-align: center + + :material-outlined:`accessibility_new;8em;sd-text-secondary` + + **New users start here!** + Topic Guides on how to install and use the key features of aiapy. + + .. grid-item-card:: Example gallery + :link: generated/gallery + :text-align: center + + :material-outlined:`palette;8em;sd-text-secondary` + + Examples including plots on accomplishing common tasks using aiapy. + + .. grid-item-card:: Reference + :link: aia_api_reference + :link-type: ref + :text-align: center + + :material-outlined:`code;8em;sd-text-secondary` + + Technical description of the inputs, outputs, and behavior of each component of aiapy. + + .. grid-item-card:: Get Help + :text-align: center + + :material-outlined:`live_help;8em;sd-text-secondary` - .. grid-item-card:: - :class-card: card + .. button-link:: https://app.element.io/#/room/#sunpy:openastronomy.org + :shadow: + :expand: + :color: warning - Getting started - ^^^^^^^^^^^^^^^ + **Join the chat** - .. toctree:: - :maxdepth: 1 + .. button-link:: https://github.com/LM-SAL/aiapy/issues + :shadow: + :expand: + :color: warning - getting_started - preparing_data - generated/gallery/index - code_ref/index + **Report an issue** - .. grid-item-card:: - :class-card: card + .. button-link:: https://community.openastronomy.org/c/sunpy/5 + :shadow: + :expand: + :color: warning - Other info - ^^^^^^^^^^ + **Post on Discourse** - .. toctree:: - :maxdepth: 1 +.. toctree:: + :maxdepth: 1 + :hidden: - citation - changelog - develop + tutorial/index + generated/gallery/index + code_ref/index + citation + changelog + develop diff --git a/docs/topic_guide/index.rst b/docs/topic_guide/index.rst new file mode 100644 index 0000000..bd2fe88 --- /dev/null +++ b/docs/topic_guide/index.rst @@ -0,0 +1,14 @@ +.. _aiapy-topic-guide-index: + +***************** +aiapy topic guide +***************** + +Welcome to the topic guide for the ``aiapy`` library. +``aiapy`` is a community-developed, free and open-source library to to convert AIA images from level 1 to level 1.5, compute point spread function deconvolution, and computing the wavelength and temperature response functions for the EUV channels. + +.. toctree:: + :maxdepth: 2 + + installation + preparing_data diff --git a/docs/topic_guide/installation.rst b/docs/topic_guide/installation.rst new file mode 100644 index 0000000..e60bc5e --- /dev/null +++ b/docs/topic_guide/installation.rst @@ -0,0 +1,102 @@ +.. _aiapy-topic-guide-installing: + +************ +Installation +************ + +This topic guide details how to get a working installation of Python and ``aiapy``. + +Installing Python +================= + +There are many ways to install Python, but even if you have Python installed somewhere on your computer we recommend following these instructions anyway. +That's because we will create a new Python environment. +As well as containing a Python installation, this environment provides an isolated place to install Python packages (like ``aiapy``) without affecting any other current Python installation. +If you already have Python and either ``conda`` or ``pip`` working you can skip the next section. + +Installing miniforge +-------------------- + +If you don't already have a Python installation then we recommend installing Python with `miniforge `__. +This will install ``conda`` and automatically configure the default channel (a channel is a remote software repository) to be ``conda-forge``, which is where ``aiapy`` is available. + +First, download the installer for your system and architecture from the links below: + +.. grid:: 3 + + .. grid-item-card:: Linux + + `x86-64 `__ + + `aarch64 `__ + + `ppc64le `__ + + .. grid-item-card:: Windows + :link: https://github.com/conda-forge/miniforge/releases/latest/download/Miniforge3-Windows-x86_64.exe + + `x86-64 `__ + + .. grid-item-card:: Mac + + `x86-64 `__ + + `arm64 (Apple + Silicon) `__ + +Then select your platform to install miniforge: + +.. tab-set:: + + .. tab-item:: Linux & Mac + :sync: platform + + Linux & Mac Run the script downloaded above, with + ``bash ``. The following should work: + + .. code-block:: console + + bash Miniforge3-$(uname)-$(uname -m).sh + + Once the installer has completed, close and reopen your terminal. + + .. tab-item:: Windows + :sync: platform + + Double click the executable file downloaded from + the links above. + + Once the installer has completed you should have a new "miniforge + Prompt" entry in your start menu. + +In a new terminal (miniforge Prompt on Windows) run ``conda list`` to test that the install has worked. + +Installing aiapy +---------------- + +To install ``aiapy``, start by launching a terminal (under a UNIX-like system) or the miniforge Prompt (under Windows). +Now we will create and activate new virtual environment to install ``aiapy`` into: + +.. code-block:: bash + + $ conda create --name aiapy + # Only run the following two lines + # if you have NOT installed miniforge or added conda-forge to your channels + $ conda config --add channels conda-forge + $ conda config --set channel_priority strict + $ conda activate aiapy + +In this case the environment is named 'aiapy'. +Feel free to change this to a different environment name. + +The benefit of using a virtual environment is that it allows you to install packages without affecting any other Python installation on your system. +This also means you can work on multiple projects (research or coding) with different package requirements without them interfering with each other. + +Now we have a fresh environment we can install ``aiapy``: + +.. code-block:: bash + + $ conda install aiapy + +This will install ``aiapy`` and all of its dependencies. +If you want to install another package later, you can run ``conda install ``. diff --git a/docs/preparing_data.rst b/docs/topic_guide/preparing_data.rst similarity index 93% rename from docs/preparing_data.rst rename to docs/topic_guide/preparing_data.rst index 8b06648..2f523ea 100644 --- a/docs/preparing_data.rst +++ b/docs/topic_guide/preparing_data.rst @@ -1,5 +1,8 @@ -Preparing AIA data -================== +.. _aiapy-topic-guide-prepping-level-1: + +============================================ +Preparing AIA data from level 1 to level 1.5 +============================================ AIA data provided by the JSOC are level 1 data products. This means that the images still include the roll angle of the satellite and diff --git a/examples/calculate_response_function.py b/examples/calculate_response_function.py index f69da19..bcc6b64 100644 --- a/examples/calculate_response_function.py +++ b/examples/calculate_response_function.py @@ -14,7 +14,7 @@ from aiapy.response import Channel -################################################## +############################################################################### # Since AIA uses narrow-band filters, other wavelengths (outside of the nominal # wavelength attributed to each filter) contribute to the image data. # Computing these response functions allow us to see which other wavelengths @@ -28,9 +28,10 @@ # this as the most recent instrument data file will # need to be downloaded from a remote server. Subsequent # calls will know that the data has been downloaded. -c = Channel(335 * u.angstrom) -################################################## +aia_335_channel = Channel(335 * u.angstrom) + +############################################################################### # From `Boerner et al. (2012) `_, # the wavelength response function is given by, # @@ -55,8 +56,8 @@ # Reflectance fig = plt.figure() ax = fig.add_subplot(221) -ax.plot(c.wavelength, c.primary_reflectance, label=r"$R_P$") -ax.plot(c.wavelength, c.secondary_reflectance, label=r"$R_S$") +ax.plot(aia_335_channel.wavelength, aia_335_channel.primary_reflectance, label=r"$R_P$") +ax.plot(aia_335_channel.wavelength, aia_335_channel.secondary_reflectance, label=r"$R_S$") ax.set_ylabel(r"Reflectance") ax.set_xlim(50, 400) ax.set_xlabel(r"$\lambda$ [Å]") @@ -64,8 +65,8 @@ # Transmittance ax = fig.add_subplot(222) -ax.plot(c.wavelength, c.entrance_filter_efficiency, label=r"$T_E$") -ax.plot(c.wavelength, c.focal_plane_filter_efficiency, label=r"$T_F$") +ax.plot(aia_335_channel.wavelength, aia_335_channel.entrance_filter_efficiency, label=r"$T_E$") +ax.plot(aia_335_channel.wavelength, aia_335_channel.focal_plane_filter_efficiency, label=r"$T_F$") ax.set_ylabel(r"Transmittance") ax.set_xlim(50, 400) ax.set_xlabel(r"$\lambda$ [Å]") @@ -73,62 +74,66 @@ # Contamination ax = fig.add_subplot(223) -ax.plot(c.wavelength, c.contamination) +ax.plot(aia_335_channel.wavelength, aia_335_channel.contamination) ax.set_ylabel(r"Contamination, $D(\lambda)$") ax.set_xlim(50, 400) ax.set_xlabel(r"$\lambda$ [Å]") # Quantumn efficiency ax = fig.add_subplot(224) -ax.plot(c.wavelength, c.quantum_efficiency) +ax.plot(aia_335_channel.wavelength, aia_335_channel.quantum_efficiency) ax.set_ylabel(r"Quantum Efficiency, $Q(\lambda)$") ax.set_xlim(50, 800) ax.set_xlabel(r"$\lambda$ [Å]") + plt.tight_layout() -plt.show() -################################################## +############################################################################### # Additionally, `aiapy.response.Channel` provides a method for calculating # the wavelength response function using the equation above, -r = c.wavelength_response() -print(r) -################################################## +wavelength_response_335 = aia_335_channel.wavelength_response() +print(wavelength_response_335) + +############################################################################### # We can then plot the response as a function of # wavelength. + fig = plt.figure() + ax = fig.gca() -ax.plot(c.wavelength, r) -ax.set_xlim((c.channel + [-10, 10] * u.angstrom).value) +ax.plot(aia_335_channel.wavelength, wavelength_response_335) +ax.set_xlim((aia_335_channel.channel + [-10, 10] * u.angstrom).value) ax.set_ylim(0, 0.03) ax.set_xlabel(r"$\lambda$ [Å]") -ax.set_ylabel(f'$R(\\lambda)$ [{r.unit.to_string("latex")}]') -plt.show() +ax.set_ylabel(f'$R(\\lambda)$ [{wavelength_response_335.unit.to_string("latex")}]') -################################################## +############################################################################### # On telescopes 1, 3, and 4, both channels are always illuminated. # This can lead to "crosstalk" contamination in a channel from the channel with # which it shares a telescope. This impacts the 94 Å and 304 Å channels # as well as the 131 Å and 335 Å channels. This effect is included # by default in the wavelength response calculation. To exclude this # effect, -r_no_cross = c.wavelength_response(include_crosstalk=False) -################################################## +wavelength_response_335_no_cross = aia_335_channel.wavelength_response(include_crosstalk=False) + +############################################################################### # If we look at the response around 131 Å (the channel with which 335 Å shares # a telescope), we can see the effect that the channel crosstalk has on the # 335 Å response function. + fig = plt.figure() + ax = fig.gca() -ax.plot(c.wavelength, r, label="crosstalk") -ax.plot(c.wavelength, r_no_cross, label="no crosstalk") +ax.plot(aia_335_channel.wavelength, wavelength_response_335, label="crosstalk") +ax.plot(aia_335_channel.wavelength, wavelength_response_335_no_cross, label="no crosstalk") ax.set_xlim(50, 350) ax.set_xlabel(r"$\lambda$ [Å]") -ax.set_ylabel(f'$R(\\lambda)$ [{r.unit.to_string("latex")}]') +ax.set_ylabel(f'$R(\\lambda)$ [{wavelength_response_335.unit.to_string("latex")}]') ax.legend(loc=1, frameon=False) -plt.show() -################################################### +############################################################################### # We can also incorporate various corrections to the # response functions, including a time-dependent # degradation correction as well as a correction based @@ -136,21 +141,24 @@ # time-dependent correction. As an example, to apply the # two aforementioned corrections given the degradation as # of 1 January 2019, + obstime = astropy.time.Time("2019-01-01T00:00:00") -r_time = c.wavelength_response(obstime=obstime) -r_eve = c.wavelength_response(obstime=obstime, include_eve_correction=True) +wavelength_response_335_time = aia_335_channel.wavelength_response(obstime=obstime) +wavelength_response_335_eve = aia_335_channel.wavelength_response(obstime=obstime, include_eve_correction=True) -#################################################### +############################################################################### # We can then compare the two corrected response # functions to the uncorrected case. + fig = plt.figure() ax = fig.gca() -ax.plot(c.wavelength, r, label="uncorrected") -ax.plot(c.wavelength, r_time, label="degradation correction") -ax.plot(c.wavelength, r_eve, label="EVE correction") -ax.set_xlim((c.channel + [-20, 20] * u.angstrom).value) +ax.plot(aia_335_channel.wavelength, wavelength_response_335, label="uncorrected") +ax.plot(aia_335_channel.wavelength, wavelength_response_335_time, label="degradation correction") +ax.plot(aia_335_channel.wavelength, wavelength_response_335_eve, label="EVE correction") +ax.set_xlim((aia_335_channel.channel + [-20, 20] * u.angstrom).value) ax.set_ylim(0, 0.03) ax.set_xlabel(r"$\lambda$ [Å]") -ax.set_ylabel(f'$R(\\lambda)$ [{r.unit.to_string("latex")}]') +ax.set_ylabel(f'$R(\\lambda)$ [{wavelength_response_335.unit.to_string("latex")}]') ax.legend(loc=2, frameon=False) + plt.show() diff --git a/examples/download_specific_data.py b/examples/download_specific_data.py index d1f7a77..a6fcea7 100644 --- a/examples/download_specific_data.py +++ b/examples/download_specific_data.py @@ -19,7 +19,7 @@ from aiapy.calibrate import correct_degradation, register, update_pointing from aiapy.calibrate.util import get_correction_table, get_pointing_table -##################################################### +############################################################################### # Exporting data from the JSOC requires registering your # email first. Please replace the text after the ``=`` # with your email address once you have registered. @@ -27,7 +27,7 @@ jsoc_email = os.environ.get("JSOC_EMAIL") -##################################################### +############################################################################### # Our goal is to request data of a recent X-class flare. # The X-class flare occurred on the 2021/07/03 at 14:30:00 UTC. # We will focus on the 5 minutes before and after this time. @@ -45,7 +45,7 @@ print(query) -##################################################### +############################################################################### # Now we will download the data and "aia prep" the # data with every feature of `aiapy` and plot the # data sequence using `sunpy`. @@ -70,8 +70,8 @@ map_cropped = map_normalized.submap(bottom_left, top_right=top_right) level_15_maps.append(map_cropped) -##################################################### -# Finally, we create a sequence of maps and animate it +############################################################################### +# Finally, we create a sequence of maps and animate it: sequence = sunpy.map.Map(level_15_maps, sequence=True) diff --git a/examples/instrument_degradation.py b/examples/instrument_degradation.py index d7d5ea6..08b56f0 100644 --- a/examples/instrument_degradation.py +++ b/examples/instrument_degradation.py @@ -17,7 +17,7 @@ from aiapy.calibrate import degradation from aiapy.calibrate.util import get_correction_table -########################################################### +############################################################################### # The sensitivity of the AIA channels degrade over time. Possible causes include # the deposition of organic molecules from the telescope structure onto the # optical elements and the decrease in detector sensitivity following (E)UV @@ -37,37 +37,46 @@ # First, fetch this correction table. It is not strictly necessary to do this explicitly, # but will significantly speed up the calculation by only fetching the table # once. + correction_table = get_correction_table() -########################################################### +############################################################################### # We want to compute the degradation for each EUV channel. -channels = [94, 131, 171, 193, 211, 304, 335] * u.angstrom -########################################################### -# We can use the `~astropy.time` subpackage to create an array of times +aia_channels = [94, 131, 171, 193, 211, 304, 335] * u.angstrom + +############################################################################### +# We can use `~astropy.time.Time` to create an array of times # between now and the start of the mission with a cadence of one week. -time_0 = astropy.time.Time("2010-03-25T00:00:00", scale="utc") + +start_time = astropy.time.Time("2010-03-25T00:00:00", scale="utc") now = astropy.time.Time.now() -time = time_0 + np.arange(0, (now - time_0).to(u.day).value, 7) * u.day +time_range = start_time + np.arange(0, (now - start_time).to(u.day).value, 7) * u.day -########################################################### +############################################################################### # Finally, we can use the `aiapy.calibrate.degradation` function to # compute the degradation for a particular channel and observation time. # This is modeled as the ratio of the effective area measured at a particular # calibration epoch over the uncorrected effective area with a polynomial # interpolation to the exact time. -deg = {c: degradation(c, time, correction_table=correction_table) for c in channels} +deg = {channel: degradation(channel, time_range, correction_table=correction_table) for channel in aia_channels} -########################################################### +############################################################################### # Plotting the different degradation curves as a function of time, we can # easily visualize how the different channels have degraded over time. -time_support(format="jyear") # This lets you pass astropy.time.Time objects directly to matplotlib + +# This lets you pass astropy.time.Time objects directly to matplotlib +time_support(format="jyear") + fig = plt.figure() ax = fig.gca() -for c in channels: - ax.plot(time, deg[c], label=f"{c.value:.0f} Å") -ax.set_xlim(time[[0, -1]]) + +for channel in aia_channels: + ax.plot(time_range, deg[channel], label=f"{channel.value:.0f} Å") + +ax.set_xlim(time_range[[0, -1]]) ax.legend(frameon=False, ncol=4, bbox_to_anchor=(0.5, 1), loc="lower center") ax.set_xlabel("Time") ax.set_ylabel("Degradation") + plt.show() diff --git a/examples/prepping_level_1_data.py b/examples/prepping_level_1_data.py index 602f435..f03d4e7 100644 --- a/examples/prepping_level_1_data.py +++ b/examples/prepping_level_1_data.py @@ -9,21 +9,22 @@ In this example, we will demonstrate how to do this with `aiapy`. This corresponds to the `aia_prep.pro` procedure as described in the `SDO Analysis Guide `__. """ - +import matplotlib.pyplot as plt import sunpy.map import aiapy.data.sample as sample_data from aiapy.calibrate import register, update_pointing -########################################################### +############################################################################### # Performing multi-wavelength analysis on level 1 data can be problematic as # each of the AIA channels have slightly different spatial scales and roll # angles. Furthermore, the estimates of the pointing keywords (``CDELT1``, ``CDELT2``, ``CRPIX1``, -# ``CRPIX2``, ``CROTA2``) may have been improved due to limb fitting procedures. The -# `Joint Science Operations Center (JSOC) `_ stores +# ``CRPIX2``, ``CROTA2``) may have been improved due to limb fitting procedures after +# the level 1 file has been created. +# The `Joint Science Operations Center (JSOC) `_ stores # AIA image data and metadata separately; when users download AIA data, these # two data types are combined to produce a FITS file. While metadata are -# continuously updated at JSOC, previously downloaded FITS files will not +# continuously updated at the JSOC, previously downloaded FITS files will not # contain the most recent information. # # Thus, before performing any multi-wavelength analyses, level 1 data @@ -33,42 +34,51 @@ # # First, let's read a level 1 94 Å AIA image from the ``aiapy`` sample data into # a `~sunpy.map.Map` object. -m = sunpy.map.Map(sample_data.AIA_094_IMAGE) -########################################################### +aia_map = sunpy.map.Map(sample_data.AIA_094_IMAGE) + +############################################################################### # The first step in this process is to update the metadata of the map to the # most recent pointing using the `aiapy.calibrate.update_pointing` function. # This function queries the JSOC for the most recent pointing information, # updates the metadata, and returns a `sunpy.map.Map` with updated metadata. -m_updated_pointing = update_pointing(m) -########################################################### +aia_map_updated_pointing = update_pointing(aia_map) + +############################################################################### # If we take a look at the plate scale and rotation matrix of the map, we # find that the scale is slightly off from the expected value of :math:`0.6''` per # pixel and that the rotation matrix has off-diagonal entries. -print(m_updated_pointing.scale) -print(m_updated_pointing.rotation_matrix) -########################################################### +print(aia_map_updated_pointing.scale) +print(aia_map_updated_pointing.rotation_matrix) + +############################################################################### # We can use the `aiapy.calibrate.register` function to scale the image to # the :math:`0.6''` per pixel and derotate the image such that the y-axis is aligned # with solar North. -m_registered = register(m_updated_pointing) -########################################################### +aia_map_registered = register(aia_map_updated_pointing) + +############################################################################### # If we look again at the plate scale and rotation matrix, we # should find that the plate scale in each direction is :math:`0.6''` # per pixel and that the rotation matrix is diagonalized. # The image in `m_registered` is now a level 1.5 data product. -print(m_registered.scale) -print(m_registered.rotation_matrix) +print(aia_map_registered.scale) +print(aia_map_registered.rotation_matrix) -########################################################### +############################################################################### # Finally, we can plot the exposure-normalized map. # Note that small negative pixel values are possible because # CCD images were taken with a pedestal set at ~100 DN. # This pedestal is then subtracted when the JSOC pipeline -# performs dark (+pedestal) subtraction and flatfielding +# performs dark (+ pedestal) subtraction and flat fielding # to generate level 1 files. -m_registered.peek(vmin=0) + +fig = plt.figure() +ax = fig.add_subplot(projection=aia_map_registered) +aia_map_registered.plot(axes=ax) + +plt.show() diff --git a/examples/replace_hot_pixels.py b/examples/replace_hot_pixels.py index c61ea9a..c64045e 100644 --- a/examples/replace_hot_pixels.py +++ b/examples/replace_hot_pixels.py @@ -14,7 +14,7 @@ import aiapy.data.sample as sample_data from aiapy.calibrate import fetch_spikes, respike -#################################################### +############################################################################### # AIA level 1 images have been corrected for hot-pixels (commonly referred to # as "spikes") using an automated correction algorithm which detects them, # removes them, and replaces the "holes" left in the image via interpolation. @@ -33,12 +33,13 @@ # `Joint Science Operations Center (JSOC) `_ as a # separate segment of the `aia.lev1_euv_12s` and `aia.lev1_uv_24s` data series -#################################################### +############################################################################### # First, let's read a level 1 193 Å AIA image from the aiapy sample data # into a `~sunpy.map.Map` object. -m = sunpy.map.Map(sample_data.AIA_193_IMAGE) -########################################################### +aia_map = sunpy.map.Map(sample_data.AIA_193_IMAGE) + +############################################################################### # The spike data are stored as separate data segments in JSOC # as a :math:`3\times N` arrays, where :math:`N` is the number of spikes # removed and the three dimensions correspond to the the 1-D pixel index @@ -52,9 +53,10 @@ # spikes to the 2D pixel full-disk pixel coordinate system given a # `~sunpy.map.Map` representing a level 1 AIA image. # -positions, values = fetch_spikes(m) -########################################################### +positions, values = fetch_spikes(aia_map) + +############################################################################### # Now we are ready to respike the level 1 AIA image. The # `aiapy.calibrate.respike` function performs the respike operation on the given # input image and returns a `~sunpy.map.Map` with the respiked image. This @@ -63,47 +65,52 @@ # # Note that explicitly specifying the spike positions and values is optional. # If they are not given, they are automatically queried from the JSOC. -m_respiked = respike(m, spikes=(positions, values)) -########################################################### +aia_map_respiked = respike(aia_map, spikes=(positions, values)) + +############################################################################### # Now let's create a cutouts of the original level 1 and "re-spiked" (i.e. # level 0.5) images for a region with hot pixels. -top_right = SkyCoord(30 * u.arcsec, 420 * u.arcsec, frame=m.coordinate_frame) -bottom_left = SkyCoord(-120 * u.arcsec, 280 * u.arcsec, frame=m.coordinate_frame) -m_cutout = m.submap(bottom_left, top_right=top_right) -m_respiked_cutout = m_respiked.submap(bottom_left, top_right=top_right) -########################################################### +top_right = SkyCoord(30 * u.arcsec, 420 * u.arcsec, frame=aia_map.coordinate_frame) +bottom_left = SkyCoord(-120 * u.arcsec, 280 * u.arcsec, frame=aia_map.coordinate_frame) +aia_map_cutout = aia_map.submap(bottom_left, top_right=top_right) +aia_map_respiked_cutout = aia_map_respiked.submap(bottom_left, top_right=top_right) + +############################################################################### # Note that we can also retrieve the positions of the spikes # as `~astropy.coordinates.SkyCoord` objects in the projected coordinate # system of the image using the `as_coords=True` keyword argument. This # gives us only those spikes in the field of view of the cutout. -spike_coords, _ = fetch_spikes(m_cutout, as_coords=True) -########################################################### +spike_coords, _ = fetch_spikes(aia_map_cutout, as_coords=True) + +############################################################################### # Finally, let's plot the two cutouts for comparison and plot # the positions of the spikes in both images, denoted by white # circles. + fig = plt.figure() -ax = fig.add_subplot(121, projection=m_cutout) +ax = fig.add_subplot(121, projection=aia_map_cutout) ax.plot_coord(spike_coords, "o", color="white", fillstyle="none", markersize=15) -m_cutout.plot(axes=ax, title='Level 1 "de-spiked" data') +aia_map_cutout.plot(axes=ax, title='Level 1 "de-spiked" data') lon, lat = ax.coords lon.set_axislabel("HPC Longitude") lat.set_axislabel("HPC Latitude") -ax = fig.add_subplot(122, projection=m_respiked_cutout) +ax = fig.add_subplot(122, projection=aia_map_respiked_cutout) ax.plot_coord(spike_coords, "o", color="white", fillstyle="none", markersize=15) -m_respiked_cutout.plot(axes=ax, annotate=False) +aia_map_respiked_cutout.plot(axes=ax, annotate=False) ax.set_title('Level 0.5 "re-spiked" data') lon, lat = ax.coords lon.set_axislabel("HPC Longitude") lat.set_axislabel(" ") lat.set_ticklabel_visible(visible=False) + plt.show() -########################################################### +############################################################################### # Lastly, let's check the metadata in both the level 1 and resulting # 0.5 images to double check that the appropriate keywords have been updated. for k in ["lvl_num", "nspikes", "comments"]: - print(f"Level 1: {k}: {m_cutout.meta.get(k)}") - print(f"Level 0.5: {k}: {m_respiked_cutout.meta.get(k)}") + print(f"Level 1: {k}: {aia_map_cutout.meta.get(k)}") + print(f"Level 0.5: {k}: {aia_map_respiked_cutout.meta.get(k)}") diff --git a/examples/skip_correct_degradation.py b/examples/skip_correct_degradation.py index 2c57e0e..1a01ea0 100644 --- a/examples/skip_correct_degradation.py +++ b/examples/skip_correct_degradation.py @@ -5,7 +5,6 @@ This example demonstrates the degradation of the filters on AIA over time. """ - import astropy.time import astropy.units as u import matplotlib.pyplot as plt @@ -15,30 +14,27 @@ from aiapy.calibrate import degradation -# These are needed to allow the use of quantities and astropy -# time objects in the plot. -time_support(format="jyear") -quantity_support() - -########################################################### +############################################################################### # The performance of the AIA telescope is unfortunately degrading over time, # leading to the resulting images becoming increasingly dim. We # can correct for this by modeling the degradation over time and # then dividing the image intensity by this correction. # # First, let's fetch some metadata for the 335 Å channel of AIA between 2010 -# and 2018 at a cadence of 30 days. We choose the 335 Å channel because it has experienced +# and 2022 at a cadence of 180 days. We choose the 335 Å channel because it has experienced # significant degradation compared to the other EUV channels. + results = Fido.search( - a.Time("2010-06-01T00:00:00", "2021-06-01T00:00:00"), - a.Sample(30 * u.day), + a.Time("2010-06-01T00:00:00", "2022-01-01T00:00:00"), + a.Sample(180 * u.day), a.jsoc.Series.aia_lev1_euv_12s, a.jsoc.Wavelength(335 * u.angstrom), ) -########################################################### +############################################################################### # We only need the date and mean intensity columns from the # metadata that was returned. We select those and nothing else. + table = results["jsoc"].show("DATE__OBS", "DATAMEAN") table["DATAMEAN"].unit = u.ct table["DATE_OBS"] = astropy.time.Time(table["DATE__OBS"], scale="utc") @@ -46,22 +42,29 @@ print(table) -########################################################### +############################################################################### # Next, we pass the date column to the `aiapy.calibrate.correct_degradation` # function. This function calculates the time-dependent correction factor # based on the time and wavelength of the observation. # We then divide the mean intensity by the correction factor to get the corrected intensity. # For more details on how the correction factor is calculated, see the documentation for the # `aiapy.calibrate.degradation` function. + correction_factor = degradation(335 * u.angstrom, table["DATE_OBS"]) # This correction can be applied to a sunpy Map as well. table["DATAMEAN_DEG"] = table["DATAMEAN"] / correction_factor -########################################################### +############################################################################### # To understand the effect of the degradation and the correction factor, we # plot the corrected and uncorrected mean intensity as a function of time. # Note that the uncorrected intensity decreases monotonically over time # while the corrected intensity recovers to pre-2011 values in 2020. + +# These are needed to allow the use of astropy quantities and +# time objects in the plot. +time_support(format="jyear") +quantity_support() + plt.plot(table["DATE_OBS"], table["DATAMEAN"], label="mean", marker="o") plt.plot(table["DATE_OBS"], table["DATAMEAN_DEG"], label="mean, corrected", marker="o") plt.title(f'{(335*u.angstrom).to_string(format="latex")} Channel Degradation') diff --git a/examples/skip_psf_deconvolution.py b/examples/skip_psf_deconvolution.py index 455b9f4..3dcfedc 100644 --- a/examples/skip_psf_deconvolution.py +++ b/examples/skip_psf_deconvolution.py @@ -15,7 +15,7 @@ import aiapy.data.sample as sample_data import aiapy.psf -######################################### +############################################################################### # AIA images are subject to convolution with the instrument point-spread # function (PSF) due to effects introduced by the filter mesh of the telescope # and the CCD, among others. This has the effect of "blurring" the image. @@ -29,39 +29,46 @@ # on the CCD grid. Once deconvolved, the image can be passed to # `aiapy.calibrate.register` # (see the :ref:`sphx_glr_generated_gallery_prepping_level_1_data.py` example). -m = sunpy.map.Map(sample_data.AIA_171_IMAGE) + +aia_map = sunpy.map.Map(sample_data.AIA_171_IMAGE) + fig = plt.figure() -ax = fig.add_subplot(111, projection=m) -m.plot( +ax = fig.add_subplot(111, projection=aia_map) +aia_map.plot( axes=ax, ) -####################################### +############################################################################### # Next, we'll calculate the PSF using `aiapy.psf.psf` for the 171 Å channel. # The PSF model accounts for several different effects, including diffraction # from the mesh grating of the filters, charge spreading, and jitter. See # `Grigis et al (2012) `_ -# for more details. Currently, this only works for -# :math:`4096\times4096` full frame images. +# for more details. Currently, this only works for :math:`4096\times4096` full frame images. # -# Note that this will be significantly faster if you have a GPU and the `cupy` +# Note that this will be significantly faster if you have a Nvidia GPU and the `cupy` # package installed. -psf = aiapy.psf.psf(m.wavelength) -############################################# +psf = aiapy.psf.psf(aia_map.wavelength) + +############################################################################### # We'll plot just a 500-by-500 pixel section centered on the center pixel. The # diffraction "arms" extending from the center pixel can often be seen in # flare observations due to the intense, small-scale brightening. + fov = 500 lc_x, lc_y = psf.shape[0] // 2 - fov // 2, psf.shape[1] // 2 - fov // 2 -plt.imshow( +fig = plt.figure() +ax = fig.add_subplot(111) +ax.imshow( psf[lc_x : lc_x + fov, lc_y : lc_y + fov], norm=ImageNormalize(vmin=1e-8, vmax=1e-3, stretch=LogStretch()), + origin="lower", ) -plt.colorbar() -plt.show() +ax.set_title("PSF") +ax.set_xlabel("Pixels") +ax.set_ylabel("Pixels") -############################################### +############################################################################### # Now that we've downloaded our image and computed the PSF, we can deconvolve # the image with the PSF using the # `Richardson-Lucy deconvolution algorithm `_. @@ -70,41 +77,54 @@ # wavelength, it is most efficient to only calculate the PSF once. # # As with `aiapy.psf.psf`, this will be much faster if you have -# a GPU and `cupy` installed. -m_deconvolved = aiapy.psf.deconvolve(m, psf=psf) +# a Nvidia GPU and `cupy` installed. + +aia_map_deconvolved = aiapy.psf.deconvolve(aia_map, psf=psf) -################################################ +############################################################################### # Let's compare the convolved and deconvolved images. -norm = ImageNormalize(vmin=0, vmax=1.5e4, stretch=AsinhStretch(0.01)) + fig = plt.figure() -ax = fig.add_subplot(121, projection=m) -m.plot(axes=ax, norm=norm) -ax = fig.add_subplot(122, projection=m_deconvolved) -m_deconvolved.plot(axes=ax, annotate=False, norm=norm) +ax = fig.add_subplot(121, projection=aia_map) +norm = ImageNormalize(vmin=0, vmax=1.5e4, stretch=AsinhStretch(0.01)) +aia_map.plot(axes=ax, norm=norm) +ax.set_title("Normal") + +ax = fig.add_subplot(122, projection=aia_map_deconvolved) +aia_map_deconvolved.plot(axes=ax, annotate=False, norm=norm) +ax.set_title("Deconvolved") ax.coords[0].set_axislabel(" ") ax.coords[1].set_axislabel(" ") ax.coords[1].set_ticklabel_visible(visible=False) -plt.show() +fig.tight_layout() -################################################# +############################################################################### # The differences become a bit more obvious when we zoom in. Note that the # deconvolution has the effect of "deblurring" the image. + left_corner = 500 * u.arcsec, -600 * u.arcsec right_corner = 1000 * u.arcsec, -100 * u.arcsec -fig = plt.figure() -m_sub = m.submap( - bottom_left=SkyCoord(*left_corner, frame=m.coordinate_frame), - top_right=SkyCoord(*right_corner, frame=m.coordinate_frame), + +aia_map_sub = aia_map.submap( + bottom_left=SkyCoord(*left_corner, frame=aia_map.coordinate_frame), + top_right=SkyCoord(*right_corner, frame=aia_map.coordinate_frame), ) -ax = fig.add_subplot(121, projection=m_sub) -m_sub.plot(axes=ax, norm=norm) -m_deconvolved_sub = m_deconvolved.submap( - bottom_left=SkyCoord(*left_corner, frame=m_deconvolved.coordinate_frame), - top_right=SkyCoord(*right_corner, frame=m_deconvolved.coordinate_frame), +aia_map_deconvolved_sub = aia_map_deconvolved.submap( + bottom_left=SkyCoord(*left_corner, frame=aia_map_deconvolved.coordinate_frame), + top_right=SkyCoord(*right_corner, frame=aia_map_deconvolved.coordinate_frame), ) -ax = fig.add_subplot(122, projection=m_deconvolved_sub) -m_deconvolved_sub.plot(axes=ax, annotate=False, norm=norm) + +fig = plt.figure() + +ax = fig.add_subplot(121, projection=aia_map_sub) +aia_map_sub.plot(axes=ax, norm=norm) +ax.set_title("Normal") + +ax = fig.add_subplot(122, projection=aia_map_deconvolved_sub) +aia_map_deconvolved_sub.plot(axes=ax, annotate=False, norm=norm) +ax.set_title("Deconvolved") ax.coords[0].set_axislabel(" ") ax.coords[1].set_axislabel(" ") ax.coords[1].set_ticklabel_visible(visible=False) + plt.show() diff --git a/examples/update_header_keywords.py b/examples/update_header_keywords.py index f3e0cef..92a7d1c 100644 --- a/examples/update_header_keywords.py +++ b/examples/update_header_keywords.py @@ -8,19 +8,19 @@ information regarding the spacecraft pointing and observer position. """ - +import matplotlib.pyplot as plt import sunpy.map import aiapy.data.sample as sample_data from aiapy.calibrate import fix_observer_location, update_pointing -########################################################### +############################################################################### # An AIA FITS header contains various pieces of # `standard `_. # metadata that are critical to the physical interpretation of the data. # These include the pointing of the spacecraft, necessary for connecting # positions on the pixel grid to physical locations on the Sun, as well as -# the observer (i.e. satellite) location. +# the observer (i.e., satellite) location. # # While this metadata is recorded in the FITS header, some values in # the headers exported by data providers (e.g. @@ -32,62 +32,74 @@ # # For this example, we will read a 171 Å image from the aiapy sample data # into a `~sunpy.map.Map` object. -m = sunpy.map.Map(sample_data.AIA_171_IMAGE) -########################################################### +aia_map = sunpy.map.Map(sample_data.AIA_171_IMAGE) + +############################################################################### # To update the pointing keywords, we can pass our `~sunpy.map.Map` to the # `aiapy.calibrate.update_pointing` function. This function will query the # JSOC, using `~sunpy`, for the most recent pointing information, update # the metadata, and then return a new `~sunpy.map.Map` with this updated # metadata. -m_updated_pointing = update_pointing(m) -############################################################ -# If we inspect the reference pixel and rotation matrix of the original map -print(m.reference_pixel) -print(m.rotation_matrix) +aia_map_updated_pointing = update_pointing(aia_map) + +############################################################################### +# If we inspect the reference pixel and rotation matrix of the original map: + +print(aia_map.reference_pixel) +print(aia_map.rotation_matrix) + +############################################################################### +# and the map with the updated pointing information: -############################################################ -# and the map with the updated pointing information -print(m_updated_pointing.reference_pixel) -print(m_updated_pointing.rotation_matrix) +print(aia_map_updated_pointing.reference_pixel) +print(aia_map_updated_pointing.rotation_matrix) -############################################################ -# we find that the relevant keywords, `CRPIX1`, `CRPIX2`, `CDELT1`, `CDELT2`, +############################################################################### +# We find that the relevant keywords, `CRPIX1`, `CRPIX2`, `CDELT1`, `CDELT2`, # and `CROTA2`, have been updated. # # Similarly, the Heliographic Stonyhurst (HGS) coordinates of the observer # location in the header are inaccurate. If we check the HGS longitude keyword # in the header, we find that it is 0 degrees which is not the HGS longitude # coordinate of SDO. -print(m_updated_pointing.meta["hgln_obs"]) -print(m_updated_pointing.meta["hglt_obs"]) -############################################################ +print(aia_map_updated_pointing.meta["hgln_obs"]) +print(aia_map_updated_pointing.meta["hglt_obs"]) + +############################################################################### # To update the HGS observer coordinates, we can use the # `aiapy.calibrate.fix_observer_location` function. This function reads the # correct observer location from Heliocentric Aries Ecliptic (HAE) coordinates # in the header, converts them to HGS, and replaces the inaccurate HGS # keywords. -m_observer_fixed = fix_observer_location(m_updated_pointing) -############################################################ +aia_map_observer_fixed = fix_observer_location(aia_map_updated_pointing) + +############################################################################### # Looking again at the HGS longitude and latitude keywords, we can see that # they have been updated. -print(m_observer_fixed.meta["hgln_obs"]) -print(m_observer_fixed.meta["hglt_obs"]) +print(aia_map_observer_fixed.meta["hgln_obs"]) +print(aia_map_observer_fixed.meta["hglt_obs"]) -############################################################ +############################################################################### # Note that in `~sunpy.map.AIAMap`, the `~sunpy.map.Map.observer_coordinate` # attribute is already derived from the HAE coordinates such that it is not # strictly necessary to apply `aiapy.calibrate.fix_observer_location`. For # example, the unfixed `~sunpy.map.Map` will still have an accurate derived # observer position -print(m_updated_pointing.observer_coordinate) -############################################################ +print(aia_map_updated_pointing.observer_coordinate) + +############################################################################### # However, we suggest that users apply this fix such that the information # stored in `~sunpy.map.Map.meta` is accurate and consistent. # # Finally, plot the fixed map. -m_observer_fixed.peek() + +fig = plt.figure() +ax = fig.add_subplot(projection=aia_map_observer_fixed) +aia_map_observer_fixed.plot(axes=ax) + +plt.show()