diff --git a/irispy/new_sji.py b/irispy/new_sji.py index 3cbb0941..7fecebc0 100644 --- a/irispy/new_sji.py +++ b/irispy/new_sji.py @@ -12,6 +12,8 @@ from sunpy.time import parse_time from ndcube import NDCube from ndcube.utils.cube import convert_extra_coords_dict_to_input_format +from ndcube import utils +from ndcube.ndcube_sequence import NDCubeSequence from irispy import iris_tools @@ -101,11 +103,14 @@ def __repr__(self): #Conversion of the end date of OBS endobs = self.meta.get("ENDOBS", None) endobs = endobs.isoformat() if endobs else None - #Conversion of the instance start of OBS - instance_start = self.extra_coords["TIME"]["value"][0] + #Conversion of the instance start and end of OBS + if isinstance(self.extra_coords["TIME"]["value"], np.ndarray): + instance_start = self.extra_coords["TIME"]["value"][0] + instance_end = self.extra_coords["TIME"]["value"][-1] + else: + instance_start = self.extra_coords["TIME"]["value"] + instance_end = self.extra_coords["TIME"]["value"] instance_start = instance_start.isoformat() if instance_start else None - #Conversion of the instance end of OBS - instance_end = self.extra_coords["TIME"]["value"][-1] instance_end = instance_end.isoformat() if instance_end else None #Representation of IRISMapCube object return ( @@ -204,14 +209,93 @@ def apply_exposure_time_correction(self, undo=False, force=False): self.missing_axis)) -def read_iris_sji_level2_fits(filename, memmap=False): +class IRISMapCubeSequence(NDCubeSequence): + """Class for holding, slicing and plotting IRIS SJI data. + + This class contains all the functionality of its super class with + some additional functionalities. + + Parameters + ---------- + data_list: `list` + List of `IRISMapCube` objects from the same OBS ID. + Must also contain the 'detector type' in its meta attribute. + + meta: `dict` or header object + Metadata associated with the sequence. + + common_axis: `int` + The axis of the NDCubes corresponding to time. + + """ + def __init__(self, data_list, meta=None, common_axis=0): + # Check that all SJI data are coming from the same OBS ID. + if np.any([cube.meta["OBSID"] != data_list[0].meta["OBSID"] for cube in data_list]): + raise ValueError("Constituent IRISMapCube objects must have same " + "value of 'OBSID' in its meta.") + # Initialize Sequence. + super().__init__(data_list, meta=meta, common_axis=common_axis) + + def __repr__(self): + #Conversion of the start date of OBS + startobs = self.meta.get("STARTOBS", None) + startobs = startobs.isoformat() if startobs else None + #Conversion of the end date of OBS + endobs = self.meta.get("ENDOBS", None) + endobs = endobs.isoformat() if endobs else None + #Conversion of the instance start of OBS + instance_start = self[0].extra_coords["TIME"]["value"] + instance_start = instance_start.isoformat() if instance_start else None + #Conversion of the instance end of OBS + instance_end = self[-1].extra_coords["TIME"]["value"] + instance_end = instance_end.isoformat() if instance_end else None + #Representation of IRISMapCube object + return """ +IRISMapCubeSequence +--------------------- +Observatory:\t\t {obs} +Instrument:\t\t {instrument} + +OBS ID:\t\t\t {obs_id} +OBS Description:\t {obs_desc} +OBS period:\t\t {obs_start} -- {obs_end} + +Sequence period:\t {inst_start} -- {inst_end} +Sequence Shape:\t\t {seq_shape} +Axis Types:\t\t {axis_types} + +""".format(obs=self.meta.get('TELESCOP', None), + instrument=self.meta.get('INSTRUME', None), + obs_id=self.meta.get("OBSID", None), + obs_desc=self.meta.get("OBS_DESC", None), + obs_start=startobs, + obs_end=endobs, + inst_start=instance_start, + inst_end=instance_end, + seq_shape=self.dimensions, + axis_types=self.world_axis_physical_types) + + def __getitem__(self, item): + return self.index_as_cube[item] + + @property + def dimensions(self): + return self.cube_like_dimensions + + @property + def world_axis_physical_types(self): + return self.cube_like_world_axis_physical_types + + +def read_iris_sji_level2_fits(filenames, memmap=False): """ Read IRIS level 2 SJI FITS from an OBS into an IRISMapCube instance Parameters ---------- - filename : `str` - File name to be read + filenames: `list` of `str` or `str` + Filename or filenames to be read. They must all be associated with the same + OBS number. memmap : `bool` Default value is `False`. @@ -222,62 +306,68 @@ def read_iris_sji_level2_fits(filename, memmap=False): result: 'irispy.sji.IRISMapCube' """ - - # Open a fits file - my_file = fits.open(filename, memmap=memmap, do_not_scale_image_data=memmap) - # Derive WCS, data and mask for NDCube from fits file. - wcs = WCS(my_file[0].header) - data = my_file[0].data - data_nan_masked = my_file[0].data - if memmap: - data_nan_masked[data == BAD_PIXEL_VALUE_UNSCALED] = 0 - mask = None - scaled = False - unit = iris_tools.DN_UNIT["SJI_UNSCALED"] - uncertainty = None - elif not memmap: - data_nan_masked[data == BAD_PIXEL_VALUE_SCALED] = np.nan - mask = data_nan_masked == BAD_PIXEL_VALUE_SCALED - scaled = True - # Derive unit and readout noise from the detector - unit = iris_tools.DN_UNIT["SJI"] - readout_noise = iris_tools.READOUT_NOISE["SJI"] - # Derive uncertainty of data for NDCube from fits file. - uncertainty = u.Quantity(np.sqrt((data_nan_masked*unit).to(u.photon).value - + readout_noise.to(u.photon).value**2), - unit=u.photon).to(unit).value - # Derive exposure time from detector. - exposure_times = my_file[1].data[:, my_file[1].header["EXPTIMES"]] - # Derive extra coordinates for NDCube from fits file. - times = np.array([parse_time(my_file[0].header["STARTOBS"]) - + timedelta(seconds=s) - for s in my_file[1].data[:, my_file[1].header["TIME"]]]) - pztx = my_file[1].data[:, my_file[1].header["PZTX"]] * u.arcsec - pzty = my_file[1].data[:, my_file[1].header["PZTY"]] * u.arcsec - xcenix = my_file[1].data[:, my_file[1].header["XCENIX"]] * u.arcsec - ycenix = my_file[1].data[:, my_file[1].header["YCENIX"]] * u.arcsec - obs_vrix = my_file[1].data[:, my_file[1].header["OBS_VRIX"]] * u.m/u.s - ophaseix = my_file[1].data[:, my_file[1].header["OPHASEIX"]] - extra_coords = [('TIME', 0, times), ("PZTX", 0, pztx), ("PZTY", 0, pzty), - ("XCENIX", 0, xcenix), ("YCENIX", 0, ycenix), - ("OBS_VRIX", 0, obs_vrix), ("OPHASEIX", 0, ophaseix), - ("EXPOSURE TIME", 0, exposure_times)] - # Extraction of meta for NDCube from fits file. - startobs = my_file[0].header.get('STARTOBS', None) - startobs = parse_time(startobs) if startobs else None - endobs = my_file[0].header.get('ENDOBS', None) - endobs = parse_time(endobs) if endobs else None - meta = {'TELESCOP': my_file[0].header.get('TELESCOP', None), - 'INSTRUME': my_file[0].header.get('INSTRUME', None), - 'TWAVE1': my_file[0].header.get('TWAVE1', None), - 'STARTOBS': startobs, - 'ENDOBS': endobs, - 'NBFRAMES': my_file[0].data.shape[0], - 'OBSID': my_file[0].header.get('OBSID', None), - 'OBS_DESC': my_file[0].header.get('OBS_DESC', None)} - - my_file.close() - - return IRISMapCube(data_nan_masked, wcs, uncertainty=uncertainty, - unit=unit, meta=meta, mask=mask, extra_coords=extra_coords, - scaled=scaled) + list_of_cubes = [] + if type(filenames) is str: + filenames = [filenames] + for filename in filenames: + # Open a fits file + hdulist = fits.open(filename, memmap=memmap, do_not_scale_image_data=memmap) + hdulist.verify('fix') + # Derive WCS, data and mask for NDCube from fits file. + wcs = WCS(hdulist[0].header) + data = hdulist[0].data + data_nan_masked = hdulist[0].data + if memmap: + data_nan_masked[data == BAD_PIXEL_VALUE_UNSCALED] = 0 + mask = None + scaled = False + unit = iris_tools.DN_UNIT["SJI_UNSCALED"] + uncertainty = None + elif not memmap: + data_nan_masked[data == BAD_PIXEL_VALUE_SCALED] = np.nan + mask = data_nan_masked == BAD_PIXEL_VALUE_SCALED + scaled = True + # Derive unit and readout noise from the detector + unit = iris_tools.DN_UNIT["SJI"] + readout_noise = iris_tools.READOUT_NOISE["SJI"] + # Derive uncertainty of data for NDCube from fits file. + uncertainty = u.Quantity(np.sqrt((data_nan_masked*unit).to(u.photon).value + + readout_noise.to(u.photon).value**2), + unit=u.photon).to(unit).value + # Derive exposure time from detector. + exposure_times = hdulist[1].data[:, hdulist[1].header["EXPTIMES"]] + # Derive extra coordinates for NDCube from fits file. + times = np.array([parse_time(hdulist[0].header["STARTOBS"]) + + timedelta(seconds=s) + for s in hdulist[1].data[:, hdulist[1].header["TIME"]]]) + pztx = hdulist[1].data[:, hdulist[1].header["PZTX"]] * u.arcsec + pzty = hdulist[1].data[:, hdulist[1].header["PZTY"]] * u.arcsec + xcenix = hdulist[1].data[:, hdulist[1].header["XCENIX"]] * u.arcsec + ycenix = hdulist[1].data[:, hdulist[1].header["YCENIX"]] * u.arcsec + obs_vrix = hdulist[1].data[:, hdulist[1].header["OBS_VRIX"]] * u.m/u.s + ophaseix = hdulist[1].data[:, hdulist[1].header["OPHASEIX"]] + extra_coords = [('TIME', 0, times), ("PZTX", 0, pztx), ("PZTY", 0, pzty), + ("XCENIX", 0, xcenix), ("YCENIX", 0, ycenix), + ("OBS_VRIX", 0, obs_vrix), ("OPHASEIX", 0, ophaseix), + ("EXPOSURE TIME", 0, exposure_times)] + # Extraction of meta for NDCube from fits file. + startobs = hdulist[0].header.get('STARTOBS', None) + startobs = parse_time(startobs) if startobs else None + endobs = hdulist[0].header.get('ENDOBS', None) + endobs = parse_time(endobs) if endobs else None + meta = {'TELESCOP': hdulist[0].header.get('TELESCOP', None), + 'INSTRUME': hdulist[0].header.get('INSTRUME', None), + 'TWAVE1': hdulist[0].header.get('TWAVE1', None), + 'STARTOBS': startobs, + 'ENDOBS': endobs, + 'NBFRAMES': hdulist[0].data.shape[0], + 'OBSID': hdulist[0].header.get('OBSID', None), + 'OBS_DESC': hdulist[0].header.get('OBS_DESC', None)} + list_of_cubes.append(IRISMapCube(data_nan_masked, wcs, uncertainty=uncertainty, + unit=unit, meta=meta, mask=mask, + extra_coords=extra_coords, scaled=scaled)) + hdulist.close() + if len(filenames) == 1: + return list_of_cubes[0] + else: + return IRISMapCubeSequence(list_of_cubes, meta=meta, common_axis=0)