From 152d8d7a8a35ec211ad218c17cacc43ff193e099 Mon Sep 17 00:00:00 2001 From: Shane Maloney Date: Sun, 9 Jun 2024 18:22:12 +0100 Subject: [PATCH] Working --- examples/imaging_demo.py | 2 +- stixpy/coordinates/tests/test_transforms.py | 39 +++++++++++++++------ stixpy/coordinates/transforms.py | 36 +++++++------------ 3 files changed, 42 insertions(+), 35 deletions(-) diff --git a/examples/imaging_demo.py b/examples/imaging_demo.py index a4f47dd..0bac220 100644 --- a/examples/imaging_demo.py +++ b/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 ) diff --git a/stixpy/coordinates/tests/test_transforms.py b/stixpy/coordinates/tests/test_transforms.py index 5b5389f..2414f86 100644 --- a/stixpy/coordinates/tests/test_transforms.py +++ b/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,11 +72,13 @@ 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( @@ -74,23 +86,28 @@ def test_get_hpc_info(): -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, :]) diff --git a/stixpy/coordinates/transforms.py b/stixpy/coordinates/transforms.py index 4bbd990..4c4aeaf 100644 --- a/stixpy/coordinates/transforms.py +++ b/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)