Feat imaging maps (#107)

* Make STIXFrame work nicely with `Map`.
* Apply suggestions from code review
* Fix docs and STIX map axis labels
* Change logs and docs

---------
Co-authored-by: DanRyanIrish <ryand5@tcd.ie>

changelog/107.bugfix.rst

@@ -0,0 +1 @@
+Update `~stixpy.calibration.energy.get_elut` to use correct science energy channels for given date.

changelog/107.feature.rst

@@ -0,0 +1 @@
+Update imaging demo to use images and construct `~sunpy.map.Map` using `~stixpy.coordinates.frames.STIXImaging` and `~sunpy.map.header_helper.make_fitswcs_header`.

examples/imaging_demo.py

@@ -16,6 +16,9 @@
 import astropy.units as u
 import matplotlib.pyplot as plt
 import numpy as np
+from astropy.time import Time
+from sunpy.coordinates import HeliographicStonyhurst, Helioprojective
+from sunpy.map import Map, make_fitswcs_header
 from xrayvision.clean import vis_clean
 from xrayvision.imaging import vis_to_image, vis_to_map
 from xrayvision.mem import mem
@@ -27,11 +30,13 @@
     create_visibility,
     get_uv_points_data,
 )
+from stixpy.coordinates.frames import STIXImaging
+from stixpy.coordinates.transforms import get_hpc_info
 from stixpy.imaging.em import em
+from stixpy.map.stix import STIXMap  # noqa
 from stixpy.product import Product
 
 logger = logging.getLogger(__name__)
-logger.setLevel("DEBUG")
 
 #############################################################################
 # Read science file as Product
@@ -114,10 +119,31 @@
 ###############################################################################
 # Make a full disk back projection (inverse transform) map
 
-fd_bp_map = vis_to_map(stix_vis, imsize, pixel_size=pixel)
-fd_bp_map.meta["rsun_obs"] = cpd_sci.meta["rsun_arc"]
+bp_image = vis_to_image(stix_vis, imsize, pixel_size=pixel)
+
+date_avg = Time('2021-09-23T15:22:30')
+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)
+header = make_fitswcs_header(bp_image, coord, telescope='STIX', observatory='Solar Orbiter', scale=[10,10]*u.arcsec/u.pix)
+fd_bp_map = Map((bp_image, header))
+
+hpc_ref = Helioprojective(pointing[0], pointing[1], observer=solo, obstime=fd_bp_map.date)
+header_hp = make_fitswcs_header(bp_image, hpc_ref, scale=[10, 10]*u.arcsec/u.pix, rotation_angle=90*u.deg+roll)
+hp_map = Map((bp_image, header_hp))
+hp_map_rotated = hp_map.rotate()
+
+fig = plt.figure(figsize=(12, 8))
+ax0 = fig.add_subplot(1, 2, 1, projection=fd_bp_map)
+ax1 = fig.add_subplot(1, 2, 2, projection=hp_map_rotated)
+fd_bp_map.plot(axes=ax0)
 fd_bp_map.draw_limb()
-fd_bp_map.plot()
+fd_bp_map.draw_grid()
+
+hp_map_rotated.plot(axes=ax1)
+hp_map_rotated.draw_limb()
+hp_map_rotated.draw_grid()
 
 ###############################################################################
 # Estimate the flare location and plot on top of back projection map.
@@ -126,14 +152,9 @@
 # because WCS axes are reverse order
 max_hpc = fd_bp_map.pixel_to_world(max_pixel[1], max_pixel[0])
 
-fig = plt.figure()
-axes = fig.add_subplot(projection=fd_bp_map)
-fd_bp_map.plot(axes=axes)
-fd_bp_map.draw_limb()
-axes.plot_coord(max_hpc, marker=".", markersize=50, fillstyle="none", color="r", markeredgewidth=2)
-axes.set_xlabel("STIX Y")
-axes.set_ylabel("STIX X")
-axes.set_title(f'STIX {" ".join(time_range_sci)} {"-".join([str(e) for e in energy_range])} keV')
+ax0.plot_coord(max_hpc, marker=".", markersize=50, fillstyle="none", color="r", markeredgewidth=2)
+ax1.plot_coord(max_hpc, marker=".", markersize=50, fillstyle="none", color="r", markeredgewidth=2)
+
 
 ################################################################################
 # Use estimated flare location to create more accurate visibilities
@@ -182,13 +203,11 @@
 # Create a back projection image with natural weighting
 
 bp_nat = vis_to_image(stix_vis1, imsize, pixel_size=pixel)
-plt.imshow(bp_nat.value)
 
 ###############################################################################
 # Create a back projection image with uniform weighting
 
 bp_uni = vis_to_image(stix_vis1, imsize, pixel_size=pixel, natural=False)
-plt.imshow(bp_uni.value)
 
 ###############################################################################
 # Create a `sunpy.map.Map` with back projection

stixpy/calibration/detector.py

@@ -56,11 +56,12 @@ def get_pixel_srm():
 
 def get_sci_channels(date):
     r"""
-    Get the science energy channel info for given date
+    Get the science energy channels for given date.
 
     Parameters
     ----------
     date : `datetime.datetime`
+        Date to lookup science energy channels.
 
     Returns
     -------

stixpy/calibration/energy.py

@@ -2,13 +2,29 @@
 
 import numpy as np
 
+from stixpy.calibration.detector import get_sci_channels
 from stixpy.io.readers import read_elut, read_elut_index
 from stixpy.product import Product
 from stixpy.utils.rebining import rebin_proportional
 
 __all__ = ['get_elut','correct_counts']
 
+
 def get_elut(date):
+    r"""
+    Get the energy lookup table (ELUT) for the given date
+
+    Combines the ELUT with the science energy channels for the same date.
+
+    Parameters
+    ----------
+    date `astropy.time.Time`
+        Date to look up the ELUT.
+
+    Returns
+    -------
+
+    """
     root = Path(__file__).parent.parent
     elut_index_file = Path(root, *["config", "data", "elut", "elut_index.csv"])
 
@@ -19,7 +35,8 @@ def get_elut(date):
     elif len(elut_info) > 1:
         raise ValueError(f"Multiple ELUTs for for date {date}")
     start_date, end_date, elut_file = list(elut_info)[0]
-    elut_table = read_elut(elut_file)
+    sci_channels = get_sci_channels(date)
+    elut_table = read_elut(elut_file, sci_channels)
 
     return elut_table
 

stixpy/coordinates/frames.py

@@ -10,7 +10,7 @@
 from stixpy.net.client import STIXClient  # noqa
 from stixpy.utils.logging import get_logger
 
-logger = get_logger(__name__, "DEBUG")
+logger = get_logger(__name__)
 
 __all__ = ["STIXImaging"]
 

stixpy/coordinates/transforms.py

@@ -1,18 +1,21 @@
 import warnings
 
+import astropy.coordinates as coord
+import astropy.units as u
 import numpy as np
-from astropy import coordinates as coord
-from astropy import units as u
 from astropy.coordinates import frame_transform_graph
-from astropy.coordinates.matrix_utilities import matrix_product, matrix_transpose, rotation_matrix
+from astropy.coordinates.matrix_utilities import matrix_transpose, rotation_matrix
 from astropy.io import fits
 from astropy.table import QTable
 from astropy.time import Time
 from sunpy.coordinates import HeliographicStonyhurst, Helioprojective
 from sunpy.net import Fido
 from sunpy.net import attrs as a
 
-from stixpy.coordinates.frames import STIX_X_OFFSET, STIX_X_SHIFT, STIX_Y_OFFSET, STIX_Y_SHIFT, STIXImaging, logger
+from stixpy.coordinates.frames import STIX_X_OFFSET, STIX_X_SHIFT, STIX_Y_OFFSET, STIX_Y_SHIFT, STIXImaging
+from stixpy.utils.logging import get_logger
+
+logger = get_logger(__name__)
 
 __all__ = ['get_hpc_info', 'stixim_to_hpc', 'hpc_to_stixim']
 
@@ -41,7 +44,7 @@ def _get_rotation_matrix_and_position(obstime):
     C = rotation_matrix(-1 * spacecraft_pointing[0], "z")
 
     # Will be applied right to left
-    rmatrix = matrix_product(A, B, C, rot_to_solo)
+    rmatrix = A @ B @ C @ rot_to_solo
 
     return rmatrix, solo_position_heeq
 
@@ -124,7 +127,8 @@ def stixim_to_hpc(stxcoord, hpcframe):
 
     # Create SOLO HPC
     solo_hpc = Helioprojective(
-        newrepr, obstime=stxcoord.obstime, observer=solo_heeq.transform_to(HeliographicStonyhurst)
+        newrepr, obstime=stxcoord.obstime,
+        observer=solo_heeq.transform_to(HeliographicStonyhurst(obstime=stxcoord.obstime))
     )
     logger.debug("SOLO HPC: %s", solo_hpc)
 
@@ -154,5 +158,5 @@ def hpc_to_stixim(hpccoord, stxframe):
     # 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)
-    logger.debug("STIX: %s", stx_coord)
+    logger.debug("STIX: %s", newrepr)
     return stx_coord

stixpy/io/readers.py

@@ -93,11 +93,7 @@ def read_elut_index(elut_index):
     return elut_it
 
 
-def read_elut(elut_file):
-    root = Path(__file__).parent.parent
-    sci_channels = read_sci_energy_channels(
-        Path(root, *["config", "data", "detector", "ScienceEnergyChannels_1000.csv"])
-    )
+def read_elut(elut_file, sci_channels):
     elut_table = Table.read(elut_file, header_start=2)
 
     elut = type("ELUT", (object,), dict())
@@ -112,7 +108,10 @@ def read_elut(elut_file):
     elut.e_actual = (elut.adc - elut.offset[..., None]) * elut.gain[..., None]
     elut.e_width_actual = elut.e_actual[..., 1:] - elut.e_actual[..., :-1]
 
+    # TODO remove after updating elut correction
     elut.e = sci_channels["Elower"]
+    elut.e_sci_low = sci_channels["Elower"]
+    elut.e_sci_high = sci_channels["Eupper"]
     elut.e_width = sci_channels["Eupper"] - sci_channels["Elower"]
 
     return elut

stixpy/io/tests/test_readers.py

@@ -33,7 +33,10 @@ def test_read_elut_index(moc_table):
 def test_read_elut():
     root = Path(__file__).parent.parent.parent
     elut_file = Path(root, *["config", "data", "elut", "elut_table_20201204.csv"])
-    elut = read_elut(elut_file)
+    sci_file = root / 'config' / 'data' / 'detector' / 'ScienceEnergyChannels_1000.csv'
+
+    sci_channels = read_sci_energy_channels(sci_file)
+    elut = read_elut(elut_file, sci_channels)
 
     # make sure a few values are as expected
     assert elut.offset[0, 0] == 897.7778

stixpy/map/stix.py

@@ -17,7 +17,7 @@ def plot(self, **kwargs):
 
         # Until this is resolved need to manually override https://github.com/astropy/astropy/issues/16339
         res = super().plot(**kwargs)
-        if self.coordinate_system == 'stiximaging':
+        if self.coordinate_frame.name == 'stiximaging':
             res.axes.coords[0].set_coord_type('longitude', coord_wrap=180*u.deg)
             res.axes.coords[1].set_coord_type('latitude')