Skip to content

Commit

Permalink
Fallback to remote data files while JSOC is unavailable (#346)
Browse files Browse the repository at this point in the history
Co-authored-by: Will Barnes <[email protected]>
  • Loading branch information
nabobalis and wtbarnes authored Jan 8, 2025
1 parent ca0cca9 commit db3c654
Show file tree
Hide file tree
Showing 36 changed files with 418 additions and 511 deletions.
4 changes: 2 additions & 2 deletions .github/workflows/ci.yml
Original file line number Diff line number Diff line change
Expand Up @@ -74,7 +74,7 @@ jobs:
.tox/sample_data/
cache-key: docs-${{ github.run_id }}
envs: |
- linux: build_docs
- linux: build-docs
online:
if: "!startsWith(github.event.ref, 'refs/tags/v')"
Expand All @@ -86,7 +86,7 @@ jobs:
coverage: codecov
toxdeps: tox-pypi-filter
envs: |
- linux: build_docs-gallery
- linux: build-docs-gallery
pytest: false
cache-path: |
docs/_build/
Expand Down
1 change: 1 addition & 0 deletions .gitignore
Original file line number Diff line number Diff line change
Expand Up @@ -74,6 +74,7 @@ instance/

# Sphinx documentation
docs/_build/
docs/generated/
# automodapi
docs/api
docs/sg_execution_times.rst
Expand Down
71 changes: 7 additions & 64 deletions aiapy/calibrate/meta.py
Original file line number Diff line number Diff line change
Expand Up @@ -8,80 +8,27 @@
import numpy as np

import astropy.units as u
from astropy.coordinates import CartesianRepresentation, HeliocentricMeanEcliptic, SkyCoord

from sunpy.map import contains_full_disk

from aiapy.calibrate.util import get_pointing_table
from aiapy.util.exceptions import AIApyUserWarning

__all__ = ["fix_observer_location", "update_pointing"]
__all__ = ["update_pointing"]


def fix_observer_location(smap):
"""
Fix inaccurate ``HGS_LON`` and ``HGS_LAT`` FITS keywords.
The heliographic Stonyhurst latitude and longitude locations in the
AIA FITS headers are incorrect. This function fixes the values of these
keywords using the heliocentric aries ecliptic keywords, ``HAEX_OBS,
HAEY_OBS, HAEZ_OBS``.
.. note::
`~sunpy.map.sources.AIAMap` already accounts for the inaccurate
HGS keywords by using the HAE keywords to construct the
derived observer location.
Parameters
----------
smap : `~sunpy.map.sources.AIAMap`
Input map.
Returns
-------
`~sunpy.map.sources.AIAMap`
"""
# Create observer coordinate from HAE coordinates
coord = SkyCoord(
x=smap.meta["haex_obs"] * u.m,
y=smap.meta["haey_obs"] * u.m,
z=smap.meta["haez_obs"] * u.m,
representation_type=CartesianRepresentation,
frame=HeliocentricMeanEcliptic,
obstime=smap.reference_date,
).heliographic_stonyhurst
# Update header
new_meta = copy.deepcopy(smap.meta)
new_meta["hgln_obs"] = coord.lon.to(u.degree).value
new_meta["hglt_obs"] = coord.lat.to(u.degree).value
new_meta["dsun_obs"] = coord.radius.to(u.m).value

return smap._new_instance(smap.data, new_meta, plot_settings=smap.plot_settings, mask=smap.mask)


def update_pointing(smap, *, pointing_table=None):
def update_pointing(smap, *, pointing_table):
"""
Update the pointing information in the input map header.
This function updates the pointing information in ``smap`` by
updating the ``CRPIX1, CRPIX2, CDELT1, CDELT2, CROTA2`` keywords
in the header using the information provided in ``pointing_table``.
If ``pointing_table`` is not specified, the 3-hour pointing
information is queried from the `JSOC <http://jsoc.stanford.edu/>`_.
.. note::
The method removes any ``PCi_j`` matrix keys in the header and
updates the ``CROTA2`` keyword.
.. note::
If correcting pointing information for a large number of images,
it is strongly recommended to query the table once for the
appropriate interval and then pass this table in rather than
executing repeated queries.
.. warning::
This function is only intended to be used for full-disk images
Expand All @@ -92,13 +39,14 @@ def update_pointing(smap, *, pointing_table=None):
----------
smap : `~sunpy.map.sources.AIAMap`
Input map.
pointing_table : `~astropy.table.QTable`, optional
Table of pointing information. If not specified, the table
will be retrieved from JSOC.
pointing_table : `~astropy.table.QTable`
Table of pointing information.
You can get this table by calling `aiapy.calibrate.util.get_pointing_table`.
Returns
-------
`~sunpy.map.sources.AIAMap`
Updated map with pointing information.
See Also
--------
Expand All @@ -111,9 +59,6 @@ def update_pointing(smap, *, pointing_table=None):
if not all(d == (s * u.pixel) for d, s in zip(smap.dimensions, shape_full_frame, strict=True)):
msg = f"Input must be at the full resolution of {shape_full_frame}"
raise ValueError(msg)
if pointing_table is None:
# Make range wide enough to get closest 3-hour pointing
pointing_table = get_pointing_table(smap.date - 12 * u.h, smap.date + 12 * u.h)
# Find row in which T_START <= T_OBS < T_STOP
# The following notes are from a private communication with J. Serafin (LMSAL)
# and are preserved here to explain the reasoning for selecting the particular
Expand Down Expand Up @@ -176,9 +121,7 @@ def update_pointing(smap, *, pointing_table=None):
)
else:
new_meta[key] = value

# sunpy map converts crota to a PCi_j matrix, so we remove it to force the
# re-conversion.
# sunpy.map.Map converts crota to a PCi_j matrix, so we remove it to force the re-conversion.
new_meta.pop("PC1_1")
new_meta.pop("PC1_2")
new_meta.pop("PC2_1")
Expand Down
59 changes: 22 additions & 37 deletions aiapy/calibrate/prep.py
Original file line number Diff line number Diff line change
Expand Up @@ -13,7 +13,7 @@
from sunpy.util.decorators import add_common_docstring

from aiapy.calibrate.transform import _rotation_function_names
from aiapy.calibrate.util import _select_epoch_from_correction_table, get_correction_table
from aiapy.calibrate.util import _select_epoch_from_correction_table
from aiapy.util import AIApyUserWarning
from aiapy.util.decorators import validate_channel

Expand Down Expand Up @@ -101,7 +101,7 @@ def register(smap, *, missing=None, order=3, method="scipy"):
missing=missing,
method=method,
)
# extract center from padded smap.rotate output
# Extract center from padded smap.rotate output
# crpix1 and crpix2 will be equal (recenter=True), as prep does not work with submaps
center = np.floor(tempmap.meta["crpix1"])
range_side = (center + np.array([-1, 1]) * smap.data.shape[0] / 2) * u.pix
Expand All @@ -115,7 +115,7 @@ def register(smap, *, missing=None, order=3, method="scipy"):
return newmap


def correct_degradation(smap, *, correction_table=None, calibration_version=None):
def correct_degradation(smap, *, correction_table):
"""
Apply time-dependent degradation correction to an AIA map.
Expand All @@ -128,22 +128,14 @@ def correct_degradation(smap, *, correction_table=None, calibration_version=None
----------
smap : `~sunpy.map.sources.AIAMap`
Map to be corrected.
correction_table : `~astropy.table.Table` or `str`, optional
Table of correction parameters or path to correction table file.
If not specified, it will be queried from JSOC. See
`aiapy.calibrate.util.get_correction_table` for more information.
If you are processing many images, it is recommended to
read the correction table once and pass it with this argument to avoid
multiple redundant network calls.
calibration_version : `int`, optional
The version of the calibration to use when calculating the degradation.
By default, this is the most recent version available from JSOC. If you
are using a specific calibration response file, you may need to specify
this according to the version in that file.
correction_table : `~astropy.table.Table`
Table of correction parameters.
You can get this table by calling `aiapy.calibrate.util.get_correction_table`.
Returns
-------
`~sunpy.map.sources.AIAMap`
Degradation-corrected map.
See Also
--------
Expand All @@ -153,7 +145,6 @@ def correct_degradation(smap, *, correction_table=None, calibration_version=None
smap.wavelength,
smap.date,
correction_table=correction_table,
calibration_version=calibration_version,
)
return smap / d

Expand All @@ -164,8 +155,7 @@ def degradation(
channel: u.angstrom,
obstime,
*,
correction_table=None,
calibration_version=None,
correction_table,
) -> u.dimensionless_unscaled:
r"""
Correction to account for time-dependent degradation of the instrument.
Expand Down Expand Up @@ -195,19 +185,17 @@ def degradation(
Parameters
----------
channel : `~astropy.units.Quantity`
Wavelength of the channel.
obstime : `~astropy.time.Time`
correction_table : `~astropy.table.Table` or `str`, optional
Table of correction parameters or path to correction table file.
If not specified, it will be queried from JSOC. See
`aiapy.calibrate.util.get_correction_table` for more information.
If you are processing many images, it is recommended to
read the correction table once and pass it with this argument to avoid
multiple redundant network calls.
calibration_version : `int`, optional
The version of the calibration to use when calculating the degradation.
By default, this is the most recent version available from JSOC. If you
are using a specific calibration response file, you may need to specify
this according to the version in that file.
Observation time.
correction_table : `~astropy.table.Table`
Table of correction parameters.
You can get this table by calling `aiapy.calibrate.util.get_correction_table`.
Returns
-------
`~astropy.units.Quantity`
Degradation correction factor.
See Also
--------
Expand All @@ -219,15 +207,12 @@ def degradation(
obstime = obstime.reshape((1,))
ratio = np.zeros(obstime.shape)
poly = np.zeros(obstime.shape)
# Do this outside of the loop to avoid repeated queries
correction_table = get_correction_table(correction_table=correction_table)
for i, t in enumerate(obstime):
table = _select_epoch_from_correction_table(channel, t, correction_table, version=calibration_version)

for idx, t in enumerate(obstime):
table = _select_epoch_from_correction_table(channel, t, correction_table)
# Time difference between obstime and start of epoch
dt = (t - table["T_START"][-1]).to(u.day).value
# Correction to most recent epoch
ratio[i] = table["EFF_AREA"][-1] / table["EFF_AREA"][0]
ratio[idx] = table["EFF_AREA"][-1] / table["EFF_AREA"][0]
# Polynomial correction to interpolate within epoch
poly[i] = table["EFFA_P1"][-1] * dt + table["EFFA_P2"][-1] * dt**2 + table["EFFA_P3"][-1] * dt**3 + 1.0
poly[idx] = table["EFFA_P1"][-1] * dt + table["EFFA_P2"][-1] * dt**2 + table["EFFA_P3"][-1] * dt**3 + 1.0
return u.Quantity(poly * ratio)
11 changes: 3 additions & 8 deletions aiapy/calibrate/spikes.py
Original file line number Diff line number Diff line change
Expand Up @@ -11,11 +11,11 @@
from astropy.io import fits
from astropy.wcs.utils import pixel_to_pixel

import drms
from sunpy.map.mapbase import PixelPair
from sunpy.map.sources.sdo import AIAMap

from aiapy.util import AIApyUserWarning
from aiapy.util.net import _get_data_from_jsoc

__all__ = ["fetch_spikes", "respike"]

Expand Down Expand Up @@ -145,13 +145,8 @@ def fetch_spikes(smap, *, as_coords=False):
array-like
Original intensity values of the spikes
"""
series = r"aia.lev1_euv_12s"
if smap.wavelength in (1600, 1700, 4500) * u.angstrom:
series = r"aia.lev1_uv_24s"
file = drms.Client().query(
f'{series}[{smap.date}/12s][WAVELNTH={smap.meta["wavelnth"]}]',
seg="spikes",
)
series = "aia.lev1_uv_24s" if smap.wavelength in (1600, 1700, 4500) * u.angstrom else "aia.lev1_euv_12s"
file = _get_data_from_jsoc(f'{series}[{smap.date}/12s][WAVELNTH={smap.meta["wavelnth"]}]', key=None, seg="spikes")
_, spikes = fits.open(f'http://jsoc.stanford.edu{file["spikes"][0]}')
# Loaded as floats, but they are actually integers
spikes = spikes.data.astype(np.int32)
Expand Down
21 changes: 8 additions & 13 deletions aiapy/calibrate/tests/test_meta.py
Original file line number Diff line number Diff line change
Expand Up @@ -6,22 +6,14 @@
from astropy.table import QTable
from astropy.time import Time, TimeDelta

from aiapy.calibrate import fix_observer_location, update_pointing
from aiapy.calibrate import update_pointing
from aiapy.calibrate.util import get_pointing_table
from aiapy.util.exceptions import AIApyUserWarning


def test_fix_observer_location(aia_171_map) -> None:
smap_fixed = fix_observer_location(aia_171_map)
# NOTE: AIAMap already fixes the .observer_coordinate property with HAE
assert smap_fixed.meta["hgln_obs"] == smap_fixed.observer_coordinate.lon.value
assert smap_fixed.meta["hglt_obs"] == smap_fixed.observer_coordinate.lat.value
assert smap_fixed.meta["dsun_obs"] == smap_fixed.observer_coordinate.radius.value


@pytest.fixture
def pointing_table(aia_171_map):
return get_pointing_table(aia_171_map.date - 6 * u.h, aia_171_map.date + 6 * u.h)
def pointing_table():
return get_pointing_table("lmsal")


@pytest.fixture
Expand Down Expand Up @@ -52,7 +44,10 @@ def test_fix_pointing(aia_171_map, pointing_table) -> None:
# Remove keys to at least test that they get set
for k in keys:
aia_171_map.meta.pop(k)
aia_map_updated = update_pointing(aia_171_map)
aia_map_updated = update_pointing(
aia_171_map,
pointing_table=pointing_table,
)
# FIXME: how do we check these values are accurate?
assert all(k in aia_map_updated.meta for k in keys)
# Check the case where we have specified the pointing
Expand Down Expand Up @@ -109,7 +104,7 @@ def test_update_pointing_no_entry_raises_exception(aia_171_map, pointing_table)
# This tests that an exception is thrown when entry corresponding to
# T_START <= T_OBS < T_END cannot be found in the pointing table.
# We explicitly set the T_OBS key
aia_171_map.meta["T_OBS"] = (aia_171_map.date + 1 * u.day).isot
aia_171_map.meta["T_OBS"] = (aia_171_map.date - 1000 * u.day).isot
with pytest.raises(IndexError, match="No valid entries for"):
update_pointing(aia_171_map, pointing_table=pointing_table)

Expand Down
Loading

0 comments on commit db3c654

Please sign in to comment.