Integrate spectral cube through a filter function

[thomasorb]

Here is a first version of a scalable integration algorithm which takes not too much RAM but takes a long time to compute. It must be optimized.

class SpectralCubePatch(SpectralCube):
    def integrate(self, filter_function):
        """
        :param filter_function: Must be an orcs.core.Filter instance
        """
        if not isinstance(filter_function, Filter):
            raise TypeError('filter_function must be an orcs.core.Filter instance')

        if (filter_function.start <= self.params.base_axis[0]
            or filter_function.end >= self.params.base_axis[-1]):
            raise ValueError('filter passband (>5%) between {} - {} out of cube band {} - {}'.format(
                filter_function.start,
                filter_function.end,
                self.params.base_axis[0],
                self.params.base_axis[-1]))

        start_pix, end_pix = orb.utils.spectrum.cm12pix(
            self.params.base_axis, [filter_function.start, filter_function.end])

        sframe = np.zeros((self.dimx, self.dimy), dtype=float)
        for iz in range(start_pix, end_pix+1):
            sframe += cube.get_data_frame(iz) * filter_function(float(self.params.base_axis[iz]))

        sframe /= np.sum(filter_function(self.params.base_axis.astype(float)))
        return sframe

cube = SpectralCubePatch('/home/thomas/M31_SN3.merged.cm1.1.0.hdf5')
hdulist = io.read_fits('/home/thomas/Astro/Python/pysynphot_datafiles/grp/hst/cdbs/comp/wfc3/wfc3uvis1_f656n_006_syn.fits',
                      return_hdu_only=True)
wav = hdulist[1].data.field(0) / 10. # nm
tra = hdulist[1].data.field(1) # transmission
intframe = cube.integrate(Filter(wav,tra))def integrate
io.write_fits('integrated_frame', intframe, overwrite=True)

[BLaunet]

Here are the improvements I suggest :

• the function now supports optional xmin, xmax, ymin, ymax arguments, to reduce the region on which the integration is performed
• instead of loading frames one by one, we load them by chunk of max 10 frames (i.e. ~32Mo max loaded in RAM, which is still acceptable)

class SpectralCubePatch(SpectralCube):
    def integrate(self, filter_function, xmin=0, xmax=None, ymin=0, ymax=None):
        """
        :param filter_function: Must be an orcs.core.Filter instance
        """
        if not isinstance(filter_function, Filter):
            raise TypeError('filter_function must be an orcs.core.Filter instance')

        if (filter_function.start <= self.params.base_axis[0]
            or filter_function.end >= self.params.base_axis[-1]):
            raise ValueError('filter passband (>5%) between {} - {} out of cube band {} - {}'.format(
                filter_function.start,
                filter_function.end,
                self.params.base_axis[0],
                self.params.base_axis[-1]))

        if xmax is None:
            xmax = self.dimx
        if ymax is None:
            ymax = self.dimy
        start_pix, end_pix = orb.utils.spectrum.cm12pix(
            self.params.base_axis, [filter_function.start, filter_function.end])

        sframe = np.zeros((self.dimx, self.dimy), dtype=float)
        zsize = end_pix-start_pix+1
        izranges = np.array_split(range(start_pix, end_pix+1), zsize//10+1) #This splits the range in zsize//10 +1 chunks (not necessarily of same size). The endpix is correctly handled in the extraction
        for izrange in izranges:
            sframe[xmin:xmax, ymin:ymax] += np.sum(self.get_data(xmin, xmax, ymin, ymax, izrange.min(), izrange.max()+1, silent=True) 
                                                   * filter_function(self.params.base_axis[izrange].astype(float)), axis=2)
        sframe /= np.sum(filter_function(self.params.base_axis.astype(float)))
        return sframe