diff --git a/.gitignore b/.gitignore index d4fadb35..9f700246 100644 --- a/.gitignore +++ b/.gitignore @@ -65,3 +65,8 @@ _docs docs/tested_blocks.md docs/md/all_blocks.rst docs/md/template.py + +## Pycharm +*.xml +.idea + diff --git a/docs/md/utils_blocks.rst b/docs/md/utils_blocks.rst index 1b9d4b4f..c1c7e897 100644 --- a/docs/md/utils_blocks.rst +++ b/docs/md/utils_blocks.rst @@ -14,4 +14,13 @@ Utils Get Calibration SortSources - CleanBadPixels \ No newline at end of file + CleanBadPixels + +.. currentmodule:: prose.blocks.visualization + +.. autosummary:: + :template: blocksum.rst + :nosignatures: + + Video + VideoPlot \ No newline at end of file diff --git a/prose/blocks/detection.py b/prose/blocks/detection.py index 6e3b2255..f9ab2861 100644 --- a/prose/blocks/detection.py +++ b/prose/blocks/detection.py @@ -8,7 +8,13 @@ from prose.console_utils import info from prose.core.source import * -__all__ = ["DAOFindStars", "SEDetection", "AutoSourceDetection", "PointSourceDetection"] +__all__ = [ + "DAOFindStars", + "SEDetection", + "AutoSourceDetection", + "PointSourceDetection", + "Peaks", +] class _SourceDetection(Block): @@ -268,17 +274,24 @@ def __init__( # TODO: document class Peaks(Block): - def __init__(self, cutout=11, **kwargs): + def __init__(self, shape=11, **kwargs): + """Computation of peak values for the detected stars (in ADUs). + + Parameters + ---------- + shape : int, optional + size of the cutout image within which the peak is calculated, by default 11 + """ super().__init__(**kwargs) - self.cutout = cutout + self.shape = shape def run(self, image, **kwargs): - idxs, cuts = cutouts(image.data, image.sources.coords, size=self.cutout) - peaks = np.ones(len(image.stars_coords)) * -1 - for j, i in enumerate(idxs): - cut = cuts[j] + idxs = np.arange(len(image.sources)) + peaks = np.ones(len(idxs)) * -1 + for j in idxs: + cut = image.cutout(image.sources.coords[j], shape=self.shape) if cut is not None: - peaks[i] = np.max(cut.data) + peaks[j] = cut.data.max() image.peaks = peaks diff --git a/prose/blocks/photometry.py b/prose/blocks/photometry.py index 75dfc05a..d16d8d72 100644 --- a/prose/blocks/photometry.py +++ b/prose/blocks/photometry.py @@ -107,6 +107,7 @@ def run(self, image: Image): annulus = image.sources.annulus(rin, rout) annulus_masks = annulus.to_mask(method="center") + annulus_area = np.pi * (rout**2 - rin**2) bkg_median = [] for mask in annulus_masks: @@ -125,4 +126,5 @@ def run(self, image: Image): "rout": rin, "median": np.array(bkg_median), "sigma": self.sigma, + "area": annulus_area, } diff --git a/prose/blocks/utils.py b/prose/blocks/utils.py index e4ec63c2..94108533 100644 --- a/prose/blocks/utils.py +++ b/prose/blocks/utils.py @@ -540,12 +540,24 @@ def get_time(im): def get_aperture(im): return im.aperture["radii"] + def get_error(im): + _area = np.pi * im.aperture["radii"] ** 2 + _signal = im.aperture["fluxes"] - ( + im.annulus["median"][:, None] * _area[None, :] + ) + # TODO : figure out the correct CCD equation for error computation + _squarred_error = _signal + _area[None, :] * ( + im.read_noise**2 + (im.gain / 2) ** 2 + im.annulus["median"][:, None] + ) + return np.sqrt(_squarred_error) + super().__init__( *args, _time=get_time, _bkg=get_bkg, _fluxes=get_fluxes, _apertures=get_aperture, + _errors=get_error, name=name, **kwargs, ) @@ -564,7 +576,11 @@ def terminate(self): data = {"bkg": np.mean(self._bkg, -1)} data.update({key: value for key, value in self.values.items() if key[0] != "_"}) self.fluxes = Fluxes( - time=time, fluxes=raw_fluxes, data=data, apertures=self._apertures + time=time, + fluxes=raw_fluxes, + data=data, + apertures=self._apertures, + errors=self._errors.T, ) diff --git a/prose/blocks/visualization.py b/prose/blocks/visualization.py index bd07b7b4..bfc110b8 100644 --- a/prose/blocks/visualization.py +++ b/prose/blocks/visualization.py @@ -1,145 +1,163 @@ import io -import shutil -import tempfile -import time import imageio import matplotlib.pyplot as plt import numpy as np -from matplotlib.backends.backend_agg import FigureCanvasAgg -from skimage.transform import resize - -from prose import Block, viz -from prose.visualization import corner_text - -__all__ = ["VideoPlot"] - - -def im_to_255(image, factor=0.25): - if factor != 1: - return ( - resize( - image.astype(float), - (np.array(np.shape(image)) * factor).astype(int), - anti_aliasing=False, - ) - * 255 - ).astype("uint8") - else: - data = image.copy().astype(float) - data = data / np.max(data) - data = data * 255 - return data.astype("uint8") - - -class _Video(Block): - """Base block to build a video""" - - def __init__(self, destination, fps=10, **kwargs): - super().__init__(**kwargs) - self.destination = destination - self.images = [] - self.fps = fps - self.checked_writer = False - def run(self, image): - if not self.checked_writer: - _ = imageio.get_writer(self.destination, mode="I") - self.checked_writer = True +from prose.core.block import Block +from prose.utils import z_scale - def terminate(self): - imageio.mimsave(self.destination, self.images, fps=self.fps) +__all__ = ["VideoPlot", "Video"] - @property - def citations(self): - return super().citations + ["imageio"] +def im_to_255(image): + data = image.copy().astype(float) + data = data / np.max(data) + data = data * 255 + return data.astype("uint8") -class RawVideo(_Video): + +class Video(Block): def __init__( - self, destination, attribute="data", fps=10, function=None, scale=1, **kwargs + self, + destination, + fps=10, + compression=None, + data_function=None, + width=None, + contrast=0.1, + name=None, ): - super().__init__(destination, fps=fps, **kwargs) - if function is None: + """ + A block to create a video from images data (using ffmpeg). - def _donothing(data): - return data + Parameters + ---------- + destination : str + The path to save the resulting video. + fps : int, optional + The frames per second of the resulting video. Default is 10. + compression : int, optional + The compression rate of the resulting video (-cr value of fmmpeg). + Default is None. + data_function : callable, optional + A function to apply to each image data before adding it to the video. + If none, a z scale is applied to the data with a contrast given by + :code:`contrast`.Default is None. + width : int, optional + The width in pixels of the resulting video. + Default is None (i.e. original image size). + contrast : float, optional + The contrast of the resulting video. Default is 0.1. + Either :code:`contrast` or :code:`data_function` must be provided. + name : str, optional + The name of the block. Default is None. + + Attributes + ---------- + citations : list of str + The citations for the block. - function = _donothing + Methods + ------- + run(image) + Adds an image to the video. + terminate() + Closes the video writer. - self.function = function - self.scale = scale - self.attribute = attribute + """ + super().__init__(name=name) + if data_function is None: + + def data_function(data): + new_data = data.copy() + new_data = z_scale(new_data, c=contrast) + return new_data + + output = [] + if compression is not None: + output += ["-crf", f"{compression}"] + if width is not None: + output += ["-vf", f"scale={width}:-1"] + self.writer = imageio.get_writer( + destination, + mode="I", + fps=fps, + output_params=output if len(output) > 0 else None, + ) + self.function = data_function def run(self, image): - super().run(image) - data = self.function(image.__dict__[self.attribute]) - self.images.append(im_to_255(data, factor=self.scale)) + data = self.function(image.data) + self.writer.append_data(im_to_255(data)) + def terminate(self): + self.writer.close() + + @property + def citations(self): + return super().citations + ["imageio"] -class VideoPlot(_Video): - def __init__(self, plot_function, destination, fps=10, name=None): - """Make a video out of a plotting function + +class VideoPlot(Video): + def __init__( + self, + plot_function, + destination, + fps=10, + compression=None, + width=None, + name=None, + ): + """ + A block to create a video from a matploltib plot (using ffmpeg). Parameters ---------- - plot_function : function - a plotting function taking an :py:class:`prose.Image` as input - destination : str or Path - destination of the video, including extension + plot_function : callable + A function that takes an image as input and produce a plot. + destination : str + The path to save the resulting video. fps : int, optional - frame per seconds, by default 10 - antialias : bool, optional - whether pyplot antialias should be used, by default False + The frames per second of the resulting video. Default is 10. + compression : int, optional + The compression rate of the resulting video (-cr value of fmmpeg). + Default is None. + width : int, optional + The width in pixels of the resulting video. + Default is None (i.e. original image size). + name : str, optional + The name of the block. Default is None. + + Attributes + ---------- + citations : list of str + The citations for the block. + + Methods + ------- + run(image) + Adds a plot to the video. + terminate() + Closes the video writer. + """ - super().__init__(destination, fps=fps, name=name) + super().__init__( + destination, + fps=fps, + compression=compression, + width=width, + name=name, + ) self.plot_function = plot_function - self.destination = destination - self._temp = tempfile.mkdtemp() - self._images = [] def run(self, image): self.plot_function(image) buf = io.BytesIO() plt.savefig(buf) - self.images.append(imageio.imread(buf)) + self.writer.append_data(imageio.imread(buf)) plt.close() - def terminate(self): - super().terminate() - shutil.rmtree(self._temp) - - -class LivePlot(Block): - def __init__(self, plot_function=None, sleep=0.0, size=None, **kwargs): - super().__init__(**kwargs) - if plot_function is None: - plot_function = lambda im: viz.show_stars( - im.data, - im.stars_coords if hasattr(im, "stars_coords") else None, - size=size, - ) - - self.plot_function = plot_function - self.sleep = sleep - self.display = None - self.size = size - self.figure_added = False - - def run(self, image): - if not self.figure_added: - from IPython import display as disp - - self.display = disp - if isinstance(self.size, tuple): - plt.figure(figsize=self.size) - self.figure_added = True - - self.plot_function(image) - self.display.clear_output(wait=True) - self.display.display(plt.gcf()) - time.sleep(self.sleep) - plt.cla() - def terminate(self): plt.close() + super().terminate() diff --git a/prose/core/image.py b/prose/core/image.py index 6e7c244d..edd14bc3 100644 --- a/prose/core/image.py +++ b/prose/core/image.py @@ -252,6 +252,16 @@ def filter(self): """Filter name""" return self.metadata["filter"] + @property + def gain(self): + """Gain of the camera""" + return self.metadata.get("gain", 1.0) + + @property + def read_noise(self): + """Read noise of the camera""" + return self.metadata.get("read_noise", 0.0) + @property def fov(self): """RA-DEC field of view of the image in degrees @@ -655,6 +665,11 @@ def label(self): ] ) + @property + def fits_header(self): + """Same as :code:`header` (backward compatibility)""" + return self.header + def str_to_astropy_unit(unit_string): return u.__dict__[unit_string] @@ -713,6 +728,8 @@ def FITSImage( "jd": header.get(telescope.keyword_jd, None), "object": header.get(telescope.keyword_object, None), "pixel_scale": telescope.pixel_scale, + "gain": telescope.gain, + "read_noise": telescope.read_noise, "overscan": telescope.trimming[::-1], "path": path, "dimensions": (header.get("NAXIS1", 1), header.get("NAXIS2", 1)), diff --git a/prose/fluxes.py b/prose/fluxes.py index 9798946e..ffdda40b 100644 --- a/prose/fluxes.py +++ b/prose/fluxes.py @@ -21,7 +21,7 @@ def weights( Parameters ---------- fluxes : np.ndarray - fluxes matrix with dimensions (star, flux) or (aperture, star, flux) + fluxes matrix with dimensions (star, flux) tolerance : float, optional the minimum standard deviation of weights difference to attain (meaning weights are stable), by default 1e-3 max_iteration : int, optional @@ -64,7 +64,7 @@ def weight_function(fluxes): ) last_weights = weights - lcs = diff(dfluxes, weights=weights) + lcs, _ = diff(dfluxes, weights=weights) i += 1 @@ -73,7 +73,7 @@ def weight_function(fluxes): return weights[0] -def diff(fluxes: np.ndarray, weights: np.ndarray = None): +def diff(fluxes: np.ndarray, weights: np.ndarray = None, errors: np.ndarray = None): """Returns differential fluxes. If weights are specified, they are used to produce an artificial light curve by which all fluxes are differentiated (see Broeg 2005). @@ -85,6 +85,8 @@ def diff(fluxes: np.ndarray, weights: np.ndarray = None): fluxes matrix with dimensions (star, flux) or (aperture, star, flux) weights :np.ndarray, optional weights matrix with dimensions (star) or (aperture, star), by default None which simply returns normalized fluxes + errors: np.ndarray, optionnal + errors matrix with the same dimensions as fluxes, (star, flux) or (aperture, star, flux) Returns ------- @@ -101,15 +103,26 @@ def diff(fluxes: np.ndarray, weights: np.ndarray = None): weights @ sub[0], -1 ) diff_fluxes = diff_fluxes / artificial_light_curve - return diff_fluxes + if errors is not None: + diff_errors = errors / np.expand_dims(np.nanmean(fluxes, -1), -1) + weighted_errors = diff_errors**2 * np.expand_dims(weights, -1) ** 2 + squarred_art_error = (sub @ weighted_errors) / np.expand_dims( + weights**2 @ sub[0], -1 + ) + diff_errors = np.sqrt(diff_errors**2 + squarred_art_error) + else: + diff_errors = None + return diff_fluxes, diff_errors -def auto_diff_1d(fluxes, i=None): - dfluxes = fluxes / np.expand_dims(np.nanmean(fluxes, -1), -1) - w = weights(dfluxes) + +def auto_diff_1d(fluxes, i=None, errors=None): + w = weights(fluxes) + + # to retain minimal set of comparison stars if i is not None: idxs = np.argsort(w)[::-1] - white_noise = utils.binned_nanstd(dfluxes) + white_noise = utils.binned_nanstd(fluxes) last_white_noise = 1e10 def best_weights(j): @@ -120,7 +133,7 @@ def best_weights(j): for j in range(w.shape[-1]): _w = best_weights(j) - _df = diff(dfluxes, _w) + _df, _ = diff(fluxes, _w) _white_noise = np.take(white_noise(_df), i, axis=-1)[0] if not np.isfinite(_white_noise): continue @@ -131,19 +144,33 @@ def best_weights(j): w = best_weights(j - 1) - df = diff(dfluxes, w) + df, derr = diff(fluxes, w, errors) + df = df.reshape(fluxes.shape) + + if derr is not None: + derr = derr.reshape(errors.shape) - return df.reshape(fluxes.shape), w + return df, w, derr -def auto_diff(fluxes: np.array, i: int = None): - if fluxes.ndim == 3: - auto_diffs = [auto_diff_1d(f, i) for f in fluxes] - w = [a[1] for a in auto_diffs] - fluxes = np.array([a[0] for a in auto_diffs]) - return fluxes, np.array(w) +def auto_diff(fluxes: np.array, i: int = None, errors=None): + if errors is not None: + if fluxes.ndim == 3: + auto_diffs = [auto_diff_1d(f, i, e) for f, e in zip(fluxes, errors)] + w = [a[1] for a in auto_diffs] + fluxes = np.array([a[0] for a in auto_diffs]) + errors = np.array([a[2] for a in auto_diffs]) + return fluxes, np.array(w), errors + else: + return auto_diff_1d(fluxes, i, errors) else: - return auto_diff_1d(fluxes, i) + if fluxes.ndim == 3: + auto_diffs = [auto_diff_1d(f, i) for f in fluxes] + w = [a[1] for a in auto_diffs] + fluxes = np.array([a[0] for a in auto_diffs]) + return fluxes, np.array(w), None + else: + return auto_diff_1d(fluxes, i) def optimal_flux(diff_fluxes, method="stddiff", sigma=4): @@ -302,10 +329,13 @@ def diff(self, comps: np.ndarray = None): else: weights = None - diff_fluxes = diff(self.fluxes, weights) _new = deepcopy(self) - _new.fluxes = diff_fluxes _new.weights = weights + + diff_fluxes, diff_errors = diff(self.fluxes, weights, self.errors) + _new.fluxes = diff_fluxes + _new.errors = diff_errors + return _new def autodiff(self): @@ -318,10 +348,13 @@ def autodiff(self): """ with warnings.catch_warnings(): warnings.simplefilter("ignore", RuntimeWarning) - diff_fluxes, weights = auto_diff(self.fluxes, self.target) + diff_fluxes, weights, diff_errors = auto_diff( + self.fluxes, self.target, self.errors + ) _new = deepcopy(self) _new.fluxes = diff_fluxes + _new.errors = diff_errors _new.weights = weights _new.aperture = _new.best_aperture_index() diff --git a/prose/io/fitsmanager.py b/prose/io/fitsmanager.py index dfbece2b..0c1cb19e 100644 --- a/prose/io/fitsmanager.py +++ b/prose/io/fitsmanager.py @@ -1,4 +1,5 @@ import sqlite3 +import sys from functools import partial from pathlib import Path @@ -242,6 +243,7 @@ def scan_files( telescope=None, verbose_new=False, verbose_os=False, + leave=False, ): """Scan files and add data to database @@ -260,6 +262,12 @@ def scan_files( FITS data unit extension where header will be parsed telescope: prose.Telescope telescope to be imposed for these files, by default None + verbose_new : bool, optional + whether to show how many files are new, by default False + verbose_os : bool, optional + whether to show OS commands, by default False + leave : bool, optional + whether to leave progress bar after completion, by default False """ if len(files) > 0: @@ -285,8 +293,14 @@ def scan_files( ) _verbose = verbose and batch_size is not False - _progress = progress(_verbose, desc="Reading fits", unit="files") + _progress = progress( + _verbose, + desc="Reading fits", + unit="files", + leave=leave, + ) + # to improve: use same progress for batch and non-batch if batch_size is not False: for batch in _progress(batches): try: diff --git a/prose/io/io.py b/prose/io/io.py index 0c149f0c..4840d6a3 100644 --- a/prose/io/io.py +++ b/prose/io/io.py @@ -91,7 +91,43 @@ def fits_to_df( hdu=0, raise_oserror=False, verbose_os=False, + leave=True, ): + """ + Convert FITS files to a pandas DataFrame. + + Parameters: + ---------- + files : list + List of FITS file paths. + telescope_kw : str, optional + Keyword for telescope name in FITS header. Default is "TELESCOP". + instrument_kw : str, optional + Keyword for instrument name in FITS header. Default is "INSTRUME". + telescope : Telescope, optional + Pre-defined Telescope object. Default is None. + verbose : bool, optional + Flag to display progress bar. Default is True. + hdu : int, optional + HDU index to read from FITS file. Default is 0. + raise_oserror : bool, optional + Flag to raise OSError if encountered. Default is False. + verbose_os : bool, optional + Flag to display OS error message. Default is False. + leave : bool, optional + Flag to leave progress bar after completion. Default is True. + + Returns: + ------- + df : pandas.DataFrame + DataFrame containing FITS file information. + + Raises: + ------ + AssertionError + If files list is empty. + + """ assert len(files) > 0, "Files not provided" last_telescope = "_" diff --git a/prose/scripts/fitsmanager.py b/prose/scripts/fitsmanager.py new file mode 100644 index 00000000..3a2ab2a0 --- /dev/null +++ b/prose/scripts/fitsmanager.py @@ -0,0 +1,96 @@ +import argparse +from pathlib import Path + +import pandas as pd + +from prose import FITSImage +from prose.io import FitsManager + + +def main(): + parser = argparse.ArgumentParser( + prog="prose FITS manager", description="Explore fits" + ) + subparsers = parser.add_subparsers(required=True) + + # parse + # ----- + def _parse(args): + fm = FitsManager(folders=args.folder, depth=args.depth, file=args.file) + observations = fm.observations() + print("Observations:") + print(observations) + + parse = subparsers.add_parser(name="parse", help="parse a fits folder") + parse.add_argument("folder", type=str, help="folder to explore", default=".") + parse.add_argument( + "-d", "--depth", type=int, help="depth of the search", default=10 + ) + parse.add_argument( + "-f", "--file", type=str, help="SQLite database file", default=None + ) + parse.add_argument( + "-m", "--max-rows", type=int, help="max number of rows to display", default=10 + ) + parse.set_defaults(func=_parse) + + # observations + # ------------ + def _observations(args): + fm = FitsManager(file=args.file) + observations = fm.observations( + telescope=args.telescope, + filter=args.filter, + date=args.date, + target=args.target, + ) + print("Observations:") + print(observations) + + observations = subparsers.add_parser( + name="observations", help="list observations from an SQLite database" + ) + observations.add_argument( + "file", type=str, help="SQLite database file", default=None + ) + observations.add_argument( + "-t", "--telescope", type=str, help="telescope name", default=None + ) + observations.add_argument( + "-f", "--filter", type=str, help="filter name", default=None + ) + observations.add_argument("-d", "--date", type=str, help="date", default=None) + observations.add_argument( + "-o", "--target", type=str, help="target name", default=None + ) + observations.add_argument( + "-m", "--max-rows", type=int, help="max number of rows to display", default=10 + ) + observations.set_defaults(func=_observations) + + # info + # ---- + def _info(args): + image = FITSImage(args.filename) + print(f"filename: {Path(args.filename).stem}") + print(f"telescope: {image.telescope.name}") + print(f"date: {image.date}") + print(f"target: {image.metadata['object']}") + print(f"filter: {image.filter}") + print(f"exposure: {image.exposure}") + print(f"dimensions: {image.shape}") + print(f"JD: {image.jd}") + print(f"RA: {image.ra}") + print(f"DEC: {image.dec}") + print(f"pixel scale: {image.pixel_scale}") + + info = subparsers.add_parser(name="info", help="FITS image information") + info.add_argument("filename", type=str, help="FITS image filename") + info.set_defaults(func=_info) + + args = parser.parse_args() + args.func(args) + + +if __name__ == "__main__": + main() diff --git a/pyproject.toml b/pyproject.toml index 139b953f..3132505e 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,6 +1,6 @@ [tool.poetry] name = "prose" -version = "3.2.10" +version = "3.3.1" description = "Modular image processing pipelines for Astronomy" authors = ["Lionel Garcia"] license = "MIT" @@ -52,3 +52,6 @@ ipywidgets = "*" [build-system] requires = ["poetry-core"] build-backend = "poetry.core.masonry.api" + +[tool.poetry.scripts] +fitsmanager = 'prose.scripts.fitsmanager:main' diff --git a/tests/test_blocks.py b/tests/test_blocks.py index 2d6e9ada..69bf16d1 100644 --- a/tests/test_blocks.py +++ b/tests/test_blocks.py @@ -101,6 +101,11 @@ def test_Get(): assert g.values == {"a": [3], "b": [6], "c": [42]} +def test_peaks(): + im = image.copy() + blocks.Peaks().run(im) + + def test_LimitSources(): from prose.core.source import PointSource, Sources @@ -240,3 +245,23 @@ def test_require_sources(): with pytest.raises(AttributeError, match="sources"): Sequence([blocks.Cutouts()]).run(im) + + +def test_Video(tmp_path): + from prose.blocks import Video + + im = image.copy() + im.sources = None + + Sequence([Video(tmp_path / "video.gif", fps=3)]).run([im, im, im]) + + +def test_VideoPlot(tmp_path): + from prose.blocks import VideoPlot + + def plot(image): + image.show() + + im = image.copy() + + Sequence([VideoPlot(plot, tmp_path / "video.gif", fps=3)]).run([im, im, im]) diff --git a/tests/test_fluxes.py b/tests/test_fluxes.py index 41de8a5d..43eeeccc 100644 --- a/tests/test_fluxes.py +++ b/tests/test_fluxes.py @@ -68,3 +68,10 @@ def test_3d(): f = Fluxes(fluxes=x, errors=x, target=0, aperture=0) f.error + + +def test_diff(): + x = np.random.uniform(0, 10000, size=(3, 2, 10)) + f = Fluxes(x) + f.target = 1 + diff = f.autodiff()