diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index ed5837a..fd06a89 100644 --- a/.pre-commit-config.yaml +++ b/.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 diff --git a/changelog/41.breaking.1.rst b/changelog/41.breaking.1.rst new file mode 100644 index 0000000..fc0a649 --- /dev/null +++ b/changelog/41.breaking.1.rst @@ -0,0 +1,2 @@ +``fitsinfo`` now renamed to ``fits_info``. +The output has been updated to be nicer. diff --git a/changelog/41.breaking.rst b/changelog/41.breaking.rst new file mode 100644 index 0000000..16bee80 --- /dev/null +++ b/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. diff --git a/docs/guide.rst b/docs/guide.rst index 3e05e0a..b83ea4a 100644 --- a/docs/guide.rst +++ b/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 () diff --git a/irispy/conftest.py b/irispy/conftest.py index 32e7952..0e09759 100644 --- a/irispy/conftest.py +++ b/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) diff --git a/irispy/io/sji.py b/irispy/io/sji.py index 8bd7481..d41f1c8 100644 --- a/irispy/io/sji.py +++ b/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 ---------- diff --git a/irispy/io/spectrograph.py b/irispy/io/spectrograph.py index ce1c378..babfcfd 100644 --- a/irispy/io/spectrograph.py +++ b/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 diff --git a/irispy/io/tests/test_utils.py b/irispy/io/tests/test_utils.py index 4afbb11..522ab08 100644 --- a/irispy/io/tests/test_utils.py +++ b/irispy/io/tests/test_utils.py @@ -1,51 +1,91 @@ -import pytest +from irispy.io.utils import fits_info, read_files, tar_extract_all +from irispy.sji import SJICube +from irispy.spectrograph import Collection -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): + fits_info(raster_file) captured = capsys.readouterr() assert raster_file in captured.out - fitsinfo(sji_1330_file) + fits_info(sji_1330_file) captured = capsys.readouterr() assert sji_1330_file 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): + output = read_files(raster_file) + assert output["raster"] is not None + assert len(output["raster"]) == 1 + assert not output["sji"] + output = read_files([raster_file]) + assert output["raster"] is not None + assert len(output["raster"]) == 1 + assert not output["sji"] -def test_read_files_raster(raster_file): - # Simple test to ensure it does not error - read_files(raster_file) - read_files([raster_file]) + output = read_files([raster_file, raster_file]) + assert output["raster"] is not None + assert len(output["raster"]) == 2 + assert not output["sji"] 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]) - - -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]) + output = read_files(sji_1330_file) + assert output["sji"] is not None + assert len(output["sji"]) == 1 + assert not output["raster"] + + output = read_files(sji_1400_file) + assert output["sji"] is not None + assert len(output["sji"]) == 1 + assert not output["raster"] + + output = read_files(sji_2796_file) + assert output["sji"] is not None + assert len(output["sji"]) == 1 + assert not output["raster"] + + output = read_files(sji_2832_file) + assert output["sji"] is not None + assert len(output["sji"]) == 1 + assert not output["raster"] + + output = read_files([sji_2832_file, sji_2796_file, sji_1400_file, sji_1330_file]) + assert output["sji"] is not None + assert len(output["sji"]) == 4 + assert not output["raster"] + + +def test_read_files_with_mix(raster_file, sji_1330_file): + output = read_files([raster_file, sji_1330_file]) + assert output["sji"] is not None + assert len(output["sji"]) == 1 + assert isinstance(output["sji"][0], SJICube) + + assert output["raster"] is not None + assert len(output["raster"]) == 1 + assert isinstance(output["raster"][0], Collection) diff --git a/irispy/io/utils.py b/irispy/io/utils.py index dc1ab49..e896a27 100644 --- a/irispy/io/utils.py +++ b/irispy/io/utils.py @@ -1,13 +1,46 @@ -import logging +import sys import tarfile +import warnings +from collections.abc import Iterable from pathlib import Path from astropy.io import fits +from astropy.table import Table +from sunpy import log -__all__ = ["fitsinfo", "read_files"] +from irispy.io.sji import read_sji_lvl2 +from irispy.io.spectrograph import read_spectrograph_lvl2 +__all__ = ["fits_info", "read_files", "tar_extract_all"] -def fitsinfo(filename): + +def tar_extract_all(filename: Path | str) -> list[Path]: + """ + Extract all files a given tarfile and return the files. + + It will extract the files in the same folder as the tarfile. + + Parameters + ---------- + filename : pathlib.Path, str + The tarfile to extract. + + Returns + ------- + list of pathlib.Path + The list of files extracted from the tarfile. + """ + filename = Path(filename) + final_location = filename.parent / filename.stem + final_location.mkdir(parents=True, exist_ok=True) + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + with tarfile.open(filename, "r") as tar: + tar.extractall(final_location) + return [final_location / file for file in tar.getnames()] + + +def fits_info(filename: str) -> None: """ Prints information about the extension of a raster or SJI level 2 data file. @@ -17,41 +50,75 @@ def fitsinfo(filename): filename : str Filename to load. """ + + def get_description(idx, idx_mod): + if idx == 0 and idx_mod == 0: + return "Primary Header (no data)" + if idx not in [num_extensions - 2, num_extensions - 1]: + text_modifier = f" ({header[f'TDET{idx+idx_mod}'][:3]})" if "SPEC" in header[f"TDET{idx+idx_mod}"] else "" + return f"{header[f'TDESC{idx+idx_mod}'].replace('_', ' ')} ({header[f'TWMIN{idx+idx_mod}']:.0f} - {header[f'TWMAX{idx+idx_mod}']:.0f} AA{text_modifier})" + return "Auxiliary data" + + results = [ + f"Filename: {Path(filename).absolute()}", + f"Observation: {fits.getval(filename, 'OBS_DESC')}", + f"OBS ID: {fits.getval(filename, 'OBSID')}", + ] + table = Table( + names=("No.", "Name", "Ver", "Type", "Cards", "Dimensions", "Format", "Description"), + dtype=("S4", "S10", "S4", "S10", "S4", "S20", "S10", "S100"), + ) with fits.open(filename) as hdulist: - hdulist.info() - hdr = hdulist[0].header - msg = f"Observation description: {hdr['OBS_DESC']}" - logging.info(msg) - modifier = "" - for i in range(hdr["NWIN"]): - msg = f"Extension No. {i+1} stores data and header of {hdr[f'TDESC{i+1}']}: " - logging.info(msg) - if "SJI" not in hdr[f"TDET{i + 1}"]: - modifier = f" ({hdr[f'TDET{i + 1}'][:3]})" - msg = f"{hdr[f'TWMIN{i + 1}']:.2f} - {hdr[f'TWMAX{i + 1}']:.2f} AA{modifier}" - logging.info(msg) - - -def read_files(filename, *, spectral_windows=None, uncertainty=False, memmap=True): + hdu_info = hdulist.info(output=False) + header = hdulist[0].header + num_extensions = len(hdulist) + for idx in range(num_extensions): + description = get_description(idx, 0 if "SPEC" in header["INSTRUME"] else 1) + hdu_info[idx] = hdu_info[idx][:-1] + (description,) + table.add_row(list(map(str, hdu_info[idx]))) + + sys.stdout.write("\n".join(results) + "\n") + table.pprint() + sys.stdout.flush() + + +def read_files( + filepath: list[Path | str] | Path | str, + *, + spectral_windows: Iterable[str] | str | None = None, + uncertainty: bool = False, + memmap: bool = True, +) -> dict[str, list]: """ A wrapper function to read a raster or SJI level 2 data file. - You can provide one SJI image or a one raster image or a list of raster images. + You can provide: + + - one SJI file + - one raster file + - list of raster file + - list of SJI files + - list of mixed raster and SJI files + + If any of them are tar files, they will be extracted. + Then they will append to the list of files to read at the end. + + This does not glob, so you cannot use wildcards, so you need to provide the + full path to the files. - If you mix raster and SJI images, the function will raise an error. + The order of the files will be the same as the input list, minus the tar files. Parameters ---------- - filename : `list of `str`, `str`, `pathlib.Path` + filepath : `list of `str` or `pathlib.Path`, `str`, `pathlib.Path` Filename(s) to load. - If given a string, will load that file. - If given a list of strings, it will check they are all raster files and load them. spectral_windows: iterable of `str` or `str` - Spectral windows to extract from files. Default=None, implies, extract all - spectral windows. + Spectral windows to extract from files. + Default is None which means it will extract all spectral windows. + Only used for raster files. uncertainty : `bool`, optional - If `True` (not the default), will compute the uncertainty for the data (slower and - uses more memory). If `memmap=True`, the uncertainty is never computed. + If `True` (not the default), will compute the uncertainty for the data (slower and uses more memory). + If you set ``memmap`` to be True, the uncertainty is never computed. memmap : `bool`, optional If `True` (the default), will not load arrays into memory, and will only read from the file into memory when needed. This option is faster and uses a @@ -60,38 +127,37 @@ def read_files(filename, *, spectral_windows=None, uncertainty=False, memmap=Tru Returns ------- - The corresponding `irispy.sji.SJICube` or `irispy.spectrogram.SpectrogramCube`. + dict + The dictionary contains the following keys: + "raster" + All loaded `irispy.spectrogram.SpectrogramCube` or `irispy.spectrogram.Collection` + "sji" + All loaded `irispy.sji.SJICube` or `irispy.spectrogram.SpectrogramCube`. """ - from irispy.io.sji import read_sji_lvl2 - from irispy.io.spectrograph import read_spectrograph_lvl2 - - if isinstance(filename, Path): - filename = str(filename) - if isinstance(filename, str | Path): - if tarfile.is_tarfile(filename): - path = Path(filename.replace(".tar.gz", "")) - path.mkdir(parents=True, exist_ok=True) - with tarfile.open(filename, "r") as tar: - tar.extractall(path) - filename = [path / file for file in tar.getnames()] + if isinstance(filepath, Path | str): + filepath = [Path(filepath)] + data_cubes = {"raster": [], "sji": []} + tar_files = [] + for file in filepath: + if tarfile.is_tarfile(file): + tar_files.append(file) + log.debug(f"Extracting {file}") + filepath.remove(file) + filepath.extend(tar_extract_all(file)) + for file in filepath: + instrument = fits.getval(file, "INSTRUME") + if instrument == "SJI": + data_cubes["sji"].append(read_sji_lvl2(file, memmap=memmap, uncertainty=uncertainty)) + elif instrument == "SPEC": + data_cubes["raster"].append( + read_spectrograph_lvl2( + file, + spectral_windows=spectral_windows, + memmap=memmap, + uncertainty=uncertainty, + ) + ) else: - filename = [filename] - intrume = fits.getval(filename[0], "INSTRUME") - all_instrume = [fits.getval(f, "INSTRUME") for f in filename] - if any(intrume != i for i in all_instrume): - msg = "You cannot mix raster and SJI files." - raise ValueError(msg) - if intrume == "SJI": - if len(filename) > 1: - msg = "You cannot load more than one SJI file at a time." + msg = f"Unsupported instrument: {instrument}" raise ValueError(msg) - return read_sji_lvl2(filename[0], memmap=memmap, uncertainty=uncertainty) - if intrume == "SPEC": - return read_spectrograph_lvl2( - filename, - spectral_windows=spectral_windows, - memmap=memmap, - uncertainty=uncertainty, - ) - msg = f"Unsupported instrument: {intrume}" - raise ValueError(msg) + return data_cubes diff --git a/irispy/sji.py b/irispy/sji.py index 60b8765..2e5465a 100644 --- a/irispy/sji.py +++ b/irispy/sji.py @@ -1,7 +1,7 @@ -import logging import textwrap import matplotlib.pyplot as plt +from sunpy import log from sunraster import SpectrogramCube from irispy.utils import calculate_dust_mask @@ -117,7 +117,7 @@ def plot(self, *args, **kwargs): try: cmap = plt.get_cmap(name=f"irissji{int(self.meta['TWAVE1'])}") except Exception as e: # NOQA: BLE001 - logging.debug(e) + log.debug(e) cmap = "viridis" kwargs["cmap"] = cmap ax = Plotter(ndcube=self).plot(*args, **kwargs) diff --git a/irispy/spectrograph.py b/irispy/spectrograph.py index b4d02f4..e0cdcab 100644 --- a/irispy/spectrograph.py +++ b/irispy/spectrograph.py @@ -1,4 +1,3 @@ -import logging import textwrap import astropy.units as u @@ -7,6 +6,7 @@ from astropy.coordinates import SkyCoord from astropy.time import Time from ndcube import NDCollection +from sunpy import log from sunpy.coordinates import Helioprojective from sunraster import SpectrogramCube as SpecCube from sunraster import SpectrogramSequence as SpecSeq @@ -134,7 +134,7 @@ def plot(self, *args, **kwargs): try: cmap = plt.get_cmap(name=f"irissji{int(self.meta.detector[:3])}") except Exception as e: # NOQA: BLE001 - logging.debug(e) + log.debug(e) cmap = "viridis" kwargs["cmap"] = cmap if len(self.dimensions) == 1: