imastack.py 12.8 KB
# -*- coding: utf-8 -*-
"""
@author: aklotz@irap.omp.eu
"""

import os
import inspect
import datetime
import numpy as np

try:
    from .filenames import FileNames
except:
    from filenames import FileNames

try:
    from .ima import Ima
except:
    from ima import Ima

try:
    from .dates import Date
except:
    from dates import Date

try:
    from .guitastrotools import GuitastroTools, GuitastroException
except:
    from guitastrotools import GuitastroTools, GuitastroException

# #####################################################################
# #####################################################################
# #####################################################################
# Class Buffer
# #####################################################################
# #####################################################################
# This class provides an Astropy wrapper
# #####################################################################

class ImaStackException(GuitastroException):
    """Exception raised for errors in the ImaStack class.
    """

    SHAPES_NOT_SAME = 0
    errors = [""]*1
    errors[SHAPES_NOT_SAME] = "Shapes of images are not the same"



class ImaStack(ImaStackException, GuitastroTools):
    """Image processing applied to a series of images.

    This class process a series of images to generate only one image.

    The class ImaSeries imports the following methods from the class FileNames. See the FileNames class documentation for these methods:

        * longitude
        * namings
        * naming
        * naming_rules
        * naming_ident
        * naming_get
        * naming_set
        * path
        * extension
        * get_night
        * fullfilename
        * basename
        * itername
        * enumname
        * genename
        * genenames
        * innames
        * inoutnames
        * outfilename
        * indexnames
        * deletenames

    """

    VERBOSE_NONE = 0
    VERBOSE_SPECIFIC = 1
    VERBOSE_DEBUG = 2
    VERBOSE_ESSENTIAL = 4
    VERBOSE_ALL = 8

    def __init__(self, *args, **kwargs):
        self._ima =  Ima()
        self._verbose_level = self.VERBOSE_NONE
        self._longiau_deg = 0.0
        self._outfilename = os.path.join(self._ima.path(), "noname"+self._ima.extension())
        # --- import FileNames methods
        fn = FileNames()
        self.longitude = fn.longitude
        self.namings = fn.namings
        self.naming = fn.naming
        self.naming_rules = fn.naming_rules
        self.naming_ident = fn.naming_ident
        self.naming_get = fn.naming_get
        self.naming_set = fn.naming_set
        self.path = fn.path
        self.extension = fn.extension
        self.get_night = fn.get_night
        self.fullfilename = fn.fullfilename
        self.basename = fn.basename
        self.itername = fn.itername
        self.enumname = fn.enumname
        self.genename = fn.genename
        self.genenames = fn.genenames
        self.innames = fn.innames
        self.inoutnames = fn.inoutnames
        self.outfilename = fn.outfilename
        self.indexnames = fn.indexnames
        self.deletenames = fn.deletenames
        self._askopenfilename = fn._askopenfilename
        self._asksaveasfilename = fn._asksaveasfilename

    def longitude(self, longiau_deg:float):
        self._ima.longitude(longiau_deg)

    # =============================================
    # Print filters
    # =============================================

    def print_level(self, level:int=""):
        if isinstance(level,type(int)) == True:
            self._verbose_level = level
        return self._verbose_level

    def print_msg(self, level:int, message:str):
        if message == "":
            self._verbose_level = level
            return
        else:
            if self._verbose_level >= level:
                print(message)

    # =============================================
    # Managing filenames
    # =============================================

    def path(self, path=""):
        return self._ima.path(path)

    def extension(self, extension=""):
        return self._ima.extension(extension)

    def outshortfilename(self, fitsname:str=""):
        outfilename = self.outfilename(fitsname)
        basename = os.path.basename(outfilename)
        shortname, file_extension = os.path.splitext(basename)
        return shortname

    def outfilename(self, fitsname:str=""):
        if fitsname == "":
            # --- return the current output name
            basename = os.path.basename(self._outfilename)
            shortname, file_extension = os.path.splitext(basename)
            res = self._outfilename
        else:
            # --- return a simulation of a output name
            genename = self._ima.genename(fitsname)
            outfilename_short = genename["genename"] + genename["file_extension"]
            outfilename = os.path.join( genename["dirname"], outfilename_short)
            res = outfilename
        return res

    def fullfilename(self,fitsname:str)->str:
        return self._ima.path(fitsname)

    def basename(self,fitsname:str)->str:
        return self._ima.basename(fitsname)

    def genename(self,fitsname:str)->str:
        return self._ima.genename(fitsname)

    def innames(self, fitsname, **kwargs):
        # ---
        genename = self._ima.genename(fitsname)
        n = len(genename["indexes"])
        # ---
        first = 1
        if "first" in kwargs.keys():
            first = kwargs["first"]
            if first > n:
                first = n
        end = n
        if "end" in kwargs.keys():
            end = kwargs["end"]
            if end < first:
                end = first
            if end > n:
                end = n
        indexes = []
        for index in genename["indexes"]:
            rank = int(index)
            if rank >= first and rank <= end:
                indexes.append(index)
        genename["indexes"] = indexes
        # ---
        return genename

    def indexnames(self, genename, index, k):
        n = len(genename["indexes"])
        filename_short = genename["genename"] + genename["sep"] + index + genename["suffix"] + genename["file_extension"]
        filename = os.path.join( genename["dirname"], filename_short)
        called_by = inspect.stack()[1][3]
        print(f"ImaStack.{called_by} {k}/{n} : {filename_short}")
        return filename_short, filename

    # =============================================
    # Image processing (pixels are modified)
    # =============================================

    def median(self, fitsname, **kwargs):
        genename = self.innames(fitsname, **kwargs)
        n = len(genename["indexes"])
        k = 1
        for index in genename["indexes"]:
            filename_short, filename = self.indexnames(genename, index, k)
            self._ima.load(filename)
            naxis1, naxis2 = imastack._ima._array.shape
            if k == 1:
                naxis10, naxis20 = naxis1, naxis2
                self.allimages = np.ndarray(shape=(n, naxis10, naxis20), dtype=float)
            else:
                if naxis1 != naxis10 or naxis2 != naxis20:
                    raise ImaStackException(ImaStackException.SHAPES_NOT_SAME)
            self.allimages[k-1] = self._ima._array
            k += 1
        array = np.median(self.allimages, axis=0)
        self._ima._array = array.T
        if "outfile" in kwargs.keys():
            fitsname = kwargs["outfile"]
        self._outfilename = self.outfilename(fitsname)
        self._ima.save(self._outfilename)
        return self._outfilename

    def add(self, fitsname, **kwargs):
        genename = self.innames(fitsname, **kwargs)
        k = 1
        for index in genename["indexes"]:
            filename_short, filename = self.indexnames(genename, index, k)
            if k == 1:
                self._ima.load(filename)
            else:
                self._ima.add(filename)
            k += 1
        if "outfile" in kwargs.keys():
            fitsname = kwargs["outfile"]
        self._outfilename = self.outfilename(fitsname)
        self._ima.save(self._outfilename)
        return self._outfilename

    def mean(self, fitsname, **kwargs):
        genename = self.innames(fitsname, **kwargs)
        n = len(genename["indexes"])
        self.add(fitsname, **kwargs)
        self._ima.mult(1./n)
        self._ima.save(self._outfilename)
        return self._outfilename

    def sextractor(self, fitsname, **kwargs):
        genename = self.innames(fitsname, **kwargs)
        k = 1
        allstars = {}
        for index in genename["indexes"]:
            filename_short, filename = self.indexnames(genename, index, k)
            self._ima.load(filename)
            date_obs = self._ima.getkwd("DATE-OBS")
            if date_obs=="":
                jd = k
            else:
                jd = Date(date_obs).jd()
            stars = self._ima.sextractor("app_pho")
            if k == 1:
                star1s = stars
                for k1 in range(len(star1s)):
                    star1 = star1s[k1]
                    x1, y1, f1, flux1, fluxerr1, flag1 = star1
                    try:
                        ra, dec = self._ima.xy2radec(x1, y1)
                    except:
                        ra, dec = (0, 0)
                    allstars[str(k1)] = [(k1, k, jd, x1, y1, flux1, fluxerr1, flag1, ra, dec)]
            else:
                dr2min = 2
                for k1 in range(len(star1s)):
                    star1 = star1s[k1]
                    x1, y1, f1, flux1, fluxerr1, flag1 = star1
                    for star in stars:
                        x, y, f, flux, fluxerr, flag = star
                        dx = x-x1
                        dy = y-y1
                        dr2 = dx*dx+dy*dy
                        if dr2<dr2min:
                            try:
                                ra, dec = self._ima.xy2radec(x, y)
                            except:
                                ra, dec = (0, 0)
                            allstars[str(k1)].append( (k1, k, jd, x, y, flux, fluxerr, flag, ra, dec) )
                            break
            k += 1
        if "outfile" in kwargs.keys():
            fitsname = kwargs["outfile"]
        self._outfilename = self.outfilename(fitsname)
        names = os.path.splitext(self._outfilename)
        self._outfilename = names[0] + ".txt"
        # --- Fichier stars de sortie
        today = datetime.datetime.utcnow()
        dateiso = today.isoformat(sep='T',timespec='milliseconds')
        g = "guitastro"
        with open(self._outfilename, "wt") as fid:
            fid.write(f"# {g} method_name Imaseries.sextractor\n")
            fid.write(f"# {g} method_version 1.0\n")
            fid.write(f"# {g} process_date {dateiso}\n")
            fid.write(f"# {g} col  1 STAR_INDEX      Star index\n")
            fid.write(f"# {g} col  2 IMAGE_INDEX     Image index\n")
            fid.write(f"# {g} col  3 JD_OBS_UTC      Julian day UTC\n")
            fid.write(f"# {g} col  4 POS_X_PIX       x position in the first image (pix)\n")
            fid.write(f"# {g} col  5 POS_Y_PIX       y position in the first image (pix)\n")
            fid.write(f"# {g} col  6 FLUX_ADU        flux (adu)\n")
            fid.write(f"# {g} col  7 FLUX_ERROR_ADU  flux error (adu)\n")
            fid.write(f"# {g} col  8 SEXTRACTOR_FLAG sextractor flag\n")
            fid.write(f"# {g} col  9 RAJ2000         R.A. J2000.0\n")
            fid.write(f"# {g} col 10 DECJ2000        Dec. J2000.0\n")
            for k1 in range(len(star1s)):
                stars = allstars[str(k1)]
                for star in stars:
                    kk, k, jd, x, y, flux, fluxerr, flag, ra, dec = star
                    msg = f"{kk:5d} {k:4d} {jd:.6f} {x:7.2f} {y:7.2f} {flux:e} {fluxerr:e} {flag:2.0f} {ra:10.6f} {dec:+10.6f}"
                    print(msg)
                    fid.write(msg+"\n")
        # ---
        return self._outfilename

# #####################################################################
# #####################################################################
# #####################################################################
# Main
# #####################################################################
# #####################################################################
# #####################################################################

if __name__ == "__main__":

    example = 2
    print("Example       = {}".format(example))

    # --- Get path_images for examples
    ima = Ima()
    path_images = ima.conf_guitastro()['path_data']
    path_products = ima.conf_guitastro()['path_products']

    if example == 1:
        imastack = ImaStack()
        path = r"C:\d\mp-g2a\photometrie_pulsante\2017-11-07"
        f = os.path.join(path,"Dypeg offset-001.fit")
        res = imastack.mean(f)

    if example == 2:
        imastack = ImaStack()
        path = r"C:\d\mp-g2a\photometrie_pulsante\2017-11-07"
        imastack.path(path)
        res = imastack.median("bb")
        imastack._ima.visu()