From 3547ab23386aa33bde4af9802e400252bbd5f8d0 Mon Sep 17 00:00:00 2001 From: Shane Maloney Date: Sat, 8 Jun 2024 21:01:28 +0100 Subject: [PATCH] Working test for get_hpc_info --- stixpy/coordinates/tests/test_transforms.py | 38 ++++++++++----------- stixpy/coordinates/transforms.py | 32 +++++++++-------- 2 files changed, 37 insertions(+), 33 deletions(-) diff --git a/stixpy/coordinates/tests/test_transforms.py b/stixpy/coordinates/tests/test_transforms.py index 9db51c1..5b5389f 100644 --- a/stixpy/coordinates/tests/test_transforms.py +++ b/stixpy/coordinates/tests/test_transforms.py @@ -5,7 +5,6 @@ import pytest from astropy.tests.helper import assert_quantity_allclose from astropy.time import Time -from numpy.testing import assert_allclose from sunpy.coordinates import HeliographicStonyhurst, Helioprojective from stixpy.coordinates.frames import STIXImaging @@ -63,20 +62,26 @@ 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:16", "2022-08-28T16:00:17.000"])) + roll, solo_heeq, stix_pointing = get_hpc_info(Time(["2022-08-28T16:00:17", "2022-08-28T16:00:18.000"])) - assert_quantity_allclose(-3.3733306 * u.deg, roll) - assert_quantity_allclose([25.61525646, 57.914266] * u.arcsec, stix_pointing) - assert_quantity_allclose([-97868816.0, 62984852.0, -5531986.5] * u.km, solo_heeq) + 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) 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.7867768e07, 6.2983744e07, -5532067.5] * u.km, solo_heeq) + assert_quantity_allclose( + roll, + -3.3732104 * u.deg, + ) + assert_quantity_allclose( + stix_pointing[0, :], + [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(-3.3619127 * u.deg, roll) - assert_quantity_allclose([-26.070198, 173.48871] * u.arcsec, stix_pointing) + 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) @@ -84,13 +89,8 @@ def test_get_hpc_info(): 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[[0, -1]]) - - roll3, solo_heeq3, stix_pointing3 = get_hpc_info(t[5]) - - assert_allclose(np.mean(roll1, axis=0), roll2, rtol=1e-6) - assert_allclose(np.mean(roll1, axis=0), roll3, rtol=5e-5) - - assert_allclose(np.mean(solo_heeq1, axis=0), solo_heeq2.value, rtol=1e-6) - assert_allclose(np.mean(solo_heeq1, axis=0), solo_heeq3[0, :], rtol=1e-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, :]) diff --git a/stixpy/coordinates/transforms.py b/stixpy/coordinates/transforms.py index c5d909e..4bbd990 100644 --- a/stixpy/coordinates/transforms.py +++ b/stixpy/coordinates/transforms.py @@ -55,7 +55,7 @@ def _get_rotation_matrix_and_position(obstime): return rmatrix, solo_position_heeq -def get_hpc_info(times): +def get_hpc_info(times, end_time=None): r""" Get STIX pointing and SO location from L2 aspect files. @@ -71,9 +71,11 @@ def get_hpc_info(times): aux = _get_aux_data(times.min(), times.max()) indices = np.argwhere((aux["time"] >= times.min()) & (aux["time"] <= times.max())) + if end_time is not None: + indices = np.argwhere((aux["time"] >= times.min()) & (aux["time"] <= end_time)) indices = indices.flatten() - if times.size == 2 and indices.size >= 2: + if end_time is not None and times.size == 1 and indices.size >= 2: # mean aux = aux[indices] @@ -84,7 +86,7 @@ def get_hpc_info(times): good_sas = aux[good_solution] if len(good_sas) == 0: - warnings.warn(f"No good SAS solution found for time range: {times[0]} to {times[1]}.") + warnings.warn(f"No good SAS solution found for time range: {times} to {end_time}.") sas_x = 0 sas_y = 0 else: @@ -97,7 +99,7 @@ def get_hpc_info(times): warnings.warn(f"Pointing unstable: StD(X) = {sigma_x}, StD(Y) = {sigma_y}.") else: if indices.size < 2: - logger.info("Only one data point found interpolating between two closest times.") + 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 @@ -113,13 +115,13 @@ def get_hpc_info(times): aux = aux[start_ind:end_ind] indices = [0, 1] - # Interpolate + # Interpolate all times x = (times - aux["time"][0]).to_value(u.s) xp = (aux["time"] - aux["time"][0]).to_value(u.s) roll = np.interp(x, xp, aux["roll_angle_rpy"][:, 0].value) << aux["roll_angle_rpy"].unit - yaw = np.interp(x, xp, aux["roll_angle_rpy"][:, 1].value) << aux["roll_angle_rpy"].unit - pitch = np.interp(x, xp, aux["roll_angle_rpy"][:, 2].value) << aux["roll_angle_rpy"].unit + pitch = np.interp(x, xp, aux["roll_angle_rpy"][:, 1].value) << aux["roll_angle_rpy"].unit + yaw = np.interp(x, xp, aux["roll_angle_rpy"][:, 2].value) << aux["roll_angle_rpy"].unit solo_heeq = ( np.vstack( @@ -177,15 +179,15 @@ def get_hpc_info(times): # 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]) + spacecraft_pointing = np.vstack([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]) + sas_pointing = np.vstack([sas_x + STIX_X_OFFSET, -1 * sas_y + STIX_Y_OFFSET]) - pointing_diff = np.linalg.norm(spacecraft_pointing - sas_pointing) + 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: - if pointing_diff < 200 * u.arcsec: + if np.max(pointing_diff) < 200 * u.arcsec: logger.info(f"Using SAS pointing: {sas_pointing}") - stix_pointing = sas_pointing + stix_pointing = np.where(pointing_diff < 200 * u.arcsec, sas_pointing, spacecraft_pointing) else: warnings.warn( f"Using spacecraft pointing: {spacecraft_pointing}" f" large difference between SAS and spacecraft." @@ -193,11 +195,11 @@ def get_hpc_info(times): else: warnings.warn(f"SAS solution not available using spacecraft pointing: {stix_pointing}.") - return roll, solo_heeq, stix_pointing + return roll, solo_heeq, stix_pointing.T @lru_cache -def _get_aux_data(start_time, end_time): +def _get_aux_data(start_time, end_time=None): r""" Search, download and read L2 pointing data. @@ -210,6 +212,8 @@ def _get_aux_data(start_time, end_time): ------- """ + if end_time is None: + end_time = start_time # Find, download, read aux file with pointing, sas and position information logger.debug(f"Searching for AUX data: {start_time} - {end_time}") query = Fido.search(