Working

examples/imaging_demo.py

@@ -129,7 +129,7 @@
 roll, solo_xyz, pointing = get_hpc_info(date_avg)
 
 solo = HeliographicStonyhurst(*solo_xyz, obstime=date_avg, representation_type="cartesian")
-coord = STIXImaging(0 * u.arcsec, 0 * u.arcsec, obstime="2021-09-23T15:22:30", observer=solo)
+coord = STIXImaging(0 * u.arcsec, 0 * u.arcsec, obstime=date_avg, observer=solo)
 header = make_fitswcs_header(
     bp_image, coord, telescope="STIX", observatory="Solar Orbiter", scale=[10, 10] * u.arcsec / u.pix
 )

stixpy/coordinates/tests/test_transforms.py

@@ -3,6 +3,7 @@
 import astropy.units as u
 import numpy as np
 import pytest
+from astropy.coordinates import SkyCoord
 from astropy.tests.helper import assert_quantity_allclose
 from astropy.time import Time
 from sunpy.coordinates import HeliographicStonyhurst, Helioprojective
@@ -33,15 +34,24 @@ def test_transforms_with_know_values():
 def test_hpc_to_stx(mock):
     # fake some data to make check easier
     obstime = Time("2021-05-22T02:52:00")
-    roll = 0 * u.deg
-    sas = [10, 15] * 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, 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)
-    assert_quantity_allclose(stix_coord.Tx, -sas[1])
-    assert_quantity_allclose(stix_coord.Ty, sas[0])
+    assert_quantity_allclose(stix_coord.Tx, -sas[0, 1])
+    assert_quantity_allclose(stix_coord.Ty, sas[0, 0])
+
+
+def test_stx_to_hpc_times():
+    times = Time("2023-01-01") + np.arange(10) * u.min
+    stix_coord = SkyCoord([0] * 10 * u.deg, [0] * 10 * u.deg, frame=STIXImaging(obstime=times))
+    hp = stix_coord.transform_to(Helioprojective(obstime=times))
+    stix_coord_rt = hp.transform_to(STIXImaging(obstime=times))
+    assert_quantity_allclose(0 * u.deg, stix_coord.separation(stix_coord_rt), atol=1e-17 * u.deg)
+    assert np.all(stix_coord.obstime.isclose(stix_coord_rt.obstime))
 
 
 @pytest.mark.remote_data
@@ -62,35 +72,42 @@ def test_get_aux_data():
 @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:17", "2022-08-28T16:00:18.000"]))
+        roll, solo_heeq, stix_pointing = get_hpc_info(
+            Time("2022-08-28T16:00:17"), end_time=Time("2022-08-28T16:00:18.000")
+        )
 
     assert_quantity_allclose(roll, -3.3733292 * u.deg)
-    assert_quantity_allclose(stix_pointing[0, :], [25.622768, 57.920166] * u.arcsec)
-    assert_quantity_allclose(solo_heeq[0, :], [-97868803.82, 62984839.12, -5531987.445] * u.km)
+    assert_quantity_allclose(stix_pointing, [25.622768, 57.920166] * u.arcsec)
+    assert_quantity_allclose(solo_heeq, [-97868803.82, 62984839.12, -5531987.445] * u.km)
 
     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(
         roll,
         -3.3732104 * u.deg,
     )
     assert_quantity_allclose(
-        stix_pointing[0, :],
+        stix_pointing,
         [25.923784, 58.179676] * u.arcsec,
     )
     assert_quantity_allclose(solo_heeq, [-9.7867768e07, 6.2983744e07, -5532067.5] * u.km)
 
     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(roll, -3.3619127 * u.deg)
-    assert_quantity_allclose(stix_pointing[0, :], [-26.070198, 173.48871] * u.arcsec)
-    assert_quantity_allclose([-9.7671984e07, 6.2774768e07, -5547166.0] * u.km, solo_heeq)
+    assert_quantity_allclose(stix_pointing, [-26.070198, 173.48871] * u.arcsec)
+    assert_quantity_allclose(solo_heeq, [-9.7671984e07, 6.2774768e07, -5547166.0] * u.km)
 
 
 @pytest.mark.remote_data
 def test_get_hpc_info_shapes():
     t = Time("2022-08-28T16:00:00") + np.arange(10) * u.min
     roll1, solo_heeq1, stix_pointing1 = get_hpc_info(t)
-    roll2, solo_heeq2, stix_pointing2 = get_hpc_info(t[5])
+    roll2, solo_heeq2, stix_pointing2 = get_hpc_info(t[5:6])
+    roll3, solo_heeq3, stix_pointing3 = get_hpc_info(t[5])
 
     assert_quantity_allclose(roll1[5], roll2)
     assert_quantity_allclose(solo_heeq1[5, :], solo_heeq2[0, :])
     assert_quantity_allclose(stix_pointing1[5, :], stix_pointing2[0, :])
+
+    assert_quantity_allclose(roll3, roll2[0])
+    assert_quantity_allclose(solo_heeq3, solo_heeq2[0, :])
+    assert_quantity_allclose(stix_pointing3, stix_pointing2[0, :])

stixpy/coordinates/transforms.py

@@ -46,8 +46,8 @@ def _get_rotation_matrix_and_position(obstime):
 
     logger.debug("Pointing: %s.", spacecraft_pointing)
     A = rotation_matrix(-1 * roll, "x")
-    B = rotation_matrix(spacecraft_pointing[1], "y")
-    C = rotation_matrix(-1 * spacecraft_pointing[0], "z")
+    B = rotation_matrix(spacecraft_pointing[..., 1], "y")
+    C = rotation_matrix(-1 * spacecraft_pointing[..., 0], "z")
 
     # Will be applied right to left
     rmatrix = A @ B @ C @ rot_to_solo
@@ -98,22 +98,8 @@ def get_hpc_info(times, end_time=None):
             if sigma_x > tolerance or sigma_y > tolerance:
                 warnings.warn(f"Pointing unstable: StD(X) = {sigma_x}, StD(Y) = {sigma_y}.")
     else:
-        if indices.size < 2:
-            logger.info("Only one data contained in time interval found interpolating between two closest times.")
-            # Times contained in one time or only contains one time
-            if times.size == 2:
-                times = times[0] + (times[1] - times[0]) * 0.5
-            diff_center = (aux["time"] - times).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
-            aux = aux[start_ind:end_ind]
-            indices = [0, 1]
+        if end_time is not None and indices.size < 2:
+            times = times + (end_time - times) * 0.5
 
         # Interpolate all times
         x = (times - aux["time"][0]).to_value(u.s)
@@ -179,9 +165,9 @@ def get_hpc_info(times, end_time=None):
     # 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.vstack([STIX_X_SHIFT + rotated_yaw, STIX_Y_SHIFT + rotated_pitch])
+    spacecraft_pointing = np.vstack([STIX_X_SHIFT + rotated_yaw, STIX_Y_SHIFT + rotated_pitch]).T
     stix_pointing = spacecraft_pointing
-    sas_pointing = np.vstack([sas_x + STIX_X_OFFSET, -1 * sas_y + STIX_Y_OFFSET])
+    sas_pointing = np.vstack([sas_x + STIX_X_OFFSET, -1 * sas_y + STIX_Y_OFFSET]).T
 
     pointing_diff = np.sqrt(np.sum((spacecraft_pointing - sas_pointing) ** 2, axis=0))
     if np.all(np.isfinite(sas_pointing)) and len(good_sas) > 0:
@@ -195,7 +181,11 @@ def get_hpc_info(times, end_time=None):
     else:
         warnings.warn(f"SAS solution not available using spacecraft pointing: {stix_pointing}.")
 
-    return roll, solo_heeq, stix_pointing.T
+    if end_time is not None or times.ndim == 0:
+        solo_heeq = solo_heeq.squeeze()
+        stix_pointing = stix_pointing.squeeze()
+
+    return roll, solo_heeq, stix_pointing
 
 
 @lru_cache
@@ -257,7 +247,7 @@ def stixim_to_hpc(stxcoord, hpcframe):
     logger.debug("STIX: %s", stxcoord)
 
     rot_matrix, solo_pos_heeq = _get_rotation_matrix_and_position(stxcoord.obstime)
-    solo_heeq = HeliographicStonyhurst(solo_pos_heeq, representation_type="cartesian", obstime=stxcoord.obstime)
+    solo_heeq = HeliographicStonyhurst(solo_pos_heeq.T, representation_type="cartesian", obstime=stxcoord.obstime)
 
     # Transform from STIX imaging to SOLO HPC
     newrepr = stxcoord.cartesian.transform(rot_matrix)
@@ -285,7 +275,7 @@ def hpc_to_stixim(hpccoord, stxframe):
     """
     logger.debug("Input HPC: %s", hpccoord)
     rmatrix, solo_pos_heeq = _get_rotation_matrix_and_position(stxframe.obstime)
-    solo_hgs = HeliographicStonyhurst(solo_pos_heeq, representation_type="cartesian", obstime=stxframe.obstime)
+    solo_hgs = HeliographicStonyhurst(solo_pos_heeq.T, representation_type="cartesian", obstime=stxframe.obstime)
 
     # Create SOLO HPC
     solo_hpc_frame = Helioprojective(obstime=stxframe.obstime, observer=solo_hgs)