apply fixes for issues #129 and #131

AegeanTools/AeRes.py

@@ -19,24 +19,60 @@
 def load_sources(filename):
     """
     Open a file, read contents, return a list of all the sources in that file.
-    @param filename:
-    @return: list of ComponentSource objects
+
+    Parameters
+    ----------
+    filename : str
+        Filename to be read
+
+    Return
+    ------
+    catalog : [`class:AegeanTools.models.ComponentSource`, ...]
+        A list of source components
     """
-    catalog = catalogs.table_to_source_list(catalogs.load_table(filename))
+    table = catalogs.load_table(filename)
+    required_cols = ['ra','dec','peak_flux','a','b','pa']
+    good = True
+    for c in required_cols:
+        if c not in table.colnames:
+            logging.error("Column {0} not found".format(c))
+            good = False
+    if not good:
+        logging.error("Some required columns missing")
+        return None
+    catalog = catalogs.table_to_source_list(table)
     logging.info("read {0} sources from {1}".format(len(catalog), filename))
     return catalog
 
 
 def make_model(sources, shape, wcshelper, mask=False, frac=None, sigma=4):
     """
+    Create a model image based on a catalogue of sources.
+
+    Parameters
+    ----------
+    sources : [`class:AegeanTools.models.ComponentSource`, ...]
+        a list of sources
+
+    shape : [float, float]
+        the shape of the input (and output) image
+
+    wcshelper : 'class:AegeanTools.wcs_helpers.WCSHelper'
+        A WCSHelper object corresponding to the input image
 
-    @param sources: a list of AegeanTools.models.SimpleSource objects
-    @param shape: the shape of the input (and output) image
-    @param wcshelper: an AegeanTools.wcs_helpers.WCSHelper object corresponding to the input image
-    @param mask: If true then mask pixels instead of subtracting the sources
-    @param frac: pixels that are brighter than frac*peak_flux for each source will be masked if mask=True
-    @param sigma: pixels that are brighter than rms*sigma be masked if mask=True
-    @return:
+    mask : bool
+        If true then mask pixels instead of subtracting or adding sources
+
+    frac : float
+        pixels that are brighter than frac*peak_flux for each source will be masked if mask=True
+
+    sigma: float
+        pixels that are brighter than rms*sigma be masked if mask=True
+
+    Returns
+    -------
+    model : np.ndarray
+        The desired model.
     """
 
     # Model array
@@ -85,9 +121,9 @@ def make_model(sources, shape, wcshelper, mask=False, frac=None, sigma=4):
             logging.debug(" flux, sx, sy: {0} {1} {2}".format(src.peak_flux, sx, sy))
 
         # positions for which we want to make the model
-        x, y = np.mgrid[xmin:xmax, ymin:ymax]
-        x = list(map(int, x.ravel()))
-        y = list(map(int, y.ravel()))
+        x, y = np.mgrid[int(xmin):int(xmax), int(ymin):int(ymax)]
+        x = x.ravel()
+        y = y.ravel()
 
         # TODO: understand why xo/yo -1 is needed
         model = fitting.elliptical_gaussian(x, y, src.peak_flux, xo-1, yo-1, sx*FWHM2CC, sy*FWHM2CC, theta)
@@ -98,27 +134,55 @@ def make_model(sources, shape, wcshelper, mask=False, frac=None, sigma=4):
                 indices = np.where(model >= (frac*src.peak_flux))
             else:
                 indices = np.where(model >= (sigma*src.local_rms))
-            model[indices] = np.nan
-
-        m[x, y] += model
+            # somehow m[x,y][indices] = np.nan doesn't assign any values
+            # so we have to do the more complicated
+            # m[x[indices],y[indices]] = np.nan
+            m[x[indices], y[indices]]= np.nan
+        else:
+            m[x, y] += model
         i_count += 1
     logging.info("modeled {0} sources".format(i_count))
     return m
 
 
 def make_residual(fitsfile, catalog, rfile, mfile=None, add=False, mask=False, frac=None, sigma=4):
     """
+    Take an input image and catalogue, make a model of the catalogue, and then add/subtract or mask the input image.
+    Saving the residual and (optionally) model files.
+
+    Parameters
+    ----------
+    fitsfile : str
+        Input fits image filename
+
+    catalog : str
+        Input catalog filename of a type supported by Aegean
+
+    rfile : str
+        Filename to write residual image
+
+    mfile : str
+        Filename to write model image. Default=None means don't write the model file.
+
+    add : bool
+        If True add the model instead of subtracting it
+
+    mask : bool
+        If true then mask pixels instead of adding or subtracting the sources
+
+    frac : float
+        pixels that are brighter than frac*peak_flux for each source will be masked if mask=True
+
+    sigma : float
+        pixels that are brighter than sigma*local_rms for each source will be masked if mask=True
 
-    @param fitsfile: Input fits image filename
-    @param catalog: Input catalog filename of a type supported by Aegean
-    @param rfile: Residual image filename
-    @param mfile: model image filename
-    @param add: add the model instead of subtracting it
-    @param mask: If true then mask pixels instead of subtracting the sources
-    @param frac: pixels that are brighter than frac*peak_flux for each source will be masked if mask=True
-    @return:
+    Return
+    ------
+    None
     """
     source_list = load_sources(catalog)
+    if source_list is None:
+        return None
     # force two axes so that we dump redundant stokes/freq axes if they are present.
     hdulist = fits.open(fitsfile, naxis=2)
     # ignore dimensions of length 1

scripts/AeRes

@@ -11,8 +11,8 @@ from optparse import OptionParser
 import sys
 
 __author__ = 'Paul Hancock'
-__version__ = 'v0.2.6'
-__date__ = '2020-07-23'
+__version__ = 'v0.2.7'
+__date__ = '2020-07-30'
 
 
 # global constants

tests/test_AerRes.py → tests/test_AeRes.py

@@ -8,7 +8,7 @@
 __author__ = 'Paul Hancock'
 
 from AegeanTools import AeRes as ar
-from AegeanTools import wcs_helpers
+from AegeanTools import wcs_helpers, catalogs
 
 from astropy.io import fits
 import numpy as np
@@ -24,8 +24,22 @@ def test_load_sources():
     return
 
 
+def test_load_sources_missing_columns():
+    filename = 'tests/test_files/1904_comp.fits'
+    table = catalogs.load_table(filename)
+    table.rename_column('ra', 'RAJ2000')
+    table.write('dlme.fits')
+    cat = ar.load_sources('dlme.fits')
+    if os.path.exists('dlme.fits'):
+        os.remove('dlme.fits')
+
+    if cat is not None:
+        raise AssertionError("Missing columns should be caught, but weren't")
+    return
+
+
 def test_make_model():
-    """Test make_modell"""
+    """Test make_model"""
     filename = 'tests/test_files/1904_comp.fits'
     sources = ar.load_sources(filename)
     hdulist = fits.open('tests/test_files/1904-66_SIN.fits')
@@ -43,21 +57,37 @@ def test_make_model():
     if not np.all(model == 0.):
         raise AssertionError("Model is *not* empty")
 
+
+def test_make_masked_model():
+    """Test make_model when a mask is being used"""
+    filename = 'tests/test_files/1904_comp.fits'
+    sources = ar.load_sources(filename)
+    hdulist = fits.open('tests/test_files/1904-66_SIN.fits')
+    wcs_helper = wcs_helpers.WCSHelper.from_header(header=hdulist[0].header)
+
     # test mask with sigma
     model = ar.make_model(sources=sources, shape=hdulist[0].data.shape,
                           wcshelper=wcs_helper, mask=True)
+
+    finite = np.where(np.isfinite(model))
     if np.all(model == 0.):
         raise AssertionError("Model is empty")
     if not np.any(np.isnan(model)):
         raise AssertionError("Model is not masked")
+    if not np.all(model[finite] == 0.):
+        raise AssertionError("Model has values that are not zero or nan")
 
     # test mask with frac
     model = ar.make_model(sources=sources, shape=hdulist[0].data.shape,
                           wcshelper=wcs_helper, mask=True, frac=0.1)
+
+    finite = np.where(np.isfinite(model))
     if np.all(model == 0.):
         raise AssertionError("Model is empty")
     if not np.any(np.isnan(model)):
         raise AssertionError("Model is not masked")
+    if not np.all(model[finite] == 0.):
+        raise AssertionError("Model has values that are not zero or nan")
 
 
 def test_make_residual():