Merge branch 'ndcubify' of https://github.com/DanRyanIrish/stixpy int…

…o ndcubify

changelog/111.feature.rst

@@ -0,0 +1 @@
+Update :func:`~stixpy.coordinates.transforms.get_hpc_info` to use mean or interpolate pointing and location data based on number of data points in given timerange. Add function cache for aux data frequenly needed during coordinate transforms.

changelog/112.doc.rst

@@ -0,0 +1 @@
+Correct calculation of the HPC coordinate of the center of the STIX pointing in the imaging demo.

examples/imaging_demo.py

@@ -129,7 +129,7 @@
 header = make_fitswcs_header(bp_image, coord, telescope='STIX', observatory='Solar Orbiter', scale=[10,10]*u.arcsec/u.pix)
 fd_bp_map = Map((bp_image, header))
 
-hpc_ref = Helioprojective(pointing[0], pointing[1], observer=solo, obstime=fd_bp_map.date)
+hpc_ref = coord.transform_to(Helioprojective(observer=solo, obstime=fd_bp_map.date))  # Center of STIX pointing in HPC
 header_hp = make_fitswcs_header(bp_image, hpc_ref, scale=[10, 10]*u.arcsec/u.pix, rotation_angle=90*u.deg+roll)
 hp_map = Map((bp_image, header_hp))
 hp_map_rotated = hp_map.rotate()

stixpy/coordinates/frames.py

@@ -14,10 +14,6 @@
 
 __all__ = ["STIXImaging"]
 
-STIX_X_SHIFT = 26.1 * u.arcsec  # fall back to this when non sas solution available
-STIX_Y_SHIFT = 58.2 * u.arcsec  # fall back to this when non sas solution available
-STIX_X_OFFSET = 60.0 * u.arcsec  # remaining offset after SAS solution
-STIX_Y_OFFSET = 8.0 * u.arcsec  # remaining offset after SAS solution
 
 STIX_X_CTYPE = "SXLN-TAN"
 STIX_Y_CTYPE = "SXLT-TAN"

stixpy/coordinates/tests/test_frames.py

@@ -9,6 +9,7 @@
 from sunpy.map.mapbase import SpatialPair
 
 from stixpy.coordinates.frames import STIXImaging, stix_frame_to_wcs, stix_wcs_to_frame
+from stixpy.coordinates.transforms import hpc_to_stixim, stixim_to_hpc  # noqa
 from stixpy.map.stix import STIXMap
 
 

stixpy/coordinates/tests/test_transforms.py

@@ -1,13 +1,13 @@
 from unittest import mock
 
 import astropy.units as u
-import numpy as np
 import pytest
 from astropy.tests.helper import assert_quantity_allclose
 from astropy.time import Time
 from sunpy.coordinates import HeliographicStonyhurst, Helioprojective
 
 from stixpy.coordinates.frames import STIXImaging
+from stixpy.coordinates.transforms import _get_aux_data, get_hpc_info
 
 
 @pytest.mark.skip(reason="Test data maybe incorrect")
@@ -28,32 +28,51 @@ def test_transforms_with_know_values():
     assert_quantity_allclose(earth_coord.Ty, stix_to_earth.Ty)
 
 
-@mock.patch("stixpy.coordinates.transforms._get_aux_data")
+@mock.patch("stixpy.coordinates.transforms.get_hpc_info")
 def test_hpc_to_stx(mock):
     # fake some data to make check easier
     obstime = Time("2021-05-22T02:52:00")
-    roll, yaw, pitch, sas_x, sas_y = [0, 0, 0, 0, 0] * u.arcsec
+    roll = 0 * u.deg
+    sas = [10, 15] * u.arcsec
     solo_pos = HeliographicStonyhurst(0 * u.deg, 0 * u.deg, 1 * u.AU, obstime=obstime)
-    mock.return_value = (roll, yaw, pitch, sas_x, sas_y, solo_pos.cartesian)
+    mock.return_value = (roll, solo_pos.cartesian, sas )
     solo_hpc_coord = Helioprojective(0 * u.arcsec, 0 * u.arcsec, observer=solo_pos, obstime=obstime)
     stix_frame = STIXImaging(0 * u.arcsec, 0 * u.arcsec, obstime=obstime)
     stix_coord = solo_hpc_coord.transform_to(stix_frame)
-    # should match the offset -8, 60
-    assert_quantity_allclose(stix_coord.Tx, -8 * u.arcsec)
-    assert_quantity_allclose(stix_coord.Ty, 60 * u.arcsec)
+    assert_quantity_allclose(stix_coord.Tx, -sas[1])
+    assert_quantity_allclose(stix_coord.Ty, sas[0])
 
 
-@mock.patch("stixpy.coordinates.transforms._get_aux_data")
-def test_hpc_to_stx_no_sas(mock):
-    # fake some data to make check easier
-    obstime = Time("2021-05-22T02:52:00")
-    roll, yaw, pitch, sas_x, sas_y = [0, 10, 10, np.nan, np.nan] * u.arcsec
-    solo_pos = HeliographicStonyhurst(0 * u.deg, 0 * u.deg, 1 * u.AU, obstime=obstime)
-    mock.return_value = (roll, yaw, pitch, sas_x, sas_y, solo_pos.cartesian)
-    solo_hpc_coord = Helioprojective(0 * u.arcsec, 0 * u.arcsec, observer=solo_pos, obstime=obstime)
-    stix_frame = STIXImaging(0 * u.arcsec, 0 * u.arcsec, obstime=obstime)
-    with pytest.warns(Warning, match="SAS solution not"):
-        stix_coord = solo_hpc_coord.transform_to(stix_frame)
-        # should match the offset -8, 60 added to yaw and pitch 10, 10
-        assert_quantity_allclose(stix_coord.Tx, (10 - 8) * u.arcsec)
-        assert_quantity_allclose(stix_coord.Ty, (10 + 60) * u.arcsec)
+@pytest.mark.remote_data
+def test_get_aux_data():
+    with pytest.raises(ValueError, match="No STIX pointing data found for time range"):
+        _get_aux_data(Time('2015-06-06'))  # Before the mission started
+
+    aux_data = _get_aux_data(Time('2022-08-28T16:02:00'))
+    assert len(aux_data) == 1341
+
+    aux_data = _get_aux_data(Time('2022-08-28T16:02:00'), end_time=Time('2022-08-28T16:04:00'))
+    assert len(aux_data) == 1341
+
+    aux_data = _get_aux_data(Time('2022-08-28T23:58:00'), end_time=Time('2022-08-29T00:02:00'))
+    assert len(aux_data) == 2691
+
+
+@pytest.mark.remote_data
+def test_get_hpc_info():
+    with pytest.warns(match='SAS solution not available.*'):
+        roll, solo_heeq, stix_pointing = get_hpc_info(Time('2022-08-28T16:00:16'), end_time=Time('2022-08-28T16:00:17.000'))
+
+    assert_quantity_allclose(-3.3733306 * u.deg, roll)
+    assert_quantity_allclose([25.61525646, 57.914266] * u.arcsec, stix_pointing)
+    assert_quantity_allclose([-97868816.,  62984852.,  -5531986.5] * u.km, solo_heeq)
+
+    roll, solo_heeq, stix_pointing = get_hpc_info(Time('2022-08-28T16:00:00'), end_time=Time('2022-08-28T16:03:59.575'))
+    assert_quantity_allclose(-3.3732104 * u.deg, roll)
+    assert_quantity_allclose([25.923784, 58.179676]*u.arcsec, stix_pointing)
+    assert_quantity_allclose([-9.7867768e+07, 6.2983744e+07, -5532067.5]*u.km, solo_heeq)
+
+    roll, solo_heeq, stix_pointing = get_hpc_info(Time('2022-08-28T21:00:00'), end_time=Time('2022-08-28T21:03:59.575'))
+    assert_quantity_allclose(-3.3619127 * u.deg, roll)
+    assert_quantity_allclose([-26.070198, 173.48871]*u.arcsec, stix_pointing)
+    assert_quantity_allclose([-9.7671984e+07, 6.2774768e+07,  -5547166.0]*u.km, solo_heeq)

stixpy/coordinates/transforms.py

@@ -1,20 +1,26 @@
 import warnings
+from functools import lru_cache
 
 import astropy.coordinates as coord
 import astropy.units as u
 import numpy as np
 from astropy.coordinates import frame_transform_graph
 from astropy.coordinates.matrix_utilities import matrix_transpose, rotation_matrix
 from astropy.io import fits
-from astropy.table import QTable
+from astropy.table import QTable, vstack
 from astropy.time import Time
 from sunpy.coordinates import HeliographicStonyhurst, Helioprojective
 from sunpy.net import Fido
 from sunpy.net import attrs as a
 
-from stixpy.coordinates.frames import STIX_X_OFFSET, STIX_X_SHIFT, STIX_Y_OFFSET, STIX_Y_SHIFT, STIXImaging
+from stixpy.coordinates.frames import STIXImaging
 from stixpy.utils.logging import get_logger
 
+STIX_X_SHIFT = 26.1 * u.arcsec  # fall back to this when non sas solution available
+STIX_Y_SHIFT = 58.2 * u.arcsec  # fall back to this when non sas solution available
+STIX_X_OFFSET = 60.0 * u.arcsec  # remaining offset after SAS solution
+STIX_Y_OFFSET = 8.0 * u.arcsec  # remaining offset after SAS solution
+
 logger = get_logger(__name__)
 
 __all__ = ['get_hpc_info', 'stixim_to_hpc', 'hpc_to_stixim']
@@ -49,67 +55,145 @@ def _get_rotation_matrix_and_position(obstime):
     return rmatrix, solo_position_heeq
 
 
-def get_hpc_info(obstime):
-    roll, pitch, yaw, sas_x, sas_y, solo_position_heeq = _get_aux_data(obstime)
-    if np.isnan(sas_x) or np.isnan(sas_y):
-        warnings.warn("SAS solution not available using spacecraft pointing and average SAS offset")
-        sas_x = yaw
-        sas_y = pitch
+def get_hpc_info(start_time, end_time=None):
+    r"""
+    Get STIX pointing and SO location from L2 aspect files.
+
+    Parameters
+    ----------
+    start_time : `astropy.time.Time`
+        Time or start of a time interval.
+    end_time : `astropy.time.Time`, optional
+        End of time interval
+
+    Returns
+    -------
+
+    """
+    aux = _get_aux_data(start_time, end_time=end_time)
+    if end_time is None:
+        end_time = start_time
+    indices = np.argwhere((aux['time'] >= start_time) & (aux['time'] <= end_time))
+    indices = indices.flatten()
+
+    if indices.size < 2:
+        logger.info('Only one data point found interpolating between two closest times.')
+        # Times contained in one time bin have to interpolate
+        time_center = start_time + (end_time - start_time) * 0.5
+        diff_center = (aux['time']-time_center).to('s')
+        closest_index = np.argmin(np.abs(diff_center))
+
+        if diff_center[closest_index] > 0:
+            start_ind = closest_index-1
+            end_ind = closest_index + 1
+        else:
+            start_ind = closest_index
+            end_ind = closest_index + 2
+
+        interp = slice(start_ind, end_ind)
+        dt = np.diff(aux[start_ind:end_ind]['time'])[0].to(u.s)
+
+        roll, pitch, yaw = (aux[start_ind:start_ind+1]['roll_angle_rpy']
+                            + (np.diff(aux[interp]['roll_angle_rpy'], axis=0) / dt) * diff_center[closest_index])[0]
+        solo_heeq = (aux[start_ind:start_ind+1]['solo_loc_heeq_zxy']
+                     + (np.diff(aux[interp]['solo_loc_heeq_zxy'], axis=0) / dt) * diff_center[closest_index])[0]
+        sas_x = (aux[start_ind:start_ind+1]['y_srf']
+                 + (np.diff(aux[interp]['y_srf']) / dt) * diff_center[closest_index])[0]
+        sas_y = (aux[start_ind:start_ind+1]['z_srf']
+                 + (np.diff(aux[interp]['z_srf']) / dt) * diff_center[closest_index])[0]
+
+        good_solution = np.where(aux[interp]['sas_ok'] == 1)
+        good_sas = aux[good_solution]
+    else:
+        # average over
+        aux = aux[indices]
+
+        roll, pitch, yaw = np.mean(aux['roll_angle_rpy'], axis=0)
+        solo_heeq = np.mean(aux['solo_loc_heeq_zxy'], axis=0)
+
+        good_solution = np.where(aux['sas_ok'] == 1)
+        good_sas = aux[good_solution]
+
+        if len(good_sas) == 0:
+            warnings.warn(f"No good SAS solution found for time range: {start_time} to {end_time}.")
+            sas_x = 0
+            sas_y = 0
+        else:
+            sas_x = np.mean(good_sas["y_srf"])
+            sas_y = np.mean(good_sas["z_srf"])
+            sigma_x = np.std(good_sas["y_srf"])
+            sigma_y = np.std(good_sas["z_srf"])
+            tolerance = 3*u.arcsec
+            if sigma_x > tolerance or sigma_y > tolerance:
+                warnings.warn(f"Pointing unstable: StD(X) = {sigma_x}, StD(Y) = {sigma_y}.")
+
     # Convert the spacecraft pointing to STIX frame
     rotated_yaw = -yaw * np.cos(roll) + pitch * np.sin(roll)
     rotated_pitch = yaw * np.sin(roll) + pitch * np.cos(roll)
     spacecraft_pointing = np.hstack([STIX_X_SHIFT + rotated_yaw, STIX_Y_SHIFT + rotated_pitch])
+    stix_pointing = spacecraft_pointing
     sas_pointing = np.hstack([sas_x + STIX_X_OFFSET, -1 * sas_y + STIX_Y_OFFSET])
+
     pointing_diff = np.linalg.norm(spacecraft_pointing - sas_pointing)
-    if pointing_diff < 200 * u.arcsec:
-        logger.debug(f"Using SAS pointing {sas_pointing}")
-        spacecraft_pointing = sas_pointing
+    if np.all(np.isfinite(sas_pointing)) and len(good_sas) > 0:
+
+        if pointing_diff < 200 * u.arcsec:
+            logger.info(f"Using SAS pointing: {sas_pointing}")
+            stix_pointing = sas_pointing
+        else:
+            warnings.warn(f"Using spacecraft pointing: {spacecraft_pointing}"
+                          f" large difference between SAS and spacecraft.")
     else:
-        logger.debug(f"Using spacecraft pointing {spacecraft_pointing}")
-        warnings.warn(
-            f"Using spacecraft pointing large difference between "
-            f"SAS {sas_pointing} and spacecraft {spacecraft_pointing}."
-        )
-    return roll, solo_position_heeq, spacecraft_pointing
+        warnings.warn(f"SAS solution not available using spacecraft pointing: {stix_pointing}.")
 
+    return roll, solo_heeq, stix_pointing
 
-def _get_aux_data(obstime):
+
+@lru_cache
+def _get_aux_data(start_time, end_time=None):
+    r"""
+    Search, download and read L2 pointing data.
+
+    Parameters
+    ----------
+    start_time : `astropy.time.Time`
+        Time or start of a time interval.
+    end_time : `astropy.time.Time`, optional
+        End of time interval
+
+    Returns
+    -------
+
+    """
     # Find, download, read aux file with pointing, sas and position information
-    logger.debug("Searching for AUX data")
+    if end_time is None:
+        end_time = start_time
+    logger.debug(f"Searching for AUX data: {start_time} - {end_time}")
     query = Fido.search(
-        a.Time(obstime, obstime), a.Instrument.stix, a.Level.l2, a.stix.DataType.aux, a.stix.DataProduct.aux_ephemeris
+        a.Time(start_time, end_time), a.Instrument.stix, a.Level.l2, a.stix.DataType.aux, a.stix.DataProduct.aux_ephemeris
     )
-    if len(query["stix"]) != 1:
-        raise ValueError("Exactly one AUX file should be found but %d were found.", len(query["stix"]))
-    logger.debug("Downloading AUX data")
-    files = Fido.fetch(query)
-    if len(files.errors) > 0:
+    if len(query['stix']) == 0:
+        raise ValueError(f'No STIX pointing data found for time range {start_time} to {end_time}.')
+
+    logger.debug(f"Downloading {len(query['stix'])} AUX files")
+    aux_files = Fido.fetch(query['stix'])
+    if len(aux_files.errors) > 0:
         raise ValueError("There were errors downloading the data.")
     # Read and extract data
     logger.debug("Loading and extracting AUX data")
-    hdu = fits.getheader(files[0], ext=0)
-    aux = QTable.read(files[0], hdu=2)
-    start_time = Time(hdu.get("DATE-BEG"))
-    rel_date = (obstime - start_time).to("s")
-    idx = np.argmin(np.abs(rel_date - aux["time"]))
-    sas_x, sas_y = aux[idx][["y_srf", "z_srf"]]
-    roll, pitch, yaw = aux[idx]["roll_angle_rpy"]
-    solo_position_heeq = aux[idx]["solo_loc_heeq_zxy"]
-    logger.debug(
-        "SAS: %s, %s, Orientation: %s, %s, %s, Position: %s",
-        sas_x,
-        sas_y,
-        roll,
-        yaw.to("arcsec"),
-        pitch.to("arcsec"),
-        solo_position_heeq,
-    )
 
-    # roll, pitch, yaw = np.mean(aux[1219:1222]['roll_angle_rpy'], axis=0)
-    # sas_x = np.mean(aux[1219:1222]['y_srf'])
-    # sas_y = np.mean(aux[1219:1222]['z_srf'])
+    aux_data = []
+    for aux_file in aux_files:
+        hdu = fits.getheader(aux_file, ext=0)
+        aux = QTable.read(aux_file, hdu=2)
+        date_beg = Time(hdu.get("DATE-BEG"))
+        aux['time'] = date_beg + aux['time'] - 32*u.s  # Shift AUX data by half a time bin (starting time vs. bin centre)
+        aux_data.append(aux)
+
+    aux = vstack(aux_data)
+    aux.sort(keys=['time'])
 
-    return roll, pitch, yaw, sas_x, sas_y, solo_position_heeq
+    return aux
 
 
 @frame_transform_graph.transform(coord.FunctionTransform, STIXImaging, Helioprojective)