diff --git a/changelog/90.feature.rst b/changelog/90.feature.rst new file mode 100644 index 0000000..341235f --- /dev/null +++ b/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. diff --git a/docs/reference/coordinates.rst b/docs/reference/coordinates.rst new file mode 100644 index 0000000..2be0f49 --- /dev/null +++ b/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 diff --git a/docs/reference/frames.rst b/docs/reference/frames.rst deleted file mode 100644 index 94ecd07..0000000 --- a/docs/reference/frames.rst +++ /dev/null @@ -1,8 +0,0 @@ -.. _frames: - -Frames ('stixpy.frames') -************************ - -The ``frames`` submodule contains STIX specific coordinate frames and transforms. - -.. automodapi:: stixpy.frames diff --git a/docs/reference/index.rst b/docs/reference/index.rst index 3007e01..f12351f 100644 --- a/docs/reference/index.rst +++ b/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 diff --git a/docs/reference/map.rst b/docs/reference/map.rst new file mode 100644 index 0000000..26a07f3 --- /dev/null +++ b/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 diff --git a/setup.cfg b/setup.cfg index 186363b..1fc17e5 100644 --- a/setup.cfg +++ b/setup.cfg @@ -68,7 +68,7 @@ filterwarnings = [isort] balanced_wrapping = true -skip = +extend_skip = docs/conf.py, stixpy/__init__.py, default_section = THIRDPARTY diff --git a/stixpy/coordinates/__init__.py b/stixpy/coordinates/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/stixpy/coordinates/frames.py b/stixpy/coordinates/frames.py new file mode 100644 index 0000000..5660f07 --- /dev/null +++ b/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) diff --git a/stixpy/coordinates/tests/__init__.py b/stixpy/coordinates/tests/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/stixpy/coordinates/tests/test_frames.py b/stixpy/coordinates/tests/test_frames.py new file mode 100644 index 0000000..1402bdb --- /dev/null +++ b/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) diff --git a/stixpy/tests/test_frames.py b/stixpy/coordinates/tests/test_transforms.py similarity index 52% rename from stixpy/tests/test_frames.py rename to stixpy/coordinates/tests/test_transforms.py index 4155c71..fd168c3 100644 --- a/stixpy/tests/test_frames.py +++ b/stixpy/coordinates/tests/test_transforms.py @@ -3,44 +3,14 @@ import astropy.units as u import numpy as np import pytest -from astropy.tests.helper import assert_quantity_allclose from astropy.time import Time -from astropy.wcs import WCS -from sunpy.coordinates.frames import HeliographicStonyhurst, Helioprojective +from sunpy.coordinates import HeliographicStonyhurst, Helioprojective -from stixpy.frames import STIXImaging, stix_frame_to_wcs, stix_wcs_to_frame +from stixpy.coordinates.frames import STIXImaging -@pytest.fixture -def stix_wcs(): - w = WCS(naxis=2) - - w.wcs.dateavg = '2024-01-01' - 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 assert_quantity_allclose(x, param): + pass @pytest.mark.skip(reason="Test data maybe incorrect") @@ -61,7 +31,7 @@ def test_transforms_with_know_values(): assert_quantity_allclose(earth_coord.Ty, stix_to_earth.Ty) -@mock.patch("stixpy.frames._get_aux_data") +@mock.patch("stixpy.coordinates.transforms._get_aux_data") def test_hpc_to_stx(mock): # fake some data to make check easier obstime = Time("2021-05-22T02:52:00") @@ -72,11 +42,11 @@ def test_hpc_to_stx(mock): stix_frame = STIXImaging(0 * u.arcsec, 0 * u.arcsec, obstime=obstime) stix_coord = solo_hpc_coord.transform_to(stix_frame) # should match the offset -8, 60 - assert_quantity_allclose(stix_coord.x, -8 * u.arcsec) - assert_quantity_allclose(stix_coord.y, 60 * u.arcsec) + assert_quantity_allclose(stix_coord.Tx, -8 * u.arcsec) + assert_quantity_allclose(stix_coord.Ty, 60 * u.arcsec) -@mock.patch("stixpy.frames._get_aux_data") +@mock.patch("stixpy.coordinates.transforms._get_aux_data") def test_hpc_to_stx_no_sas(mock): # fake some data to make check easier obstime = Time("2021-05-22T02:52:00") @@ -88,37 +58,5 @@ def test_hpc_to_stx_no_sas(mock): with pytest.warns(Warning, match="SAS solution not"): stix_coord = solo_hpc_coord.transform_to(stix_frame) # should match the offset -8, 60 added to yaw and pitch 10, 10 - assert_quantity_allclose(stix_coord.x, (10 - 8) * u.arcsec) - assert_quantity_allclose(stix_coord.y, (10 + 60) * u.arcsec) - - -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 + assert_quantity_allclose(stix_coord.Tx, (10 - 8) * u.arcsec) + assert_quantity_allclose(stix_coord.Ty, (10 + 60) * u.arcsec) diff --git a/stixpy/coordinates/transforms.py b/stixpy/coordinates/transforms.py new file mode 100644 index 0000000..b075263 --- /dev/null +++ b/stixpy/coordinates/transforms.py @@ -0,0 +1,157 @@ +import warnings + +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.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 + +__all__ = ['get_hpc_info', 'stixim_to_hpc', 'hpc_to_stixim'] + +def _get_rotation_matrix_and_position(obstime): + r""" + Return rotation matrix STIX Imaging to SOLO HPC and position of SOLO in HEEQ. + + Parameters + ---------- + obstime : `astropy.time.Time` + Time + Returns + ------- + """ + roll, solo_position_heeq, spacecraft_pointing = get_hpc_info(obstime) + + # Generate the rotation matrix using the x-convention (see Goldstein) + # Rotate +90 clockwise around the first axis + # STIX is mounted on the -Y panel (SOLO_SRF) need to rotate to from STIX frame to SOLO_SRF + rot_to_solo = np.array([[1, 0, 0], [0, 0, -1], [0, 1, 0]]) + + logger.debug("Pointing: %s.", spacecraft_pointing) + A = rotation_matrix(-1 * roll, "x") + B = rotation_matrix(spacecraft_pointing[1], "y") + C = rotation_matrix(-1 * spacecraft_pointing[0], "z") + + # Will be applied right to left + rmatrix = matrix_product(A, B, C, rot_to_solo) + + return rmatrix, solo_position_heeq + + +def get_hpc_info(obstime): + roll, pitch, yaw, sas_x, sas_y, solo_position_heeq = _get_aux_data(obstime) + if np.isnan(sas_x) or np.isnan(sas_y): + warnings.warn("SAS solution not available using spacecraft pointing and average SAS offset") + sas_x = yaw + sas_y = pitch + # 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) + spacecraft_pointing = np.hstack([STIX_X_SHIFT + rotated_yaw, STIX_Y_SHIFT + rotated_pitch]) + sas_pointing = np.hstack([sas_x + STIX_X_OFFSET, -1 * sas_y + STIX_Y_OFFSET]) + pointing_diff = np.linalg.norm(spacecraft_pointing - sas_pointing) + if pointing_diff < 200 * u.arcsec: + logger.debug(f"Using SAS pointing {sas_pointing}") + spacecraft_pointing = sas_pointing + else: + logger.debug(f"Using spacecraft pointing {spacecraft_pointing}") + warnings.warn( + f"Using spacecraft pointing large difference between " + f"SAS {sas_pointing} and spacecraft {spacecraft_pointing}." + ) + return roll, solo_position_heeq, spacecraft_pointing + + +def _get_aux_data(obstime): + # Find, download, read aux file with pointing, sas and position information + logger.debug("Searching for AUX data") + query = Fido.search( + a.Time(obstime, obstime), a.Instrument.stix, a.Level.l2, a.stix.DataType.aux, a.stix.DataProduct.aux_ephemeris + ) + if len(query["stix"]) != 1: + raise ValueError("Exactly one AUX file should be found but %d were found.", len(query["stix"])) + logger.debug("Downloading AUX data") + files = Fido.fetch(query) + if len(files.errors) > 0: + raise ValueError("There were errors downloading the data.") + # Read and extract data + logger.debug("Loading and extracting AUX data") + hdu = fits.getheader(files[0], ext=0) + aux = QTable.read(files[0], hdu=2) + start_time = Time(hdu.get("DATE-BEG")) + rel_date = (obstime - start_time).to("s") + idx = np.argmin(np.abs(rel_date - aux["time"])) + sas_x, sas_y = aux[idx][["y_srf", "z_srf"]] + roll, pitch, yaw = aux[idx]["roll_angle_rpy"] + solo_position_heeq = aux[idx]["solo_loc_heeq_zxy"] + logger.debug( + "SAS: %s, %s, Orientation: %s, %s, %s, Position: %s", + sas_x, + sas_y, + roll, + yaw.to("arcsec"), + pitch.to("arcsec"), + solo_position_heeq, + ) + + # roll, pitch, yaw = np.mean(aux[1219:1222]['roll_angle_rpy'], axis=0) + # sas_x = np.mean(aux[1219:1222]['y_srf']) + # sas_y = np.mean(aux[1219:1222]['z_srf']) + + return roll, pitch, yaw, sas_x, sas_y, solo_position_heeq + + +@frame_transform_graph.transform(coord.FunctionTransform, STIXImaging, Helioprojective) +def stixim_to_hpc(stxcoord, hpcframe): + r""" + Transform STIX Imaging coordinates to given HPC frame + """ + logger.debug("STIX: %s", stxcoord) + + rot_matrix, solo_pos_heeq = _get_rotation_matrix_and_position(stxcoord.obstime) + solo_heeq = HeliographicStonyhurst(solo_pos_heeq, representation_type="cartesian", obstime=stxcoord.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) + ) + logger.debug("SOLO HPC: %s", solo_hpc) + + # Transform from SOLO HPC to input HPC + if hpcframe.observer is None: + hpcframe = type(hpcframe)(obstime=hpcframe.obstime, observer=solo_heeq) + hpc = solo_hpc.transform_to(hpcframe) + logger.debug("Target HPC: %s", hpc) + return hpc + + +@frame_transform_graph.transform(coord.FunctionTransform, Helioprojective, STIXImaging) +def hpc_to_stixim(hpccoord, stxframe): + r""" + 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, representation_type="cartesian", obstime=stxframe.obstime) + + # Create SOLO HPC + solo_hpc_frame = Helioprojective(obstime=stxframe.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) + logger.debug("STIX: %s", stx_coord) + return stx_coord diff --git a/stixpy/frames.py b/stixpy/frames.py deleted file mode 100644 index 974880e..0000000 --- a/stixpy/frames.py +++ /dev/null @@ -1,365 +0,0 @@ -import warnings - -import astropy -import astropy.coordinates as coord -import astropy.units as u -import numpy as np -from astropy.coordinates import QuantityAttribute, frame_transform_graph -from astropy.coordinates.matrix_utilities import matrix_product, matrix_transpose, rotation_matrix -from astropy.io import fits -from astropy.table import QTable -from astropy.time import Time -from astropy.wcs import WCS -from sunpy.coordinates.frameattributes import ObserverCoordinateAttribute -from sunpy.coordinates.frames import HeliographicStonyhurst, Helioprojective, SunPyBaseCoordinateFrame -from sunpy.net import Fido -from sunpy.net import attrs as a -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", "_get_rotation_matrix_and_position", "hpc_to_stixim", "stixim_to_hpc"] - -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", "x", u.arcsec), - coord.RepresentationMapping("lat", "y", u.arcsec), - coord.RepresentationMapping("distance", "distance"), - ] - } - - def __init__(self, *args, **kwargs): - super().__init__(*args, **kwargs) - - # def make_3d(self): - # """ - # This method calculates the third coordinate of the Helioprojective - # frame. It assumes that the coordinate point is on the surface of the Sun. - # If a point in the frame is off limb then NaN will be returned. - # Returns - # ------- - # new_frame : `~sunpy.coordinates.frames.Helioprojective` - # A new frame instance with all the attributes of the original but - # now with a third coordinate. - # """ - # # Skip if we already are 3D - # if not self._is_2d: - # return self - # - # if not isinstance(self.observer, BaseCoordinateFrame): - # raise ConvertError("Cannot calculate distance to the Sun " - # f"for observer '{self.observer}' " - # "without `obstime` being specified.") - # - # rep = self.represent_as(UnitSphericalRepresentation) - # lat, lon = rep.lat, rep.lon - # - # # Check for the use of floats with lower precision than the native Python float - # if not set([lon.dtype.type, lat.dtype.type]).issubset([float, np.float64, np.longdouble]): - # warn_user("The Helioprojective component values appear to be lower " - # "precision than the native Python float: " - # f"Tx is {lon.dtype.name}, and Ty is {lat.dtype.name}. " - # "To minimize precision loss, you may want to cast the values to " - # "`float` or `numpy.float64` via the NumPy method `.astype()`.") - # - # # Calculate the distance to the surface of the Sun using the law of cosines - # cos_alpha = np.cos(lat) * np.cos(lon) - # c = self.observer.radius ** 2 - self.rsun ** 2 - # b = -2 * self.observer.radius * cos_alpha - # # Ignore sqrt of NaNs - # with np.errstate(invalid='ignore'): - # d = ((-1 * b) - np.sqrt(b ** 2 - 4 * c)) / 2 # use the "near" solution - # - # if self._spherical_screen: - # sphere_center = self._spherical_screen['center'].transform_to(self).cartesian - # c = sphere_center.norm() ** 2 - self._spherical_screen['radius'] ** 2 - # b = -2 * sphere_center.dot(rep) - # # Ignore sqrt of NaNs - # with np.errstate(invalid='ignore'): - # dd = ((-1 * b) + np.sqrt(b ** 2 - 4 * c)) / 2 # use the "far" solution - # - # d = np.fmin(d, dd) if self._spherical_screen['only_off_disk'] else dd - # - # # This warning can be triggered in specific draw calls when plt.show() is called - # # we can not easily prevent this, so we check the specific function is being called - # # within the stack trace. - # stack_trace = traceback.format_stack() - # matching_string = 'wcsaxes.*_draw_grid' - # bypass = any([re.search(matching_string, string) for string in stack_trace]) - # if not bypass and np.all(np.isnan(d)) and np.any(np.isfinite(cos_alpha)): - # warn_user("The conversion of these 2D helioprojective coordinates to 3D is all NaNs " - # "because off-disk coordinates need an additional assumption to be mapped to " - # "calculate distance from the observer. Consider using the context manager " - # "`Helioprojective.assume_spherical_screen()`.") - # - # return self.realize_frame(SphericalRepresentation(lon=lon, - # lat=lat, - # distance=d)) - - -def _get_rotation_matrix_and_position(obstime): - r""" - Return rotation matrix STIX Imaging to SOLO HPC and position of SOLO in HEEQ. - - Parameters - ---------- - obstime : `astropy.time.Time` - Time - Returns - ------- - """ - roll, solo_position_heeq, spacecraft_pointing = get_hpc_info(obstime) - - # Generate the rotation matrix using the x-convention (see Goldstein) - # Rotate +90 clockwise around the first axis - # STIX is mounted on the -Y panel (SOLO_SRF) need to rotate to from STIX frame to SOLO_SRF - rot_to_solo = np.array([[1, 0, 0], [0, 0, -1], [0, 1, 0]]) - - logger.debug("Pointing: %s.", spacecraft_pointing) - A = rotation_matrix(-1 * roll, "x") - B = rotation_matrix(spacecraft_pointing[1], "y") - C = rotation_matrix(-1 * spacecraft_pointing[0], "z") - - # Will be applied right to left - rmatrix = matrix_product(A, B, C, rot_to_solo) - - return rmatrix, solo_position_heeq - - -def get_hpc_info(obstime): - roll, pitch, yaw, sas_x, sas_y, solo_position_heeq = _get_aux_data(obstime) - if np.isnan(sas_x) or np.isnan(sas_y): - warnings.warn("SAS solution not available using spacecraft pointing and average SAS offset") - sas_x = yaw - sas_y = pitch - # 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) - spacecraft_pointing = np.hstack([STIX_X_SHIFT + rotated_yaw, STIX_Y_SHIFT + rotated_pitch]) - sas_pointing = np.hstack([sas_x + STIX_X_OFFSET, -1 * sas_y + STIX_Y_OFFSET]) - pointing_diff = np.linalg.norm(spacecraft_pointing - sas_pointing) - if pointing_diff < 200 * u.arcsec: - logger.debug(f"Using SAS pointing {sas_pointing}") - spacecraft_pointing = sas_pointing - else: - logger.debug(f"Using spacecraft pointing {spacecraft_pointing}") - warnings.warn( - f"Using spacecraft pointing large difference between " - f"SAS {sas_pointing} and spacecraft {spacecraft_pointing}." - ) - return roll, solo_position_heeq, spacecraft_pointing - - -def _get_aux_data(obstime): - # Find, download, read aux file with pointing, sas and position information - logger.debug("Searching for AUX data") - query = Fido.search( - a.Time(obstime, obstime), a.Instrument.stix, a.Level.l2, a.stix.DataType.aux, a.stix.DataProduct.aux_ephemeris - ) - if len(query["stix"]) != 1: - raise ValueError("Exactly one AUX file should be found but %d were found.", len(query["stix"])) - logger.debug("Downloading AUX data") - files = Fido.fetch(query) - if len(files.errors) > 0: - raise ValueError("There were errors downloading the data.") - # Read and extract data - logger.debug("Loading and extracting AUX data") - hdu = fits.getheader(files[0], ext=0) - aux = QTable.read(files[0], hdu=2) - start_time = Time(hdu.get("DATE-BEG")) - rel_date = (obstime - start_time).to("s") - idx = np.argmin(np.abs(rel_date - aux["time"])) - sas_x, sas_y = aux[idx][["y_srf", "z_srf"]] - roll, pitch, yaw = aux[idx]["roll_angle_rpy"] - solo_position_heeq = aux[idx]["solo_loc_heeq_zxy"] - logger.debug( - "SAS: %s, %s, Orientation: %s, %s, %s, Position: %s", - sas_x, - sas_y, - roll, - yaw.to("arcsec"), - pitch.to("arcsec"), - solo_position_heeq, - ) - - # roll, pitch, yaw = np.mean(aux[1219:1222]['roll_angle_rpy'], axis=0) - # sas_x = np.mean(aux[1219:1222]['y_srf']) - # sas_y = np.mean(aux[1219:1222]['z_srf']) - - return roll, pitch, yaw, sas_x, sas_y, solo_position_heeq - - -@frame_transform_graph.transform(coord.FunctionTransform, STIXImaging, Helioprojective) -def stixim_to_hpc(stxcoord, hpcframe): - r""" - Transform STIX Imaging coordinates to given HPC frame - """ - logger.debug("STIX: %s", stxcoord) - - rot_matrix, solo_pos_heeq = _get_rotation_matrix_and_position(stxcoord.obstime) - solo_heeq = HeliographicStonyhurst(solo_pos_heeq, representation_type="cartesian", obstime=stxcoord.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) - ) - logger.debug("SOLO HPC: %s", solo_hpc) - - # Transform from SOLO HPC to input HPC - if hpcframe.observer is None: - hpcframe = type(hpcframe)(obstime=hpcframe.obstime, observer=solo_heeq) - hpc = solo_hpc.transform_to(hpcframe) - logger.debug("Target HPC: %s", hpc) - return hpc - - -@frame_transform_graph.transform(coord.FunctionTransform, Helioprojective, STIXImaging) -def hpc_to_stixim(hpccoord, stxframe): - r""" - 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, representation_type="cartesian", obstime=stxframe.obstime) - - # Create SOLO HPC - solo_hpc_frame = Helioprojective(obstime=stxframe.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) - logger.debug("STIX: %s", stx_coord) - return stx_coord - - -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 - """ - - # Not a STIX wcs bail out early - if set(wcs.wcs.ctype) != {'SXLN-TAN', 'SXLT-TAN'}: - return None - - dateobs = wcs.wcs.dateavg - - 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): - 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 - - -astropy.wcs.utils.WCS_FRAME_MAPPINGS.append([stix_wcs_to_frame]) -astropy.wcs.utils.FRAME_WCS_MAPPINGS.append([stix_frame_to_wcs]) diff --git a/stixpy/map/__init__.py b/stixpy/map/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/stixpy/map/stix.py b/stixpy/map/stix.py new file mode 100644 index 0000000..9816575 --- /dev/null +++ b/stixpy/map/stix.py @@ -0,0 +1,79 @@ +import astropy.units as u +from sunpy.map import GenericMap + +__all__ = ['STIXMap'] + + +class STIXMap(GenericMap): + r""" + A class to represent a STIX Map. + """ + + @classmethod + def is_datasource_for(cls, data, header, **kwargs): + return str(header.get('TELESCOP', '')).endswith('STIX') + + # @property + # def coordinate_system(self): + # """ + # Coordinate system used + # + # Overrides the values in the header which are not understood by Astropy WCS + # """ + # return SpatialPair("SXLN-TAN", "SXLT-TAN") + # + # @property + # def spatial_units(self): + # return SpatialPair(u.arcsec, u.arcsec) + + def plot(self, **kwargs): + + # Until this is resolved need to manually override https://github.com/astropy/astropy/issues/16339 + res = super().plot(**kwargs) + + res.axes.coords[0].set_coord_type('longitude', coord_wrap=180*u.deg) + res.axes.coords[1].set_coord_type('latitude') + + res.axes.coords[0].set_format_unit(u.arcsec) + res.axes.coords[1].set_format_unit(u.arcsec) + + return res + + # @property + # def wcs(self): + # if not self.coordinate_system.axis1 == 'SXLN-TAN' and self.coordinate_system.axis2 == 'SXLT-TAN': + # return super().wcs + # + # w2 = astropy.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.dateobs = self.date.isot + # w2.wcs.aux.rsun_ref = self.rsun_meters.to_value(u.m) + # + # obs_frame = self.observer_coordinate + # if hasattr(obs_frame, 'frame'): + # obs_frame = obs_frame.frame + # + # w2.wcs.aux.hgln_obs = obs_frame.lon.to_value(u.deg) + # w2.wcs.aux.hglt_obs = obs_frame.lat.to_value(u.deg) + # w2.wcs.aux.dsun_obs = obs_frame.radius.to_value(u.m) + # + # w2.wcs.dateobs = self.date.utc.iso + # + # # Set the shape of the data array + # w2.array_shape = self.data.shape + # + # return w2