Source code for aiapy.calibrate.util

"""
Utilities for computing intensity corrections.
"""
import os
import pathlib
import warnings
from urllib.parse import urljoin

import astropy.io.ascii
import astropy.units as u
import drms
import numpy as np
from astropy.table import QTable
from astropy.time import Time
from erfa.core import ErfaWarning
from sunpy.data import manager
from sunpy.net import attrs, jsoc

from aiapy import _SSW_MIRRORS
from aiapy.util.decorators import validate_channel
from aiapy.util.exceptions import AiapyUserWarning

__all__ = ["get_correction_table", "get_pointing_table", "get_error_table"]

# Default version of the degradation calibration curve to use.
# This needs to be incremented as the calibration is updated in JSOC.
CALIBRATION_VERSION = 10
# Error table filename available from SSW
AIA_ERROR_FILE = "sdo/aia/response/aia_V{}_error_table.txt"
# Most recent version number for error tables; increment as new versions become available
ERROR_VERSION = 3
# URLs and SHA-256 hashes for each version of the error tables
URL_HASH = {
    2: (
        [urljoin(mirror, AIA_ERROR_FILE.format(2)) for mirror in _SSW_MIRRORS],
        "ac97ccc48057809723c27e3ef290c7d78ee35791d9054b2188baecfb5c290d0a",
    ),
    3: (
        [urljoin(mirror, AIA_ERROR_FILE.format(3)) for mirror in _SSW_MIRRORS],
        "66ff034923bb0fd1ad20e8f30c7d909e1a80745063957dd6010f81331acaf894",
    ),
}


[docs] def get_correction_table(*, correction_table=None): """ Return table of degradation correction factors. This function returns a table of parameters for estimating the time-dependent degradation of the instrument. By default, this table is queried from `aia.response` series in `JSOC <http://jsoc.stanford.edu/>`_. The correction table can also be read from a file by passing a filepath to `correction_table`. These files are typically included in the SDO tree of an SSW installation in `$SSW/sdo/aia/response/` with filenames like `aia_V*_response_table.txt`. Parameters ---------- correction_table: `str` or `~astropy.table.QTable`, optional Path to correction table file or an existing correction table. If None, the table will be queried from JSOC. Returns ------- `~astropy.table.QTable` See Also -------- aiapy.calibrate.degradation """ if isinstance(correction_table, astropy.table.QTable): return correction_table if correction_table is not None: if isinstance(correction_table, (str, pathlib.Path)): table = QTable(astropy.io.ascii.read(correction_table)) else: raise ValueError("correction_table must be a file path, an existing table, or None.") else: # NOTE: the [!1=1!] disables the drms PrimeKey logic and enables # the query to find records that are ordinarily considered # identical because the PrimeKeys for this series are WAVE_STR # and T_START. Without the !1=1! the query only returns the # latest record for each unique combination of those keywords. table = drms.Client().query("aia.response[][!1=1!]", key="**ALL**") table = QTable.from_pandas(table) selected_cols = [ "DATE", "VER_NUM", "WAVE_STR", "WAVELNTH", "T_START", "T_STOP", "EFFA_P1", "EFFA_P2", "EFFA_P3", "EFF_AREA", "EFF_WVLN", ] table = table[selected_cols] table["T_START"] = Time(table["T_START"], scale="utc") # NOTE: The warning from erfa here is due to the fact that dates in # this table include at least one date from 2030 and converting this # date to UTC is ambiguous as the UTC conversion is not well defined # at this date. with warnings.catch_warnings(): warnings.simplefilter("ignore", category=ErfaWarning) table["T_STOP"] = Time(table["T_STOP"], scale="utc") table["WAVELNTH"].unit = "Angstrom" table["EFF_WVLN"].unit = "Angstrom" table["EFF_AREA"].unit = "cm2" return table
@u.quantity_input @validate_channel("channel") def _select_epoch_from_correction_table(channel: u.angstrom, obstime, table, *, version=None): """ Return correction table with only the first epoch and the epoch in which `obstime` falls and for only one given calibration version. Parameters ---------- channel : `~astropy.units.Quantity` obstime : `~astropy.time.Time` table : `~astropy.table.QTable` version : `int` """ version = CALIBRATION_VERSION if version is None else version # Select only this channel # NOTE: The WAVE_STR prime keys for the aia.response JSOC series for the # non-EUV channels do not have a thick/thin designation thin = "_THIN" if channel not in (1600, 1700, 4500) * u.angstrom else "" wave = channel.to(u.angstrom).value table = table[table["WAVE_STR"] == f"{wave:.0f}{thin}"] table = table[table["VER_NUM"] == version] table.sort("DATE") # Newest entries will be last if len(table) == 0: extra_msg = " Max version is 3." if channel == 4500 * u.AA else "" raise ValueError( f"Correction table does not contain calibration for version {version} for {channel}." + extra_msg, ) # Select the epoch for the given observation time obstime_in_epoch = np.logical_and(obstime >= table["T_START"], obstime < table["T_STOP"]) if not obstime_in_epoch.any(): raise ValueError(f"No valid calibration epoch for {obstime}") # NOTE: In some cases, there may be multiple entries for a single epoch. We want to # use the most up-to-date one. i_epoch = np.where(obstime_in_epoch)[0] if i_epoch.shape[0] > 1: warnings.warn( f"Multiple valid epochs for {obstime}. Using the most recent one", AiapyUserWarning, stacklevel=3, ) # Create new table with only first and obstime epochs return QTable(table[[0, i_epoch[-1]]])
[docs] def get_pointing_table(start, end): """ Retrieve 3-hourly master pointing table from the JSOC. This function queries `JSOC <http://jsoc.stanford.edu/>`__ for the 3-hourly master pointing table (MPT) in the interval defined by `start` and `end`. The 3-hourly MPT entries are computed from limb fits of images with ``T_OBS`` between ``T_START`` and ``T_STOP``. .. note:: A MPT entry covers the interval ``[T_START:T_STOP)``; that is, the interval includes ``T_START`` and excludes ``T_STOP``. .. note:: While it is generally true that ``TSTOP = T_START + 3 hours``, there are edge cases where ``T_STOP`` is more than 3 hours after ``T_START`` because of a calibration, an eclipse, or other reasons, but the fits are still calculated based on images from ``T_START`` to ``T_START + 3 hours``. Pointing is not stable during these periods, so the question of which MPT entry to use is not relevant. Parameters ---------- start : `~astropy.time.Time` end : `~astropy.time.Time` Returns ------- `~astropy.table.QTable` See Also -------- aiapy.calibrate.update_pointing """ q = jsoc.JSOCClient().search( attrs.Time(start, end=end), attrs.jsoc.Series.aia_master_pointing3h, ) table = QTable(q) if len(table.columns) == 0: # If there's no pointing information available between these times, # JSOC will raise a cryptic KeyError # (see https://github.com/LM-SAL/aiapy/issues/71) raise RuntimeError(f"Could not find any pointing information between {start} and {end}") table["T_START"] = Time(table["T_START"], scale="utc") table["T_STOP"] = Time(table["T_STOP"], scale="utc") for c in table.colnames: if "X0" in c or "Y0" in c: table[c].unit = "pixel" if "IMSCALE" in c: table[c].unit = "arcsecond / pixel" if "INSTROT" in c: table[c].unit = "degree" # Remove masking on columns with pointing parameters for c in table.colnames: if any(n in c for n in ["X0", "Y0", "IMSCALE", "INSTROT"]): if hasattr(table[c], "mask"): table[c] = table[c].filled(np.nan) return table
[docs] def get_error_table(error_table=None): if error_table is None: # This is to work around a parfive bug # https://github.com/Cadair/parfive/issues/121 os.environ["PARFIVE_DISABLE_RANGE"] = "1" error_table = fetch_error_table() os.environ.pop("PARFIVE_DISABLE_RANGE") if isinstance(error_table, (str, pathlib.Path)): table = astropy.io.ascii.read(error_table) elif isinstance(error_table, QTable): table = error_table else: raise ValueError("error_table must be a file path, an existing table, or None.") table = QTable(table) table["DATE"] = Time(table["DATE"], scale="utc") table["T_START"] = Time(table["T_START"], scale="utc") # NOTE: The warning from erfa here is due to the fact that dates in # this table include at least one date from 2030 and converting this # date to UTC is ambiguous as the UTC conversion is not well defined # at this date. with warnings.catch_warnings(): warnings.simplefilter("ignore", category=ErfaWarning) table["T_STOP"] = Time(table["T_STOP"], scale="utc") table["WAVELNTH"] = u.Quantity(table["WAVELNTH"], "Angstrom") table["DNPERPHT"] = u.Quantity(table["DNPERPHT"], "ct photon-1") return table
@manager.require("error_table", *URL_HASH[ERROR_VERSION]) def fetch_error_table(): return manager.get("error_table")