DataBlock.py

#
# A class to pack the content of a FITS file (IMAGE) and related data.
# --------------------------------------------------------------------
#
# Michel Caillat
# Observatoire de Paris - LERMA
# 27/06/2019
#
from astropy.io import fits
from astropy import wcs
from astropy.coordinates import SkyCoord
import numpy as np
import math
import json
from io import StringIO
from io import BytesIO
import matplotlib as mpl
import png
from PIL import Image, PngImagePlugin
import dask
import dask.array as da
import dask.delayed as dd

import traceback
from datetime import datetime, timedelta
import os
import re
from functools import reduce
from matplotlib import pyplot as pp

import yt
from spectral_cube import SpectralCube

import cv2

def cmap2palette (palette_name):
    x = pp.get_cmap(palette_name)
    return [tuple(map(lambda x: int(round(x*255)), x(i)[0:3])) for i in range(x.N)]

class DataBlock:
    # What is the FITS files root directory
    __FITSFilePrefix = ""

    # What is the PNG files root directory
    __PNGFilePrefix  = ""

    # What is the PNG files root directory
    __OBJFilePrefix  = ""

    # What are the palette names
    __palette_names = ["Greys", "RdYlBu", "hsv", "gist_ncar", "gist_rainbow", "gist_gray", "Spectral", "jet", "plasma", "inferno", "magma", "afmhot", "gist_heat"]

    __palette = {palette_name: cmap2palette(palette_name) for palette_name in __palette_names}

    __default_palette_name = "gist_rainbow"

    # What are the transformation names ?
    __transformation_names = ["minmax", "percent98"]
    __default_transformation_name = "minmax"

    # What are the video mode names ?
    __video_mode_names = ["direct", "inverse"]
    __default_video_mode_name = "direct"

    # What are the rendering capabilities ?
    __renderingCapabilities = {
        "itts": __transformation_names,
        "default_itt_index": 0,
        "luts": __palette_names,
        "default_lut_index": 2,
        "vmodes": __video_mode_names,
        "default_vmode_index": 0
    }

    # date time format
    __dateTimeFormat = "%m/%d/%Y - %H:%M:%S"

    # A regulare expression to catch the NAXISn
    __naxisn_rexp = re.compile("NAXIS[1-9]")

    # Maximum idleness duration for an instance of DataBlock
    # retrieved from the env var YAFITSS_MAXIDLE
    __max_idle = float(os.environ["YAFITSS_MAXIDLE"])

    #===========================================================================
    # Class methods - Setters and getters
    @classmethod
    def getPaletteFromName(cls, name):
        return DataBlock.__palette[name] if name in DataBlock.__palette_names else  DataBlock.__palette[DataBlock.__default_palette_name]

    @classmethod
    def getDefaultPaletteName(cls):
        return DataBlock.__default_palette_name

    @classmethod
    def getDefaults(cls):
        return (DataBlock.__default_palette_name, DataBlock.__default_transformation_name, DataBlock.__default_video_mode_name)

    @classmethod
    def setFITSFilePrefix(cls, FITSFilePrefix):
        DataBlock.__FITSFilePrefix = FITSFilePrefix

    @classmethod
    def setPNGFilePrefix(cls, PNGFilePrefix):
        DataBlock.__PNGFilePrefix = PNGFilePrefix

    @classmethod
    def setOBJFilePrefix(cls, OBJFilePrefix):
        DataBlock.__OBJFilePrefix = OBJFilePrefix

    @classmethod
    def getFITSFilePrefix(cls):
        return DataBlock.__FITSFilePrefix

    @classmethod
    def getPNGFilePrefix(cls):
        return DataBlock.__PNGFilePrefix

    @classmethod
    def getOBJFilePrefix(cls):
        return DataBlock.__OBJFilePrefix

    @classmethod
    def getRenderingCapabilities(cls):
        return DataBlock.__renderingCapabilities

    @classmethod
    def getDefaultITTName(cls):
        return DataBlock.__default_transformation_name

    @classmethod
    def getDefaultLUTName(cls):
        return DataBlock.__default_palette_name

    @classmethod
    def getDefaultVMName(cls):
        return DataBlock.__default_video_mode_name

    @classmethod
    def convert_size(cls, sizeInBytes):
        if sizeInBytes == 0:
            return "0B"
        size_name = ("B", "KB", "MB", "GB", "TB", "PB", "EB", "ZB", "YB")
        i = int(math.floor(math.log(sizeInBytes, 1024)))
        p = math.pow(1024, i)
        s = round(sizeInBytes / p, 2)
        return "%s %s" % (s, size_name[i])

    @classmethod
    def getDataBlockInfoNames(cls):
        return ["path", "creationTime", "lastAccess", "naxisn", "size", "idle"]

    @classmethod
    def getMaxIdle(cls):
        return DataBlock.__max_idle

    @classmethod
    def getDateTimeFormat(cls):
        return DataBlock.__dateTimeFormat

    #===========================================================================
    # CTOR
    #

    def __init__(self, logger):
        self.__logger = logger
        self.__relFITSFilePath = None
        self.__data = None      # The data part of the FITS file
        self.__header = None      # The FITS file header as a dictionary
        self.__sliceMargin = None
        self.__wcs2 = None      # WCS object
        self.__wcs3 = None      # WCS object
        self.__convert = None
        self.__cdelt = None
        self.__lastAccess = None    # the last time a method was called on self
        self.__creationTime = None      # the date of birth
        self.__sizeInBytes = None
        self.__naxisn = None
        self.__logger.debug("A DataBlock has been just built")

        self.__infod = {
            "path": lambda: self.__relFITSFilePath,
            "creationTime": lambda: self.__creationTime.strftime(DataBlock.getDateTimeFormat()),
            "lastAccess": lambda: self.__lastAccess.strftime(DataBlock.getDateTimeFormat()),
            "naxisn": lambda: self.__naxisn,
            "size": lambda: DataBlock.convert_size(self.__sizeInBytes),
            "idle": None
        }

        self.__statistics ={}

    #===========================================================================
    #
    # Properties
    #
    @property
    def relFITSFilePath(self):
        return self.__relFITSFilePath

    @property
    def header(self):
        return self.__header

    @property
    def lastAccess(self):
        return self.__lastAccess

    @property
    def creationTime(self):
        return self.__creationTime

    @property
    def sizeInBytes(self):
        return self.__sizeInBytes

    @property
    def naxisn(self):
        return self.__naxisn
#
#===========================================================================
#
# Public setters
#

    # Set a lastAccess on self
    def setLastAccess(self):
        self.__lastAccess = datetime.now()

    #
    # Given a FITS file path:
    # * read the content of the FITS file in memory
    # * populates the instance variables accordingly
    #
    def setData(self, relFITSFilePath):
        self.__logger.debug("setData: entering")
        self.__logger.debug( "Setting data: %s" % relFITSFilePath)
        absFITSFilePath = DataBlock.getFITSFilePrefix() + relFITSFilePath
        self.__logger.debug("Full path is '%s'" % absFITSFilePath)

        try:
            hdu_list = fits.open(absFITSFilePath)
            self.__logger.debug("Opened {absFITSFilePath}")
            # MUSE detecttion.
            data_index = "PRIMARY"
            primary_header = hdu_list[data_index].header
            if "EXTEND" in primary_header and \
                primary_header["EXTEND"] and \
                "INSTRUME" in primary_header and \
                primary_header["INSTRUME"].strip() == "MUSE" and \
                "DATA" in hdu_list :
                    data_index = "DATA"
                    self.__header = hdu_list[data_index].header
                    self.__header["ORIGIN"] = primary_header["ORIGIN"].strip()
                    self.__header["INSTRUME"] = primary_header["INSTRUME"].strip()
                    self.__header["TELESCOP"] = primary_header["TELESCOP"].strip()
                    self.__header["OBSERVER"] = primary_header["OBSERVER"].strip()
                    self.__header["CDELT1"] = self.__header["CD1_1"];
                    self.__header["CDELT2"] = self.__header["CD2_2"];
                    if "CD3_3" in self.__header:
                        self.__header["CDELT3"] = self.__header["CD3_3"];
                    self.__header["RADESYS"] = primary_header["RADECSYS"].strip();
                    self.__header["DATE-OBS"] = primary_header["DATE"].strip();
                    self.__header["DATE"] = primary_header["DATE"].strip();
                    self.__logger.debug("We are in MUSE")
            # This is NOT a MUSE observation
            else:
                self.__header = hdu_list[data_index].header

            numDims = len(hdu_list[data_index].data.shape)
            if numDims > 4 :
                self.__logger.debug("Unacceptable value of numDims: %s" % numDims)
            if numDims == 4:
                self.__data = da.from_array(hdu_list[data_index].data[0], chunks = (hdu_list[data_index].data[0].shape[0], 128, 128))
            elif numDims == 3:
                self.__data = da.from_array(hdu_list[data_index].data, chunks = (hdu_list[data_index].data.shape[0], 128, 128))
            elif numDims == 2:
                self.__data = da.from_array(hdu_list[data_index].data, chunks = (128, 128))
            else:
                self.__logger.debug("Unacceptable value of numDims: %s" % numDims)

            shp = self.__data.shape
            if numDims > 2 :
                if shp[1] > shp[2]:
                    self.__sliceMargin = da.from_array(np.nan * np.ones(shape=(shp[1], shp[1] - shp[2])), chunks=(128, 128))
                elif shp[1] < shp[2]:
                    self.__sliceMargin = da.from_array(np.nan * np.ones(shape=(shp[2] - shp[1], shp[2])), chunks=(128, 128))
                else:
                    pass
            else :
                if shp[0] > shp[1]:
                    self.__sliceMargin = da.from_array(np.nan * np.ones(shape=(shp[0], shp[0] - shp[1])), chunks=(128, 128))
                elif shp[0] < shp[1]:
                    self.__sliceMargin = da.from_array(np.nan * np.ones(shape=(shp[1] - shp[0], shp[1])), chunks=(128, 128))
                else:
                    pass


            #
            # Header "normalization"
            #
            if "BUNIT" in self.__header and self.__header["BUNIT"] == "JY/BEAM":
                self.__header["BUNIT"] = "Jy/beam"

            if "TELESCOP" not in self.__header:
                self.__header["TELESCOP"] = "Unknown"

            if "INSTRUME" not in self.__header:
                self.__header["INSTRUME"] = "Unknown"

            if "OBSERVER" not in self.__header:
                self.__header["OBSERVER"] = "John Doe"

            if "ORIGIN" not in self.__header:
                self.__header["ORIGIN"] = "Unknown"

            if "SPECSYS" not in self.__header:
                self.__header["SPECSYS"] = "Unknown"

            if "RESTFRQ" not in self.__header:
                if "RESTFREQ" in self.__header:
                    self.__header["RESTFRQ"] = self.__header["RESTFREQ"]
                else:
                    self.__header["RESTFRQ"] = 0.0

            if self.__header["ORIGIN"] == "GILDAS Consortium":
                if "CUNIT3" not in self.__header:
                    self.__header["CUNIT3"] = "km/s"
                if "RADESYS" not in self.__header:
                    self.__header["RADESYS"] = "ICRS"
                if "DATE-OBS" not in self.__header:
                    self.__header["DATE-OBS"] = self.__header["DATE"]
                if "TELESCOP" not in self.__header:
                    self.__header["TELESCOP"] = "NOEMA"

            if "BUNIT" not in self.__header:
                self.__header["BUNIT"] = "Undefined"
            elif self.__header["INSTRUME"] == "SITELLE" \
                and self.__header["BUNIT"] == "FLUX":
                self.__header["BUNIT"] = "erg/s/cm^2/A/arcsec^2"

            self.__logger.debug("Working on '%s'" % (relFITSFilePath))

            # *None* prevents throwing error
            self.__header.pop("HISTORY", None)
            self.__header.pop("COMMENT", None)
            self.__header.pop("", None)
            try:
                self.__logger.debug("Got RESTFRQ %f" % self.__header["RESTFRQ"])
            except KeyError:
                self.__header["RESTFRQ"] = self.__header["RESTFREQ"]

            self.__wcs2 = wcs.WCS(self.__header, naxis=2)
            if self.__header["INSTRUME"] != "SITELLE" and self.__header["NAXIS"] > 2:
                self.__wcs3 = wcs.WCS(self.__header, naxis=3)

            pi180 = math.pi / 180 / 4.86e-6
            try:
                bmaj = self.__header["BMAJ"] * pi180
                bmin = self.__header["BMIN"] * pi180
                self.__convert = math.pi * bmaj * bmin
                self.__cdelt =  4 * math.log(2) * math.fabs(self.__header["CDELT1"] * pi180) * math.fabs(self.__header["CDELT2"] * pi180)
            except KeyError:
                self.__convert = 1.0
                self.__cdelt = 1.0

            self.__relFITSFilePath = relFITSFilePath

            # date of birth == lastAccess == now
            self.__lastAccess = self.__creationTime = datetime.now()
            naxisnValues = [self.__header[key] for key in self.__header if DataBlock.__naxisn_rexp.match(key)]
            self.__sizeInBytes = reduce(lambda x, y: x*y, naxisnValues) * abs(self.__header["BITPIX"]) / 8
            self.__naxisn = reduce(lambda x, y: "%sx%s"%(x,y), naxisnValues)

            result = self.getHeader()
        except Exception as e:
            message = f"Error while opening file {relFITSFilePath}: {e}"
            self.__logger.debug(message)
            self.__logger.debug(f"Traceback is : {traceback.format_exc()}")
            result = {"status": False, "message": message, "result": None}

        self.__logger.debug("setData: Exiting")
        return result


    #===========================================================================
    #
    # Private methods
    #
    def __getPixelValueAtiFreqiRAiDEC(self, iFreq, iRA, iDEC):

        self.__logger.debug("__getPixelValuAtiFreqiRAiDEC: entering")
        self.__logger.debug("iFreq = {:d}, iRA = {:d}, iDEC = {:d}".format(iFreq, iRA, iDEC) )
        result = self.__data[iFreq, iDEC, iRA].compute()
        if np.isnan(result):
            result = None
        else:
            result = float(result)
        self.__logger.debug("__getPixelValuAtiFreqiRAiDEC: exiting")
        return result

    def __getFreqsAtiRAiDEC(self, iRA, iDEC):
        self.__logger.debug("__getFreqsAtiRAiDEC: entering")
        x = [[iRA, iDEC, iFreq] for iFreq in range(self.__header["NAXIS3"])]
        result = [crdnn[2] for crdnn in self.__wcs3.all_pix2world(x, 0)]
        self.__logger.debug("__getFreqsAtiRAiDEC: exiting")
        return result

    def __getSumOverSliceRectArea_0(self, iFREQ, iRA0=None, iRA1=None, iDEC0=None, iDEC1=None):
        self.__logger.debug("__getSumOverSliceRectArea_0: entering")

        if self.__header["NAXIS"] > 2:
            result = (dask.array.nansum(self.__data[iFREQ, iDEC0:iDEC1, iRA0:iRA1]) / self.__convert * self.__cdelt).compute().tolist()
        else :
            result = (dask.array.nansum(self.__data[iDEC0:iDEC1, iRA0:iRA1]) / self.__convert * self.__cdelt).compute().tolist()

        self.__logger.debug("__getSumOverSliceRectArea__0: exiting")

        return result

    def __getAverage_0(self, iFREQ0=None, iFREQ1=None, iDEC0=None, iDEC1=None, iRA0=None, iRA1=None):
        self.__logger.debug("__getAverage_0: entering")
        self.__logger.debug("shape length: %d" % len(self.__data.shape))

        if iRA0 == None:
            iRA0 = 0
        if iRA1 == None:
            iRA1 = self.__data.shape[2]
        if iDEC0 == None:
            iDEC0 = 0
        if iDEC1 == None:
            iDEC1 = self.__data.shape[1]
        if iFREQ0 == None:
            iFREQ0 = 0
        if iFREQ1 == None:
            iFREQ1 = self.__data.shape[0]

        self.__logger.debug(f"{iFREQ0}, {iFREQ1}, {iDEC0}, {iDEC1}, {iRA0}, {iRA1}")

        if (self.__header["ORIGIN"]=="GILDAS Consortium"):
            self.__logger.debug("Gildas data")
            cdelt = math.fabs(self.__header["CDELT3"])/1000.
            self.__logger.debug("CDELT = %f"%cdelt)
        elif self.__header["ORIGIN"].startswith("CASA"):
            self.__logger.debug("ALMA data")
            cunit3 = self.__header["CUNIT3"]
            if cunit3 == "m/s":
                cdelt = math.fabs(self.__header["CDELT3"]) / 1000.0
            elif cunit3 == "Hz":
                cdelt = 3e+5 * math.fabs(self.__header["CDELT3"] / self.__header["RESTFRQ"])
            else:
                cdelt = 0.0
        elif self.__header["INSTRUME"] == "SITELLE":
            cdelt = 1.0
        elif self.__header["INSTRUME"] == "MUSE":
            cdelt = 1.0
        else:
            cdelt = 0.0

        self.__logger.debug("slice sum begin")

        if self.__header["NAXIS"] > 2:
            result = (dask.array.nansum(self.__data[iFREQ0:iFREQ1, iDEC0:iDEC1, iRA0:iRA1], 0)*cdelt).compute()
        else:
            result = self.__data[iDEC0:iDEC1, iRA0:iRA1]*cdelt.compute()

        self.__collectStatistics("%d-%d"%(iFREQ0, iFREQ1), result)


        self.__logger.debug("slice sum end")
        self.__logger.debug(f"result={result}")
        self.__logger.debug("__getAverage_0: exiting")
        return result

    def __getPercentile(self, a, percent):
        return np.nanpercentile(a, percent)

    def __convertOneSlice2PNG(self, PNGPath, iFREQ, sliceData,  transformation_name, palette_name, video_mode_name):
        self.__logger.debug("__convertOneSlice2PNG: entering")
        self.__logger.debug("__convertOneSlice2PNG: exiting")
        return self.__convertOneSlice2PNG_0(PNGPath, iFREQ, sliceData, transformation_name, palette_name, video_mode_name)

    def __convertOneSlice2PNG_0 (self, PNGPath, sliceData, transformation_name, palette_name, video_mode_name):
        self.__logger.debug("__convertOneSlice2PNG_0: entering")

        img_f = np.zeros((sliceData.shape[1], sliceData.shape[0]), dtype=np.float32)
        img_i = np.zeros((sliceData.shape[1], sliceData.shape[0]), dtype=np.uint8)

        data_steps = {}
        palette = DataBlock.getPaletteFromName(palette_name)
        N = len(palette)
        if transformation_name == "percent98":
            self.__logger.debug("A '%s' transformation will be applied " % transformation_name)
            img_f = np.array(sliceData)
            p_high = self.__getPercentile(sliceData, 99.9)
            p_low = self.__getPercentile(sliceData, 0.1)

            data_min, data_max = self.__min_max(sliceData)
            shape = sliceData.shape

            n_low_steps = 3
            n_high_steps = 3
            n_mid_steps = N - n_low_steps - n_high_steps

            data_step = (p_low - data_min) / n_low_steps
            for i in range (0, n_low_steps):
                data_steps["%d_%d_%d"%(palette[i][0], palette[i][1], palette[i][2])] = data_min + i * data_step

            data_step = (p_high - p_low) / n_mid_steps
            for i in range(0, n_mid_steps):
                j = n_low_steps + i
                data_steps["%d_%d_%d"%(palette[j][0], palette[j][1], palette[j][2])] = p_low + i * data_step

            data_step = (data_max - p_high) / n_high_steps
            for i in range(0, n_high_steps):
                j = n_low_steps + n_mid_steps + i
                data_steps["%d_%d_%d"%(palette[j][0], palette[j][1], palette[j][2])] = p_high + i * data_step

            img_f_flat = img_f.flatten()
            img_f_flat_N = np.zeros(img_f_flat.shape)

            low_mask = img_f_flat < p_low
            img_f_flat_N[low_mask] = n_low_steps * (img_f_flat[low_mask] - data_min) / (p_low - data_min)

            mid_mask = np.logical_and((img_f_flat >= p_low), (img_f_flat < p_high))
            img_f_flat_N[mid_mask] = n_low_steps + n_mid_steps * (img_f_flat[mid_mask] - p_low) / ( p_high - p_low)

            high_mask = img_f_flat >= p_high
            img_f_flat_N[high_mask] = n_low_steps + n_mid_steps + n_high_steps * (img_f_flat[high_mask] - p_high) / (data_max - p_high) - 1

            img_f = img_f_flat_N.reshape(shape)

        else:
            self.__logger.debug("A '%s' transformation will be applied " % "minmax")
            data_min, data_max = self.__min_max(sliceData)

            data_step = (data_max - data_min) / (N - 1)
            for i in range(0, N):
                data_steps["%d_%d_%d"%(palette[i][0], palette[i][1], palette[i][2])] = data_min + i * data_step

            data_range = data_max - data_min
            img_f = ((N - 1) * (np.array(sliceData) - data_min) / data_range)

        img_i = img_f.astype(np.uint8)

        self.__logger.debug("video_mode_name = %s" % video_mode_name)
        if (video_mode_name == "inverse"):
            self.__logger.debug("Display in inverse mode ")
            img_i = (N - 1) - img_i
        else:
            self.__logger.debug("Display in direct mode")

        for irow in range(img_i.shape[1]//2):
            img_i[irow,:], img_i[img_i.shape[1] - 1 - irow,:] = img_i[img_i.shape[1] - 1 - irow,:], img_i[irow,:].copy()

        f = open(PNGPath, 'wb')
        w = png.Writer(sliceData.shape[1], sliceData.shape[0],palette=palette)
        w.write(f, img_i.tolist())
        f.close()

        self.__logger.debug("__convertOneSlice2PNG_0: exiting")
        return data_steps

    def __convertSummedSliceRange2PNG(self, iFREQ0, iFREQ1, sliceData, transformation_name, palette_name, video_mode_name):
        self.__logger.debug("__convertSummedSliceRange2PNG: entering")
        PNGPath = "%s/%d-%d.%s.%s.%s.png" % (self.__PNGDir, iFREQ0, iFREQ1, transformation_name, palette_name, video_mode_name)

        return self.__convertOneSlice2PNG_0(PNGPath, sliceData, transformation_name, palette_name, video_mode_name)
        self.__logger.debug("__convertSummedSliceRange2PNG: exiting")
        return

    def __min_max(self, dataArray):
        data_max = np.nanmax(dataArray)
        data_min = np.nanmin(dataArray)
        return (data_min.astype(float), data_max.astype(float))

    def __getForcedTransformationName(self,  transformation_name):
        self.__logger.debug("__getForcedTransformationName: entering")
        result = transformation_name
        self.__logger.debug("__getForcedTransformationName: exiting")
        return result

    def __squareSliceData (self, sliceData):
        sliceShape = sliceData.shape
        if sliceShape[0] == sliceShape[1]:
            squaredSliceData = sliceData
        else:
            tmparr = [sliceData, self.__sliceMargin]
            if sliceShape[0] > sliceShape[1]:
                axis = 1
            else:
                axis = 0

            squaredSliceData = da.concatenate(tmparr, axis=axis)
        return squaredSliceData

    def __addWCStoPNG(self, absPNGFilePath, headerInfos):
        self.__logger.debug("__addWCStoPNG: entering")
        im = Image.open(absPNGFilePath)
        info = PngImagePlugin.PngInfo()
        naxis2 = max(headerInfos["NAXIS1"], headerInfos["NAXIS2"])
        naxis1 = naxis2
        info.add_text("zTXtcomment",
            "SIMPLE  = T "+
            "\n BITPIX = "+str(headerInfos["BITPIX"])+
            "\n NAXIS = 2 "+
            "\n NAXIS1  = "+str(naxis1)+
            "\n NAXIS2  = "+str(naxis2)+
            "\n CRPIX1  = "+str(headerInfos["CRPIX1"])+
            "\n CRPIX2  = "+str(headerInfos["CRPIX2"])+
            "\n EQUINOX = 2000.0" +
            "\n CRVAL1  = "+str(headerInfos["CRVAL1"])+
            "\n CRVAL2  = "+str(headerInfos["CRVAL2"])+
            "\n CTYPE1  = "+str(headerInfos["CTYPE1"])+
            "\n CTYPE2  = "+str(headerInfos["CTYPE2"])+
            "\n RADESYS= "+str(headerInfos["RADESYS"])+
            "\n CD1_1   = "+str(headerInfos["CDELT1"])+
            "\n CD1_2   = -0.0 " +
            "\n CD2_1   = -0.0 " +
            "\n CD2_2 = " +str(headerInfos["CDELT2"]) )

        im.save(absPNGFilePath, "PNG", pnginfo=info)
        self.__logger.debug("__addWCStoPNG: exiting")

    def __getSpectrumAtiRAiDEC(self, iRA, iDEC):
        self.__logger.debug("__getSpectrumAtiRAiDEC: entering")
        result = self.__data[:, iDEC, iRA]
        self.__logger.debug("__getSpectrumAtiRAiDEC: exiting")
        return result

    def __createFITSSpectrumFromData(self, crds, crdsUnit, spectrum, spectrumUnit):
        self.__logger.debug("__createFITSSpectrumFromData: entering")
        header = self.__header
        #outmp = StringIO()
        outmp = BytesIO()
        #create table
        xcolumn = fits.Column(name = 'Frequency', format='E', unit = crdsUnit, array=crds)
        ycolumn = fits.Column(name = 'Flux', format='E', unit = spectrumUnit, array=spectrum)
        hdu = fits.TableHDU.from_columns(fits.ColDefs([xcolumn, ycolumn]))

        #create header and primary hdu
        hdr = fits.Header()
        hdr['DATE'] = header["DATE"]
        hdr['DATE-OBS'] = header["DATE-OBS"]
        hdr['TELESCOP'] = header["TELESCOP"]
        hdr['OBSERVER'] = header["OBSERVER"]
        hdr['RESTFRQ'] = header["RESTFRQ"]
        hdr['SPECSYS'] = header["SPECSYS"]
        hdr['COMMENT'] = 'Spectrum of a single pixel of the cube'
        primary_hdu = fits.PrimaryHDU(header = hdr)
        hdulist = fits.HDUList([primary_hdu,hdu])
        self.__logger.debug("type of hdlist = %s" % type(hdulist))
        hdulist.writeto(outmp,overwrite=True)
        fits_content = outmp.getvalue().decode("utf-8")
        outmp.close()
        self.__logger.debug("__createFITSSpectrumFromData: exiting")
        return fits_content

    def __createFITS0(self, iRA, iDEC):
        self.__logger.debug("__createFITS0: entering")
        result = {"status": False, "message": "", "result": None}
        try:
            header = self.__header
            #
            # Retrieve the data to be published via SAMP
            #
            # frequencies ( or velocities )
            crds = self.__getFreqsAtiRAiDEC(iRA, iDEC)
            spectrum = self.__getSpectrumAtiRAiDEC(iRA, iDEC)


            fits_content = self.__createFITSSpectrumFromData( crds, header["CUNIT3"], spectrum, header["BUNIT"])
            result = {"status": True, "message": "", "result": fits_content}
        except Exception as e:
            result["message"] = "An exception occurred with message '%s'" % str(e)
            self.__logger.debug("An exception occurred with message '%s'" % str(e))
            self.__logger.debug("Traceback: %s"  % traceback.format_exc())
        self.__logger.debug("__createFITS0: exiting")
        return result

    #===========================================================================
    # Public getters and accessors.
    #
    def getlastAccess(self):
        return {"status": True, "message": "", "result": self.__lastAccess}

    def getHeader(self):
        return {"status": True, "message": "", "result": json.dumps(dict(self.__header))}

    def getDimensions(self):
        return {"status": True, "message": "", "result": self.__data.shape}

    def __collectStatistics(self, key, data): # data is expected to be a numpy array
        self.__logger.debug("__collectStatistics : entering")
        if key in self.__statistics:
            self.__logger.debug("Statistics are already collected for key %s" % key)
        else:
            self.__logger.debug("Statistics have to be collected for key %s" % key)
            results = []
            t = dd(np.nanmin)(data); results.append(t)
            t = dd(np.nanmax)(data); results.append(t)
            t = dd(np.nanmean)(data); results.append(t)
            t = dd(np.nanstd)(data); results.append(t)
            t = dd(np.histogram)(data[~np.isnan(data)],bins=100 ); results.append(t)

            results = dask.compute(*results)
            self.__statistics[key]={}
            self.__statistics[key]["min"] = results[0].item()
            self.__statistics[key]["max"] = results[1].item()
            self.__statistics[key]["mean"] = results[2].item()
            self.__statistics[key]["std"] = results[3].item()
            self.__statistics[key]["histogram"] = [x.tolist() for x in results[4]]
            population = self.__statistics[key]["histogram"][0]
            bins = self.__statistics[key]["histogram"][1]
            aux = [0 for x in range(len(bins))]
            for i in range (1, len(aux)):
                aux[i] = aux[i-1]+population[i-1]
            normfactor = np.count_nonzero(~np.isnan(data)) #reduce(lambda a,b : a*b, self.__data.shape);
            self.__statistics[key]["cumuldist"] = list(map(lambda x: x/normfactor, aux));
        self.__logger.debug(self.__statistics[key]);
        self.__logger.debug("__collectStatistics : exiting")


    def __getSlice(self, iFREQ, step=1):
        self.__logger.debug("__getSlice: entering")

        result = None
        numDimensions = len(self.__data.shape)

        if numDimensions == 2:
            result = self.__data[0:self.__data.shape[0]:step, 0:self.__data.shape[1]:step]
        elif numDimensions == 3:
            if iFREQ == None:
                iFREQ = self.__data.shape[0] / 2
            result = self.__data[iFREQ, 0:self.__data.shape[1]:step, 0:self.__data.shape[2]:step]

        #
        # Collect some statistics if not already done
        #
        self.__collectStatistics("%d"%iFREQ, result.compute())
        self.__logger.debug("__getSlice: exiting")
        return result # A DASK array

    def getSlice(self, iFREQ, step=1):
        self.__logger.debug("getSlice: entering")
        result = dict()
        numDimensions = len(self.__data.shape)

        try :
            if numDimensions > 3 or numDimensions < 2:
                message = f"Can't process data with such a shape: {self.__data.shape}"
                result["status"]=False
                result["message"]= message
            else:
                x = self.__getSlice(iFREQ, step)
                self.__logger.debug(f"type self.__data = {type(self.__data)}, type x = {type(x)}")
                self.__logger.debug("Avant isnan")
                self.__logger.debug(x.shape)
                x = x.compute()
                b = np.isnan(x)
                self.__logger.debug(b.shape)
                x[b]=None
                self.__logger.debug("Apres isnan")
                result["status"]=True
                result["message"]=""
                result["result"]={}
                result["result"]["slice"]=x.tolist()
                result["result"]["statistics"]=self.__statistics["%d"%iFREQ]
        except Exception as e:
            result["message"]=f"An exception occurred with message {e}"
            result["status"]= False
            result["result"]=None

        self.__logger.debug("getSlice: exiting")
        return result

    def getSpectrum(self, iRA=None, iDEC=None, iFREQ0=None, iFREQ1=None):
        self.__logger.debug( "getSpectrum: entering")
        if iRA == None:
            iRA = 0
        if iDEC == None:
            iDEC = 0
        if iFREQ0 == None:
            iFREQ0 = 0
        if iFREQ1 == None:
            iFREQ1 = self.__data.shape[0]

        nonanArray = np.nan_to_num(self.__data[iFREQ0:iFREQ1, iDEC, iRA].compute())
        result = nonanArray.tolist()
        self.__logger.debug( "getSpectrum: exiting")

        return {"status": True, "message": "", "result": result}

    def getPixelValueAtiFreqiRAiDEC(self, iFreq, iRA, iDEC):
        self.__logger.debug("getPixelValueAtiFreqiRAiDEC: entering")
        try:
            data = self.__data
            shape=data.shape
            if (len(shape) == 3 ) and (0 <= iFreq < shape[0]) and ( 0 <= iRA < shape[1]) and (0 <= iDEC < shape[2]):
                result = {"status": True, "message": "", "result": self.__getPixelValueAtiFreqiRAiDEC(iFreq, iRA, iDEC)}
            else:
                result = {"status": False, "message": "Invalid coordinates '{0}, {1}, {2}' or data with inappropriate shape '{3}' ".format(iFreq, iRA, iDEC, shape), "result": None}
        except Exception as e:
            result = {"status": False, "message": "{0} - {1}".format(type(e), e.args), "result": None}

        self.__logger.debug("getPixelValueAtiFreqiRAiDEC: exiting")
        return result

    def getAverageSpectrum(self, iDEC0=None, iDEC1=None, iRA0=None, iRA1=None, retFITS=False):
        self.__logger.debug("getAverageSpectrum: entering")

        if iRA0 == None:
            iRA0 = 0
        if iRA1 == None:
            iRA1 = self.__data.shape[2]
        if iDEC0 == None:
            iDEC0 = 0
        if iDEC1 == None:
            iDEC1 = self.__data.shape[1]

        pi180 = math.pi / 180 / 4.86e-6

        averageSpectrum = None
        with_dask=True
        if with_dask:
            averageSpectrum = (dask.array.nansum(self.__data[:, iDEC0:iDEC1, iRA0:iRA1], (1,2)) / self.__convert * self.__cdelt).compute().tolist()
        else:
            averageSpectrum = np.nansum(self.__data[:, iDEC0:iDEC1, iRA0:iRA1], (1,2)) / self.__convert * self.__cdelt
        if retFITS:
            crds = self.__getFreqsAtiRAiDEC(iRA0, iDEC0)
            averageSpectrumFits = self.__createFITSSpectrumFromData(crds,self.__header["CUNIT3"], averageSpectrum, "Jy")
            result = {"averageSpectrum": averageSpectrum, "averageSpectrumFits": averageSpectrumFits}
        else:
            result = {"averageSpectrum": averageSpectrum, "averageSpectrumFits": None}

        result = {"status": True, "message": "", "result": result}
        self.__logger.debug("getAverageSpectrum: exiting")
        return result

    def getSumOverSliceRectArea(self, iFREQ, iRA0=None, iRA1=None, iDEC0=None, iDEC1=None):
        self.__logger.debug("getSumOnSliceRectArea: entering")
        if (len(self.__data.shape) == 2 and iFREQ >= 1) or (len(self.__data.shape) == 3 and iFREQ >= self.__data.shape[0]):
            result = {"status": False, "message": "Invalid slice index '%d'." % (iFREQ), "result": None }
        else:
            result = {"status": True, "message": "", "result": self.__getSumOverSliceRectArea_0(iFREQ, iRA0, iRA1, iDEC0, iDEC1)}

        self.__logger.debug("getSumOnSliceRectArea: exiting")
        return result

    def getAverage(self, iFREQ0=None, iFREQ1=None, iDEC0=None, iDEC1=None, iRA0=None, iRA1=None, retFITS=False):
        self.__logger.debug("getAverage: entering")
        result = self.__getAverage_0(iFREQ0, iFREQ1, iDEC0, iDEC1, iRA0, iRA1)
        self.__logger.debug("getAverage: exiting")
        return {"status": True, "message": "", "result": result.tolist()}

    def getOneSliceAsPNG (self, iFREQ,  ittName=__default_transformation_name, lutName=__default_palette_name, vmName=__default_video_mode_name):
        self.__logger.debug("getOneSliceAsPNG: entering.")
        relPNGFileDir = ('.').join(self.__relFITSFilePath.split('.')[:-1])
        absPNGFileDir = DataBlock.getPNGFilePrefix()+"/"+relPNGFileDir

        ittName = self.__getForcedTransformationName(ittName)

        status = None
        sliceData = None

        data = self.__data
        shape = data.shape
        self.__logger.debug(shape)

        if ( iFREQ >= shape[0] ):
            self.__logger.debug("No such slice index '%d' . Max. possible value is '%d'" %(iFREQ, shape[1]))
        elif ( len(shape) > 4 ):
            self.__logger.debug("Can't process data with more than 4 dimensions")
        else:
            sliceData = self.__getSlice(iFREQ)

            relPNGFilePath = ("%s/%d.%s.%s.%s.png" % (relPNGFileDir, iFREQ, ittName, lutName, vmName))
            absPNGFilePath = DataBlock.getPNGFilePrefix() + "/" + relPNGFilePath

            #if not os.path.exists(absPNGFilePath):
            if True:
                self.__logger.debug("File '%s' does not exists; it must be created..." % absPNGFilePath)

                try:
                    self.__PNGDir = absPNGFileDir
                    if not os.path.exists(self.__PNGDir):
                        os.makedirs(self.__PNGDir)
                except Exception as e:
                    self.__logger.debug("Failed to open '%s'. Message was '%s'" % (self.__PNGDir, e))
                    self.__logger.debug("Traceback: %s" %traceback.format_exc())
                    return {"status": False, "message": "Problem while creating the PNG file: '%s'" % e }

                self.__logger.debug("... in cache directory is %s" % (self.__PNGDir))
                try:
                    squaredSliceData = self.__squareSliceData( sliceData)
                    data_steps = self.__convertOneSlice2PNG_0(absPNGFilePath, squaredSliceData, ittName, lutName, vmName)

                    h = self.__header
                    if "CDELT1" in h and "CDELT2" in h:
                        x = {"BITPIX":   h["BITPIX"],
                            "NAXIS1":   squaredSliceData.shape[0],
                            "NAXIS2":   squaredSliceData.shape[0],
                            "CRPIX1":   h["CRPIX1"],
                            "CRPIX2":   h["CRPIX2"],
                            "CRVAL1":   h["CRVAL1"],
                            "CRVAL2":   h["CRVAL2"],
                            "CRTYPE1":  h["CTYPE1"],
                            "CRTYPE2":  h["CTYPE2"],
                            "RADESYS":  h["RADESYS"],
                            "CD1_1":    h["CDELT1"],
                            "CD2_2":    h["CDELT2"]}
                        self.__addWCStoPNG(absPNGFilePath, h)
                    status = True
                    self.__logger.debug("File '%s' has been created" % absPNGFilePath)

                except Exception as e:
                        self.__logger.debug("Failed to convert slice #%d. Error was '%r'" % (iFREQ, e))
                        self.__logger.debug("Traceback: %s" %traceback.format_exc())
                        return {"status": False, "message": "Problem while creating the PNG file: '%s'" % e }
            else:
                status = True
                self.__logger.debug("Using cached file '%s'" % absPNGFilePath)

        result = {"data_steps": data_steps, "path_to_png": relPNGFilePath, "statistics": self.__statistics["%d"%iFREQ]}
        self.__logger.debug("getOneSliceAsPNG: exiting.")

        return {"status": True, "message": "", "result": result}

    def getSummedSliceRangeAsPNG( self, iFREQ0, iFREQ1, ittName=__default_transformation_name, lutName=__default_palette_name, vmName=__default_video_mode_name):
        self.__logger.debug("getSummedSliceRangeAsPNG: entering.")

        relPNGFileDir = ('.').join(self.__relFITSFilePath.split('.')[:-1])
        absPNGFileDir = DataBlock.getPNGFilePrefix()+"/"+relPNGFileDir

        ittName = self.__getForcedTransformationName(ittName)

        relPNGFilePath = ("%s/%d-%d.%s.%s.%s.png" % (relPNGFileDir, iFREQ0, iFREQ1, ittName, lutName, vmName))
        absPNGFilePath = DataBlock.getPNGFilePrefix() + "/" + relPNGFilePath
        summedSliceRangeData = self.__getAverage_0(iFREQ0, iFREQ1)

        status = None
        #if not os.path.exists(absPNGFilePath):
        if True:
            self.__logger.debug("File '%s' does not exists; it must be created..." % absPNGFilePath)
            try:
                self.__PNGDir = absPNGFileDir
                if not os.path.exists(self.__PNGDir):
                    os.makedirs(self.__PNGDir)
                status = True
            except Exception as e:
                self.__logger.debug("Failed to open '%s'. Message was '%s'" % (self.__PNGDir, e))
                self.__logger.debug("Traceback: %s" %traceback.format_exc())
                return {"status": False, "message": "Problem while creating the PNG file: '%s'" % e }


            self.__logger.debug("... in cache directory is %s" % (self.__PNGDir))
            try:
                squaredData = self.__squareSliceData(summedSliceRangeData)

                data_steps = self.__convertSummedSliceRange2PNG(iFREQ0, iFREQ1, squaredData, ittName, lutName, vmName)
                h = self.__header
                if "CDELT1" in h and "CDELT2" in h:
                    x = {"BITPIX":   h["BITPIX"],
                    "NAXIS1":   squaredData.shape[0],
                    "NAXIS2":   squaredData.shape[0],
                    "CRPIX1":   h["CRPIX1"],
                    "CRPIX2":   h["CRPIX2"],
                    "CRVAL1":   h["CRVAL1"],
                    "CRVAL2":   h["CRVAL2"],
                    "CRTYPE1":  h["CTYPE1"],
                    "CRTYPE2":  h["CTYPE2"],
                    "RADESYS":  h["RADESYS"],
                    "CD1_1":    h["CDELT1"],
                    "CD2_2":    h["CDELT2"]}
                    self.__addWCStoPNG(absPNGFilePath, h)
                status = True
            except Exception as e:
                self.__logger.debug("Problem while creating the PNG file: '%s'" % e )
                self.__logger.debug("Traceback: %s" %traceback.format_exc())
                return {"status": False, "message": "Problem while creating the PNG file: '%s'" % e }
        else:
            status = True