Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Creation of IRISMapCubeSequence class #83

Merged
merged 8 commits into from
May 16, 2018
Merged

Creation of IRISMapCubeSequence class #83

Changes from 1 commit
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
718b467
First version of IRISMapCubeSequence
BaptistePellorceAstro May 7, 2018
945c9b0
Return IRISMapCube or IRISMapCubeSequence
BaptistePellorceAstro May 7, 2018
e613257
Some new properties
BaptistePellorceAstro May 8, 2018
2b29ab3
Rewriting of properties
BaptistePellorceAstro May 14, 2018
d36a565
Removing commented lines
BaptistePellorceAstro May 14, 2018
ccb98fc
Removing trailing spaces
BaptistePellorceAstro May 14, 2018
fa0c541
Bugs fixed on getitem and dimensions
BaptistePellorceAstro May 15, 2018
757d07c
Method name change
BaptistePellorceAstro May 16, 2018
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
197 changes: 135 additions & 62 deletions irispy/new_sji.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,7 @@
from sunpy.time import parse_time
from ndcube import NDCube
from ndcube.utils.cube import convert_extra_coords_dict_to_input_format
from ndcube.ndcube_sequence import NDCubeSequence

from irispy import iris_tools

Expand Down Expand Up @@ -204,14 +205,82 @@ 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 len(np.unique([cube.meta["OBSID"] for cube in data_list])) != 1:
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This could be rewritten as:

if any([cube.meta["OBSID"] for cube in data_list] != data_list[0].meta["OBSID"]):

The list might have to be an array. I'm not sure. Check this out.

Copy link
Contributor Author

@BaptistePellorceAstro BaptistePellorceAstro May 7, 2018 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

With this line, I get an error that I solve by doing :

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"][0]
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"][-1]
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 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 of filenames to be read. They must all be associated with the same
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Typo: "of" should be "or"

OBS number.

memmap : `bool`
Default value is `False`.
Expand All @@ -222,62 +291,66 @@ 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()

return IRISMapCubeSequence(list_of_cubes, meta=meta, common_axis=0)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think an IRISMapCubeSequence should only be returned if there is more than one file input. Otherwise, an IRISMapCube should be returned as before.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ok, this is boring to use it with only one file (I am working with). I will modify it.