Working end to end

changelog/90.feature.rst

@@ -0,0 +1,2 @@
+Add a custom coordinate frame :class:`stixpy.coordinate.frames.STIXImaging` and the associated transformations (`~stixpy.coordinates.transforms.stixim_to_hpc`, ~stixpy.coordinates.transforms.hpc_to_stixim') to and from `~sunpy.coordinates.frames.Helioprojective`.
+Also add a `stixpy.map.stix.STIXMap` for `~sunpy.map.Map` source which properly displays the :class:`stixpy.coordinate.frames.STIXImaging` coordinate frame.

docs/reference/coordinates.rst

@@ -0,0 +1,12 @@
+.. _coordinates:
+
+Coordinates ('stixpy.coordinates')
+**********************************
+
+The ``coordinates`` submodule contains STIX specific coordinate frames and transforms.
+
+.. automodapi:: stixpy.coordinates
+
+.. automodapi:: stixpy.coordinates.frames
+
+.. automodapi:: stixpy.coordinates.transforms

docs/reference/index.rst

@@ -7,12 +7,15 @@ Reference
 .. toctree::
    :maxdepth: 2
 
+   coordinates
    data
    calibration
    timeseries
    science
+   map
    product
+   science
+   timeseries
    visualisation
-   frames
 
    ../whatsnew/index

docs/reference/map.rst

@@ -0,0 +1,10 @@
+.. _map:
+
+Map ('stixpy.map')
+******************
+
+The ``map`` submodule contains stix specific map source for `~sunpy.map.Map`
+
+.. automodapi:: stixpy.map
+
+.. automodapi:: stixpy.map.stix

setup.cfg

@@ -68,7 +68,7 @@ filterwarnings =
 
 [isort]
 balanced_wrapping = true
-skip =
+extend_skip =
     docs/conf.py,
     stixpy/__init__.py,
 default_section = THIRDPARTY

stixpy/coordinates/frames.py

@@ -0,0 +1,160 @@
+import astropy
+import astropy.coordinates as coord
+import astropy.units as u
+from astropy.coordinates import QuantityAttribute
+from astropy.wcs import WCS
+from sunpy.coordinates.frameattributes import ObserverCoordinateAttribute
+from sunpy.coordinates.frames import HeliographicStonyhurst, SunPyBaseCoordinateFrame
+from sunpy.sun.constants import radius as _RSUN
+
+from stixpy.net.client import STIXClient  # noqa
+from stixpy.utils.logging import get_logger
+
+logger = get_logger(__name__, "DEBUG")
+
+__all__ = ["STIXImaging"]
+
+STIX_X_SHIFT = 26.1 * u.arcsec  # fall back to this when non sas solution available
+STIX_Y_SHIFT = 58.2 * u.arcsec  # fall back to this when non sas solution available
+STIX_X_OFFSET = 60.0 * u.arcsec  # remaining offset after SAS solution
+STIX_Y_OFFSET = 8.0 * u.arcsec  # remaining offset after SAS solution
+
+
+class STIXImaging(SunPyBaseCoordinateFrame):
+    r"""
+    STIX Imaging Frame
+
+    - ``x`` (aka "theta_x") along the pixel rows (e.g. 0, 1, 2, 3; 4, 5, 6, 8).
+    - ``y`` (aka "theta_y") along the pixel columns (e.g. A, B, C, D).
+    - ``distance`` is the Sun-observer distance.
+
+    Aligned with SIIX 'pixels' +X corresponds direction along pixel row toward pixels
+    0, 4 and +Y corresponds direction along columns towards pixels 0, 1, 2, 3.
+
+    .. code-block:: text
+
+      Col\Row
+            ---------
+           | 3  | 7  |
+        D  |   _|_   |
+           |  |11 |  |
+            ---------
+           | 2  | 6  |
+        C  |   _|_   |
+           |  |10 |  |
+            ---------
+           | 1  | 5  |
+        B  |   _|_   |
+           |  | 9 |  |
+            ---------
+           | 0  | 4  |
+        A  |   _|_   |
+           |  | 8 |  |
+            ---------
+
+    """
+    observer = ObserverCoordinateAttribute(HeliographicStonyhurst)
+    rsun = QuantityAttribute(default=_RSUN, unit=u.km)
+
+    frame_specific_representation_info = {
+        coord.SphericalRepresentation: [
+            coord.RepresentationMapping("lon", "Tx", u.arcsec),
+            coord.RepresentationMapping("lat", "Ty", u.arcsec),
+            coord.RepresentationMapping("distance", "distance"),
+        ],
+        coord.UnitSphericalRepresentation: [coord.RepresentationMapping('lon', 'Tx', u.arcsec),
+                                            coord.RepresentationMapping('lat', 'Ty', u.arcsec)]
+    }
+
+
+def stix_wcs_to_frame(wcs):
+    r"""
+    This function registers the coordinate frames to their FITS-WCS coordinate
+    type values in the `astropy.wcs.utils.wcs_to_celestial_frame` registry.
+
+    Parameters
+    ----------
+    wcs : astropy.wcs.WCS
+
+    Returns
+    -------
+    astropy.coordinates.BaseCoordinateFrame
+    """
+    if hasattr(wcs, "coordinate_frame"):
+        return wcs.coordinate_frame
+
+    # Not a STIX wcs bail out early
+    if set(wcs.wcs.ctype) != {'SXLN-TAN', 'SXLT-TAN'}:
+        return None
+
+    dateobs = wcs.wcs.dateobs
+
+    rsun = wcs.wcs.aux.rsun_ref
+    if rsun is not None:
+        rsun *= u.m
+
+    hgs_longitude = wcs.wcs.aux.hgln_obs
+    hgs_latitude = wcs.wcs.aux.hglt_obs
+    hgs_distance = wcs.wcs.aux.dsun_obs
+
+    observer = HeliographicStonyhurst(hgs_latitude * u.deg,
+                                      hgs_longitude * u.deg,
+                                      hgs_distance * u.m,
+                                      obstime=dateobs,
+                                      rsun=rsun)
+
+    frame_args = {'obstime': dateobs,
+                  'observer': observer,
+                  'rsun': rsun}
+
+    return STIXImaging(**frame_args)
+
+
+def stix_frame_to_wcs(frame, projection='TAN'):
+    r"""
+    For a given frame, this function returns the corresponding WCS object.
+    It registers the WCS coordinates types from their associated frame in the
+    `astropy.wcs.utils.celestial_frame_to_wcs` registry.
+
+    Parameters
+    ----------
+    frame : astropy.coordinates.BaseCoordinateFrame
+    projection : str, optional
+
+    Returns
+    -------
+    astropy.wcs.WCS
+    """
+    # Bail out early if not STIXImaging frame
+    if not isinstance(frame, STIXImaging):
+        return None
+
+    wcs = WCS(naxis=2)
+    wcs.wcs.aux.rsun_ref = frame.rsun.to_value(u.m)
+
+    # Sometimes obs_coord can be a SkyCoord, so convert down to a frame
+    obs_frame = frame.observer
+    if hasattr(obs_frame, 'frame'):
+        obs_frame = frame.observer.frame
+
+    wcs.wcs.aux.hgln_obs = obs_frame.lon.to_value(u.deg)
+    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.cunit = ['arcsec', 'arcsec']
+    wcs.wcs.ctype = ['SXLN-TAN', 'SXLT-TAN']
+
+    return wcs
+
+
+# Remove once min version of sunpy has https://github.com/sunpy/sunpy/pull/7594
+astropy.wcs.utils.WCS_FRAME_MAPPINGS.insert(1, [stix_wcs_to_frame])
+astropy.wcs.utils.FRAME_WCS_MAPPINGS.insert(1, [stix_frame_to_wcs])
+
+STIX_CTYPE_TO_UCD1 = {
+    "SXLT": "custom:pos.stiximaging.lat",
+    "SXLN": "custom:pos.stiximaging.lon",
+    "SXRZ": "custom:pos.stiximaging.z"
+}
+astropy.wcs.wcsapi.fitswcs.CTYPE_TO_UCD1.update(STIX_CTYPE_TO_UCD1)

stixpy/coordinates/tests/test_frames.py

@@ -0,0 +1,82 @@
+import astropy.units as u
+import numpy as np
+import pytest
+from astropy.coordinates import SkyCoord
+from astropy.wcs import WCS
+from sunpy.coordinates.frames import HeliographicStonyhurst, Helioprojective
+from sunpy.map import Map, make_fitswcs_header
+
+from stixpy.coordinates.frames import STIXImaging, stix_frame_to_wcs, stix_wcs_to_frame
+from stixpy.map.stix import STIXMap
+
+
+@pytest.fixture
+def stix_wcs():
+    w = WCS(naxis=2)
+
+    w.wcs.dateobs = '2024-01-01 00:00:00.000'
+    w.wcs.crpix = [10, 20]
+    w.wcs.cdelt = np.array([2, 2])
+    w.wcs.crval = [0, 0]
+    w.wcs.ctype = ["SXLN-TAN", "SXLT-TAN"]
+
+    w.wcs.aux.hgln_obs = 10
+    w.wcs.aux.hglt_obs = 20
+    w.wcs.aux.dsun_obs = 1.5e11
+
+    return w
+
+
+@pytest.fixture
+def stix_frame():
+    obstime = '2024-01-01'
+    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 = STIXImaging(**frame_args)
+    return frame
+
+def test_stix_wcs_to_frame(stix_wcs):
+    frame = stix_wcs_to_frame(stix_wcs)
+    assert isinstance(frame, STIXImaging)
+
+
+def test_stix_wcs_to_frame_none():
+    w = WCS(naxis=2)
+    w.wcs.ctype = ['ham', 'cheese']
+    frame = stix_wcs_to_frame(w)
+
+    assert frame is None
+
+
+def test_stix_frame_to_wcs(stix_frame):
+    wcs = 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.aux.rsun_ref == stix_frame.rsun.to_value(u.m)
+    assert wcs.wcs.aux.dsun_obs == 1.5e11
+    assert wcs.wcs.aux.hgln_obs == 10
+    assert wcs.wcs.aux.hglt_obs == 20
+
+
+def test_stix_frame_to_wcs_none():
+    wcs = stix_frame_to_wcs(Helioprojective())
+    assert wcs is None
+
+def test_stix_frame_map():
+    data = np.random.rand(10, 10)
+    coord = SkyCoord(0 * u.arcsec, 0 * u.arcsec, obstime='now', observer='earth', frame=STIXImaging)
+    header = make_fitswcs_header(data, coord, scale=[2, 2]*u.arcsec/u.pix, telescope='STIX',
+                                 instrument='STIX', observatory='Solar Orbiter')
+    smap = Map((data, header))
+    assert isinstance(smap, STIXMap)