Add generic imaging frame. (#76)

* Add generic projective frame.

* Update to allow observer information to be missing

---------

Co-authored-by: Shane Maloney <[email protected]>

changelog/76.feature.rst

@@ -0,0 +1 @@
+Add `~xrayvision.coordinates.frames.Projective` coordinate frame to represent generic observer based projective coordinate system.

xrayvision/coordinates/frames.py

@@ -0,0 +1,142 @@
+import astropy.coordinates
+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 HeliographicCarrington, HeliographicStonyhurst, SunPyBaseCoordinateFrame
+from sunpy.sun.constants import radius as _RSUN
+
+__all__ = ["Projective"]
+
+
+X_CTYPE = "PJLN"
+Y_CTYPE = "PJLT"
+
+
+class Projective(SunPyBaseCoordinateFrame):
+    """A generic projective coordinate frame for an arbitrary observer."""
+
+    observer = ObserverCoordinateAttribute(HeliographicStonyhurst)
+    rsun = QuantityAttribute(default=_RSUN, unit=u.km)
+    frame_specific_representation_info = {
+        astropy.coordinates.SphericalRepresentation: [
+            astropy.coordinates.RepresentationMapping("lon", "Tx", u.arcsec),
+            astropy.coordinates.RepresentationMapping("lat", "Ty", u.arcsec),
+            astropy.coordinates.RepresentationMapping("distance", "distance"),
+        ],
+        astropy.coordinates.UnitSphericalRepresentation: [
+            astropy.coordinates.RepresentationMapping("lon", "Tx", u.arcsec),
+            astropy.coordinates.RepresentationMapping("lat", "Ty", u.arcsec),
+        ],
+    }
+
+
+def projective_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
+    -------
+    : `Projective`
+    """
+    if hasattr(wcs, "coordinate_frame"):
+        return wcs.coordinate_frame
+
+    # Not a lat,lon coordinate system bail out early
+    if X_CTYPE not in wcs.wcs.ctype[0] or Y_CTYPE not in wcs.wcs.ctype[1]:
+        return None
+
+    dateavg = wcs.wcs.dateobs
+
+    # Get observer coordinate from the WCS auxiliary information
+    # Note: the order of the entries is important, as it determines which set
+    # of header keys is given priority below. Stonyhurst should usually be
+    # prioritized, as it is defined more consistently across implementations,
+    # and so it should occur before Carrington here.
+    required_attrs = {
+        HeliographicStonyhurst: ["hgln_obs", "hglt_obs", "dsun_obs"],
+        HeliographicCarrington: ["crln_obs", "hglt_obs", "dsun_obs"],
+    }
+
+    # Get rsun from the WCS auxiliary information
+    rsun = wcs.wcs.aux.rsun_ref
+    if rsun is not None:
+        rsun *= u.m
+
+    observer = None
+    for frame, attr_names in required_attrs.items():
+        attrs = [getattr(wcs.wcs.aux, attr_name) for attr_name in attr_names]
+        if all([attr is not None for attr in attrs]):
+            kwargs = {"obstime": dateavg}
+            if rsun is not None:
+                kwargs["rsun"] = rsun
+            if issubclass(frame, HeliographicCarrington):
+                kwargs["observer"] = "self"
+
+            observer = frame(attrs[0] * u.deg, attrs[1] * u.deg, attrs[2] * u.m, **kwargs)
+            break
+
+    frame_args = {"obstime": dateavg, "observer": observer, "rsun": rsun}
+
+    return Projective(**frame_args)
+
+
+def projective_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 : `Projective`
+    projection : `str`, optional
+
+    Returns
+    -------
+    `astropy.wcs.WCS`
+    """
+    # Bail out early if not Projective frame
+    if not isinstance(frame, Projective):
+        return None
+    else:
+        conjunction = "-"
+        ctype = [conjunction.join([X_CTYPE, projection]), conjunction.join([Y_CTYPE, projection])]
+        cunit = ["arcsec"] * 2
+
+    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
+
+    if obs_frame is not None:
+        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 = cunit
+    wcs.wcs.ctype = ctype
+
+    return wcs
+
+
+# Remove once min version of sunpy has https://github.com/sunpy/sunpy/pull/7594
+astropy.wcs.utils.WCS_FRAME_MAPPINGS.insert(1, [projective_wcs_to_frame])
+astropy.wcs.utils.FRAME_WCS_MAPPINGS.insert(1, [projective_frame_to_wcs])
+
+PROJECTIVE_CTYPE_TO_UCD1 = {
+    "PJLT": "custom:pos.projective.lat",
+    "PJLN": "custom:pos.projective.lon",
+    "PJRZ": "custom:pos.projective.z",
+}
+astropy.wcs.wcsapi.fitswcs.CTYPE_TO_UCD1.update(PROJECTIVE_CTYPE_TO_UCD1)

xrayvision/coordinates/tests/test_frames.py

@@ -0,0 +1,112 @@
+import astropy.units as u
+import numpy as np
+import pytest
+from astropy.wcs import WCS
+from sunpy.coordinates import HeliographicStonyhurst
+
+from xrayvision.coordinates.frames import Projective, projective_frame_to_wcs, projective_wcs_to_frame
+
+
+@pytest.fixture
+def projective_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 = ["PJLN-TAN", "PJLT-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 projective_wcs_no_observer():
+    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 = ["PJLN-TAN", "PJLT-TAN"]
+
+    return w
+
+
+@pytest.fixture
+def projective_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 = Projective(**frame_args)
+    return frame
+
+
+@pytest.fixture
+def projective_frame_no_obsever():
+    obstime = "2024-01-01"
+    frame_args = {"obstime": obstime, "rsun": 695_700_000 * u.m}
+    frame = Projective(**frame_args)
+    return frame
+
+
+def test_projective_wcs_to_frame(projective_wcs):
+    frame = projective_wcs_to_frame(projective_wcs)
+    assert isinstance(frame, Projective)
+
+    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"
+    )
+
+
+def test_projective_wcs_no_observer_to_frame(projective_wcs_no_observer):
+    frame = projective_wcs_to_frame(projective_wcs_no_observer)
+    assert isinstance(frame, Projective)
+
+    assert frame.obstime.isot == "2024-01-01T00:00:00.000"
+    assert frame.rsun == 695700 * u.km
+    assert frame.observer is None
+
+
+def test_projective_wcs_to_frame_none():
+    w = WCS(naxis=2)
+    w.wcs.ctype = ["ham", "cheese"]
+    frame = projective_wcs_to_frame(w)
+
+    assert frame is None
+
+
+def test_projective_frame_to_wcs(projective_frame):
+    wcs = projective_frame_to_wcs(projective_frame)
+
+    assert isinstance(wcs, WCS)
+    assert wcs.wcs.ctype[0] == "PJLN-TAN"
+    assert wcs.wcs.cunit[0] == "arcsec"
+    assert wcs.wcs.dateobs == "2024-01-01 00:00:00.000"
+
+    assert wcs.wcs.aux.rsun_ref == projective_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_projective_frame_no_oberver(projective_frame_no_obsever):
+    wcs = projective_frame_to_wcs(projective_frame_no_obsever)
+
+    assert isinstance(wcs, WCS)
+    assert wcs.wcs.ctype[0] == "PJLN-TAN"
+    assert wcs.wcs.cunit[0] == "arcsec"
+    assert wcs.wcs.dateobs == "2024-01-01 00:00:00.000"
+
+
+def test_projective_frame_to_wcs_none():
+    wcs = projective_frame_to_wcs(HeliographicStonyhurst())
+    assert wcs is None