Commit 47931c2f authored by Simon Conseil's avatar Simon Conseil

Allow to correct zero offset in mpdaf_merging_sigma_clipping

parent 14b1b6c4
......@@ -40,6 +40,7 @@ from astropy import units as u
from astropy.table import Table
from astropy.utils.console import ProgressBar
from ctypes import c_char_p
from datetime import datetime
from numpy import allclose, array_equal
from .cube import Cube
......@@ -139,7 +140,7 @@ def _pycombine(self, nmax=2, nclip=5.0, var='propagate', nstop=2, nl=None,
if self.flux_scales is None and self.flux_offsets is None:
for l in range(nl):
if l % 100 == 0:
info('%d/%d', l, nl)
info('%d/%d %s', l, nl, datetime.now())
arr.fill(np.nan)
for i, f in enumerate(data):
x, y, x2, y2 = pos[i]
......@@ -158,17 +159,17 @@ def _pycombine(self, nmax=2, nclip=5.0, var='propagate', nstop=2, nl=None,
scales = np.ones(self.nfiles)
else:
scales = np.asarray(self.flux_scales)
self._logger.info('Using scales: %s', scales)
self._logger.info('Using scales')
if self.flux_offsets is None:
offsets = np.zeros(self.nfiles)
else:
offsets = np.asarray(self.flux_offsets)
self._logger.info('Using offsets: %s', offsets)
self._logger.info('Using offsets')
for l in range(nl):
if l % 100 == 0:
info('%d/%d', l, nl)
info('%d/%d %s', l, nl, datetime.now())
arr.fill(np.nan)
for i, f in enumerate(data):
x, y, x2, y2 = pos[i]
......@@ -335,6 +336,11 @@ class CubeList:
str(c.wcs.wcs.wcs.crpix), str(c.wcs.wcs.wcs.crval))
for c in self.cubes]
t = Table(rows=rows, names=('filename', 'shape', 'crpix', 'crval'))
if self.flux_scales is not None:
t['scale'] = self.flux_scales
if self.flux_offsets is not None:
t['offset'] = self.flux_offsets
for line in t.pformat():
self._logger.info(line)
......@@ -492,12 +498,19 @@ class CubeList:
files = files.encode('utf8')
if self.flux_scales is None:
scales = np.ones(self.nfiles)
scale = np.ones(self.nfiles)
else:
scales = np.asarray(self.flux_scales)
scale = np.asarray(self.flux_scales)
self._logger.info('Using scales')
if self.flux_offsets is None:
offset = np.zeros(self.nfiles)
else:
offset = np.asarray(self.flux_offsets)
self._logger.info('Using offsets')
ctools.mpdaf_merging_sigma_clipping(
c_char_p(files), data, vardata, expmap, scales, select_pix,
c_char_p(files), data, vardata, expmap, scale, offset, select_pix,
valid_pix, nmax, np.float64(nclip_low), np.float64(nclip_up),
nstop, np.int32(var_mean), np.int32(mad))
......
......@@ -167,8 +167,11 @@ class TestCubeList(unittest.TestCase):
@pytest.mark.skipif(not HAS_CFITSIO, reason="requires cfitsio")
def test_combine_scale(self):
clist = CubeList(self.cubenames, scalelist=[2.] * self.ncubes)
combined_cube = np.full(self.shape, 2 * 2, dtype=float)
clist = CubeList(self.cubenames, scalelist=[2.] * self.ncubes,
offsetlist=[2.] * self.ncubes)
combined_cube = np.full(self.shape,
np.mean((np.array(self.cubevals) + 2) * 2),
dtype=float)
cube, expmap, stat_pix = clist.combine(header={'FOO': 'BAR'})
assert_array_equal(cube.data, combined_cube)
......@@ -182,7 +185,9 @@ class TestCubeList(unittest.TestCase):
clist = CubeList(self.cubenames, scalelist=[2.] * self.ncubes,
offsetlist=[0.5] * self.ncubes)
combined_cube = np.full(self.shape, 5, dtype=float)
combined_cube = np.full(self.shape,
np.mean((np.array(self.cubevals) + 0.5) * 2),
dtype=float)
cube2, expmap2, _, _ = clist.pycombine(header={'FOO': 'BAR'})
assert_array_equal(cube2.data, combined_cube)
......
......@@ -81,6 +81,7 @@ else:
array_1d_double, # double* var
array_1d_int, # int* expmap
array_1d_double, # double* scale
array_1d_double, # double* offset
array_1d_int, # int* selected_pix
array_1d_int, # int* valid_pix
ctypes.c_int, # int nmax
......
......@@ -289,7 +289,7 @@ int mpdaf_merging_median(char* input, double* data, int* expmap, int* valid_pix)
// var=0: 'propagate'
// var=1: 'stat_mean'
// var=2: 'stat_one'
int mpdaf_merging_sigma_clipping(char* input, double* data, double* var, int* expmap, double* scale, int* selected_pix, int* valid_pix, int nmax, double nclip_low, double nclip_up, int nstop, int typ_var, int mad)
int mpdaf_merging_sigma_clipping(char* input, double* data, double* var, int* expmap, double* scale, double*offset, int* selected_pix, int* valid_pix, int nmax, double nclip_low, double nclip_up, int nstop, int typ_var, int mad)
{
char* filenames[MAX_FILES];
char buffer[80], begin[80];
......@@ -418,7 +418,7 @@ int mpdaf_merging_sigma_clipping(char* input, double* data, double* var, int* ex
{
if (!isnan(pix[i][ii]))
{
wdata[n] = pix[i][ii]*scale[i];
wdata[n] = (offset[i] + pix[i][ii]) * scale[i];
files_id[n] = i;
indx[n] = n;
if (typ_var==0)
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment