#147 Add option to have weight and badpix images from a different file than the main image.

piff/input.py

@@ -208,6 +208,10 @@ class InputFiles(Input):
                         weight == 0. [default: None]
         :noise:         Rather than a weight image, provide the noise variance in the image.
                         (Useful for simulations where this is a known value.) [default: None]
+        :weight_file_name: If the weight image is in a different file than the main image,
+                        this gives the name of the file to use. [default: None]
+        :badpix_file_name: If the badpix image is in a different file than the main image,
+                        this gives the name of the file to use. [default: None]
 
         :cat_hdu:       The hdu to use in the catalog files. [default: 1]
         :x_col:         The name of the X column in the input catalogs. [default: 'x']
@@ -291,6 +295,8 @@ def __init__(self, config, logger=None):
                 'image_hdu' : int,
                 'weight_hdu' : int,
                 'badpix_hdu' : int,
+                'weight_file_name' : str,
+                'badpix_file_name' : str,
                 'sky_file_name' : str,
                 'sky_hdu': int,
                 'cat_hdu' : int,
@@ -450,7 +456,9 @@ def __init__(self, config, logger=None):
             # Read the image
             image_file_name = params['image_file_name']
             image_hdu = params.get('image_hdu', None)
+            weight_file_name = params.get('weight_file_name', None)
             weight_hdu = params.get('weight_hdu', None)
+            badpix_file_name = params.get('badpix_file_name', None)
             badpix_hdu = params.get('badpix_hdu', None)
             sky_file_name = params.get('sky_file_name', None)
             sky_hdu = params.get('sky_hdu', None)
@@ -460,7 +468,9 @@ def __init__(self, config, logger=None):
             self.image_kwargs.append({
                     'image_file_name' : image_file_name,
                     'image_hdu' : image_hdu,
+                    'weight_file_name' : weight_file_name,
                     'weight_hdu' : weight_hdu,
+                    'badpix_file_name' : badpix_file_name,
                     'badpix_hdu' : badpix_hdu,
                     'sky_file_name' : sky_file_name,
                     'sky_hdu' : sky_hdu,
@@ -730,13 +740,15 @@ def _removeSignalFromWeight(image, weight, gain=None):
         return newweight, gain
 
     @staticmethod
-    def readImage(image_file_name, image_hdu, weight_hdu, badpix_hdu,
-                  sky_file_name, sky_hdu, noise, logger):
+    def readImage(image_file_name, image_hdu, weight_file_name, weight_hdu,
+                  badpix_file_name, badpix_hdu, sky_file_name, sky_hdu, noise, logger):
         """Read in the image and weight map (or make one if no weight information is given
 
         :param image_file_name: The name of the file to read.
         :param image_hdu:       The hdu of the main image.
+        :param weight_file_name: The name if the file with the weight image (if any).
         :param weight_hdu:      The hdu of the weight image (if any).
+        :param badpix_file_name: The name of the file with the bad pixel mask (if any).
         :param badpix_hdu:      The hdu of the bad pixel mask (if any).
         :param sky_file_name:   A file to use for a sky background to subtract from the image
                                 (if any).
@@ -758,9 +770,11 @@ def readImage(image_file_name, image_hdu, weight_hdu, badpix_hdu,
             image -= sky
 
         # Either read in the weight image, or build a dummy one
-        if weight_hdu is not None:
-            logger.info("Reading weight image from hdu %d.", weight_hdu)
-            weight = galsim.fits.read(image_file_name, hdu=weight_hdu)
+        if weight_file_name is not None or weight_hdu is not None:
+            weight_file_name = weight_file_name or image_file_name
+            weight_hdu = weight_hdu or (1 if weight_file_name.endswith('.fz') else 0)
+            logger.info("Reading weight image from %s, hdu %d.", weight_file_name, weight_hdu)
+            weight = galsim.fits.read(weight_file_name, hdu=weight_hdu)
             if np.all(weight.array == 0):
                 logger.error("According to the weight mask in %s, all pixels have zero weight!",
                              image_file_name)
@@ -776,9 +790,11 @@ def readImage(image_file_name, image_hdu, weight_hdu, badpix_hdu,
             weight = galsim.ImageF(image.bounds, init_value=1)
 
         # If requested, set wt=0 for any bad pixels
-        if badpix_hdu is not None:
-            logger.info("Reading badpix image from hdu %d.", badpix_hdu)
-            badpix = galsim.fits.read(image_file_name, hdu=badpix_hdu)
+        if badpix_file_name is not None or badpix_hdu is not None:
+            badpix_file_name = badpix_file_name or image_file_name
+            badpix_hdu = badpix_hdu or (1 if badpix_file_name.endswith('.fz') else 0)
+            logger.info("Reading badpix image from %s, hdu %d.", badpix_file_name, badpix_hdu)
+            badpix = galsim.fits.read(badpix_file_name, hdu=badpix_hdu)
             # The badpix image may be offset by 32768 from the true value.
             # If so, subtract it off.
             if np.any(badpix.array > 32767):  # pragma: no cover

tests/input/make_input.yaml

@@ -50,7 +50,6 @@ output:
         hdu: 2
 
     truth:
-        dir: 'input'
         file_name:
             type: FormattedStr
             format: 'test_input_cat_%02d.fits'

tests/test_input.py

@@ -18,14 +18,18 @@
 import numpy as np
 import fitsio
 import piff
+import pytest
 
 from piff_test_helper import get_script_name, timer, CaptureLog
 
-@timer
+@pytest.fixture(scope="module", autouse=True)
 def setup():
     """Make sure the images and catalogs that we'll use throughout this module are done first.
     """
-    gs_config = galsim.config.ReadConfig(os.path.join('input','make_input.yaml'))[0]
+    dir = os.path.dirname(os.path.realpath(__file__))
+    config_file_name = os.path.join(dir, 'input', 'make_input.yaml')
+    gs_config = galsim.config.ReadConfig(config_file_name)[0]
+    gs_config['output']['dir'] = os.path.join(dir, 'input')
     galsim.config.BuildFiles(2, galsim.config.CopyConfig(gs_config))
 
     # For the third file, add in a real wcs and ra, dec to the output catalog.
@@ -42,7 +46,7 @@ def setup():
     gs_config['output']['truth']['columns']['dec'] = '$wcs.toWorld(image_pos).dec / galsim.degrees'
     galsim.config.BuildFiles(1, galsim.config.CopyConfig(gs_config), file_num=2)
 
-    cat_file_name = os.path.join('input', 'test_input_cat_00.fits')
+    cat_file_name = os.path.join(dir, 'input', 'test_input_cat_00.fits')
     data = fitsio.read(cat_file_name)
     sky = np.mean(data['sky'])
     gain = np.mean(data['gain'])
@@ -51,7 +55,7 @@ def setup():
 
     # Write a sky image file
     sky_im = galsim.Image(1024,1024, init_value=sky)
-    sky_im.write(os.path.join('input', 'test_input_sky_00.fits.fz'))
+    sky_im.write(os.path.join(dir, 'input', 'test_input_sky_00.fits.fz'))
 
     # Add two color columns to first catalog file
     rng1 = np.random.default_rng(123)
@@ -65,7 +69,7 @@ def setup():
 
     # Add some header values to the first one.
     # Also add some alternate weight and badpix maps to enable some edge-case tests
-    image_file = os.path.join('input','test_input_image_00.fits')
+    image_file = os.path.join(dir, 'input', 'test_input_image_00.fits')
     with fitsio.FITS(image_file, 'rw') as f:
         f[0].write_key('SKYLEVEL', sky, 'sky level')
         f[0].write_key('GAIN_A', gain, 'gain')
@@ -91,10 +95,21 @@ def setup():
         bp[:,:] = 0
         bp[1::2,:] = 16 # Odd cols
         f.write(bp) # hdu = 9
-        wt = f[1].read().copy()
-        var = 1/wt + (f[0].read() - sky) / gain
+        wtx = f[1].read().copy()
+        var = 1/wtx + (f[0].read() - sky) / gain
         f.write(var) # hdu = 10
 
+    weight_file = os.path.join(dir, 'input', 'test_input_weight_00.fits')
+    with fitsio.FITS(weight_file, 'rw', clobber=True) as f:
+        # hdu 0 here is equivalent to hdu 10 above.
+        f.write(var)
+        # hdu 1 here is equivalent to hdu 8 above.
+        f.write(wt)
+
+    badpix_file = os.path.join(dir, 'input', 'test_input_badpix_00.fits')
+    with fitsio.FITS(badpix_file, 'rw', clobber=True) as f:
+        # hdu 0 here is equivalent to hdu 9 above.
+        f.write(bp)
 
 @timer
 def test_basic():
@@ -1047,7 +1062,20 @@ def test_weight():
     assert weight.array.shape == (1024, 1024)
     np.testing.assert_almost_equal(weight.array, 0.)
 
+    # Use a different file for badpix and weight.  These are hdu 2,1 in this file.
+    config['weight_file_name'] = 'input/test_input_weight_00.fits'
+    config['weight_hdu'] = 1
+    config['badpix_file_name'] = 'input/test_input_badpix_00.fits'
+    config['badpix_hdu'] = 0
+    input = piff.InputFiles(config, logger=logger)
+    assert input.nimages == 1
+    _, weight, _, _ = input.getRawImageData(0)
+    assert weight.array.shape == (1024, 1024)
+    np.testing.assert_almost_equal(weight.array, 0.)
+
     # Negative valued weights are invalid
+    del config['weight_file_name']
+    del config['badpix_file_name']
     config['weight_hdu'] = 4
     input = piff.InputFiles(config)
     with CaptureLog() as cl:
@@ -1127,6 +1155,29 @@ def test_lsst_weight():
     print('var = ',1./weight.array)
     np.testing.assert_allclose(weight.array, expected_noise**-1, rtol=1.e-5)
 
+    # Finally, check using a different file for the weight image.
+    config = {
+                'image_file_name' : 'input/test_input_image_00.fits',
+                'cat_file_name' : 'input/test_input_cat_00.fits',
+                'weight_file_name' : 'input/test_input_weight_00.fits',
+                'gain' : 'GAIN_A',
+                'invert_weight' : True,
+                'remove_signal_from_weight' : True,
+             }
+    input = piff.InputFiles(config, logger=logger)
+    assert input.nimages == 1
+    image, weight1, image_pos, props = input.getRawImageData(0)
+    assert len(image_pos) == 100
+    assert image.array.shape == (1024, 1024)
+    assert weight.array.shape == (1024, 1024)
+    np.testing.assert_array_equal(weight1, weight)  # From previous test using hdu 10
+    read_noise = 10
+    expected_noise = read_noise**2 / gain**2
+    print('expected noise = ',expected_noise)
+    print('var = ',1./weight.array)
+    np.testing.assert_allclose(weight.array, expected_noise**-1, rtol=1.e-5)
+
+
 
 @timer
 def test_stars():

tests/test_meanify.py

@@ -20,6 +20,7 @@
 import os
 import fitsio
 import glob
+import pytest
 
 from piff_test_helper import get_script_name, timer
 
@@ -123,6 +124,7 @@ def make_average(coord=None, gp=True):
 
     return params
 
+@pytest.fixture(scope="module", autouse=True)
 def setup():
     np.random.seed(42)
     if __name__ == '__main__':
@@ -133,7 +135,8 @@ def setup():
         stars_per_image = 50
 
     # Delete any existing image files
-    for image_file in glob.glob(os.path.join('output','test_mean_image_*.fits')):
+    dir = os.path.dirname(os.path.realpath(__file__))
+    for image_file in glob.glob(os.path.join(dir, 'output', 'test_mean_image_*.fits')):
         os.remove(image_file)
 
     for k in range(nimages):
@@ -163,19 +166,19 @@ def setup():
             offset = galsim.PositionD( x-int(x)-0.5 , y-int(y)-0.5 )
             psf.drawImage(image=image[bounds], method='no_pixel', offset=offset)
 
-        image_file = os.path.join('output','test_mean_image_%02i.fits'%k)
+        image_file = os.path.join(dir, 'output', 'test_mean_image_%02i.fits'%k)
         image.write(image_file)
 
         dtype = [ ('x','f8'), ('y','f8') ]
         data = np.empty(len(x_list), dtype=dtype)
         data['x'] = x_list
         data['y'] = y_list
-        cat_file = os.path.join('output','test_mean_cat_%02i.fits'%k)
+        cat_file = os.path.join(dir, 'output', 'test_mean_cat_%02i.fits'%k)
         fitsio.write(cat_file, data, clobber=True)
 
-        image_file = os.path.join('output','test_mean_image_%02i.fits'%k)
-        cat_file = os.path.join('output','test_mean_cat_%02i.fits'%k)
-        psf_file = os.path.join('output','test_mean_%02i.piff'%k)
+        image_file = os.path.join(dir, 'output', 'test_mean_image_%02i.fits'%k)
+        cat_file = os.path.join(dir, 'output', 'test_mean_cat_%02i.fits'%k)
+        psf_file = os.path.join(dir, 'output', 'test_mean_%02i.piff'%k)
 
         config = {
             'input' : {

tests/test_stats.py

@@ -20,6 +20,7 @@
 import fitsio
 import yaml
 import subprocess
+import pytest
 
 from piff_test_helper import get_script_name, timer, CaptureLog
 
@@ -222,7 +223,7 @@ def generate_starlist(n_samples=500):
 
     return star_list, model
 
-@timer
+@pytest.fixture(scope="module", autouse=True)
 def setup():
     """Build an input image and catalog used by a few tests below.
     """
@@ -254,15 +255,16 @@ def setup():
     image.addNoise(galsim.GaussianNoise(rng=galsim.BaseDeviate(1234), sigma=1e-6))
 
     # Write out the image to a file
-    image_file = os.path.join('output','test_stats_image.fits')
+    dir = os.path.dirname(os.path.realpath(__file__))
+    image_file = os.path.join(dir, 'output', 'test_stats_image.fits')
     image.write(image_file)
 
     # Write out the catalog to a file
     dtype = [ ('x','f8'), ('y','f8') ]
     data = np.empty(len(x_list), dtype=dtype)
     data['x'] = x_list
     data['y'] = y_list
-    cat_file = os.path.join('output','test_stats_cat.fits')
+    cat_file = os.path.join(dir, 'output', 'test_stats_cat.fits')
     fitsio.write(cat_file, data, clobber=True)
 
 @timer