From 2e95ade224264def51fa288b333b8e112577d485 Mon Sep 17 00:00:00 2001 From: Roy Smart Date: Sat, 10 Oct 2026 10:53:55 -0600 Subject: [PATCH 1/3] Correct the pointing of XRT images with the co-alignment databases of SolarSoft, by default from the cross-correlation with AIA. Ports xrt_read_coaldb.pro (with and without /aia_cc) and xrt_rollangle.pro as the `coalign` option of `hinode.xrt.open` and `Filtergram`: "aia" (the default) uses coaldb_c and falls back to coaldb_u, "ufss" uses coaldb_u, and None keeps the headers. `Filtergram.coalignment` records the CALIBRATION_TYPE of each image. The Level 1 files already carry the UFSS pointing, so "aia" moves the 2019-09-30 18:08:00 image by -0.34", +0.68", and the roll is recomputed, which changes the files made in 2022 from the pre-anomaly roll by more than a degree. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_011VGPHLGJi8acbuLYSfBHTk --- docs/index.rst | 3 + docs/refs.bib | 11 + hinode/xrt/_coalign.py | 538 ++++++++++++++++++++++++++++++++ hinode/xrt/_coalign_test.py | 333 ++++++++++++++++++++ hinode/xrt/_data.py | 154 ++++++--- hinode/xrt/_filtergrams.py | 93 +++++- hinode/xrt/_filtergrams_test.py | 32 +- hinode/xrt/_xrt.py | 10 + hinode/xrt/_xrt_test.py | 11 + pyproject.toml | 1 + 10 files changed, 1141 insertions(+), 45 deletions(-) create mode 100644 hinode/xrt/_coalign.py create mode 100644 hinode/xrt/_coalign_test.py diff --git a/docs/index.rst b/docs/index.rst index 7054cbd..35dcd6d 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -17,6 +17,9 @@ the visible light leaking into each pixel, :attr:`hinode.xrt.Filtergram.leak`, and the uncertainty of each pixel, :meth:`hinode.xrt.Filtergram.uncertainty`. +The pointing of each image is corrected with the co-alignment databases of +SolarSoft, so that it lines up with AIA, +:attr:`hinode.xrt.Filtergram.coalignment`. Installation ============ diff --git a/docs/refs.bib b/docs/refs.bib index fe47555..a92dbff 100644 --- a/docs/refs.bib +++ b/docs/refs.bib @@ -63,3 +63,14 @@ @article{Takeda2016 year = {2016}, doi = {10.1007/s11207-015-0823-8}, } + +@article{Yoshimura2015, + author = {Yoshimura, Keiji and McKenzie, David E.}, + title = {Calibration of Hinode/XRT for Coalignment}, + journal = {Solar Physics}, + volume = {290}, + number = {8}, + pages = {2355--2372}, + year = {2015}, + doi = {10.1007/s11207-015-0746-4}, +} diff --git a/hinode/xrt/_coalign.py b/hinode/xrt/_coalign.py new file mode 100644 index 0000000..8b49887 --- /dev/null +++ b/hinode/xrt/_coalign.py @@ -0,0 +1,538 @@ +from typing import Sequence, Literal +import re +import time +import pathlib +import numpy as np +import hinode +from ._data import _get, _path_local, _write, _download_file + +__all__ = [] + +_url_coalign = "https://sohoftp.nascom.nasa.gov/sdb/hinode/xrt/xrt_msu_coalign/" +""" +The directory of the SolarSoft database which holds the co-alignment +databases of XRT and the database of its roll angle. +""" + +_directories_coalign = { + "aia": "coaldb_c", + "ufss": "coaldb_u", +} +""" +The directory of each co-alignment database: +``coaldb_c``, from the cross-correlation of the images with AIA 335 Å, +and ``coaldb_u``, from the Ultra Fine Sun Sensors (UFSS) of Hinode. +""" + +_pattern_coalign = re.compile(r'href="(\d{8}_\d{4})\.geny"') +""" +The name of a file of a co-alignment database in a listing of its directory, +which is the start of the time it covers, as in ``20190930_1503.geny``. +""" + +_name_listing = "listing.txt" +""" +The file in which the names of the files of a co-alignment database are kept +once its directory has been listed, +since the listing of ``coaldb_c`` is several megabytes. +""" + +_signature_genx = b"SR\x00\x04" +""" +The first bytes of every IDL save file, +which is what a GENX file of SolarSoft is. +""" + +_age_refresh = 86400 +""" +How old, in seconds, a downloaded listing or file of a database must be +before it is downloaded again to look for images it may not have had yet, +since the databases are updated about once a week. +""" + +_tolerance = 0.1 +""" +The largest difference, in seconds, between the start of the exposure of an +image and the time of an entry of a database for the entry to be that of the +image, as in ``xrt_read_coaldb.pro``. +""" + +_scale_xray = 1.0286 +"""The plate scale of the X-ray images, in arcseconds per unsummed pixel.""" + +_scale_gband = 1.0302 +"""The plate scale of the G-band images, in arcseconds per unsummed pixel.""" + +_position_gband = 3 +"""The position of the G-band filter in the second filter wheel, ``EC_FW2``.""" + +_epoch = np.datetime64("1979-01-01T00:00:00", "us") +"""The time from which ``anytim`` in SolarSoft counts seconds.""" + +_times_roll = np.array( + [ + "2012-03-10", + "2012-06-08", + "2013-01-21", + "2021-12-27T13:00", + "2021-12-27T13:28", + "2021-12-27T14:40", + "2021-12-27T15:10", + ], + dtype="datetime64[us]", +) +""" +The times at which ``xrt_rollangle.pro`` changes from one model of the roll +to the next, +the last four of them during the attitude anomaly of 2021 December 27. +""" + + +def _seconds(time: Sequence[str] | np.ndarray) -> np.ndarray: + """ + The seconds since 1979 January 1 of times in the ISOT format, + counted without leap seconds, as ``anytim`` in SolarSoft counts them. + + Parameters + ---------- + time + The times, in UTC. + """ + return (np.asarray(time, dtype="datetime64[us]") - _epoch) / np.timedelta64(1, "s") + + +def _start(names: Sequence[str]) -> np.ndarray: + """ + The start of the time covered by each file of a co-alignment database, + in seconds since 1979 January 1, from its name. + + Parameters + ---------- + names + The names of the files, without their extension, + such as ``20190930_1503``. + """ + isot = [f"{n[:4]}-{n[4:6]}-{n[6:8]}T{n[9:11]}:{n[11:13]}" for n in names] + return _seconds(isot) + + +def _stale(path: pathlib.Path) -> bool: + """ + Whether a downloaded file is old enough to be downloaded again. + + Parameters + ---------- + path + The file. + """ + return time.time() - path.stat().st_mtime > _age_refresh + + +def _restore(path: pathlib.Path) -> dict: + """ + The variables of an IDL save file, such as a GENX file of SolarSoft. + + Parameters + ---------- + path + The file. + """ + # Imported here, since only the pointing needs it + import scipy.io + + return scipy.io.readsav(str(path)) + + +def _interpol( + x: np.ndarray, + xp: np.ndarray, + fp: np.ndarray, +) -> np.ndarray: + """ + Linear interpolation which extrapolates beyond the ends from the first and + the last interval, as ``interpol`` in IDL does. + + Parameters + ---------- + x + Where to interpolate. + xp + The increasing coordinates of the samples. + fp + The value of each sample. + """ + result = np.interp(x, xp, fp) + below = x < xp[0] + above = x > xp[-1] + slope_below = (fp[1] - fp[0]) / (xp[1] - xp[0]) + slope_above = (fp[-1] - fp[-2]) / (xp[-1] - xp[-2]) + result[below] = fp[0] + slope_below * (x[below] - xp[0]) + result[above] = fp[-1] + slope_above * (x[above] - xp[-1]) + return result + + +def _database_roll( + latest: float, + directory: pathlib.Path, + num_retry: int = 5, +) -> tuple[np.ndarray, np.ndarray]: + """ + The database of the roll of XRT relative to SDO since the attitude anomaly + of 2021 December 27, which ``xrt_rollangle.pro`` interpolates: + the time of each measurement, in seconds since 1979 January 1, + and the roll, in degrees, opposite in sign to ``CROTA2``. + + It is downloaded again if an image is later than its last measurement, + unless it was downloaded less than :data:`_age_refresh` ago. + + Parameters + ---------- + latest + The start of the latest image which needs the roll, + in seconds since 1979 January 1. + directory + The directory to place the database in. + num_retry + The number of times to try to connect to the server. + """ + url = _url_coalign + "xrt_rollangle_db.geny" + path = _path_local(url, directory) + + refresh = False + if path.is_file() and _stale(path): + refresh = latest > _restore(path)["p0"][-1] + + path = _download_file( + url=url, + directory=directory, + overwrite=refresh, + num_retry=num_retry, + signature=_signature_genx, + kind="an IDL save file", + ) + + variables = _restore(path) + return variables["p0"].astype(float), variables["p1"].astype(float) + + +def _roll( + seconds: np.ndarray, + directory: pathlib.Path, + num_retry: int = 5, +) -> np.ndarray: + """ + The roll of each image, ``CROTA2``, in degrees, + as ``xrt_rollangle.pro`` in SolarSoft finds it: + from fits to the roll of XRT relative to SDO before the attitude anomaly + of 2021 December 27, + from a model of the roll during the anomaly, + and from the database of the roll after it, interpolated linearly. + + The models and their ranges are those of ``xrt_rollangle.pro``, + including its gap at exactly 2012 June 8 00:00, + when the roll is zero, + and its linear extrapolation beyond the last measurement of the + database. + + Parameters + ---------- + seconds + The start of each image, in seconds since 1979 January 1, + as ``anytim`` counts them. + directory + The directory to place the database of the roll in. + num_retry + The number of times to try to connect to the server. + """ + t = np.asarray(seconds, dtype=float) + t1, t2, t3, t4, t5, t6, t7 = _seconds(_times_roll) + + roll = np.zeros(t.shape) + + # Before 2012 March 10 + where = t <= t1 + roll[where] = 0.314067 + 0.0850355 * np.sin(1.99892e-07 * t[where] + 3.45231) + + # From 2012 March 10 to June 8, and from 2013 January 21 to the anomaly + where = ((t1 < t) & (t < t2)) | ((t3 < t) & (t <= t4)) + roll[where] = 0.314457 - 0.0712843 * np.sin(1.99474e-07 * t[where] + 0.445648) + + # From 2012 June 8 to 2013 January 21, with a variation over each day + where = (t2 < t) & (t <= t3) + roll[where] = ( + 0.314457 + - 0.0712843 * np.sin(1.99474e-07 * t[where] + 0.445648) + + 0.00833730 + + 0.0472532 * np.sin(7.27225e-05 * t[where] + 0.661514) + ) + + # During the anomaly + d = t - t4 + where = (t4 < t) & (t <= t5) + roll[where] = 5.2946095e-07 * d[where] ** 2 + 0.0038027661 * d[where] + 0.26089384 + where = (t5 < t) & (t <= t6) + roll[where] = 0.0015707323 * d[where] + 5.6717680 + where = (t6 < t) & (t <= t7) + roll[where] = 0.0062752709 * d[where] - 22.536388 + + # After the anomaly + where = t > t7 + if np.any(where): + times, angles = _database_roll(t[where].max(), directory, num_retry) + roll[where] = _interpol(t[where], times, angles) + + return -roll + + +def _listing( + database: Literal["aia", "ufss"], + latest: float, + directory: pathlib.Path, + num_retry: int = 5, +) -> list[str]: + """ + The names of the files of a co-alignment database, + without their extension, in order of time. + + The listing of the directory is kept in :data:`_name_listing`, + and the directory is listed again if an image is later than the start of + the last file, + unless it was listed less than :data:`_age_refresh` ago. + + Parameters + ---------- + database + The database, ``"aia"`` or ``"ufss"``. + latest + The start of the latest image to be looked up, + in seconds since 1979 January 1. + directory + The directory to place the listing in. + num_retry + The number of times to try to connect to the server. + + Raises + ------ + ValueError + If the listing of the directory has no files of the database. + """ + url = f"{_url_coalign}{_directories_coalign[database]}/" + path = _path_local(url + _name_listing, directory) + + if path.is_file(): + names = path.read_text().split() + if names and not (latest >= _start(names[-1:])[0] and _stale(path)): + return names + + names = sorted(set(_pattern_coalign.findall(_get(url, num_retry).text))) + if not names: + raise ValueError(f"The listing of {url} has no files of the database.") + + _write(path, "\n".join(names).encode(), overwrite=True) + + return names + + +def _lookup( + database: Literal["aia", "ufss"], + seconds: np.ndarray, + directory: pathlib.Path, + num_retry: int = 5, +) -> tuple[np.ndarray, np.ndarray, np.ndarray]: + """ + The center of the field of view of each image, in arcseconds, + and the type of its reference, + from the entry of a co-alignment database which began within + :data:`_tolerance` of it, as ``xrt_read_coaldb.pro`` finds it, + or NaN and -1 for an image the database has no entry for. + + Each file of the database covers the time from its start to the start of + the next one. + A file is downloaded again if an image later than its last entry has no + entry, since it may have been downloaded before the database reached the + image, + unless it was downloaded less than :data:`_age_refresh` ago. + An image before the last entry may have no entry of its own, + as an image the cross-correlation with AIA was not made for has none in + ``coaldb_c``. + + Parameters + ---------- + database + The database, ``"aia"`` or ``"ufss"``. + seconds + The start of each image, in seconds since 1979 January 1. + directory + The directory to place the files of the database in. + num_retry + The number of times to try to connect to the server. + """ + num = len(seconds) + xcen = np.full(num, np.nan) + ycen = np.full(num, np.nan) + ref_type = np.full(num, -1, dtype=int) + + names = _listing(database, float(np.max(seconds)), directory, num_retry) + starts = _start(names) + + # The files which cover the time within the tolerance of each image + first = np.searchsorted(starts, seconds - _tolerance, side="right") - 1 + last = np.searchsorted(starts, seconds + _tolerance, side="right") - 1 + + url_directory = f"{_url_coalign}{_directories_coalign[database]}/" + + for i in range(max(int(first.min()), 0), int(last.max()) + 1): + images = (first <= i) & (i <= last) + if not np.any(images): + continue + + url = f"{url_directory}{names[i]}.geny" + + refresh = False + while True: + path = _download_file( + url=url, + directory=directory, + overwrite=refresh, + num_retry=num_retry, + signature=_signature_genx, + kind="an IDL save file", + ) + entries = _restore(path)["p0"] + times = _seconds([t.decode() for t in entries["DATE_OBS"]]) + + for j in np.flatnonzero(images & (ref_type == -1)): + difference = np.abs(times - seconds[j]) + k = int(np.argmin(difference)) + if difference[k] <= _tolerance: + xcen[j] = entries["XCEN"][k] + ycen[j] = entries["YCEN"][k] + ref_type[j] = entries["REF_TYPE"][k] + + missing = np.any(images & (ref_type == -1) & (seconds > times.max())) + if refresh or not missing or not _stale(path): + break + refresh = True + + return xcen, ycen, ref_type + + +def _calibration(ref_type: np.ndarray) -> np.ndarray: + """ + The type of calibration of the pointing of each image, + ``CALIBRATION_TYPE`` of ``xrt_read_coaldb.pro``, + from the type of the reference of its entry in a co-alignment database. + + Parameters + ---------- + ref_type + The type of the reference, ``REF_TYPE``, or -1 for no entry. + """ + result = ref_type.copy() + result[ref_type >= 10] = 1 + result[ref_type == 1] = 2 + result[ref_type == 5] = 4 + return result + + +def _coalign( + time: Sequence[str], + filter_2: Sequence[int], + chip_sum: Sequence[int], + num_x: Sequence[int], + num_y: Sequence[int], + coalign: Literal["aia", "ufss"], + directory: None | pathlib.Path = None, + num_retry: int = 5, +) -> tuple[list[dict[str, float | str]], np.ndarray]: + """ + The keywords to change in the header of each Level 1 image to correct its + pointing as ``xrt_read_coaldb.pro`` in SolarSoft corrects it, + and the type of calibration of the pointing of each image. + + The center of the field of view of an image found in a co-alignment + database is placed at the center of the image, + with the plate scale of ``xrt_read_coaldb.pro``, + and every image is given the roll of ``xrt_rollangle.pro``. + An image the databases have no entry for keeps the center and the plate + scale of its header. + + Parameters + ---------- + time + The start of the exposure of each image, ``DATE_OBS``, + in UTC in the ISOT format. + filter_2 + The position of the second filter wheel of each image, ``EC_FW2``. + chip_sum + The number of pixels of the CCD summed along each axis into each + pixel of each image, ``CHIP_SUM``. + num_x + The number of columns of each image, ``NAXIS1``. + num_y + The number of rows of each image, ``NAXIS2``. + coalign + Which co-alignment database to find the center in: + ``"aia"`` for the cross-correlation with AIA, + or the UFSS database for an image it has no entry for, + as ``xrt_read_coaldb.pro`` does with ``/aia_cc``, + or ``"ufss"`` for the UFSS database alone, + as it does by default. + directory + The directory to place the databases in. + If :obj:`None` (the default), :data:`hinode.directory_default` is used. + num_retry + The number of times to try to connect to the server. + """ + if coalign not in _directories_coalign: + raise ValueError( + f"`coalign` must be one of {list(_directories_coalign)} or None, " + f"got {coalign!r}." + ) + + root: pathlib.Path = hinode.directory_default if directory is None else directory + + num = len(time) + seconds = _seconds(time) + + databases: list[Literal["aia", "ufss"]] = ["ufss"] + if coalign == "aia": + databases = ["aia", "ufss"] + + xcen = np.full(num, np.nan) + ycen = np.full(num, np.nan) + calibration = np.full(num, -1, dtype=int) + for database in databases: + missing = calibration == -1 + if not np.any(missing): + break + x, y, ref_type = _lookup(database, seconds[missing], root, num_retry) + xcen[missing] = x + ycen[missing] = y + calibration[missing] = _calibration(ref_type) + + roll = _roll(seconds, root, num_retry) + + keywords: list[dict[str, float | str]] = [] + for i in range(num): + result: dict[str, float | str] = { + "CROTA1": float(roll[i]), + "CROTA2": float(roll[i]), + } + if calibration[i] != -1: + scale = _scale_xray + if filter_2[i] == _position_gband: + scale = _scale_gband + scale = scale * chip_sum[i] + for axis, center, size in ( + ("1", xcen[i], num_x[i]), + ("2", ycen[i], num_y[i]), + ): + result[f"CRVAL{axis}"] = float(center) + # Integer division, as in ``xrt_read_coaldb.pro`` + result[f"CRPIX{axis}"] = size // 2 + 0.5 + result[f"CDELT{axis}"] = scale + result[f"CUNIT{axis}"] = "arcsec" + keywords.append(result) + + return keywords, calibration diff --git a/hinode/xrt/_coalign_test.py b/hinode/xrt/_coalign_test.py new file mode 100644 index 0000000..571f5fc --- /dev/null +++ b/hinode/xrt/_coalign_test.py @@ -0,0 +1,333 @@ +import os +import time +import pathlib +import datetime +import pytest +import requests +import numpy as np +import astropy.io.fits +import hinode +import hinode.xrt._coalign +from hinode.xrt._coalign import ( + _seconds, + _start, + _interpol, + _roll, + _listing, + _lookup, + _calibration, + _coalign, + _age_refresh, + _name_listing, + _url_coalign, +) + +_date_obs = "2019-09-30T18:08:00.577" +"""The start of the Al_poly image which the tutorial inverts for a DEM.""" + + +def _header(time_start: str, time_stop: str, filter: str) -> astropy.io.fits.Header: + """The header of the one file through a filter which began during a time range.""" + urls = hinode.xrt.urls(time_start, time_stop, filter=filter) + (path,) = hinode.xrt.download(urls).ndarray + return astropy.io.fits.getheader(path) + + +def _corrected( + header: astropy.io.fits.Header, + coalign: str, +) -> tuple[astropy.io.fits.Header, int]: + """ + A copy of a header with the keywords :func:`_coalign` gives it, + and the type of calibration of its pointing. + """ + (keywords,), calibration = _coalign( + time=[header["DATE_OBS"]], + filter_2=[int(header["EC_FW2"])], + chip_sum=[int(header["CHIP_SUM"])], + num_x=[int(header["NAXIS1"])], + num_y=[int(header["NAXIS2"])], + coalign=coalign, # type: ignore[arg-type] + ) + result = header.copy() + result.update(keywords) + return result, int(calibration[0]) + + +def _response(text: str) -> requests.Response: + """A response with the given text, made without connecting to a server.""" + response = requests.Response() + response.status_code = 200 + response._content = text.encode() + return response + + +@pytest.mark.parametrize( + argnames="isot,expected", + argvalues=[ + ("1979-01-01T00:00:00", 0), + ("1979-01-02T00:00:01.5", 86401.5), + ( + _date_obs, + ( + datetime.datetime(2019, 9, 30, 18, 8, 0, 577000) + - datetime.datetime(1979, 1, 1) + ).total_seconds(), + ), + ], +) +def test_seconds(isot: str, expected: float) -> None: + """Seconds since 1979 are counted without leap seconds, as ``anytim`` counts them.""" + assert _seconds([isot])[0] == expected + + +def test_start() -> None: + result = _start(["20190930_1503", "20190930_2100"]) + assert np.all(result == _seconds(["2019-09-30T15:03", "2019-09-30T21:00"])) + + +def test_interpol() -> None: + """Values beyond the samples are extrapolated from the first and the last interval.""" + xp = np.array([0.0, 1.0, 3.0]) + fp = np.array([0.0, 2.0, 3.0]) + x = np.array([-1.0, 0.5, 2.0, 5.0]) + assert np.allclose(_interpol(x, xp, fp), [-2.0, 1.0, 2.5, 4.0]) + + +@pytest.mark.parametrize( + argnames="isot,expected", + argvalues=[ + # Before 2012 March 10 + ("2007-06-01T12:43:00.331", -0.353642642498), + # From 2012 March 10 to June 8 + ("2012-06-01T12:28:28.295", -0.336584538221), + # From 2012 June 8 to 2013 January 21 + ("2012-09-01T12:00:02.749", -0.347388714552), + # From 2013 January 21 to the attitude anomaly of 2021 December 27 + ("2015-07-01T12:03:06.961", -0.368229240179), + (_date_obs, -0.35857245326), + # After the anomaly, from the database + ("2025-06-01T12:35:26.993", 1.22919750214), + ("2026-09-01T12:30:03.762", -0.392326414585), + ], +) +def test_roll(isot: str, expected: float) -> None: + """ + The roll is the one ``xrt_prep`` wrote to the headers of these files, + which it took from ``xrt_rollangle.pro``. + + The headers are rounded to single precision, + and between 2012 June and 2013 January differ by up to 5e-5 degrees. + """ + result = _roll(_seconds([isot]), hinode.directory_default) + assert np.isclose(result[0], expected, rtol=0, atol=1e-4) + + +def test_roll_anomaly() -> None: + """ + During the attitude anomaly the roll follows the model of + ``xrt_rollangle.pro``, which ends at -23 degrees, + rather than the roll of about -0.26 degrees in the headers of the files + ``xrt_prep`` made before the model was written. + """ + isot = ["2021-12-27T13:00:13.873", "2021-12-27T15:01:41.946"] + result = _roll(_seconds(isot), hinode.directory_default) + assert np.allclose(result, [-0.31375, -23.2853], rtol=0, atol=1e-4) + + +@pytest.mark.parametrize( + argnames="ref_type,expected", + argvalues=[ + (-1, -1), + (1, 2), + (3, 3), + (5, 4), + (6, 6), + (10, 1), + (12, 1), + ], +) +def test_calibration(ref_type: int, expected: int) -> None: + assert _calibration(np.array([ref_type]))[0] == expected + + +def test_listing(tmp_path: pathlib.Path, monkeypatch: pytest.MonkeyPatch) -> None: + """ + The directory is listed again only if an image is later than the start of + the last file and the listing is older than a day. + """ + calls = [] + + def get(url: str, num_retry: int = 5) -> requests.Response: + calls.append(url) + return _response(' ') + + monkeypatch.setattr(hinode.xrt._coalign, "_get", get) + + early = _seconds(["2019-09-30T12:00"])[0] + late = _seconds(["2019-10-01T00:00"])[0] + + assert _listing("aia", early, tmp_path) == ["20190930_0631", "20190930_1503"] + assert calls == [f"{_url_coalign}coaldb_c/"] + + # Kept, since the image is in a file before the last one + _listing("aia", early, tmp_path) + # Kept, since the listing is new + _listing("aia", late, tmp_path) + assert len(calls) == 1 + + path = tmp_path / "sdb/hinode/xrt/xrt_msu_coalign/coaldb_c" / _name_listing + old = time.time() - 2 * _age_refresh + os.utime(path, (old, old)) + _listing("aia", early, tmp_path) + assert len(calls) == 1 + _listing("aia", late, tmp_path) + assert len(calls) == 2 + + +def test_listing_empty(tmp_path: pathlib.Path, monkeypatch: pytest.MonkeyPatch) -> None: + monkeypatch.setattr(hinode.xrt._coalign, "_get", lambda url, n=5: _response("")) + with pytest.raises(ValueError, match="has no files of the database"): + _listing("ufss", 0, tmp_path) + + +@pytest.mark.parametrize( + argnames="database,expected", + argvalues=[ + ("aia", (5.99925, -16.8217, 10)), + ("ufss", (6.334633, -17.498896, 3)), + ], +) +def test_lookup(database: str, expected: tuple[float, float, int]) -> None: + """ + An image within the tolerance of an entry is given it, + and one further from every entry is given none. + """ + isot = [_date_obs, "2019-09-30T18:08:00.667", "2019-09-30T18:08:00.877"] + xcen, ycen, ref_type = _lookup(database, _seconds(isot), hinode.directory_default) + assert np.allclose(xcen[:2], expected[0]) + assert np.allclose(ycen[:2], expected[1]) + assert np.all(ref_type[:2] == expected[2]) + assert np.isnan(xcen[2]) and np.isnan(ycen[2]) and ref_type[2] == -1 + + +def test_lookup_before() -> None: + """An image before the first file of a database has no entry.""" + xcen, ycen, ref_type = _lookup( + "ufss", _seconds(["2006-01-01T00:00"]), hinode.directory_default + ) + assert np.isnan(xcen[0]) and ref_type[0] == -1 + + +@pytest.mark.parametrize( + argnames="coalign,expected", + argvalues=[ + ("aia", (5.99925, -16.8217, 1)), + ("ufss", (6.334633, -17.498896, 3)), + ], +) +def test_coalign(coalign: str, expected: tuple[float, float, int]) -> None: + header = _header("2019-09-30T18:07:59", "2019-09-30T18:08:01", "Al_poly") + assert header["DATE_OBS"] == _date_obs + + result, calibration = _corrected(header, coalign) + + assert calibration == expected[2] + assert np.isclose(result["CRVAL1"], expected[0]) + assert np.isclose(result["CRVAL2"], expected[1]) + assert result["CRPIX1"] == result["CRPIX2"] == 192.5 + assert result["CDELT1"] == result["CDELT2"] == 1.0286 + assert np.isclose(result["CROTA2"], header["CROTA2"], rtol=0, atol=1e-6) + assert result["CROTA1"] == result["CROTA2"] + + # The header given is not changed + assert header["CRVAL1"] == 6.33463287354 + + +def test_coalign_gband() -> None: + """The images through G-band have a plate scale of their own.""" + header = _header("2019-09-30T18:09:32", "2019-09-30T18:09:33", "Gband") + result, calibration = _corrected(header, "aia") + assert calibration == 1 + assert result["CDELT1"] == result["CDELT2"] == 1.0302 * header["CHIP_SUM"] + + +def test_coalign_no_entry() -> None: + """ + An image no database has an entry for keeps the center and the plate + scale of its header, but is given the roll of ``xrt_rollangle.pro``. + """ + header = _header("2019-09-30T18:07:59", "2019-09-30T18:08:01", "Al_poly") + header["DATE_OBS"] = "2006-01-01T00:00:00.000" + header["CDELT1"] = 1.03 + + result, calibration = _corrected(header, "aia") + + assert calibration == -1 + for key in ["CRVAL1", "CRVAL2", "CRPIX1", "CRPIX2", "CDELT1", "CDELT2"]: + assert result[key] == header[key] + expected = _roll(_seconds(["2006-01-01T00:00:00.000"]), hinode.directory_default) + assert result["CROTA2"] == expected[0] + + +def test_coalign_invalid() -> None: + with pytest.raises(ValueError, match="`coalign` must be one of"): + _coalign([], [], [], [], [], "aia_cc") # type: ignore[arg-type] + + +@pytest.mark.parametrize( + argnames="isot,stale,num_downloads", + argvalues=[ + # After the last entry of the file, which was downloaded long ago + ("2019-09-30T16:00:00.000", True, 2), + # After the last entry, but the file is new + ("2019-09-30T16:00:00.000", False, 1), + # Before the last entry, with no entry of its own + ("2019-09-30T15:30:00.000", True, 1), + ], +) +def test_lookup_refresh( + monkeypatch: pytest.MonkeyPatch, + isot: str, + stale: bool, + num_downloads: int, +) -> None: + """ + A file is downloaded again only if an image later than its last entry has + no entry, and the file is old. + """ + downloads = [] + + def download_file( + url: str, directory: pathlib.Path, overwrite: bool, **kwargs: object + ) -> pathlib.Path: + downloads.append(overwrite) + return directory / "file.geny" + + def restore(path: pathlib.Path) -> dict: + times = ["2019-09-30T15:10:00.000", "2019-09-30T15:50:00.000"] + if len(downloads) > 1: + times.append("2019-09-30T16:00:00.000") + num = len(times) + entries = np.rec.fromarrays( + [ + np.array([t.encode() for t in times]), + np.ones(num), + np.ones(num), + np.full(num, 10), + ], + names=["DATE_OBS", "XCEN", "YCEN", "REF_TYPE"], + ) + return {"p0": entries} + + monkeypatch.setattr( + hinode.xrt._coalign, "_listing", lambda *args: ["20190930_1503"] + ) + monkeypatch.setattr(hinode.xrt._coalign, "_download_file", download_file) + monkeypatch.setattr(hinode.xrt._coalign, "_restore", restore) + monkeypatch.setattr(hinode.xrt._coalign, "_stale", lambda path: stale) + + _, _, ref_type = _lookup("aia", _seconds([isot]), pathlib.Path("cache")) + + assert downloads == [False, True][:num_downloads] + assert ref_type[0] == (10 if num_downloads == 2 else -1) diff --git a/hinode/xrt/_data.py b/hinode/xrt/_data.py index ef69913..6ea5ff8 100644 --- a/hinode/xrt/_data.py +++ b/hinode/xrt/_data.py @@ -472,6 +472,120 @@ def _replace( os.replace(source, destination) +def _path_local( + url: str, + directory: pathlib.Path, +) -> pathlib.Path: + """ + Where a file on a server is placed under a directory, + at the same path as on the server. + + Parameters + ---------- + url + The URL of the file. + directory + The directory to place the file in. + + Raises + ------ + ValueError + If the URL is not that of a file. + """ + relative = "/".join(url.split("/")[3:]) + if not relative or relative.endswith("/"): + raise ValueError(f"{url} is not the URL of a file.") + return directory / relative + + +def _write( + path: pathlib.Path, + content: bytes, + overwrite: bool, +) -> None: + """ + Write a file, so that an interrupted write never leaves it half written. + + Parameters + ---------- + path + The file to write. + content + What to write to it. + overwrite + Whether the new file must replace one which another process has open, + rather than keeping that one. + """ + path.parent.mkdir(parents=True, exist_ok=True) + + # Written next to its final place, under a name no other download + # uses, and then moved, so an interrupted download never looks + # finished. + # It is made with :func:`open` rather than :mod:`tempfile`, + # so that it has the permissions the umask gives a new file. + part = path.with_name(f"{path.name}.{uuid.uuid4().hex}.part") + try: + with open(part, "xb") as file: + file.write(content) + try: + _replace(part, path) + except PermissionError: + # Another process may have downloaded the file first and + # still have it open, in which case its copy is kept, + # unless a new copy was asked for. + if overwrite or not path.is_file(): + raise + finally: + part.unlink(missing_ok=True) + + +def _download_file( + url: str, + directory: pathlib.Path, + overwrite: bool = False, + num_retry: int = 5, + signature: bytes = _signature_fits, + kind: str = "a FITS file", +) -> pathlib.Path: + """ + Download a file to a directory, unless it has been downloaded already, + and check that it starts the way a file of its kind does. + + Parameters + ---------- + url + The URL to download. + directory + The directory to place the file in, + under the same path as on the server. + overwrite + Boolean flag controlling whether to download the file if it is + already in `directory`. + num_retry + The number of times to try to connect to the server. + signature + The first bytes of every file of its kind. + kind + The kind of file, for the error message. + + Raises + ------ + ValueError + If the file does not start with `signature`. + """ + path = _path_local(url, directory) + + if overwrite or not path.is_file(): + content = _get(url, num_retry).content + + if not content.startswith(signature): + raise ValueError(f"{url} is not {kind}, it starts with {content[:40]!r}.") + + _write(path, content, overwrite) + + return path + + def download( urls: na.AbstractScalarArray, directory: None | pathlib.Path = None, @@ -517,48 +631,12 @@ def download( hinode.xrt.download(urls) """ - if directory is None: - directory = hinode.directory_default + root: pathlib.Path = hinode.directory_default if directory is None else directory ndarray = np.asarray(urls.ndarray) def get(url: str) -> str: - relative = "/".join(url.split("/")[3:]) - if not relative or relative.endswith("/"): - raise ValueError(f"{url} is not the URL of a file.") - path = directory / relative - - if overwrite or not path.is_file(): - content = _get(url, num_retry).content - - if not content.startswith(_signature_fits): - raise ValueError( - f"{url} is not a FITS file, it starts with {content[:40]!r}." - ) - - path.parent.mkdir(parents=True, exist_ok=True) - - # Written next to its final place, under a name no other download - # uses, and then moved, so an interrupted download never looks - # finished. - # It is made with :func:`open` rather than :mod:`tempfile`, - # so that it has the permissions the umask gives a new file. - part = path.with_name(f"{path.name}.{uuid.uuid4().hex}.part") - try: - with open(part, "xb") as file: - file.write(content) - try: - _replace(part, path) - except PermissionError: - # Another process may have downloaded the file first and - # still have it open, in which case its copy is kept, - # unless a new copy was asked for. - if overwrite or not path.is_file(): - raise - finally: - part.unlink(missing_ok=True) - - return str(path) + return str(_download_file(url, root, overwrite, num_retry)) # Each URL is downloaded once, even if it is given more than once. unique = list(dict.fromkeys(str(url) for url in ndarray.flat)) diff --git a/hinode/xrt/_filtergrams.py b/hinode/xrt/_filtergrams.py index b26d051..10ee754 100644 --- a/hinode/xrt/_filtergrams.py +++ b/hinode/xrt/_filtergrams.py @@ -1,4 +1,4 @@ -from typing import Sequence +from typing import Sequence, Literal import typing import os import re @@ -14,6 +14,7 @@ from ._leak import _leak, _history_leak from ._uncertainty import _quality, _error_jpeg, _error_dark from ._vignetting import _vignetting, _error_vignetting +from ._coalign import _coalign __all__ = [ "Filtergram", @@ -119,6 +120,11 @@ class Filtergram( Each image is one frame along :attr:`axis_time`, and its coordinates come from its own header, since the pointing changes from one image to the next. + By default the pointing is corrected with the co-alignment database of + SolarSoft from the cross-correlation of XRT with AIA + :cite:p:`Yoshimura2015`, + so that the images line up with those of AIA, + as :attr:`coalignment` records. The time of each image, ``inputs.time``, is the start of its exposure, ``DATE_OBS``, as in the other sun-data packages, @@ -272,6 +278,20 @@ class Filtergram( named as in :func:`hinode.xrt.urls`. """ + coalignment: None | na.ScalarArray = None + """ + How the pointing of each image was found, + numbered as ``CALIBRATION_TYPE`` of ``xrt_read_coaldb.pro`` in SolarSoft: + 1 for the cross-correlation with AIA, + 2 for a fit to the limb of the Sun in G-band, + 3 for the Ultra Fine Sun Sensors (UFSS), + 4 for a fit to the limb in X-rays, + 6 for another method, + and -1 for an image no database had an entry for, + which keeps the pointing of its header, + or :obj:`None` if the pointing was not corrected. + """ + axis_time: str = "time" """The logical axis corresponding to changes in time.""" @@ -295,6 +315,7 @@ def from_time_range( num_retry: int = 5, leak: bool = False, uncertainty: bool = False, + coalign: None | Literal["aia", "ufss"] = "aia", ) -> "Filtergram": """ Download the Level 1 images which began during a given time range @@ -330,6 +351,17 @@ def from_time_range( uncertainty Whether to load :attr:`vignetting` and :attr:`uncertainty_map`, which :meth:`uncertainty` needs. + coalign + Which co-alignment database of SolarSoft to correct the pointing + of the images with, as ``xrt_read_coaldb.pro`` does, + or :obj:`None` to keep the pointing of the headers. + ``"aia"`` (the default) uses the cross-correlation of each image + with AIA 335 Å :cite:p:`Yoshimura2015`, + and the Ultra Fine Sun Sensors (UFSS) of Hinode for an image it + has no entry for, + as ``xrt_read_coaldb.pro`` does with ``/aia_cc``. + ``"ufss"`` uses the UFSS alone, as it does by default. + See :attr:`coalignment`. """ urls = hinode.xrt.urls( time_start=time_start, @@ -358,7 +390,9 @@ def from_time_range( axis_detector_y=axis_detector_y, leak=leak, uncertainty=uncertainty, + coalign=coalign, directory=directory, + num_retry=num_retry, ) @classmethod @@ -370,7 +404,9 @@ def from_fits( axis_detector_y: str = "detector_y", leak: bool = False, uncertainty: bool = False, + coalign: None | Literal["aia", "ufss"] = "aia", directory: None | pathlib.Path = None, + num_retry: int = 5, ) -> "Filtergram": """ Load one or more Level 1 XRT files taken through the same filter. @@ -399,9 +435,24 @@ def from_fits( which :meth:`uncertainty` needs. They are not loaded by default, since each takes as much memory as the images. + coalign + Which co-alignment database of SolarSoft to correct the pointing + of the images with, as ``xrt_read_coaldb.pro`` does, + or :obj:`None` to keep the pointing of the headers. + ``"aia"`` (the default) uses the cross-correlation of each image + with AIA 335 Å :cite:p:`Yoshimura2015`, + and the Ultra Fine Sun Sensors (UFSS) of Hinode for an image it + has no entry for, + as ``xrt_read_coaldb.pro`` does with ``/aia_cc``. + ``"ufss"`` uses the UFSS alone, as it does by default. + See :attr:`coalignment`. directory - The directory to place the downloaded images of the leak in. + The directory to place the downloaded images of the leak and the + co-alignment databases in. If :obj:`None` (the default), :data:`hinode.directory_default` is used. + num_retry + The number of times to try to connect to the server + for the co-alignment databases. Raises ------ @@ -412,8 +463,9 @@ def from_fits( ``xrt_prep`` set the saturated pixels to, if `leak` is :obj:`True` and SolarSoft has no image of the leak for the filter and the time of an image, - or if `uncertainty` is :obj:`True` and the compression of an - image is not known. + if `uncertainty` is :obj:`True` and the compression of an + image is not known, + or if `coalign` is not ``"aia"``, ``"ufss"``, or :obj:`None`. Notes ----- @@ -423,6 +475,22 @@ def from_fits( The headers give the roll of each image as ``CROTA2`` rather than as a ``PC`` matrix, and it is converted to one the way :cite:t:`Calabretta2002` describe for that older keyword. + + ``xrt_prep`` has already corrected the pointing of the Level 1 files + with the UFSS, so the center of an image moves only with ``"aia"``, + the default. + Either database also gives each image the roll ``xrt_rollangle.pro`` + gives it now. + Since the attitude anomaly of 2021 December 27, + the roll is interpolated in a database which grows as the roll is + measured, + and some of the files ``xrt_prep`` made before it existed have a roll + in their header which is more than a degree from it, + and up to 23 degrees during the anomaly. + Like ``xrt_read_coaldb.pro``, + the correction places the center the database gives at pixel + ``NAXIS // 2 + 0.5``, counted from one, + with the plate scale of the X-ray or the G-band images. """ if isinstance(path, na.AbstractScalarArray): path = [str(p) for p in np.ravel(np.asarray(path.ndarray))] @@ -443,6 +511,22 @@ def from_fits( _check_normalized(header, p) levels = [_level_saturation(header, p) for p, header in zip(path, headers)] + coalignment = None + if coalign is not None: + keywords, calibration = _coalign( + time=[h["DATE_OBS"] for h in headers], + filter_2=[int(h["EC_FW2"]) for h in headers], + chip_sum=[int(h["CHIP_SUM"]) for h in headers], + num_x=[int(h["NAXIS1"]) for h in headers], + num_y=[int(h["NAXIS2"]) for h in headers], + coalign=coalign, + directory=directory, + num_retry=num_retry, + ) + for header, keywords_header in zip(headers, keywords): + header.update(keywords_header) + coalignment = na.ScalarArray(calibration, axes=axis_time) + shape = { axis_time: len(headers), axis_detector_y: max(int(h["NAXIS2"]) for h in headers), @@ -666,6 +750,7 @@ def angle(key: str, key_unit: str) -> na.ScalarArray: uncertainty_map=uncertainty_map, leak=leaks, filter=filters[0], + coalignment=coalignment, axis_time=axis_time, axis_detector_x=axis_detector_x, axis_detector_y=axis_detector_y, diff --git a/hinode/xrt/_filtergrams_test.py b/hinode/xrt/_filtergrams_test.py index 7286a2f..514e5d9 100644 --- a/hinode/xrt/_filtergrams_test.py +++ b/hinode/xrt/_filtergrams_test.py @@ -10,6 +10,7 @@ import named_arrays as na import hinode from hinode.xrt._vignetting_test import _vignetting_numpy +from hinode.xrt._coalign_test import _corrected _time_start = astropy.time.Time("2019-09-30T18:08:00") _time_stop = astropy.time.Time("2019-09-30T18:09:00") @@ -118,6 +119,12 @@ def _num_saturated(header: astropy.io.fits.Header) -> int: axis_detector_y="y", leak=True, uncertainty=True, + coalign="ufss", + ), + hinode.xrt.Filtergram.from_time_range( + time_start=_time_start, + time_stop=_time_stop, + coalign=None, ), ], ) @@ -166,6 +173,21 @@ def test_saturated(self, array: hinode.xrt.Filtergram) -> None: def test_filter(self, array: hinode.xrt.Filtergram) -> None: assert array.filter == "Al_poly" + def test_coalignment(self, array: hinode.xrt.Filtergram) -> None: + coalignment = array.coalignment + if coalignment is None: + return + assert isinstance(coalignment, na.ScalarArray) + assert coalignment.shape == {array.axis_time: array.shape[array.axis_time]} + + # The cross-correlation with AIA has every one of these images, + # and the UFSS database alone gives the pointing of the headers. + assert np.all(coalignment.ndarray == 1) or np.all(coalignment.ndarray == 3) + + image = array[{array.axis_time: 0}] + assert image.coalignment is not None + assert image.coalignment.shape == {} + def test_vignetting(self, array: hinode.xrt.Filtergram) -> None: vignetting = array.vignetting if vignetting is None: @@ -563,17 +585,21 @@ def test_from_fits_history( hinode.xrt.Filtergram.from_fits(file) -def test_inputs_against_astropy_wcs() -> None: +@pytest.mark.parametrize("coalign", [None, "aia"]) +def test_inputs_against_astropy_wcs(coalign: None | str) -> None: """ The coordinates of each pixel must be the ones :mod:`astropy.wcs` - makes of the header, which gives the roll as ``CROTA2``. + makes of the header, which gives the roll as ``CROTA2``, + once the pointing of the header is corrected. """ path = _path("2019-09-30T18:08:30", "2019-09-30T18:08:40") - result = hinode.xrt.Filtergram.from_fits(path) + result = hinode.xrt.Filtergram.from_fits(path, coalign=coalign) header = astropy.io.fits.getheader(path) assert header["CROTA2"] != 0 + if coalign is not None: + header, _ = _corrected(header, coalign) # Without the projection, # which :class:`named_arrays.AbstractWcsVector` does not apply. diff --git a/hinode/xrt/_xrt.py b/hinode/xrt/_xrt.py index b9ce3d9..e6e3b9a 100644 --- a/hinode/xrt/_xrt.py +++ b/hinode/xrt/_xrt.py @@ -1,3 +1,4 @@ +from typing import Literal import pathlib import astropy.time from ._filtergrams import Filtergram @@ -19,6 +20,7 @@ def open( num_retry: int = 5, leak: bool = False, uncertainty: bool = False, + coalign: None | Literal["aia", "ufss"] = "aia", ) -> Filtergram: """ Download the Level 1 XRT images which began during a given time range @@ -58,6 +60,13 @@ def open( Whether to load :attr:`~hinode.xrt.Filtergram.vignetting` and :attr:`~hinode.xrt.Filtergram.uncertainty_map`, which :meth:`~hinode.xrt.Filtergram.uncertainty` needs. + coalign + Which co-alignment database of SolarSoft to correct the pointing of + the images with, + ``"aia"`` (the default) for the cross-correlation with AIA, + ``"ufss"`` for the Ultra Fine Sun Sensors, + or :obj:`None` to keep the pointing of the headers, + as in :meth:`~hinode.xrt.Filtergram.from_fits`. """ return Filtergram.from_time_range( time_start=time_start, @@ -71,4 +80,5 @@ def open( num_retry=num_retry, leak=leak, uncertainty=uncertainty, + coalign=coalign, ) diff --git a/hinode/xrt/_xrt_test.py b/hinode/xrt/_xrt_test.py index e093bfb..c69f565 100644 --- a/hinode/xrt/_xrt_test.py +++ b/hinode/xrt/_xrt_test.py @@ -13,6 +13,17 @@ def test_open() -> None: assert result.filter == "Al_poly" assert result.outputs.unit == u.DN / u.s assert result.leak is None + assert result.coalignment is not None + assert result.coalignment.ndarray.tolist() == [1, 1] + + +def test_open_without_coalign() -> None: + result = hinode.xrt.open( + time_start="2019-09-30T18:08:30", + time_stop="2019-09-30T18:09:00", + coalign=None, + ) + assert result.coalignment is None def test_open_leak() -> None: diff --git a/pyproject.toml b/pyproject.toml index 5c42673..1a32d68 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -20,6 +20,7 @@ dependencies = [ "urllib3>=2", "joblib", "astropy", + "scipy", "named-arrays~=2.14", "xrtpy~=0.5.1", ] From 400a2f25d4acd9b3fe8ffbbbefb8a3060f1cb793 Mon Sep 17 00:00:00 2001 From: Roy Smart Date: Sat, 10 Oct 2026 11:14:02 -0600 Subject: [PATCH 2/3] Test the refresh of the roll database and the files skipped between images. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_011VGPHLGJi8acbuLYSfBHTk --- hinode/xrt/_coalign_test.py | 86 +++++++++++++++++++++++++++++++++++++ 1 file changed, 86 insertions(+) diff --git a/hinode/xrt/_coalign_test.py b/hinode/xrt/_coalign_test.py index 571f5fc..0724d0e 100644 --- a/hinode/xrt/_coalign_test.py +++ b/hinode/xrt/_coalign_test.py @@ -13,6 +13,7 @@ _start, _interpol, _roll, + _database_roll, _listing, _lookup, _calibration, @@ -21,6 +22,7 @@ _name_listing, _url_coalign, ) +from hinode.xrt._data import _path_local _date_obs = "2019-09-30T18:08:00.577" """The start of the Al_poly image which the tutorial inverts for a DEM.""" @@ -331,3 +333,87 @@ def restore(path: pathlib.Path) -> dict: assert downloads == [False, True][:num_downloads] assert ref_type[0] == (10 if num_downloads == 2 else -1) + + +@pytest.mark.parametrize( + argnames="latest,stale,overwrite", + argvalues=[ + # After the last measurement, in a database downloaded long ago + (20.0, True, True), + # Within the database + (5.0, True, False), + # After the last measurement, in a new database + (20.0, False, False), + ], +) +def test_database_roll_refresh( + tmp_path: pathlib.Path, + monkeypatch: pytest.MonkeyPatch, + latest: float, + stale: bool, + overwrite: bool, +) -> None: + """ + The database of the roll is downloaded again only for an image after its + last measurement, and only if it is old. + """ + path = _path_local(_url_coalign + "xrt_rollangle_db.geny", tmp_path) + path.parent.mkdir(parents=True) + path.touch() + + downloads = [] + + def download_file( + url: str, directory: pathlib.Path, overwrite: bool, **kwargs: object + ) -> pathlib.Path: + downloads.append(overwrite) + return path + + def restore(path: pathlib.Path) -> dict: + return {"p0": np.array([0.0, 10.0]), "p1": np.array([1.0, 2.0])} + + monkeypatch.setattr(hinode.xrt._coalign, "_download_file", download_file) + monkeypatch.setattr(hinode.xrt._coalign, "_restore", restore) + monkeypatch.setattr(hinode.xrt._coalign, "_stale", lambda path: stale) + + times, angles = _database_roll(latest, tmp_path) + + assert downloads == [overwrite] + assert np.all(times == [0, 10]) and np.all(angles == [1, 2]) + + +def test_lookup_gap(monkeypatch: pytest.MonkeyPatch) -> None: + """ + Only the files which cover an image are downloaded, + and not the ones between them. + """ + isot = ["2019-09-30T00:30:00.000", "2019-09-30T02:30:00.000"] + downloads = [] + + def download_file(url: str, directory: pathlib.Path, **kwargs: object) -> str: + downloads.append(url.split("/")[-1]) + return url + + def restore(path: str) -> dict: + times = [t for t in isot if t[11:13] == path.split("_")[-1][:2]] + num = len(times) + entries = np.rec.fromarrays( + [ + np.array([t.encode() for t in times]), + np.ones(num), + np.ones(num), + np.full(num, 3), + ], + names=["DATE_OBS", "XCEN", "YCEN", "REF_TYPE"], + ) + return {"p0": entries} + + names = ["20190930_0000", "20190930_0100", "20190930_0200"] + monkeypatch.setattr(hinode.xrt._coalign, "_listing", lambda *args: names) + monkeypatch.setattr(hinode.xrt._coalign, "_download_file", download_file) + monkeypatch.setattr(hinode.xrt._coalign, "_restore", restore) + + _, _, ref_type = _lookup("ufss", _seconds(isot), pathlib.Path("cache")) + + assert downloads == ["20190930_0000.geny", "20190930_0200.geny"] + assert ref_type.tolist() == [3, 3] From dfb4f18e840be9d7b25c5b1e0a46de673495f2d1 Mon Sep 17 00:00:00 2001 From: Roy Smart Date: Sat, 10 Oct 2026 15:39:12 -0600 Subject: [PATCH 3/3] Download a database file again only when the listing shows it changed, use cached copies offline, and warn when the roll is extrapolated. The listing now keeps when each file last changed, so a complete file before a gap is no longer downloaded every day, and files the server rewrites are noticed, with the listings checked again after a month for any image. A download which fails keeps the copy from before with a warning, and the roll past the last measurement is extrapolated as in xrt_rollangle.pro but with a warning. A wrong coalign is caught before anything is downloaded. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_011VGPHLGJi8acbuLYSfBHTk --- hinode/xrt/_coalign.py | 357 ++++++++++++++++++++++++-------- hinode/xrt/_coalign_test.py | 342 +++++++++++++++++++++++------- hinode/xrt/_filtergrams.py | 38 +++- hinode/xrt/_filtergrams_test.py | 66 +++++- 4 files changed, 630 insertions(+), 173 deletions(-) diff --git a/hinode/xrt/_coalign.py b/hinode/xrt/_coalign.py index 8b49887..ce19127 100644 --- a/hinode/xrt/_coalign.py +++ b/hinode/xrt/_coalign.py @@ -2,6 +2,8 @@ import re import time import pathlib +import warnings +import requests import numpy as np import hinode from ._data import _get, _path_local, _write, _download_file @@ -24,10 +26,13 @@ and ``coaldb_u``, from the Ultra Fine Sun Sensors (UFSS) of Hinode. """ -_pattern_coalign = re.compile(r'href="(\d{8}_\d{4})\.geny"') +_pattern_coalign = re.compile( + r'href="(\d{8}_\d{4})\.geny"(?:[^\n]*?(\d{4}-\d{2}-\d{2} \d{2}:\d{2}))?' +) """ The name of a file of a co-alignment database in a listing of its directory, -which is the start of the time it covers, as in ``20190930_1503.geny``. +which is the start of the time it covers, as in ``20190930_1503.geny``, +and when the file was last changed, if the listing shows it on the same line. """ _name_listing = "listing.txt" @@ -45,9 +50,27 @@ _age_refresh = 86400 """ -How old, in seconds, a downloaded listing or file of a database must be -before it is downloaded again to look for images it may not have had yet, -since the databases are updated about once a week. +How old, in seconds, a listing of a database or the database of the roll +must be before it is downloaded again to look for images it may not have had +yet, since the databases are updated about once a week. +""" + +_age_check = 30 * 86400 +""" +How old, in seconds, a listing of a database or the database of the roll +must be before it is downloaded again whatever the images, +since the server sometimes changes files it has had for years, +as it did for more than a thousand files of each database on one day in +2015 and again in 2020. +""" + +_errors_network = ( + ConnectionError, + requests.exceptions.RequestException, +) +""" +The errors of a download which mean the server could not be reached or +refused it, so that a copy downloaded before can be used instead. """ _tolerance = 0.1 @@ -116,7 +139,24 @@ def _start(names: Sequence[str]) -> np.ndarray: return _seconds(isot) -def _stale(path: pathlib.Path) -> bool: +def _isot(seconds: float) -> str: + """ + A time in seconds since 1979 January 1 in the ISOT format, + the inverse of :func:`_seconds`. + + Parameters + ---------- + seconds + The time, counted as ``anytim`` counts it. + """ + microseconds = np.timedelta64(round(seconds * 1e6), "us") + return str(np.datetime_as_string(_epoch + microseconds, unit="s")) + + +def _stale( + path: pathlib.Path, + age: float = _age_refresh, +) -> bool: """ Whether a downloaded file is old enough to be downloaded again. @@ -124,8 +164,55 @@ def _stale(path: pathlib.Path) -> bool: ---------- path The file. + age + How old, in seconds, the file must be. + """ + return time.time() - path.stat().st_mtime > age + + +def _download_database( + url: str, + directory: pathlib.Path, + overwrite: bool = False, + num_retry: int = 5, +) -> pathlib.Path: """ - return time.time() - path.stat().st_mtime > _age_refresh + Download a file of a database to a directory, + unless it has been downloaded already, + and keep the copy downloaded before, with a warning, + if it cannot be downloaded again. + + Parameters + ---------- + url + The URL of the file. + directory + The directory to place the file in, + under the same path as on the server. + overwrite + Whether to download the file again if it is already in `directory`. + num_retry + The number of times to try to connect to the server. + """ + path = _path_local(url, directory) + try: + return _download_file( + url=url, + directory=directory, + overwrite=overwrite, + num_retry=num_retry, + signature=_signature_genx, + kind="an IDL save file", + ) + except _errors_network as e: + if not path.is_file(): + raise + warnings.warn( + f"Could not download {url} again, " + f"so the copy downloaded before is used: {e!r}", + stacklevel=2, + ) + return path def _restore(path: pathlib.Path) -> dict: @@ -182,8 +269,9 @@ def _database_roll( the time of each measurement, in seconds since 1979 January 1, and the roll, in degrees, opposite in sign to ``CROTA2``. - It is downloaded again if an image is later than its last measurement, - unless it was downloaded less than :data:`_age_refresh` ago. + It is downloaded again if an image is later than its last measurement + and it was downloaded more than :data:`_age_refresh` ago, + or if it was downloaded more than :data:`_age_check` ago. Parameters ---------- @@ -198,18 +286,14 @@ def _database_roll( url = _url_coalign + "xrt_rollangle_db.geny" path = _path_local(url, directory) - refresh = False - if path.is_file() and _stale(path): - refresh = latest > _restore(path)["p0"][-1] - - path = _download_file( - url=url, - directory=directory, - overwrite=refresh, - num_retry=num_retry, - signature=_signature_genx, - kind="an IDL save file", - ) + downloaded = path.is_file() + if downloaded: + variables = _restore(path) + after = latest > variables["p0"][-1] + if not ((after and _stale(path)) or _stale(path, _age_check)): + return variables["p0"].astype(float), variables["p1"].astype(float) + + path = _download_database(url, directory, downloaded, num_retry) variables = _restore(path) return variables["p0"].astype(float), variables["p1"].astype(float) @@ -231,8 +315,9 @@ def _roll( The models and their ranges are those of ``xrt_rollangle.pro``, including its gap at exactly 2012 June 8 00:00, when the roll is zero, - and its linear extrapolation beyond the last measurement of the - database. + and its linear extrapolation beyond the first and the last measurement + of the database, which is done with a warning, + since the measurements reach the server days after the images. Parameters ---------- @@ -280,6 +365,16 @@ def _roll( if np.any(where): times, angles = _database_roll(t[where].max(), directory, num_retry) roll[where] = _interpol(t[where], times, angles) + outside = (t[where] < times[0]) | (t[where] > times[-1]) + if np.any(outside): + warnings.warn( + f"{np.count_nonzero(outside)} of {t.size} images began " + f"outside the measurements of the roll, " + f"from {_isot(times[0])} to {_isot(times[-1])}, " + f"so their roll is extrapolated, " + f"as xrt_rollangle.pro extrapolates it, and may be wrong.", + stacklevel=3, + ) return -roll @@ -289,15 +384,23 @@ def _listing( latest: float, directory: pathlib.Path, num_retry: int = 5, -) -> list[str]: +) -> dict[str, float]: """ The names of the files of a co-alignment database, - without their extension, in order of time. + without their extension, in order of time, + and when the listing first showed the version of each file the server + has now, from :func:`time.time`. The listing of the directory is kept in :data:`_name_listing`, and the directory is listed again if an image is later than the start of - the last file, - unless it was listed less than :data:`_age_refresh` ago. + the last file and it was listed more than :data:`_age_refresh` ago, + or if it was listed more than :data:`_age_check` ago. + A file whose time of change in the listing changes, + such as the last file as the database grows, + is given the time it was listed again, + so that a copy downloaded before then is downloaded again. + If the directory cannot be listed again, + the listing from before is used, with a warning. Parameters ---------- @@ -319,18 +422,49 @@ def _listing( url = f"{_url_coalign}{_directories_coalign[database]}/" path = _path_local(url + _name_listing, directory) + # The name of each file, when it was last changed, and when the listing + # first showed that change + before: dict[str, tuple[str, float]] = {} if path.is_file(): - names = path.read_text().split() - if names and not (latest >= _start(names[-1:])[0] and _stale(path)): - return names + for line in path.read_text().splitlines(): + name, changed, listed = line.split() + before[name] = (changed, float(listed)) + + if before: + last = _start(list(before)[-1:])[0] + if not ((latest >= last and _stale(path)) or _stale(path, _age_check)): + return {name: listed for name, (_, listed) in before.items()} + + try: + text = _get(url, num_retry).text + except _errors_network as e: + if not before: + raise + warnings.warn( + f"Could not list {url} again, so the listing from before is used: {e!r}", + stacklevel=2, + ) + return {name: listed for name, (_, listed) in before.items()} - names = sorted(set(_pattern_coalign.findall(_get(url, num_retry).text))) - if not names: + # A time of change the listing does not show is kept as "-" + found = {n: c.replace(" ", "T") or "-" for n, c in _pattern_coalign.findall(text)} + if not found: raise ValueError(f"The listing of {url} has no files of the database.") - _write(path, "\n".join(names).encode(), overwrite=True) + now = time.time() + result: dict[str, float] = {} + lines = [] + for name in sorted(found): + changed = found[name] + changed_before, listed = before.get(name, ("", now)) + if changed != changed_before: + listed = now + result[name] = listed + lines.append(f"{name} {changed} {listed!r}") + + _write(path, "\n".join(lines).encode(), overwrite=True) - return names + return result def _lookup( @@ -347,12 +481,12 @@ def _lookup( or NaN and -1 for an image the database has no entry for. Each file of the database covers the time from its start to the start of - the next one. - A file is downloaded again if an image later than its last entry has no - entry, since it may have been downloaded before the database reached the - image, - unless it was downloaded less than :data:`_age_refresh` ago. - An image before the last entry may have no entry of its own, + the next one, + and only the files which cover an image without an entry yet are read. + A file is downloaded again if the copy was downloaded before the listing + of :func:`_listing` showed the version the server has now, + as for the last file once the database has grown. + An image may have no entry, as an image the cross-correlation with AIA was not made for has none in ``coaldb_c``. @@ -372,7 +506,8 @@ def _lookup( ycen = np.full(num, np.nan) ref_type = np.full(num, -1, dtype=int) - names = _listing(database, float(np.max(seconds)), directory, num_retry) + listed = _listing(database, float(np.max(seconds)), directory, num_retry) + names = list(listed) starts = _start(names) # The files which cover the time within the tolerance of each image @@ -381,42 +516,58 @@ def _lookup( url_directory = f"{_url_coalign}{_directories_coalign[database]}/" - for i in range(max(int(first.min()), 0), int(last.max()) + 1): - images = (first <= i) & (i <= last) - if not np.any(images): + for i in np.unique(np.concatenate([first, last])): + images = np.flatnonzero((first <= i) & (i <= last) & (ref_type == -1)) + if i < 0 or images.size == 0: continue - url = f"{url_directory}{names[i]}.geny" - - refresh = False - while True: - path = _download_file( - url=url, - directory=directory, - overwrite=refresh, - num_retry=num_retry, - signature=_signature_genx, - kind="an IDL save file", - ) - entries = _restore(path)["p0"] - times = _seconds([t.decode() for t in entries["DATE_OBS"]]) - - for j in np.flatnonzero(images & (ref_type == -1)): - difference = np.abs(times - seconds[j]) - k = int(np.argmin(difference)) - if difference[k] <= _tolerance: - xcen[j] = entries["XCEN"][k] - ycen[j] = entries["YCEN"][k] - ref_type[j] = entries["REF_TYPE"][k] - - missing = np.any(images & (ref_type == -1) & (seconds > times.max())) - if refresh or not missing or not _stale(path): - break - refresh = True + name = names[i] + url = f"{url_directory}{name}.geny" + path = _path_local(url, directory) + overwrite = path.is_file() and path.stat().st_mtime < listed[name] + path = _download_database(url, directory, overwrite, num_retry) + + entries = _restore(path)["p0"] + times = _seconds([t.decode() for t in entries["DATE_OBS"]]) + + index, difference = _nearest(times, seconds[images]) + match = difference <= _tolerance + images = images[match] + index = index[match] + xcen[images] = entries["XCEN"][index] + ycen[images] = entries["YCEN"][index] + ref_type[images] = entries["REF_TYPE"][index] return xcen, ycen, ref_type +def _nearest( + times: np.ndarray, + seconds: np.ndarray, +) -> tuple[np.ndarray, np.ndarray]: + """ + The index of the entry nearest in time to each image, + the earlier one if two are as near, + and how far it is, in seconds. + + Parameters + ---------- + times + The time of each entry, in any order. + seconds + The start of each image. + """ + order = np.argsort(times, kind="stable") + ordered = times[order] + position = np.searchsorted(ordered, seconds) + before = np.clip(position - 1, 0, ordered.size - 1) + after = np.clip(position, 0, ordered.size - 1) + difference_before = np.abs(ordered[before] - seconds) + difference_after = np.abs(ordered[after] - seconds) + nearest = np.where(difference_after < difference_before, after, before) + return order[nearest], np.minimum(difference_before, difference_after) + + def _calibration(ref_type: np.ndarray) -> np.ndarray: """ The type of calibration of the pointing of each image, @@ -435,6 +586,29 @@ def _calibration(ref_type: np.ndarray) -> np.ndarray: return result +def _check_coalign(coalign: object) -> None: + """ + Check that a co-alignment database is one of those of + :data:`_directories_coalign`, + so that a wrong one is found before any image is downloaded. + + Parameters + ---------- + coalign + The database, which :obj:`None` is not. + + Raises + ------ + ValueError + If `coalign` is not ``"aia"`` or ``"ufss"``. + """ + if coalign not in tuple(_directories_coalign): + raise ValueError( + f"`coalign` must be one of {list(_directories_coalign)} or None, " + f"got {coalign!r}." + ) + + def _coalign( time: Sequence[str], filter_2: Sequence[int], @@ -483,12 +657,16 @@ def _coalign( If :obj:`None` (the default), :data:`hinode.directory_default` is used. num_retry The number of times to try to connect to the server. + + Raises + ------ + ValueError + If `coalign` is not ``"aia"`` or ``"ufss"``. + ConnectionError + If a file of the databases which has not been downloaded before + cannot be downloaded. """ - if coalign not in _directories_coalign: - raise ValueError( - f"`coalign` must be one of {list(_directories_coalign)} or None, " - f"got {coalign!r}." - ) + _check_coalign(coalign) root: pathlib.Path = hinode.directory_default if directory is None else directory @@ -502,16 +680,23 @@ def _coalign( xcen = np.full(num, np.nan) ycen = np.full(num, np.nan) calibration = np.full(num, -1, dtype=int) - for database in databases: - missing = calibration == -1 - if not np.any(missing): - break - x, y, ref_type = _lookup(database, seconds[missing], root, num_retry) - xcen[missing] = x - ycen[missing] = y - calibration[missing] = _calibration(ref_type) - - roll = _roll(seconds, root, num_retry) + try: + for database in databases: + missing = calibration == -1 + if not np.any(missing): + break + x, y, ref_type = _lookup(database, seconds[missing], root, num_retry) + xcen[missing] = x + ycen[missing] = y + calibration[missing] = _calibration(ref_type) + + roll = _roll(seconds, root, num_retry) + except _errors_network as e: + raise ConnectionError( + f"Could not download the co-alignment databases of SolarSoft " + f"from {_url_coalign}. " + f"Pass coalign=None to keep the pointing of the headers." + ) from e keywords: list[dict[str, float | str]] = [] for i in range(num): diff --git a/hinode/xrt/_coalign_test.py b/hinode/xrt/_coalign_test.py index 0724d0e..68e4bee 100644 --- a/hinode/xrt/_coalign_test.py +++ b/hinode/xrt/_coalign_test.py @@ -2,23 +2,30 @@ import time import pathlib import datetime +import warnings import pytest import requests import numpy as np import astropy.io.fits import hinode +import hinode.xrt._data import hinode.xrt._coalign from hinode.xrt._coalign import ( _seconds, + _isot, _start, _interpol, _roll, _database_roll, _listing, _lookup, + _nearest, _calibration, + _check_coalign, _coalign, + _download_database, _age_refresh, + _age_check, _name_listing, _url_coalign, ) @@ -88,6 +95,10 @@ def test_start() -> None: assert np.all(result == _seconds(["2019-09-30T15:03", "2019-09-30T21:00"])) +def test_isot() -> None: + assert _isot(_seconds(["2026-10-04T05:01:43"])[0]) == "2026-10-04T05:01:43" + + def test_interpol() -> None: """Values beyond the samples are extrapolated from the first and the last interval.""" xp = np.array([0.0, 1.0, 3.0]) @@ -153,10 +164,25 @@ def test_calibration(ref_type: int, expected: int) -> None: assert _calibration(np.array([ref_type]))[0] == expected +def _row(name: str, changed: str) -> str: + """A row of a listing of a directory of a database, as the server writes it.""" + return ( + f'{name}.geny' + f'{changed} \n' + ) + + +def _age(path: pathlib.Path, seconds: float) -> None: + """Make a file look as if it was written some seconds ago.""" + old = time.time() - seconds + os.utime(path, (old, old)) + + def test_listing(tmp_path: pathlib.Path, monkeypatch: pytest.MonkeyPatch) -> None: """ The directory is listed again only if an image is later than the start of - the last file and the listing is older than a day. + the last file and the listing is older than a day, + or if the listing is older than a month. """ calls = [] @@ -169,23 +195,88 @@ def get(url: str, num_retry: int = 5) -> requests.Response: early = _seconds(["2019-09-30T12:00"])[0] late = _seconds(["2019-10-01T00:00"])[0] - assert _listing("aia", early, tmp_path) == ["20190930_0631", "20190930_1503"] + result = _listing("aia", early, tmp_path) + assert list(result) == ["20190930_0631", "20190930_1503"] assert calls == [f"{_url_coalign}coaldb_c/"] # Kept, since the image is in a file before the last one - _listing("aia", early, tmp_path) + assert _listing("aia", early, tmp_path) == result # Kept, since the listing is new _listing("aia", late, tmp_path) assert len(calls) == 1 path = tmp_path / "sdb/hinode/xrt/xrt_msu_coalign/coaldb_c" / _name_listing - old = time.time() - 2 * _age_refresh - os.utime(path, (old, old)) + _age(path, 2 * _age_refresh) _listing("aia", early, tmp_path) assert len(calls) == 1 _listing("aia", late, tmp_path) assert len(calls) == 2 + # Listed again for any image once the listing is a month old + _age(path, 2 * _age_check) + _listing("aia", early, tmp_path) + assert len(calls) == 3 + + +def test_listing_changed( + tmp_path: pathlib.Path, + monkeypatch: pytest.MonkeyPatch, +) -> None: + """ + A file is given the time it was listed again only if the listing shows + that it has changed, or if it is new. + """ + text = ( + _row("20190930_0631", "2019-10-01 10:00") + + _row("20190930_1503", "2019-10-08 10:00") + + _row("20191001_0000", "2019-10-08 10:00") + ) + monkeypatch.setattr(hinode.xrt._coalign, "_get", lambda url, n=5: _response(text)) + + path = tmp_path / "sdb/hinode/xrt/xrt_msu_coalign/coaldb_c" / _name_listing + path.parent.mkdir(parents=True) + path.write_text( + "20190930_0631 2019-10-01T10:00 1.0\n20190930_1503 2019-10-01T10:00 1.0" + ) + _age(path, 2 * _age_check) + + start = time.time() + result = _listing("aia", 0, tmp_path) + + assert result["20190930_0631"] == 1.0 + assert result["20190930_1503"] >= start + assert result["20191001_0000"] >= start + + # The times are kept for the next listing + assert _listing("aia", 0, tmp_path) == result + assert "20190930_1503 2019-10-08T10:00 " in path.read_text() + + +def test_listing_unreachable( + tmp_path: pathlib.Path, + monkeypatch: pytest.MonkeyPatch, +) -> None: + """ + If the directory cannot be listed again, the listing from before is used, + and if there is none, the error is raised. + """ + + def get(url: str, num_retry: int = 5) -> requests.Response: + raise ConnectionError("The server cannot be reached.") + + monkeypatch.setattr(hinode.xrt._coalign, "_get", get) + + with pytest.raises(ConnectionError): + _listing("ufss", 0, tmp_path) + + path = tmp_path / "sdb/hinode/xrt/xrt_msu_coalign/coaldb_u" / _name_listing + path.parent.mkdir(parents=True) + path.write_text("20190930_0631 - 1.0") + _age(path, 2 * _age_check) + + with pytest.warns(UserWarning, match="the listing from before is used"): + assert _listing("ufss", 0, tmp_path) == {"20190930_0631": 1.0} + def test_listing_empty(tmp_path: pathlib.Path, monkeypatch: pytest.MonkeyPatch) -> None: monkeypatch.setattr(hinode.xrt._coalign, "_get", lambda url, n=5: _response("")) @@ -277,94 +368,166 @@ def test_coalign_invalid() -> None: _coalign([], [], [], [], [], "aia_cc") # type: ignore[arg-type] +def test_check_coalign() -> None: + _check_coalign("aia") + _check_coalign("ufss") + for coalign in ["AIA", None, True]: + with pytest.raises(ValueError, match="`coalign` must be one of"): + _check_coalign(coalign) + + +def test_coalign_unreachable(monkeypatch: pytest.MonkeyPatch) -> None: + """A database which cannot be downloaded says how to do without it.""" + + def lookup(*args: object) -> None: + raise ConnectionError("The server cannot be reached.") + + monkeypatch.setattr(hinode.xrt._coalign, "_lookup", lookup) + + with pytest.raises(ConnectionError, match="coalign=None"): + _coalign([_date_obs], [1], [1], [384], [384], "aia") + + +def test_download_database_unreachable( + tmp_path: pathlib.Path, + monkeypatch: pytest.MonkeyPatch, +) -> None: + """ + A file which cannot be downloaded again is kept, with a warning, + and one which has never been downloaded raises the error. + """ + + def get(url: str, num_retry: int = 5) -> requests.Response: + raise ConnectionError("The server cannot be reached.") + + monkeypatch.setattr(hinode.xrt._data, "_get", get) + + url = f"{_url_coalign}coaldb_u/20190930_0631.geny" + + with pytest.raises(ConnectionError): + _download_database(url, tmp_path) + + path = _path_local(url, tmp_path) + path.parent.mkdir(parents=True) + path.write_bytes(b"SR\x00\x04") + + assert _download_database(url, tmp_path) == path + with pytest.warns(UserWarning, match="the copy downloaded before is used"): + assert _download_database(url, tmp_path, overwrite=True) == path + + +def _entries(times: list[str], ref_type: int) -> np.recarray: + """The entries of a file of a co-alignment database at the given times.""" + num = len(times) + return np.rec.fromarrays( + [ + np.array([t.encode() for t in times]), + np.ones(num), + np.ones(num), + np.full(num, ref_type), + ], + names=["DATE_OBS", "XCEN", "YCEN", "REF_TYPE"], + ) + + @pytest.mark.parametrize( - argnames="isot,stale,num_downloads", + argnames="cached,age,overwrite", argvalues=[ - # After the last entry of the file, which was downloaded long ago - ("2019-09-30T16:00:00.000", True, 2), - # After the last entry, but the file is new - ("2019-09-30T16:00:00.000", False, 1), - # Before the last entry, with no entry of its own - ("2019-09-30T15:30:00.000", True, 1), + # Never downloaded + (False, 0, False), + # Downloaded before the listing showed the version on the server + (True, 10, True), + # Downloaded after it + (True, -10, False), ], ) def test_lookup_refresh( + tmp_path: pathlib.Path, monkeypatch: pytest.MonkeyPatch, - isot: str, - stale: bool, - num_downloads: int, + cached: bool, + age: float, + overwrite: bool, ) -> None: """ - A file is downloaded again only if an image later than its last entry has - no entry, and the file is old. + A file is downloaded again only if the copy is older than the version + the listing shows. """ + listed = time.time() - 100 + url = f"{_url_coalign}coaldb_c/20190930_1503.geny" + path = _path_local(url, tmp_path) + if cached: + path.parent.mkdir(parents=True) + path.touch() + os.utime(path, (listed - age, listed - age)) + downloads = [] - def download_file( - url: str, directory: pathlib.Path, overwrite: bool, **kwargs: object + def download_database( + url: str, + directory: pathlib.Path, + overwrite: bool = False, + num_retry: int = 5, ) -> pathlib.Path: downloads.append(overwrite) - return directory / "file.geny" + return path - def restore(path: pathlib.Path) -> dict: - times = ["2019-09-30T15:10:00.000", "2019-09-30T15:50:00.000"] - if len(downloads) > 1: - times.append("2019-09-30T16:00:00.000") - num = len(times) - entries = np.rec.fromarrays( - [ - np.array([t.encode() for t in times]), - np.ones(num), - np.ones(num), - np.full(num, 10), - ], - names=["DATE_OBS", "XCEN", "YCEN", "REF_TYPE"], - ) - return {"p0": entries} + isot = "2019-09-30T15:30:00.000" monkeypatch.setattr( - hinode.xrt._coalign, "_listing", lambda *args: ["20190930_1503"] + hinode.xrt._coalign, "_listing", lambda *args: {"20190930_1503": listed} + ) + monkeypatch.setattr(hinode.xrt._coalign, "_download_database", download_database) + monkeypatch.setattr( + hinode.xrt._coalign, "_restore", lambda path: {"p0": _entries([isot], 10)} ) - monkeypatch.setattr(hinode.xrt._coalign, "_download_file", download_file) - monkeypatch.setattr(hinode.xrt._coalign, "_restore", restore) - monkeypatch.setattr(hinode.xrt._coalign, "_stale", lambda path: stale) - _, _, ref_type = _lookup("aia", _seconds([isot]), pathlib.Path("cache")) + _, _, ref_type = _lookup("aia", _seconds([isot]), tmp_path) - assert downloads == [False, True][:num_downloads] - assert ref_type[0] == (10 if num_downloads == 2 else -1) + assert downloads == [overwrite] + assert ref_type[0] == 10 @pytest.mark.parametrize( - argnames="latest,stale,overwrite", + argnames="cached,latest,age,downloads_expected", argvalues=[ + # Never downloaded + (False, 5.0, 0, [False]), # After the last measurement, in a database downloaded long ago - (20.0, True, True), + (True, 20.0, 2 * _age_refresh, [True]), # Within the database - (5.0, True, False), + (True, 5.0, 2 * _age_refresh, []), # After the last measurement, in a new database - (20.0, False, False), + (True, 20.0, 0, []), + # Within the database, but checked again after a month + (True, 5.0, 2 * _age_check, [True]), ], ) def test_database_roll_refresh( tmp_path: pathlib.Path, monkeypatch: pytest.MonkeyPatch, + cached: bool, latest: float, - stale: bool, - overwrite: bool, + age: float, + downloads_expected: list[bool], ) -> None: """ The database of the roll is downloaded again only for an image after its - last measurement, and only if it is old. + last measurement, if it is a day old, or for any image, if it is a month + old. """ path = _path_local(_url_coalign + "xrt_rollangle_db.geny", tmp_path) - path.parent.mkdir(parents=True) - path.touch() + if cached: + path.parent.mkdir(parents=True) + path.touch() + _age(path, age) downloads = [] - def download_file( - url: str, directory: pathlib.Path, overwrite: bool, **kwargs: object + def download_database( + url: str, + directory: pathlib.Path, + overwrite: bool = False, + num_retry: int = 5, ) -> pathlib.Path: downloads.append(overwrite) return path @@ -372,48 +535,75 @@ def download_file( def restore(path: pathlib.Path) -> dict: return {"p0": np.array([0.0, 10.0]), "p1": np.array([1.0, 2.0])} - monkeypatch.setattr(hinode.xrt._coalign, "_download_file", download_file) + monkeypatch.setattr(hinode.xrt._coalign, "_download_database", download_database) monkeypatch.setattr(hinode.xrt._coalign, "_restore", restore) - monkeypatch.setattr(hinode.xrt._coalign, "_stale", lambda path: stale) times, angles = _database_roll(latest, tmp_path) - assert downloads == [overwrite] + assert downloads == downloads_expected assert np.all(times == [0, 10]) and np.all(angles == [1, 2]) +def test_roll_extrapolated(monkeypatch: pytest.MonkeyPatch) -> None: + """ + The roll of an image after the last measurement of the database is + extrapolated, with a warning. + """ + end = _seconds(["2021-12-27T15:10"])[0] + times = end + np.array([100.0, 200.0]) + monkeypatch.setattr( + hinode.xrt._coalign, + "_database_roll", + lambda *args: (times, np.array([1.0, 2.0])), + ) + + with warnings.catch_warnings(): + warnings.simplefilter("error") + assert np.allclose(_roll(end + np.array([150.0]), pathlib.Path()), -1.5) + + with pytest.warns(UserWarning, match="1 of 2 images began outside"): + result = _roll(end + np.array([150.0, 300.0]), pathlib.Path()) + assert np.allclose(result, [-1.5, -3.0]) + + +def test_nearest() -> None: + """The entries may be in any order, and a tie goes to the earlier one.""" + times = np.array([30.0, 10.0, 20.0]) + index, difference = _nearest(times, np.array([0.0, 14.0, 16.0, 25.0, 40.0])) + assert index.tolist() == [1, 1, 2, 2, 0] + assert difference.tolist() == [10, 4, 4, 5, 10] + + def test_lookup_gap(monkeypatch: pytest.MonkeyPatch) -> None: """ - Only the files which cover an image are downloaded, - and not the ones between them. + Only the files which cover an image without an entry yet are downloaded, + and not the ones between them, + or the next one for an image near its start which the file before has + an entry for. """ - isot = ["2019-09-30T00:30:00.000", "2019-09-30T02:30:00.000"] + isot = [ + "2019-09-30T00:30:00.000", + "2019-09-30T02:30:00.000", + "2019-09-30T02:59:59.950", + ] downloads = [] - def download_file(url: str, directory: pathlib.Path, **kwargs: object) -> str: + def download_database(url: str, *args: object) -> str: downloads.append(url.split("/")[-1]) return url def restore(path: str) -> dict: times = [t for t in isot if t[11:13] == path.split("_")[-1][:2]] - num = len(times) - entries = np.rec.fromarrays( - [ - np.array([t.encode() for t in times]), - np.ones(num), - np.ones(num), - np.full(num, 3), - ], - names=["DATE_OBS", "XCEN", "YCEN", "REF_TYPE"], - ) - return {"p0": entries} - - names = ["20190930_0000", "20190930_0100", "20190930_0200"] - monkeypatch.setattr(hinode.xrt._coalign, "_listing", lambda *args: names) - monkeypatch.setattr(hinode.xrt._coalign, "_download_file", download_file) + return {"p0": _entries(times, 3)} + + names = ["20190930_0000", "20190930_0100", "20190930_0200", "20190930_0300"] + monkeypatch.setattr( + hinode.xrt._coalign, "_listing", lambda *args: dict.fromkeys(names, 0.0) + ) + monkeypatch.setattr(hinode.xrt._coalign, "_download_database", download_database) monkeypatch.setattr(hinode.xrt._coalign, "_restore", restore) _, _, ref_type = _lookup("ufss", _seconds(isot), pathlib.Path("cache")) assert downloads == ["20190930_0000.geny", "20190930_0200.geny"] - assert ref_type.tolist() == [3, 3] + assert ref_type.tolist() == [3, 3, 3] diff --git a/hinode/xrt/_filtergrams.py b/hinode/xrt/_filtergrams.py index 10ee754..6092d6c 100644 --- a/hinode/xrt/_filtergrams.py +++ b/hinode/xrt/_filtergrams.py @@ -14,7 +14,7 @@ from ._leak import _leak, _history_leak from ._uncertainty import _quality, _error_jpeg, _error_dark from ._vignetting import _vignetting, _error_vignetting -from ._coalign import _coalign +from ._coalign import _check_coalign, _coalign __all__ = [ "Filtergram", @@ -288,7 +288,8 @@ class Filtergram( 4 for a fit to the limb in X-rays, 6 for another method, and -1 for an image no database had an entry for, - which keeps the pointing of its header, + which keeps the center and the plate scale of its header + but is still given the roll of ``xrt_rollangle.pro``, or :obj:`None` if the pointing was not corrected. """ @@ -362,7 +363,19 @@ def from_time_range( as ``xrt_read_coaldb.pro`` does with ``/aia_cc``. ``"ufss"`` uses the UFSS alone, as it does by default. See :attr:`coalignment`. + + Raises + ------ + ValueError + If `coalign` is not ``"aia"``, ``"ufss"``, or :obj:`None`, + which is checked before any image is downloaded, + or for the reasons :meth:`from_fits` gives. + ConnectionError + For the reasons :meth:`from_fits` gives. """ + if coalign is not None: + _check_coalign(coalign) + urls = hinode.xrt.urls( time_start=time_start, time_stop=time_stop, @@ -466,6 +479,10 @@ def from_fits( if `uncertainty` is :obj:`True` and the compression of an image is not known, or if `coalign` is not ``"aia"``, ``"ufss"``, or :obj:`None`. + ConnectionError + If a file of the co-alignment databases an image needs has not + been downloaded before and cannot be downloaded, + in which case ``coalign=None`` keeps the pointing of the headers. Notes ----- @@ -491,7 +508,24 @@ def from_fits( the correction places the center the database gives at pixel ``NAXIS // 2 + 0.5``, counted from one, with the plate scale of the X-ray or the G-band images. + + The files of the databases are downloaded to `directory` the first + time an image needs them. + They are downloaded again when the listing of the server shows they + have changed, + and the listings and the database of the roll are checked again after + a day for an image later than they reach, + and after a month for any image. + If the server cannot be reached, + the copies downloaded before are used, with a warning. + A warning is also given for an image the database of the roll does + not reach yet, + since it is extrapolated, as ``xrt_rollangle.pro`` extrapolates it, + and the measurements reach the server days after the images. """ + if coalign is not None: + _check_coalign(coalign) + if isinstance(path, na.AbstractScalarArray): path = [str(p) for p in np.ravel(np.asarray(path.ndarray))] elif isinstance(path, (str, os.PathLike)): diff --git a/hinode/xrt/_filtergrams_test.py b/hinode/xrt/_filtergrams_test.py index 514e5d9..7346dcf 100644 --- a/hinode/xrt/_filtergrams_test.py +++ b/hinode/xrt/_filtergrams_test.py @@ -1,6 +1,7 @@ import re import pathlib import dataclasses +import unittest.mock import pytest import numpy as np import astropy.units as u @@ -180,14 +181,6 @@ def test_coalignment(self, array: hinode.xrt.Filtergram) -> None: assert isinstance(coalignment, na.ScalarArray) assert coalignment.shape == {array.axis_time: array.shape[array.axis_time]} - # The cross-correlation with AIA has every one of these images, - # and the UFSS database alone gives the pointing of the headers. - assert np.all(coalignment.ndarray == 1) or np.all(coalignment.ndarray == 3) - - image = array[{array.axis_time: 0}] - assert image.coalignment is not None - assert image.coalignment.shape == {} - def test_vignetting(self, array: hinode.xrt.Filtergram) -> None: vignetting = array.vignetting if vignetting is None: @@ -235,7 +228,7 @@ def test_getitem(self, array: hinode.xrt.Filtergram) -> None: assert result.saturated is not None assert array.saturated is not None assert np.all(result.saturated == array.saturated[index]) - for name in ["vignetting", "uncertainty_map", "leak"]: + for name in ["vignetting", "uncertainty_map", "leak", "coalignment"]: field = getattr(array, name) if field is not None: assert np.all(getattr(result, name) == field[index]) @@ -585,6 +578,61 @@ def test_from_fits_history( hinode.xrt.Filtergram.from_fits(file) +@pytest.mark.parametrize( + argnames="coalign,expected", + argvalues=[ + ("aia", 1), + ("ufss", 3), + (None, None), + ], +) +def test_from_fits_coalign(coalign: None | str, expected: None | int) -> None: + """ + The cross-correlation with AIA has an entry for this image, + and the UFSS database alone gives the pointing of its header. + """ + path = _path("2019-09-30T18:08:30", "2019-09-30T18:08:40") + result = hinode.xrt.Filtergram.from_fits(path, coalign=coalign) # type: ignore[arg-type] + if expected is None: + assert result.coalignment is None + else: + assert result.coalignment is not None + assert result.coalignment.ndarray.tolist() == [expected] + + +def test_from_fits_coalign_directory(tmp_path: pathlib.Path) -> None: + """The co-alignment databases are placed in the directory given.""" + path = _path("2019-09-30T18:08:30", "2019-09-30T18:08:40") + hinode.xrt.Filtergram.from_fits(path, coalign="ufss", directory=tmp_path) + database = tmp_path / "sdb/hinode/xrt/xrt_msu_coalign/coaldb_u" + assert (database / "listing.txt").is_file() + assert len(list(database.glob("*.geny"))) == 1 + + +def test_from_time_range_coalign_invalid(monkeypatch: pytest.MonkeyPatch) -> None: + """A wrong database is found before the archive is searched.""" + urls = unittest.mock.Mock() + monkeypatch.setattr(hinode.xrt, "urls", urls) + + with pytest.raises(ValueError, match="`coalign` must be one of"): + hinode.xrt.Filtergram.from_time_range( + time_start=_time_start, + time_stop=_time_stop, + coalign="AIA", # type: ignore[arg-type] + ) + + urls.assert_not_called() + + +def test_from_fits_coalign_invalid(tmp_path: pathlib.Path) -> None: + """A wrong database is found before the files are read.""" + with pytest.raises(ValueError, match="`coalign` must be one of"): + hinode.xrt.Filtergram.from_fits( + tmp_path / "missing.fits", + coalign="AIA", # type: ignore[arg-type] + ) + + @pytest.mark.parametrize("coalign", [None, "aia"]) def test_inputs_against_astropy_wcs(coalign: None | str) -> None: """