read_files now handles more types and returns a dict

fits_info renamed and improved

.pre-commit-config.yaml

@@ -1,3 +1,4 @@
+exclude: ".*(.fits|.fts|.fit|.txt|.csv)$"
 repos:
   - repo: https://github.com/myint/docformatter
     rev: v1.7.5
@@ -14,7 +15,6 @@ repos:
             "--remove-all-unused-imports",
             "--remove-unused-variable",
           ]
-        exclude: ".*(.fits|.fts|.fit|.txt|tca.*|extern.*|.rst|.md|docs/conf.py)$"
   - repo: https://github.com/charliermarsh/ruff-pre-commit
     rev: "v0.2.1"
     hooks:
@@ -27,11 +27,8 @@ repos:
       - id: check-ast
       - id: check-case-conflict
       - id: trailing-whitespace
-        exclude: ".*(.fits|.fts|.fit|.txt|.csv)$"
       - id: mixed-line-ending
-        exclude: ".*(.fits|.fts|.fit|.txt|.csv)$"
       - id: end-of-file-fixer
-        exclude: ".*(.fits|.fts|.fit|.txt|.csv)$"
       - id: check-yaml
       - id: debug-statements
   - repo: https://github.com/codespell-project/codespell
@@ -40,7 +37,3 @@ repos:
       - id: codespell
         additional_dependencies:
           - tomli
-  - repo: https://github.com/pre-commit/mirrors-prettier
-    rev: v4.0.0-alpha.8
-    hooks:
-      - id: prettier

changelog/41.breaking.1.rst

@@ -0,0 +1,2 @@
+``fitsinfo`` now renamed to ``fits_info``.
+The output has been updated to be nicer.

changelog/41.breaking.rst

@@ -0,0 +1 @@
+``read_files`` now returns a dictionary with keys "raster" and "sji" where the values are a list of loaded IRIS classes.

docs/guide.rst

@@ -279,8 +279,8 @@ If you would like a bit more information, we have a similar function within `iri
 
 .. code-block:: python
 
-    >>> from irispy.io import fitsinfo  # doctest: +REMOTE_DATA
-    >>> fitsinfo(sample_data.SJI_1330)  # doctest: +REMOTE_DATA
+    >>> from irispy.io import fits_info  # doctest: +REMOTE_DATA
+    >>> fits_info(sample_data.SJI_1330)  # doctest: +REMOTE_DATA
     Filename: ...iris_l2_20211001_060925_3683602040_SJI_1330_t000.fits.gz
     No.    Name      Ver    Type      Cards   Dimensions   Format
       0  PRIMARY       1 PrimaryHDU     162   (555, 548, 20)   int16 (rescales to float32)
@@ -289,7 +289,7 @@ If you would like a bit more information, we have a similar function within `iri
 
 .. code-block:: python
 
-    >>> fitsinfo("iris_l2_20211001_060925_3683602040_raster_t000_r00000.fits") # doctest: +SKIP
+    >>> fits_info("iris_l2_20211001_060925_3683602040_raster_t000_r00000.fits") # doctest: +SKIP
     Filename: iris_l2_20211001_060925_3683602040_raster_t000_r00000.fits
     No.    Name      Ver    Type      Cards   Dimensions   Format
       0  PRIMARY       1 PrimaryHDU     215   ()

irispy/conftest.py

@@ -1,4 +1,6 @@
 import os
+import tarfile
+from pathlib import Path
 
 import numpy as np
 import pytest
@@ -10,84 +12,84 @@
 from irispy.utils import get_iris_response
 
 
-@pytest.fixture()
+@pytest.fixture(scope="session")
 def raster_file():
     return get_test_filepath("iris_l2_20210905_001833_3620258102_raster_t000_r00000_test.fits")
 
 
-@pytest.fixture()
+@pytest.fixture(scope="session")
 def sji_1330_file():
     return get_test_filepath("iris_l2_20210905_001833_3620258102_SJI_1330_t000_test.fits")
 
 
-@pytest.fixture()
+@pytest.fixture(scope="session")
 def sji_1400_file():
     return get_test_filepath("iris_l2_20210905_001833_3620258102_SJI_1400_t000_test.fits")
 
 
-@pytest.fixture()
+@pytest.fixture(scope="session")
 def sji_2796_file():
     return get_test_filepath("iris_l2_20210905_001833_3620258102_SJI_2796_t000_test.fits")
 
 
-@pytest.fixture()
+@pytest.fixture(scope="session")
 def sji_2832_file():
     return get_test_filepath("iris_l2_20210905_001833_3620258102_SJI_2832_t000_test.fits")
 
 
-@pytest.fixture()
+@pytest.fixture(scope="session")
 def SJICube_1330():
     from irispy.io.sji import read_sji_lvl2
 
     return read_sji_lvl2(get_test_filepath("iris_l2_20210905_001833_3620258102_SJI_1330_t000_test.fits"))
 
 
-@pytest.fixture()
+@pytest.fixture(scope="session")
 def SJICube_1400():
     return read_sji_lvl2(get_test_filepath("iris_l2_20210905_001833_3620258102_SJI_1400_t000_test.fits"))
 
 
-@pytest.fixture()
+@pytest.fixture(scope="session")
 def SJICube_2796():
     return read_sji_lvl2(get_test_filepath("iris_l2_20210905_001833_3620258102_SJI_2796_t000_test.fits"))
 
 
-@pytest.fixture()
+@pytest.fixture(scope="session")
 def SJICube_2832():
     return read_sji_lvl2(get_test_filepath("iris_l2_20210905_001833_3620258102_SJI_2832_t000_test.fits"))
 
 
-@pytest.fixture()
+@pytest.fixture(scope="session")
 def iris_response_v1():
     return get_iris_response(time_obs=parse_time("2013-09-03"), response_version=1)
 
 
-@pytest.fixture()
+@pytest.fixture(scope="session")
 def iris_response_v2():
     return get_iris_response(time_obs=parse_time("2013-09-03"), response_version=2)
 
 
-@pytest.fixture()
+@pytest.fixture(scope="session")
 def iris_response_v3():
     return get_iris_response(time_obs=parse_time("2013-09-03"), response_version=3)
 
 
-@pytest.fixture()
+@pytest.fixture(scope="session")
 def iris_response_v4():
     return get_iris_response(time_obs=parse_time("2013-09-03"), response_version=4)
 
 
-@pytest.fixture()
+@pytest.fixture(scope="session")
 def iris_response_v5():
     return get_iris_response(time_obs=parse_time("2013-09-03"), response_version=5)
 
 
-@pytest.fixture()
+@pytest.fixture(scope="session")
 def iris_response_v6():
     return get_iris_response(time_obs=parse_time("2013-09-03"), response_version=6)
 
 
-@pytest.fixture()
+@pytest.fixture(scope="session")
 def filelist():
     return [
         get_test_filepath("iris_l2_20210905_001833_3620258102_SJI_1330_t000_test.fits"),
@@ -110,3 +112,13 @@ def fake_long_obs(tmp_path_factory):
     fits_file = os.fspath(temp_dir.joinpath("iris_l2_20210905_001833_3620258102_SJI_2832_t000_test.fits"))
     hdu.writeto(fits_file)
     return [fits_file]
+
+
+@pytest.fixture(scope="session")
+def fake_tar(tmp_path_factory, filelist):
+    temp_dir = tmp_path_factory.mktemp("IRIS")
+    tar_file = os.fspath(temp_dir.joinpath("iris_l2_20210905_001833_3620258102_SJI_2832_t000_test.tar"))
+    with tarfile.open(tar_file, "w") as tar:
+        for file in filelist:
+            tar.add(file, arcname=Path(file).name)
+    return Path(tar_file)

irispy/io/sji.py

@@ -118,7 +118,7 @@ def _create_wcs(hdulist):
 
 def read_sji_lvl2(filename, *, uncertainty=False, memmap=False):
     """
-    Reads a level 2 SJI FITS.
+    Reads one level 2 SJI FITS.
 
     Parameters
     ----------

irispy/io/spectrograph.py

@@ -1,23 +1,125 @@
-import logging
 import tarfile
 from copy import copy
 from pathlib import Path
 
+import astropy.modeling.models as m
 import astropy.units as u
+import gwcs
+import gwcs.coordinate_frames as cf
 import numpy as np
 from astropy.coordinates import SkyCoord
 from astropy.io import fits
 from astropy.time import Time, TimeDelta
 from astropy.wcs import WCS
+from dkist.wcs.models import CoupledCompoundModel, VaryingCelestialTransform
+from sunpy import log
 from sunpy.coordinates import Helioprojective
+from sunpy.coordinates.ephemeris import get_body_heliographic_stonyhurst
 
 from irispy.spectrograph import Collection, SGMeta, SpectrogramCube, SpectrogramCubeSequence
-from irispy.utils import calculate_uncertainty
+from irispy.utils import calculate_uncertainty, tar_extract_all
 from irispy.utils.constants import DN_UNIT, READOUT_NOISE, SLIT_WIDTH
 
 __all__ = ["read_spectrograph_lvl2"]
 
 
+def _create_gwcs(hdulist: fits.HDUList) -> gwcs.WCS:
+    """
+    Creates the GWCS object for the SJI file.
+
+    Parameters
+    ----------
+    hdulist : `astropy.io.fits.HDUList`
+        The HDU list of the SJI file.
+
+    Returns
+    -------
+    `gwcs.WCS`
+        GWCS object for the SJI file.
+    """
+    pc_table = hdulist[1].data[:, hdulist[1].header["PC1_1IX"] : hdulist[1].header["PC2_2IX"] + 1].reshape(-1, 2, 2)
+    crval_table = hdulist[1].data[:, hdulist[1].header["XCENIX"] : hdulist[1].header["YCENIX"] + 1]
+    crpix = [hdulist[0].header["CRPIX1"], hdulist[0].header["CRPIX2"]]
+    cdelt = [hdulist[0].header["CDELT1"], hdulist[0].header["CDELT2"]]
+    celestial = VaryingCelestialTransform(
+        crpix=crpix * u.pixel,
+        cdelt=cdelt * u.arcsec / u.pixel,
+        pc_table=pc_table * u.pixel,
+        crval_table=crval_table * u.arcsec,
+    )
+    base_time = Time(hdulist[0].header["STARTOBS"], format="isot", scale="utc")
+    times = hdulist[1].data[:, hdulist[1].header["TIME"]] * u.s
+    # We need to account for a non-zero time delta.
+    base_time += times[0]
+    times -= times[0]
+    temporal = m.Tabular1D(
+        np.arange(hdulist[1].data.shape[0]) * u.pix,
+        lookup_table=times,
+        fill_value=np.nan,
+        bounds_error=False,
+        method="linear",
+    )
+    forward_transform = CoupledCompoundModel("&", left=celestial, right=temporal)
+    celestial_frame = cf.CelestialFrame(
+        axes_order=(0, 1),
+        unit=(u.arcsec, u.arcsec),
+        reference_frame=Helioprojective(observer="earth", obstime=base_time),
+        axis_physical_types=[
+            "custom:pos.helioprojective.lon",
+            "custom:pos.helioprojective.lat",
+        ],
+        axes_names=("Longitude", "Latitude"),
+    )
+    temporal_frame = cf.TemporalFrame(Time(base_time), unit=(u.s,), axes_order=(2,), axes_names=("Time (UTC)",))
+    output_frame = cf.CompositeFrame([celestial_frame, temporal_frame])
+    input_frame = cf.CoordinateFrame(
+        axes_order=(0, 1, 2),
+        naxes=3,
+        axes_type=["PIXEL", "PIXEL", "PIXEL"],
+        unit=(u.pix, u.pix, u.pix),
+    )
+    return gwcs.WCS(forward_transform, input_frame=input_frame, output_frame=output_frame)
+
+
+def _create_wcs(hdulist):
+    """
+    This is required as occasionally we need a normal WCS instead of a gWCS due
+    to compatibility issues.
+
+    This has been set to have an Earth Observer at the time of the
+    observation.
+    """
+    wcses = []
+    base_time = Time(hdulist[0].header["STARTOBS"], format="isot", scale="utc")
+    times = hdulist[1].data[:, hdulist[1].header["TIME"]] * u.s
+    # We need to account for a non-zero time delta.
+    base_time += times[0]
+    times -= times[0]
+    for i in range(hdulist[0].header["NAXIS3"]):
+        header = copy(hdulist[0].header)
+        header.pop("NAXIS3")
+        header.pop("PC3_1")
+        header.pop("PC3_2")
+        header.pop("CTYPE3")
+        header.pop("CUNIT3")
+        header.pop("CRVAL3")
+        header.pop("CRPIX3")
+        header.pop("CDELT3")
+        header["NAXIS"] = 2
+        header["CRVAL1"] = hdulist[1].data[i, hdulist[1].header["XCENIX"]]
+        header["CRVAL2"] = hdulist[1].data[i, hdulist[1].header["YCENIX"]]
+        header["PC1_1"] = hdulist[1].data[0, hdulist[1].header["PC1_1IX"]]
+        header["PC1_2"] = hdulist[1].data[0, hdulist[1].header["PC1_2IX"]]
+        header["PC2_1"] = hdulist[1].data[0, hdulist[1].header["PC2_1IX"]]
+        header["PC2_2"] = hdulist[1].data[0, hdulist[1].header["PC2_2IX"]]
+        header["DATE_OBS"] = (base_time + times[i]).isot
+        location = get_body_heliographic_stonyhurst("Earth", header["DATE_OBS"])
+        header["HGLN_OBS"] = location.lon.value
+        header["HGLT_OBS"] = location.lat.value
+        wcses.append(WCS(header))
+    return wcses
+
+
 def _pc_matrix(lam, angle_1, angle_2):
     return angle_1, -1 * lam * angle_2, 1 / lam * angle_2, angle_1
 
@@ -34,12 +136,16 @@ def read_spectrograph_lvl2(
     Reads IRIS level 2 spectrograph FITS from an OBS into an
     `.IRISSpectrograph` instance.
 
+    .. warning::
+
+        This function does not handle tar files.
+        Either extract them manually or use `irispy.io.read_files`.
+
     Parameters
     ----------
     filenames: `list` of `str` or `str`
         Filename of filenames to be read. They must all be associated with the same
         OBS number.
-        If you provide a tar file, it will be extracted at the same location.
     spectral_windows: iterable of `str` or `str`
         Spectral windows to extract from files. Default=None, implies, extract all
         spectral windows.
@@ -59,16 +165,6 @@ def read_spectrograph_lvl2(
     -------
     `ndcube.NDCollection`
     """
-    if isinstance(filenames, str):
-        if tarfile.is_tarfile(filenames):
-            parent = Path(filenames.replace(".tar.gz", "")).mkdir(parents=True, exist_ok=True)
-            with tarfile.open(filenames, "r") as tar:
-                tar.extractall(parent)
-                filenames = [parent / file for file in tar.getnames()]
-        else:
-            filenames = [filenames]
-
-    # Collecting the window observations
     with fits.open(filenames[0], memmap=memmap, do_not_scale_image_data=memmap) as hdulist:
         v34 = bool(hdulist[0].header["OBSID"].startswith("34"))
         hdulist.verify("silentfix")
@@ -151,7 +247,7 @@ def read_spectrograph_lvl2(
                         "The loading will continue but this will be missing in the final cube. "
                         f"Spectral window: {window_name}, step {i} in file: {filename}"
                     )
-                    logging.warning(msg)
+                    log.warning(msg)
                     continue
                 out_uncertainty = None
                 data_mask = None