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

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/spectrograph.py

@@ -1,15 +1,20 @@
-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
@@ -18,6 +23,103 @@
 __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
 
@@ -151,7 +253,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

irispy/io/tests/test_utils.py

@@ -1,51 +1,56 @@
-import pytest
+from irispy.io.utils import fits_info, read_files, tar_extract_all
 
-from irispy.io.utils import fitsinfo, read_files
 
+def test_tar_extract_all(fake_tar):
+    # fake_tar is a Path
+    extracted_files = tar_extract_all(fake_tar)
+    assert len(extracted_files) == 4
+    for file in extracted_files:
+        assert file.exists()
+        assert file.is_file()
+    assert all(file.suffix == ".fits" for file in extracted_files)
+    assert all(file.parent == fake_tar.parent / fake_tar.stem for file in extracted_files)
 
-def test_fitsinfo(capsys, raster_file, sji_1330_file, sji_1400_file, sji_2796_file, sji_2832_file):
-    fitsinfo(raster_file)
+
+def test_fits_info(capsys, raster_file, sji_1330_file, sji_1400_file, sji_2796_file, sji_2832_file):
+    RASTER = "\nObservation: Medium sit-and-stare 0.3x60 1s  C II   Si IV   Mg II h/k   Mg II w s\nOBS ID: 3620258102\nNo.   Name  Ver ...          Format                 Description       \n--- ------- --- ... ------------------------ -------------------------\n  0 PRIMARY   1 ...                  float32 SJI 1330 (1310 - 1350 AA)\n  1           1 ...                  float64            Auxiliary data\n  2           1 ... [A10, A10, A4, A66, A63]            Auxiliary data\n' = CaptureResult(out='Filename: /home/nabil/Git/irispy-lmsal/irispy/data/test/iris_l2_20210905_001833_3620258102_SJI_1330_t000_test.fits\nObservation: Medium sit-and-stare 0.3x60 1s  C II   Si IV   Mg II h/k   Mg II w s\nOBS ID: 3620258102\nNo.   Name  Ver ...          Format                 Description       \n--- ------- --- ... ------------------------ -------------------------\n  0 PRIMARY   1 ...                  float32 SJI 1330 (1310 - 1350 AA)\n  1           1 ...                  float64            Auxiliary data\n  2           1 ... [A10, A10, A4, A66, A63]            Auxiliary data\n'"
+    fits_info(raster_file)
     captured = capsys.readouterr()
     assert raster_file in captured.out
+    assert RASTER in captured.out
 
-    fitsinfo(sji_1330_file)
+    SJI = ""
+    fits_info(sji_1330_file)
     captured = capsys.readouterr()
-    assert sji_1330_file in captured.out
+    assert sji_1400_file in captured.out
+    assert SJI in captured.out
 
-    fitsinfo(sji_1400_file)
+    fits_info(sji_1400_file)
     captured = capsys.readouterr()
     assert sji_1400_file in captured.out
 
-    fitsinfo(sji_2796_file)
+    fits_info(sji_2796_file)
     captured = capsys.readouterr()
     assert sji_2796_file in captured.out
 
-    fitsinfo(sji_2832_file)
+    fits_info(sji_2832_file)
     captured = capsys.readouterr()
     assert sji_2832_file in captured.out
 
 
-def test_read_files_errors_with_mix(raster_file, sji_1330_file):
-    with pytest.raises(ValueError, match="You cannot mix raster and SJI files."):
-        read_files([raster_file, sji_1330_file])
-
-
 def test_read_files_raster(raster_file):
-    # Simple test to ensure it does not error
     read_files(raster_file)
     read_files([raster_file])
+    read_files([raster_file, raster_file])
 
 
 def test_read_files_sji(sji_1330_file, sji_1400_file, sji_2796_file, sji_2832_file):
-    # Simple test to ensure it does not error
     read_files(sji_1330_file)
     read_files(sji_1400_file)
     read_files(sji_2796_file)
     read_files(sji_2832_file)
-    read_files([sji_2832_file])
+    read_files([sji_2832_file, sji_2796_file, sji_1400_file, sji_1330_file])
 
 
-def test_read_files_sji_more_than_one(sji_1330_file):
-    # Don't allow more than one SJI file for now.
-    with pytest.raises(ValueError, match="You cannot load more than one SJI file at a time."):
-        read_files([sji_1330_file, sji_1330_file])
+def test_read_files_with_mix(raster_file, sji_1330_file):
+    read_files([raster_file, sji_1330_file])