diff --git a/changelog/76.feature.rst b/changelog/76.feature.rst new file mode 100644 index 0000000..6f1c88d --- /dev/null +++ b/changelog/76.feature.rst @@ -0,0 +1 @@ +Add `~xrayvision.coordinates.frames.Projective` coordinate frame. This is intended to represent the projective coordinate system of images with unknown pointing. diff --git a/xrayvision/coordinates/__init__.py b/xrayvision/coordinates/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/xrayvision/coordinates/frames.py b/xrayvision/coordinates/frames.py new file mode 100644 index 0000000..2db23d5 --- /dev/null +++ b/xrayvision/coordinates/frames.py @@ -0,0 +1,124 @@ +import astropy.coordinates +import astropy.units as u +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 + +__all__ = ["Projective"] + + +X_CTYPE = "PJLN" +Y_CTYPE = "PJLT" + + +class Projective(SunPyBaseCoordinateFrame): + """A generic projective coordinate frame for an image taken by an arbitrary imager.""" + + observer = ObserverCoordinateAttribute(HeliographicStonyhurst) + rsun = astropy.coordinates.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 + + 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( + lat=hgs_latitude * u.deg, lon=hgs_longitude * u.deg, radius=hgs_distance * u.m, obstime=dateavg, rsun=rsun + ) + + 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 + + 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) diff --git a/xrayvision/coordinates/tests/__init__.py b/xrayvision/coordinates/tests/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/xrayvision/coordinates/tests/test_frames.py b/xrayvision/coordinates/tests/test_frames.py new file mode 100644 index 0000000..81c89d1 --- /dev/null +++ b/xrayvision/coordinates/tests/test_frames.py @@ -0,0 +1,73 @@ +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_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 + + +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_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_to_wcs_none(): + wcs = projective_frame_to_wcs(HeliographicStonyhurst()) + assert wcs is None