Both mean and instantaneous versions working

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.

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]
 

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"
     )

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"):

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,17 +210,22 @@ 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)
 
     # 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

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

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
     <stixpy.timeseries.quicklook.HKMaxi ...>
     """

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"])