From ae5c8b704a684cdb785e3d5e0e116293ee3ccfe1 Mon Sep 17 00:00:00 2001 From: Shane Maloney Date: Wed, 12 Jun 2024 10:56:37 +0100 Subject: [PATCH] Both mean and instantaneous versions working --- changelog/128.bugfix.rst | 2 +- stixpy/coordinates/frames.py | 30 ++++++++-- stixpy/coordinates/tests/test_frames.py | 17 ++++-- stixpy/coordinates/tests/test_transforms.py | 13 +++++ stixpy/coordinates/transforms.py | 65 ++++++--------------- stixpy/map/stix.py | 40 +++++++++++++ stixpy/timeseries/quicklook.py | 3 +- stixpy/timeseries/tests/test_quicklook.py | 27 ++++----- 8 files changed, 123 insertions(+), 74 deletions(-) diff --git a/changelog/128.bugfix.rst b/changelog/128.bugfix.rst index 23d6a61..4c55884 100644 --- a/changelog/128.bugfix.rst +++ b/changelog/128.bugfix.rst @@ -1 +1 @@ -Fix bug where coordinate transforms did not support arrays of coordinates or times. +Fix bug where coordinate transforms did not support arrays of coordinates or times. To so support both integrate (e.g. images) and instantaneous (e.g. coarse flare flag) a new property `obstime_end` was added to `~stixpy.coordinates.frames.STIXImaging` when set if possible mean pointing and position information are use for transformations other wise they are interpolated between the two nearest data points. diff --git a/stixpy/coordinates/frames.py b/stixpy/coordinates/frames.py index abee101..568207c 100644 --- a/stixpy/coordinates/frames.py +++ b/stixpy/coordinates/frames.py @@ -1,9 +1,11 @@ import astropy import astropy.coordinates as coord import astropy.units as u +import numpy as np from astropy.coordinates import QuantityAttribute +from astropy.time import Time from astropy.wcs import WCS -from sunpy.coordinates.frameattributes import ObserverCoordinateAttribute +from sunpy.coordinates.frameattributes import ObserverCoordinateAttribute, TimeFrameAttributeSunPy from sunpy.coordinates.frames import HeliographicStonyhurst, SunPyBaseCoordinateFrame from sunpy.sun.constants import radius as _RSUN @@ -53,6 +55,8 @@ class STIXImaging(SunPyBaseCoordinateFrame): """ + obstime_end = TimeFrameAttributeSunPy() + observer = ObserverCoordinateAttribute(HeliographicStonyhurst) rsun = QuantityAttribute(default=_RSUN, unit=u.km) @@ -68,6 +72,17 @@ class STIXImaging(SunPyBaseCoordinateFrame): ], } + @property + def obstime_avg(self): + r"""Average time of the observation 'mean(obstime, obstime_end)'.""" + return np.mean(Time([self.obstime, self.obstime_end])) + + def __init__(self, *args, **kwargs): + super().__init__(*args, **kwargs) + if self.obstime is not None and self.obstime_end is not None: + if self.obstime.shape != self.obstime_end.shape: + raise ValueError("Both obstime and obstime_end must be either scaler or 1d array os same size.") + def stix_wcs_to_frame(wcs): r""" @@ -89,7 +104,9 @@ def stix_wcs_to_frame(wcs): if set(wcs.wcs.ctype) != {STIX_X_CTYPE, STIX_Y_CTYPE}: return None - dateobs = wcs.wcs.dateobs + datebeg = wcs.wcs.datebeg + dateavg = wcs.wcs.dateavg + dateend = wcs.wcs.dateend rsun = wcs.wcs.aux.rsun_ref if rsun is not None: @@ -100,10 +117,10 @@ def stix_wcs_to_frame(wcs): hgs_distance = wcs.wcs.aux.dsun_obs observer = HeliographicStonyhurst( - lat=hgs_latitude * u.deg, lon=hgs_longitude * u.deg, radius=hgs_distance * u.m, obstime=dateobs, rsun=rsun + lat=hgs_latitude * u.deg, lon=hgs_longitude * u.deg, radius=hgs_distance * u.m, obstime=dateavg, rsun=rsun ) - frame_args = {"obstime": dateobs, "observer": observer, "rsun": rsun} + frame_args = {"obstime": datebeg, "obstime_end": dateend, "observer": observer, "rsun": rsun} return STIXImaging(**frame_args) @@ -140,7 +157,10 @@ def stix_frame_to_wcs(frame, projection="TAN"): wcs.wcs.aux.hglt_obs = obs_frame.lat.to_value(u.deg) wcs.wcs.aux.dsun_obs = obs_frame.radius.to_value(u.m) - wcs.wcs.dateobs = frame.obstime.utc.iso + wcs.wcs.datebeg = frame.obstime.utc.iso + wcs.wcs.dateend = frame.obstime_end.utc.iso + wcs.wcs.dateavg = frame.obstime_avg.utc.iso + wcs.wcs.dateobs = frame.obstime_avg.utc.iso wcs.wcs.cunit = ["arcsec", "arcsec"] wcs.wcs.ctype = [STIX_X_CTYPE, STIX_Y_CTYPE] diff --git a/stixpy/coordinates/tests/test_frames.py b/stixpy/coordinates/tests/test_frames.py index af0ed48..1dcc42e 100644 --- a/stixpy/coordinates/tests/test_frames.py +++ b/stixpy/coordinates/tests/test_frames.py @@ -17,7 +17,9 @@ def stix_wcs(): w = WCS(naxis=2) - w.wcs.dateobs = "2024-01-01 00:00:00.000" + w.wcs.datebeg = "2024-01-01 00:00:00.000" + w.wcs.dateavg = "2024-01-01 00:00:01.000" + w.wcs.dateend = "2024-01-01 00:00:02.000" w.wcs.crpix = [10, 20] w.wcs.cdelt = np.array([2, 2]) w.wcs.crval = [0, 0] @@ -33,9 +35,10 @@ def stix_wcs(): @pytest.fixture def stix_frame(): obstime = "2024-01-01" + obstime_end = "2024-01-01 00:00:02.000" observer = HeliographicStonyhurst(10 * u.deg, 20 * u.deg, 1.5e11 * u.m, obstime=obstime) - frame_args = {"obstime": obstime, "observer": observer, "rsun": 695_700_000 * u.m} + frame_args = {"obstime": obstime, "observer": observer, "obstime_end": obstime_end, "rsun": 695_700_000 * u.m} frame = STIXImaging(**frame_args) return frame @@ -48,7 +51,7 @@ def test_stix_wcs_to_frame(stix_wcs): assert frame.obstime.isot == "2024-01-01T00:00:00.000" assert frame.rsun == 695700 * u.km assert frame.observer == HeliographicStonyhurst( - 10 * u.deg, 20 * u.deg, 1.5e11 * u.m, obstime="2024-01-01T00:00:00.000" + 10 * u.deg, 20 * u.deg, 1.5e11 * u.m, obstime="2024-01-01T00:00:01.000" ) @@ -66,7 +69,9 @@ def test_stix_frame_to_wcs(stix_frame): assert isinstance(wcs, WCS) assert wcs.wcs.ctype[0] == "SXLN-TAN" assert wcs.wcs.cunit[0] == "arcsec" - assert wcs.wcs.dateobs == "2024-01-01 00:00:00.000" + assert wcs.wcs.datebeg == "2024-01-01 00:00:00.000" + assert wcs.wcs.dateavg == "2024-01-01 00:00:01.000" + assert wcs.wcs.dateend == "2024-01-01 00:00:02.000" assert wcs.wcs.aux.rsun_ref == stix_frame.rsun.to_value(u.m) assert wcs.wcs.aux.dsun_obs == 1.5e11 @@ -84,7 +89,9 @@ def test_stix_frame_map(): data = np.random.rand(512, 512) obstime = "2023-01-01 12:00:00" solo = get_horizons_coord("solo", time=obstime) - coord = SkyCoord(0 * u.arcsec, 0 * u.arcsec, obstime=obstime, observer=solo, frame=STIXImaging) + coord = SkyCoord( + 0 * u.arcsec, 0 * u.arcsec, obstime=obstime, obstime_end="2023-01-01 12:02:00", observer=solo, frame=STIXImaging + ) header = make_fitswcs_header( data, coord, scale=[8, 8] * u.arcsec / u.pix, telescope="STIX", instrument="STIX", observatory="Solar Orbiter" ) diff --git a/stixpy/coordinates/tests/test_transforms.py b/stixpy/coordinates/tests/test_transforms.py index 2414f86..7c0072a 100644 --- a/stixpy/coordinates/tests/test_transforms.py +++ b/stixpy/coordinates/tests/test_transforms.py @@ -54,6 +54,19 @@ def test_stx_to_hpc_times(): assert np.all(stix_coord.obstime.isclose(stix_coord_rt.obstime)) +def test_stx_to_hpc_obstime_end(): + times = Time("2023-01-01") + [0, 2] * u.min + time_avg = np.mean(times) + stix_coord = SkyCoord(0 * u.deg, 0 * u.deg, frame=STIXImaging(obstime=times[0], obstime_end=times[-1])) + hp = stix_coord.transform_to(Helioprojective(obstime=time_avg)) + stix_coord_rt = hp.transform_to(STIXImaging(obstime=times[0], obstime_end=times[-1])) + + stix_coord_rt_interp = hp.transform_to(STIXImaging(obstime=times[1])) # noqa: F841 + 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)) + assert np.all(stix_coord.obstime_end.isclose(stix_coord_rt.obstime_end)) + + @pytest.mark.remote_data def test_get_aux_data(): with pytest.raises(ValueError, match="No STIX pointing data found for time range"): diff --git a/stixpy/coordinates/transforms.py b/stixpy/coordinates/transforms.py index 4c4aeaf..b3d02d0 100644 --- a/stixpy/coordinates/transforms.py +++ b/stixpy/coordinates/transforms.py @@ -26,7 +26,7 @@ __all__ = ["get_hpc_info", "stixim_to_hpc", "hpc_to_stixim"] -def _get_rotation_matrix_and_position(obstime): +def _get_rotation_matrix_and_position(obstime, obstime_end=None): r""" Return rotation matrix STIX Imaging to SOLO HPC and position of SOLO in HEEQ. @@ -37,7 +37,7 @@ def _get_rotation_matrix_and_position(obstime): Returns ------- """ - roll, solo_position_heeq, spacecraft_pointing = get_hpc_info(obstime) + roll, solo_position_heeq, spacecraft_pointing = get_hpc_info(obstime, obstime_end) # Generate the rotation matrix using the x-convention (see Goldstein) # Rotate +90 clockwise around the first axis @@ -126,42 +126,6 @@ def get_hpc_info(times, end_time=None): good_solution = np.where(aux[indices]["sas_ok"] == 1) good_sas = aux[good_solution] - # 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: - # 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) @@ -246,8 +210,13 @@ 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.T, representation_type="cartesian", obstime=stxcoord.obstime) + rot_matrix, solo_pos_heeq = _get_rotation_matrix_and_position(stxcoord.obstime, stxcoord.obstime_end) + + obstime = stxcoord.obstime + if stxcoord.obstime_end is not None: + obstime = np.mean(Time([stxcoord.obstime, stxcoord.obstime_end])) + + solo_heeq = HeliographicStonyhurst(solo_pos_heeq.T, representation_type="cartesian", obstime=obstime) # Transform from STIX imaging to SOLO HPC newrepr = stxcoord.cartesian.transform(rot_matrix) @@ -255,8 +224,8 @@ def stixim_to_hpc(stxcoord, hpcframe): # Create SOLO HPC solo_hpc = Helioprojective( newrepr, - obstime=stxcoord.obstime, - observer=solo_heeq.transform_to(HeliographicStonyhurst(obstime=stxcoord.obstime)), + obstime=obstime, + observer=solo_heeq.transform_to(HeliographicStonyhurst(obstime=obstime)), ) logger.debug("SOLO HPC: %s", solo_hpc) @@ -274,17 +243,21 @@ def hpc_to_stixim(hpccoord, stxframe): Transform HPC coordinate to STIX Imaging frame. """ logger.debug("Input HPC: %s", hpccoord) - rmatrix, solo_pos_heeq = _get_rotation_matrix_and_position(stxframe.obstime) - solo_hgs = HeliographicStonyhurst(solo_pos_heeq.T, representation_type="cartesian", obstime=stxframe.obstime) + rmatrix, solo_pos_heeq = _get_rotation_matrix_and_position(stxframe.obstime, stxframe.obstime_end) + + obstime = stxframe.obstime + if stxframe.obstime_end is not None: + obstime = np.mean(Time([stxframe.obstime, stxframe.obstime_end])) + solo_hgs = HeliographicStonyhurst(solo_pos_heeq.T, representation_type="cartesian", obstime=obstime) # Create SOLO HPC - solo_hpc_frame = Helioprojective(obstime=stxframe.obstime, observer=solo_hgs) + solo_hpc_frame = Helioprojective(obstime=obstime, observer=solo_hgs) # Transform input to SOLO HPC solo_hpc_coord = hpccoord.transform_to(solo_hpc_frame) logger.debug("SOLO HPC: %s", solo_hpc_coord) # Transform from SOLO HPC to STIX imaging newrepr = solo_hpc_coord.cartesian.transform(matrix_transpose(rmatrix)) - stx_coord = stxframe.realize_frame(newrepr, obstime=solo_hpc_frame.obstime) + stx_coord = stxframe.realize_frame(newrepr) logger.debug("STIX: %s", newrepr) return stx_coord diff --git a/stixpy/map/stix.py b/stixpy/map/stix.py index db3e553..37f1766 100644 --- a/stixpy/map/stix.py +++ b/stixpy/map/stix.py @@ -1,4 +1,5 @@ import astropy.units as u +from astropy import wcs from sunpy.map import GenericMap __all__ = ["STIXMap"] @@ -27,3 +28,42 @@ def plot(self, **kwargs): res.axes.coords[1].set_axislabel(r"STIX $\Theta_{Y}$") return res + + @property + def wcs(self): + """ + The `~astropy.wcs.WCS` property of the map. + """ + w2 = wcs.WCS(naxis=2) + + # Add one to go from zero-based to one-based indexing + w2.wcs.crpix = u.Quantity(self.reference_pixel) + 1 * u.pix + # Make these a quantity array to prevent the numpy setting element of + # array with sequence error. + # Explicitly call ``.to()`` to check that scale is in the correct units + w2.wcs.cdelt = u.Quantity( + [self.scale[0].to(self.spatial_units[0] / u.pix), self.scale[1].to(self.spatial_units[1] / u.pix)] + ) + w2.wcs.crval = u.Quantity([self._reference_longitude, self._reference_latitude]) + w2.wcs.ctype = self.coordinate_system + w2.wcs.pc = self.rotation_matrix + w2.wcs.set_pv(self._pv_values) + # FITS standard doesn't allow both PC_ij *and* CROTA keywords + w2.wcs.crota = (0, 0) + w2.wcs.cunit = self.spatial_units + + w2.wcs.datebeg = self.date_start.isot + w2.wcs.dateavg = self.date_average.isot + w2.wcs.dateend = self.date_end.isot + + w2.wcs.aux.rsun_ref = self.rsun_meters.to_value(u.m) + w2.wcs.aux.hgln_obs = self.observer_coordinate.lon.to_value(u.deg) + w2.wcs.aux.hglt_obs = self.observer_coordinate.lat.to_value(u.deg) + w2.wcs.aux.dsun_obs = self.observer_coordinate.radius.to_value(u.m) + + # Set the shape of the data array + w2.array_shape = self.data.shape + + # Validate the WCS here. + w2.wcs.set() + return w2 diff --git a/stixpy/timeseries/quicklook.py b/stixpy/timeseries/quicklook.py index b837693..373dd84 100644 --- a/stixpy/timeseries/quicklook.py +++ b/stixpy/timeseries/quicklook.py @@ -519,8 +519,7 @@ class HKMaxi(GenericTimeSeries): >>> from stixpy.data import test >>> from sunpy.timeseries import TimeSeries >>> from stixpy.timeseries.quicklook import HKMaxi - >>> hk = TimeSeries("https://pub099.cs.technik.fhnw.ch/fits/L1/2020/05/06/HK/" - ... "solo_L1_stix-hk-maxi_20200506_V02U.fits") # doctest: +REMOTE_DATA + >>> hk = TimeSeries("https://pub099.cs.technik.fhnw.ch/fits/L1/2020/05/06/HK/solo_L1_stix-hk-maxi_20200506_V02.fits") # doctest: +REMOTE_DATA >>> hk # doctest: +SKIP """ diff --git a/stixpy/timeseries/tests/test_quicklook.py b/stixpy/timeseries/tests/test_quicklook.py index 993236b..f10bb23 100644 --- a/stixpy/timeseries/tests/test_quicklook.py +++ b/stixpy/timeseries/tests/test_quicklook.py @@ -32,13 +32,13 @@ def ql_variance(): return ql_var -# @pytest.mark.remote_data -# @pytest.fixture() -# def hk_maxi(): -# hk_maxi = TimeSeries( -# r"https://pub099.cs.technik.fhnw.ch/fits/L1/2020/05/06/HK/solo_L1_stix-hk-maxi_20200506_V02U.fits" -# ) -# return hk_maxi +@pytest.mark.remote_data +@pytest.fixture() +def hk_maxi(): + hk_maxi = TimeSeries( + r"https://pub099.cs.technik.fhnw.ch/fits/L1/2020/05/06/HK/solo_L1_stix-hk-maxi_20200506_V02.fits" + ) + return hk_maxi @pytest.mark.remote_data @@ -78,15 +78,12 @@ def test_ql_variance_plot(ql_variance): @pytest.mark.remote_data -def test_hk_maxi(): - hk_maxi = TimeSeries( - "https://pub099.cs.technik.fhnw.ch/fits/L1/2020/05/06/HK/solo_L1_stix-hk-maxi_20200506_V02U.fits" - ) +def test_hk_maxi(hk_maxi): assert isinstance(hk_maxi, HKMaxi) assert hk_maxi.quantity("hk_att_c").unit == u.mA -# @pytest.mark.remote_data -# def test_hk_maxi_plot(hk_maxi): -# hk_maxi.plot() -# hk_maxi.plot(columns=['hk_asp_photoa0_v']) +@pytest.mark.remote_data +def test_hk_maxi_plot(hk_maxi): + hk_maxi.plot() + hk_maxi.plot(columns=["hk_asp_photoa0_v"])