1664

pending completion

Build # 1664

Build Type

Pull #257

travis-ci

Committed by

web-flow

Commit Message

<a href="https://github.com/LCOGT/banzai/commit/<a class=hub.com/LCOGT/banzai/commit/1f06634c885c93fc3036fc6dee40a1b52be1590f">1f06634c8<a href="https://github.com/LCOGT/banzai/commit/1f06634c885c93fc3036fc6dee40a1b52be1590f">">Merge </a><a class="double-link" href="https://github.com/LCOGT/banzai/commit/<a class="double-link" href="https://github.com/LCOGT/banzai/commit/370048bffc9bd12ed45e88a24e8d42a54c73ff5d">370048bff</a>">370048bff</a><a href="https://github.com/LCOGT/banzai/commit/1f06634c885c93fc3036fc6dee40a1b52be1590f"> into 1bf5aec56">1bf5aec56</a>

Pull Request Pull Request #257: Seek cursor to beginning of buffer before returning.

Run Details

18 of 18 new or added lines in 6 files covered. (100.0%)

1769 of 2814 relevant lines covered (62.86%)

0.63 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

89.41

/banzai/data.py

import abc
import tempfile
from typing import Union, Type

import numpy as np
from astropy.io import fits
from astropy.table import Table

from banzai.utils.image_utils import Section
from banzai.utils import fits_utils, stats


class Data(metaclass=abc.ABCMeta):
    _file_handles = []

    def __init__(self, data: Union[np.array, Table], meta: Union[dict, fits.Header],
                 mask: np.array = None, name: str = '', memmap=True):
        self.memmap = memmap
        self.data = self._init_array(data)
        self.meta = meta.copy()
        self._validate_mask(mask)
        self.mask = self._init_array(mask, dtype=np.uint8)
        self.name = name

    def _validate_mask(self, mask):
        if mask is not None:
            if mask.shape != self.data.shape:
                raise ValueError('Mask must have the same dimensions as the data')

    def _init_array(self, array: np.array = None, dtype: Type = None):
        if not self.memmap:
            return array
        file_handle = tempfile.NamedTemporaryFile('w+b')
        if array is None:
            shape = self.data.shape
            if dtype is None:
                dtype = self.data.dtype
            array = np.zeros(shape, dtype=dtype)
        if array.size > 0:
            memory_mapped_array = np.memmap(file_handle, shape=array.shape, dtype=array.dtype, mode='readwrite')
            memory_mapped_array.ravel()[:] = array.ravel()[:]
            self._file_handles.append(file_handle)
        else:
            memory_mapped_array = array
        return memory_mapped_array

    def add_mask(self, mask: np.array):
        self._validate_mask(mask)
        self.mask = self._init_array(mask)

    def __del__(self):
        for handle in self._file_handles:
            handle.close()
        del self.data
        del self.mask

    @property
    def extension_name(self):
        return self.meta.get('EXTNAME')

    @extension_name.setter
    def extension_name(self, value):
        self.meta['EXTNAME'] = value

    @classmethod
    def from_fits(cls, hdu: Union[fits.ImageHDU, fits.TableHDU, fits.BinTableHDU]):
        return cls(hdu.data, hdu.header, name=hdu.header.get('EXTNAME'))

    @abc.abstractmethod
    def to_fits(self, context) -> Union[fits.HDUList, list]:
        pass


class HeaderOnly(Data):
    def __init__(self, meta: Union[dict, fits.Header]):
        super().__init__(data=np.zeros(0), meta=meta, memmap=False)

    def to_fits(self, context):
        return fits.HDUList([fits.ImageHDU(data=None, header=self.meta)])


class DataTable(Data):
    def __init__(self, data, name, meta=None):
        self.data = data
        if meta is None:
            self.meta = fits.BinTableHDU(data).header
        else:
            self.meta=meta
        self.name = name

    def to_fits(self, context) -> Union[fits.HDUList, list]:
        hdu = fits.BinTableHDU(self.data)
        hdu.name = self.name
        # Put in the description keywords
        for k in self.meta.keys():
            if 'TTYPE' in k:
                column_name = self.meta[k].lower()
                description = self.data[column_name].description
                hdu.header[k] = (column_name.upper(), description)
                # Get the value of n in TTYPEn
                n = k[5:]
                hdu.header['TCOMM{0}'.format(n)] = description

        return [hdu]


class ArrayData(Data):
    def __init__(self, data, meta, name, memmap=True):
        super().__init__(data, meta, name=name, memmap=memmap)

    def to_fits(self, context) -> Union[fits.HDUList, list]:
        return [fits.ImageHDU(data=self.data, header=self.meta)]


class CCDData(Data):
    def __init__(self, data: Union[np.array, Table], meta: fits.Header,
                 mask: np.array = None, name: str = '', uncertainty: np.array = None, memmap=True):
        super().__init__(data=data, meta=meta, mask=mask, name=name, memmap=memmap)
        if uncertainty is None:
            uncertainty = self.read_noise * np.ones(data.shape, dtype=data.dtype) / self.gain
        self.uncertainty = self._init_array(uncertainty)
        self._detector_section = Section.parse_region_keyword(self.meta.get('DETSEC'))
        self._data_section = Section.parse_region_keyword(self.meta.get('DATASEC'))
        self._background = None

    def __getitem__(self, section):
        """
        Return a new CCDData object with the given section of data
        :param section: needs to be  in data coords
        :return:
        """
        return self.trim(trim_section=section)

    def __imul__(self, value):
        self.data *= value
        self.uncertainty *= value
        self.meta['SATURATE'] *= value
        self.meta['GAIN'] *= value
        self.meta['MAXLIN'] *= value
        return self

    def __itruediv__(self, value):
        if isinstance(value, CCDData):
            self.uncertainty = np.abs(self.data / value.data) * \
                               np.sqrt((self.uncertainty / self.data) ** 2 + (value.uncertainty / value.data) ** 2)
            self.data = np.divide(self.data, value.data)
            self.meta['SATURATE'] /= value.meta['SATURATE']
            self.meta['GAIN'] /= value.meta['GAIN']
            self.meta['MAXLIN'] /= value.meta['MAXLIN']
            self.mask |= value.mask
        else:
            self.__imul__(1.0 / value)
        return self

    def to_fits(self, context):
        data_hdu = fits.ImageHDU(data=self.data, header=fits.Header(self.meta))
        bpm_hdu = fits_utils.to_fits_image_extension(self.mask, self.extension_name, 'BPM', context,
                                                     extension_version=self.meta.get('EXTVER'))
        uncertainty_hdu = fits_utils.to_fits_image_extension(self.uncertainty, self.extension_name, 'ERR', context,
                                                             extension_version=self.meta.get('EXTVER'))
        hdulist = fits.HDUList([data_hdu, bpm_hdu, uncertainty_hdu])
        return hdulist

    def __del__(self):
        super().__del__()
        del self.uncertainty

    def __isub__(self, value):
        if isinstance(value, CCDData):
            self.data -= value.data
            self.uncertainty = np.sqrt(value.uncertainty * value.uncertainty + self.uncertainty * self.uncertainty)
            self.mask |= value.mask
        else:
            self.data -= value
        return self

    def __sub__(self, other):
        uncertainty = np.sqrt(self.uncertainty * self.uncertainty + other.uncertainty * other.uncertainty)
        return type(self)(data=self.data - other.data, meta=self.meta, mask=self.mask|other.mask,
                          uncertainty=uncertainty)

    def signal_to_noise(self):
        return np.abs(self.data) / self.uncertainty

    def get_overscan_region(self):
        return Section.parse_region_keyword(self.meta.get('BIASSEC', 'N/A'))

    def get_inner_image_section(self, inner_edge_width=0.25):
        """
         Extract the inner section of the image with dimensions:
         ny * inner_edge_width * 2.0 x nx * inner_edge_width * 2.0

         Parameters
         ----------

         inner_edge_width: float
                           Size of inner edge as fraction of total image size

         Returns
         -------
         inner_section: array
                        Inner section of image
         """
        inner_nx = round(self.data.shape[1] * inner_edge_width)
        inner_ny = round(self.data.shape[0] * inner_edge_width)
        return self.data[inner_ny: -inner_ny, inner_nx: -inner_nx]

    def trim(self, trim_section=None):
        """
        :param trim_section: Always in data coords
        :return:
        """
        if trim_section is None:
            trim_section = Section.parse_region_keyword(self.meta.get('TRIMSEC', 'N/A'))

        trimmed_image = type(self)(data=self.data[trim_section.to_slice()], meta=self.meta,
                                   mask=self.mask[trim_section.to_slice()], name=self.name,
                                   uncertainty=self.uncertainty[trim_section.to_slice()], memmap=self.memmap)
        trimmed_image.detector_section = self.data_to_detector_section(trim_section)
        trimmed_image.data_section = Section(x_start=1, y_start=1,
                                             x_stop=trimmed_image.data.shape[1],
                                             y_stop=trimmed_image.data.shape[0])
        return trimmed_image

    @property
    def dtype(self):
        return self.data.dtype

    @property
    def shape(self):
        return self.data.shape

    @property
    def gain(self):
        return self.meta.get('GAIN', 1.0)

    @gain.setter
    def gain(self, value):
        self.meta['GAIN'] = value

    @property
    def saturate(self):
        return self.meta.get('SATURATE')

    @saturate.setter
    def saturate(self, value):
        self.meta['SATURATE'] = value

    @property
    def max_linearity(self):
        return self.meta.get('MAXLIN')

    @max_linearity.setter
    def max_linearity(self, value):
        self.meta['MAXLIN'] = value

    @property
    def read_noise(self):
        return self.meta.get('RDNOISE', 0.0)

    @property
    def binning(self):
        return [int(b) for b in self.meta.get('CCDSUM', '1 1').split(' ')]

    @binning.setter
    def binning(self, value):
        x_binning, y_binning = value
        self.meta['CCDSUM'] = f'{x_binning} {y_binning}'

    @property
    def detector_section(self):
        return self._detector_section

    @detector_section.setter
    def detector_section(self, section):
        self.meta['DETSEC'] = section.to_region_keyword()
        self._detector_section = section

    @property
    def data_section(self):
        return self._data_section

    @data_section.setter
    def data_section(self, section):
        self.meta['DATASEC'] = section.to_region_keyword()
        self._data_section = section

    def rebin(self, binning):
        # TODO: Implement me
        return self

    def get_overlap(self, detector_section):
        return self.detector_section.overlap(detector_section)

    def detector_to_data_section_oned(self, axis, section):
        binning_indices = {'x': 0, 'y': 1}
        sign = np.sign(getattr(self.detector_section, f'{axis}_stop') - getattr(self.detector_section, f'{axis}_start'))
        sign *= np.sign(getattr(self.data_section, f'{axis}_stop') - getattr(self.data_section, f'{axis}_start'))

        start = sign * (getattr(section, f'{axis}_start') - getattr(self.detector_section, f'{axis}_start'))
        start //= self.binning[binning_indices[axis]]
        start += getattr(self.data_section, f'{axis}_start')

        stop = sign * (getattr(section, f'{axis}_stop') - getattr(self.detector_section, f'{axis}_start'))
        stop //= self.binning[binning_indices[axis]]
        stop += getattr(self.data_section, f'{axis}_start')
        return start, stop

    def detector_to_data_section(self, section):
        """Given a detector region, figure out the corresponding data region
        Note the + and - 1 factors cancel
        Really this is just doing the same thing as a CD matrix calculation
        r_data - data_0 = M (r_det - det_0) where M is a transformation matrix and everything else is a 2d vector
        M is either the 1/binning * identity or the 1/binning * negative Identity depending on if datasec and detsec
        are in the same ordering (both increasing or both decreasing gives positive)
        Note if you want to add a rotation, M can contain that as well.
        """
        x_start, x_stop = self.detector_to_data_section_oned('x', section)
        y_start, y_stop = self.detector_to_data_section_oned('y', section)

        return Section(x_start, x_stop, y_start, y_stop)

    def data_to_detector_section(self, section):
        """Given a data region, get the detector section that this covers.
        This is the inverse of get_data section"""
        x_start, x_stop = self.data_to_detector_section_oned('x', section)
        y_start, y_stop = self.data_to_detector_section_oned('y', section)

        return Section(x_start, x_stop, y_start, y_stop)

    def data_to_detector_section_oned(self, axis, section):
        binning_indices = {'x': 0, 'y': 1}
        sign = np.sign(getattr(self.detector_section, f'{axis}_stop') - getattr(self.detector_section, f'{axis}_start'))
        sign *= np.sign(getattr(self.data_section, f'{axis}_stop') - getattr(self.data_section, f'{axis}_start'))

        start = sign * (getattr(section, f'{axis}_start') - getattr(self.data_section, f'{axis}_start'))
        start *= self.binning[binning_indices[axis]]
        start += getattr(self.detector_section, f'{axis}_start')

        stop = sign * (getattr(section, f'{axis}_stop') - getattr(self.data_section, f'{axis}_start'))
        stop *= self.binning[binning_indices[axis]]
        stop += getattr(self.detector_section, f'{axis}_start')

        # when converting from binned to unbinned coordinates, we need to adjust our stopping pixel
        # by (binning - 1)
        stop += np.sign(stop - start) * (self.binning[binning_indices[axis]] - 1)

        return start, stop

    def copy_in(self, data):
        """
        Copy in the data from another CCDData object based on the detector sections

        :param data_to_copy:
        :return:
        """
        overlap_section = self.get_overlap(data.detector_section)
        data_to_copy = data.trim(trim_section=data.detector_to_data_section(overlap_section))
        data_to_copy = data_to_copy.rebin(self.binning)
        for array_name_to_copy in ['data', 'mask', 'uncertainty']:
            array_to_copy = getattr(data_to_copy, array_name_to_copy)
            my_overlap = self.detector_to_data_section(overlap_section).to_slice()
            getattr(self, array_name_to_copy)[my_overlap][:] = array_to_copy[:]

    def init_poisson_uncertainties(self):
        self.uncertainty = np.sqrt(self.uncertainty ** 2.0 + np.abs(self.data))

    @property
    def background(self):
        return self._background

    @background.setter
    def background(self, value):
        if self._background is not None:
            self.data += self._background
        self._background = value
        self.data -= self._background


def stack(data_to_stack, nsigma_reject) -> CCDData:
    """
    """
    shape3d = [len(data_to_stack)] + list(data_to_stack[0].shape)
    a = np.zeros(shape3d, dtype=data_to_stack[0].dtype)
    uncertainties = np.zeros(shape3d, dtype=data_to_stack[0].dtype)
    mask = np.zeros(shape3d, dtype=np.uint8)

    for i, data in enumerate(data_to_stack):
        a[i, :, :] = data.data[:, :]
        mask[i, :, :] = data.mask[:, :]
        uncertainties[i, :, :] = data.uncertainty[:, :]

    abs_deviation = stats.absolute_deviation(a, axis=0, mask=mask)

    robust_std = stats.robust_standard_deviation(a, axis=0, abs_deviation=abs_deviation, mask=mask)

    robust_std = np.expand_dims(robust_std, axis=0)

    # Mask values that are N sigma from the median
    sigma_mask = abs_deviation > (nsigma_reject * robust_std)

    mask3d = np.logical_or(sigma_mask, mask > 0)
    n_good_pixels = np.logical_not(mask3d).sum(axis=0)

    stacked_mask = np.zeros(n_good_pixels.shape, dtype=np.uint8)

    # If a pixel is bad in all images, make sure we don't divide by zero
    bad_pixels = n_good_pixels == 0

    # If pixel is bad in all of the images, we take the logical or of all of the bits to go in the final mask
    stacked_mask[bad_pixels] = np.bitwise_or.reduce(mask, axis=0)[bad_pixels]

    # If a pixel is bad in all images, fill that pixel with the mean from the images
    n_good_pixels[bad_pixels] = len(data_to_stack)
    mask3d[:, bad_pixels] = False

    a[mask3d] = 0.0
    stacked_data = a.sum(axis=0) / n_good_pixels

    # Again if a pixel is bad in all images, fill the uncertainties with the quadrature sum / N images
    uncertainties[mask3d] = 0.0
    uncertainties *= uncertainties
    stacked_uncertainty = np.sqrt(uncertainties.sum(axis=0) / (n_good_pixels ** 2.0))

    return CCDData(data=stacked_data, meta=data_to_stack[0].meta, uncertainty=stacked_uncertainty, mask=stacked_mask)

1	import abc	1✔
2	import tempfile	1✔
3	from typing import Union, Type	1✔
4
5	import numpy as np	1✔
6	from astropy.io import fits	1✔
7	from astropy.table import Table	1✔
8
9	from banzai.utils.image_utils import Section	1✔
10	from banzai.utils import fits_utils, stats	1✔
11
12
13	class Data(metaclass=abc.ABCMeta):	1✔
14	_file_handles = []	1✔
15
16	def __init__(self, data: Union[np.array, Table], meta: Union[dict, fits.Header],	1✔
17	mask: np.array = None, name: str = '', memmap=True):
18	self.memmap = memmap	1✔
19	self.data = self._init_array(data)	1✔
20	self.meta = meta.copy()	1✔
21	self._validate_mask(mask)	1✔
22	self.mask = self._init_array(mask, dtype=np.uint8)	1✔
23	self.name = name	1✔
24
25	def _validate_mask(self, mask):	1✔
26	if mask is not None:	1✔
27	if mask.shape != self.data.shape:	1✔
28	raise ValueError('Mask must have the same dimensions as the data')	1✔
29
30	def _init_array(self, array: np.array = None, dtype: Type = None):	1✔
31	if not self.memmap:	1✔
32	return array	1✔
33	file_handle = tempfile.NamedTemporaryFile('w+b')	1✔
34	if array is None:	1✔
35	shape = self.data.shape	1✔
36	if dtype is None:	1✔
37	dtype = self.data.dtype	×
38	array = np.zeros(shape, dtype=dtype)	1✔
39	if array.size > 0:	1✔
40	memory_mapped_array = np.memmap(file_handle, shape=array.shape, dtype=array.dtype, mode='readwrite')	1✔
41	memory_mapped_array.ravel()[:] = array.ravel()[:]	1✔
42	self._file_handles.append(file_handle)	1✔
43	else:
44	memory_mapped_array = array	×
45	return memory_mapped_array	1✔
46
47	def add_mask(self, mask: np.array):	1✔
48	self._validate_mask(mask)	1✔
49	self.mask = self._init_array(mask)	1✔
50
51	def __del__(self):	1✔
52	for handle in self._file_handles:	1✔
53	handle.close()	1✔
54	del self.data	1✔
55	del self.mask	1✔
56
57	@property	1✔
58	def extension_name(self):
59	return self.meta.get('EXTNAME')	1✔
60
61	@extension_name.setter	1✔
62	def extension_name(self, value):
63	self.meta['EXTNAME'] = value	1✔
64
65	@classmethod	1✔
66	def from_fits(cls, hdu: Union[fits.ImageHDU, fits.TableHDU, fits.BinTableHDU]):	1✔
67	return cls(hdu.data, hdu.header, name=hdu.header.get('EXTNAME'))	×
68
69	@abc.abstractmethod	1✔
70	def to_fits(self, context) -> Union[fits.HDUList, list]:	1✔
71	pass	×
72
73
74	class HeaderOnly(Data):	1✔
75	def __init__(self, meta: Union[dict, fits.Header]):	1✔
76	super().__init__(data=np.zeros(0), meta=meta, memmap=False)	1✔
77
78	def to_fits(self, context):	1✔
79	return fits.HDUList([fits.ImageHDU(data=None, header=self.meta)])	×
80
81
82	class DataTable(Data):	1✔
83	def __init__(self, data, name, meta=None):	1✔
84	self.data = data	×
85	if meta is None:	×
86	self.meta = fits.BinTableHDU(data).header	×
87	else:
88	self.meta=meta	×
89	self.name = name	×
90
91	def to_fits(self, context) -> Union[fits.HDUList, list]:	1✔
92	hdu = fits.BinTableHDU(self.data)	×
93	hdu.name = self.name	×
94	# Put in the description keywords
95	for k in self.meta.keys():	×
96	if 'TTYPE' in k:	×
97	column_name = self.meta[k].lower()	×
98	description = self.data[column_name].description	×
99	hdu.header[k] = (column_name.upper(), description)	×
100	# Get the value of n in TTYPEn
101	n = k[5:]	×
102	hdu.header['TCOMM{0}'.format(n)] = description	×
103
104	return [hdu]	×
105
106
107	class ArrayData(Data):	1✔
108	def __init__(self, data, meta, name, memmap=True):	1✔
109	super().__init__(data, meta, name=name, memmap=memmap)	×
110
111	def to_fits(self, context) -> Union[fits.HDUList, list]:	1✔
112	return [fits.ImageHDU(data=self.data, header=self.meta)]	×
113
114
115	class CCDData(Data):	1✔
116	def __init__(self, data: Union[np.array, Table], meta: fits.Header,	1✔
117	mask: np.array = None, name: str = '', uncertainty: np.array = None, memmap=True):
118	super().__init__(data=data, meta=meta, mask=mask, name=name, memmap=memmap)	1✔
119	if uncertainty is None:	1✔
120	uncertainty = self.read_noise * np.ones(data.shape, dtype=data.dtype) / self.gain	1✔
121	self.uncertainty = self._init_array(uncertainty)	1✔
122	self._detector_section = Section.parse_region_keyword(self.meta.get('DETSEC'))	1✔
123	self._data_section = Section.parse_region_keyword(self.meta.get('DATASEC'))	1✔
124	self._background = None	1✔
125
126	def __getitem__(self, section):	1✔
127	"""
128	Return a new CCDData object with the given section of data
129	:param section: needs to be in data coords
130	:return:
131	"""
132	return self.trim(trim_section=section)	1✔
133
134	def __imul__(self, value):	1✔
135	self.data *= value	1✔
136	self.uncertainty *= value	1✔
137	self.meta['SATURATE'] *= value	1✔
138	self.meta['GAIN'] *= value	1✔
139	self.meta['MAXLIN'] *= value	1✔
140	return self	1✔
141
142	def __itruediv__(self, value):	1✔
143	if isinstance(value, CCDData):	1✔
144	self.uncertainty = np.abs(self.data / value.data) * \	1✔
145	np.sqrt((self.uncertainty / self.data) 2 + (value.uncertainty / value.data) 2)
146	self.data = np.divide(self.data, value.data)	1✔
147	self.meta['SATURATE'] /= value.meta['SATURATE']	1✔
148	self.meta['GAIN'] /= value.meta['GAIN']	1✔
149	self.meta['MAXLIN'] /= value.meta['MAXLIN']	1✔
150	self.mask \|= value.mask	1✔
151	else:
152	self.__imul__(1.0 / value)	1✔
153	return self	1✔
154
155	def to_fits(self, context):	1✔
156	data_hdu = fits.ImageHDU(data=self.data, header=fits.Header(self.meta))	1✔
157	bpm_hdu = fits_utils.to_fits_image_extension(self.mask, self.extension_name, 'BPM', context,	1✔
158	extension_version=self.meta.get('EXTVER'))
159	uncertainty_hdu = fits_utils.to_fits_image_extension(self.uncertainty, self.extension_name, 'ERR', context,	1✔
160	extension_version=self.meta.get('EXTVER'))
161	hdulist = fits.HDUList([data_hdu, bpm_hdu, uncertainty_hdu])	1✔
162	return hdulist	1✔
163
164	def __del__(self):	1✔
165	super().__del__()	1✔
166	del self.uncertainty	1✔
167
168	def __isub__(self, value):	1✔
169	if isinstance(value, CCDData):	1✔
170	self.data -= value.data	1✔
171	self.uncertainty = np.sqrt(value.uncertainty * value.uncertainty + self.uncertainty * self.uncertainty)	1✔
172	self.mask \|= value.mask	1✔
173	else:
174	self.data -= value	1✔
175	return self	1✔
176
177	def __sub__(self, other):	1✔
178	uncertainty = np.sqrt(self.uncertainty * self.uncertainty + other.uncertainty * other.uncertainty)	1✔
179	return type(self)(data=self.data - other.data, meta=self.meta, mask=self.mask\|other.mask,	1✔
180	uncertainty=uncertainty)
181
182	def signal_to_noise(self):	1✔
183	return np.abs(self.data) / self.uncertainty	1✔
184
185	def get_overscan_region(self):	1✔
186	return Section.parse_region_keyword(self.meta.get('BIASSEC', 'N/A'))	1✔
187
188	def get_inner_image_section(self, inner_edge_width=0.25):	1✔
189	"""
190	Extract the inner section of the image with dimensions:
191	ny * inner_edge_width * 2.0 x nx * inner_edge_width * 2.0
192
193	Parameters
194	----------
195
196	inner_edge_width: float
197	Size of inner edge as fraction of total image size
198
199	Returns
200	-------
201	inner_section: array
202	Inner section of image
203	"""
204	inner_nx = round(self.data.shape[1] * inner_edge_width)	1✔
205	inner_ny = round(self.data.shape[0] * inner_edge_width)	1✔
206	return self.data[inner_ny: -inner_ny, inner_nx: -inner_nx]	1✔
207
208	def trim(self, trim_section=None):	1✔
209	"""
210	:param trim_section: Always in data coords
211	:return:
212	"""
213	if trim_section is None:	1✔
214	trim_section = Section.parse_region_keyword(self.meta.get('TRIMSEC', 'N/A'))	1✔
215
216	trimmed_image = type(self)(data=self.data[trim_section.to_slice()], meta=self.meta,	1✔
217	mask=self.mask[trim_section.to_slice()], name=self.name,
218	uncertainty=self.uncertainty[trim_section.to_slice()], memmap=self.memmap)
219	trimmed_image.detector_section = self.data_to_detector_section(trim_section)	1✔
220	trimmed_image.data_section = Section(x_start=1, y_start=1,	1✔
221	x_stop=trimmed_image.data.shape[1],
222	y_stop=trimmed_image.data.shape[0])
223	return trimmed_image	1✔
224
225	@property	1✔
226	def dtype(self):
227	return self.data.dtype	1✔
228
229	@property	1✔
230	def shape(self):
231	return self.data.shape	1✔
232
233	@property	1✔
234	def gain(self):
235	return self.meta.get('GAIN', 1.0)	1✔
236
237	@gain.setter	1✔
238	def gain(self, value):
239	self.meta['GAIN'] = value	1✔
240
241	@property	1✔
242	def saturate(self):
243	return self.meta.get('SATURATE')	1✔
244
245	@saturate.setter	1✔
246	def saturate(self, value):
247	self.meta['SATURATE'] = value	1✔
248
249	@property	1✔
250	def max_linearity(self):
251	return self.meta.get('MAXLIN')	1✔
252
253	@max_linearity.setter	1✔
254	def max_linearity(self, value):
255	self.meta['MAXLIN'] = value	1✔
256
257	@property	1✔
258	def read_noise(self):
259	return self.meta.get('RDNOISE', 0.0)	1✔
260
261	@property	1✔
262	def binning(self):
263	return [int(b) for b in self.meta.get('CCDSUM', '1 1').split(' ')]	1✔
264
265	@binning.setter	1✔
266	def binning(self, value):
267	x_binning, y_binning = value	1✔
268	self.meta['CCDSUM'] = f'{x_binning} {y_binning}'	1✔
269
270	@property	1✔
271	def detector_section(self):
272	return self._detector_section	1✔
273
274	@detector_section.setter	1✔
275	def detector_section(self, section):
276	self.meta['DETSEC'] = section.to_region_keyword()	1✔
277	self._detector_section = section	1✔
278
279	@property	1✔
280	def data_section(self):
281	return self._data_section	1✔
282
283	@data_section.setter	1✔
284	def data_section(self, section):
285	self.meta['DATASEC'] = section.to_region_keyword()	1✔
286	self._data_section = section	1✔
287
288	def rebin(self, binning):	1✔
289	# TODO: Implement me
290	return self	1✔
291
292	def get_overlap(self, detector_section):	1✔
293	return self.detector_section.overlap(detector_section)	1✔
294
295	def detector_to_data_section_oned(self, axis, section):	1✔
296	binning_indices = {'x': 0, 'y': 1}	1✔
297	sign = np.sign(getattr(self.detector_section, f'{axis}_stop') - getattr(self.detector_section, f'{axis}_start'))	1✔
298	sign *= np.sign(getattr(self.data_section, f'{axis}_stop') - getattr(self.data_section, f'{axis}_start'))	1✔
299
300	start = sign * (getattr(section, f'{axis}_start') - getattr(self.detector_section, f'{axis}_start'))	1✔
301	start //= self.binning[binning_indices[axis]]	1✔
302	start += getattr(self.data_section, f'{axis}_start')	1✔
303
304	stop = sign * (getattr(section, f'{axis}_stop') - getattr(self.detector_section, f'{axis}_start'))	1✔
305	stop //= self.binning[binning_indices[axis]]	1✔
306	stop += getattr(self.data_section, f'{axis}_start')	1✔
307	return start, stop	1✔
308
309	def detector_to_data_section(self, section):	1✔
310	"""Given a detector region, figure out the corresponding data region
311	Note the + and - 1 factors cancel
312	Really this is just doing the same thing as a CD matrix calculation
313	r_data - data_0 = M (r_det - det_0) where M is a transformation matrix and everything else is a 2d vector
314	M is either the 1/binning * identity or the 1/binning * negative Identity depending on if datasec and detsec
315	are in the same ordering (both increasing or both decreasing gives positive)
316	Note if you want to add a rotation, M can contain that as well.
317	"""
318	x_start, x_stop = self.detector_to_data_section_oned('x', section)	1✔
319	y_start, y_stop = self.detector_to_data_section_oned('y', section)	1✔
320
321	return Section(x_start, x_stop, y_start, y_stop)	1✔
322
323	def data_to_detector_section(self, section):	1✔
324	"""Given a data region, get the detector section that this covers.
325	This is the inverse of get_data section"""
326	x_start, x_stop = self.data_to_detector_section_oned('x', section)	1✔
327	y_start, y_stop = self.data_to_detector_section_oned('y', section)	1✔
328
329	return Section(x_start, x_stop, y_start, y_stop)	1✔
330
331	def data_to_detector_section_oned(self, axis, section):	1✔
332	binning_indices = {'x': 0, 'y': 1}	1✔
333	sign = np.sign(getattr(self.detector_section, f'{axis}_stop') - getattr(self.detector_section, f'{axis}_start'))	1✔
334	sign *= np.sign(getattr(self.data_section, f'{axis}_stop') - getattr(self.data_section, f'{axis}_start'))	1✔
335
336	start = sign * (getattr(section, f'{axis}_start') - getattr(self.data_section, f'{axis}_start'))	1✔
337	start *= self.binning[binning_indices[axis]]	1✔
338	start += getattr(self.detector_section, f'{axis}_start')	1✔
339
340	stop = sign * (getattr(section, f'{axis}_stop') - getattr(self.data_section, f'{axis}_start'))	1✔
341	stop *= self.binning[binning_indices[axis]]	1✔
342	stop += getattr(self.detector_section, f'{axis}_start')	1✔
343
344	# when converting from binned to unbinned coordinates, we need to adjust our stopping pixel
345	# by (binning - 1)
346	stop += np.sign(stop - start) * (self.binning[binning_indices[axis]] - 1)	1✔
347
348	return start, stop	1✔
349
350	def copy_in(self, data):	1✔
351	"""
352	Copy in the data from another CCDData object based on the detector sections
353
354	:param data_to_copy:
355	:return:
356	"""
357	overlap_section = self.get_overlap(data.detector_section)	1✔
358	data_to_copy = data.trim(trim_section=data.detector_to_data_section(overlap_section))	1✔
359	data_to_copy = data_to_copy.rebin(self.binning)	1✔
360	for array_name_to_copy in ['data', 'mask', 'uncertainty']:	1✔
361	array_to_copy = getattr(data_to_copy, array_name_to_copy)	1✔
362	my_overlap = self.detector_to_data_section(overlap_section).to_slice()	1✔
363	getattr(self, array_name_to_copy)[my_overlap][:] = array_to_copy[:]	1✔
364
365	def init_poisson_uncertainties(self):	1✔
366	self.uncertainty = np.sqrt(self.uncertainty ** 2.0 + np.abs(self.data))	1✔
367
368	@property	1✔
369	def background(self):
370	return self._background	×
371
372	@background.setter	1✔
373	def background(self, value):
374	if self._background is not None:	×
375	self.data += self._background	×
376	self._background = value	×
377	self.data -= self._background	×
378
379
380	def stack(data_to_stack, nsigma_reject) -> CCDData:	1✔
381	"""
382	"""
383	shape3d = [len(data_to_stack)] + list(data_to_stack[0].shape)	1✔
384	a = np.zeros(shape3d, dtype=data_to_stack[0].dtype)	1✔
385	uncertainties = np.zeros(shape3d, dtype=data_to_stack[0].dtype)	1✔
386	mask = np.zeros(shape3d, dtype=np.uint8)	1✔
387
388	for i, data in enumerate(data_to_stack):	1✔
389	a[i, :, :] = data.data[:, :]	1✔
390	mask[i, :, :] = data.mask[:, :]	1✔
391	uncertainties[i, :, :] = data.uncertainty[:, :]	1✔
392
393	abs_deviation = stats.absolute_deviation(a, axis=0, mask=mask)	1✔
394
395	robust_std = stats.robust_standard_deviation(a, axis=0, abs_deviation=abs_deviation, mask=mask)	1✔
396
397	robust_std = np.expand_dims(robust_std, axis=0)	1✔
398
399	# Mask values that are N sigma from the median
400	sigma_mask = abs_deviation > (nsigma_reject * robust_std)	1✔
401
402	mask3d = np.logical_or(sigma_mask, mask > 0)	1✔
403	n_good_pixels = np.logical_not(mask3d).sum(axis=0)	1✔
404
405	stacked_mask = np.zeros(n_good_pixels.shape, dtype=np.uint8)	1✔
406
407	# If a pixel is bad in all images, make sure we don't divide by zero
408	bad_pixels = n_good_pixels == 0	1✔
409
410	# If pixel is bad in all of the images, we take the logical or of all of the bits to go in the final mask
411	stacked_mask[bad_pixels] = np.bitwise_or.reduce(mask, axis=0)[bad_pixels]	1✔
412
413	# If a pixel is bad in all images, fill that pixel with the mean from the images
414	n_good_pixels[bad_pixels] = len(data_to_stack)	1✔
415	mask3d[:, bad_pixels] = False	1✔
416
417	a[mask3d] = 0.0	1✔
418	stacked_data = a.sum(axis=0) / n_good_pixels	1✔
419
420	# Again if a pixel is bad in all images, fill the uncertainties with the quadrature sum / N images
421	uncertainties[mask3d] = 0.0	1✔
422	uncertainties *= uncertainties	1✔
423	stacked_uncertainty = np.sqrt(uncertainties.sum(axis=0) / (n_good_pixels ** 2.0))	1✔
424
425	return CCDData(data=stacked_data, meta=data_to_stack[0].meta, uncertainty=stacked_uncertainty, mask=stacked_mask)	1✔

LCOGT / banzai / 1664

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous