adding tests

astrocut/__init__.py

@@ -24,4 +24,4 @@ class UnsupportedPythonError(Exception):
 if not _ASTROPY_SETUP_:  # noqa
     from .make_cube import CubeFactory  # noqa
     from .cube_cut import CutoutFactory  # noqa
-    from .fits_cut import fits_cut  # noqa
+    from .cutouts import fits_cut  # noqa

astrocut/fits_cut.py → astrocut/cutouts.py

@@ -53,7 +53,7 @@ def _get_cutout_limits(img_wcs, center_coord, cutout_size):
     except wcs.NoConvergence:  # If wcs can't converge, center coordinate is far from the footprint
         raise InvalidQueryError("Cutout location is not in image footprint!")
 
-    print(f"Center pixel: {center_pixel}")
+    #print(f"Center pixel: {center_pixel}")
 
     lims = np.zeros((2, 2), dtype=int)
 
@@ -74,7 +74,7 @@ def _get_cutout_limits(img_wcs, center_coord, cutout_size):
         # The case where the requested area is so small it rounds to zero
         if lims[axis, 0] == lims[axis, 1]:
             lims[axis, 0] = int(np.floor(center_pixel[axis] - 1))
-            lims[axis, 1] = int(np.ceil(center_pixel[axis] - 1))
+            lims[axis, 1] = lims[axis, 0] + 1 #int(np.ceil(center_pixel[axis] - 1))
 
     return lims
 
@@ -198,13 +198,13 @@ def _hducut(img_hdu, center_coord, cutout_size, correct_wcs=False, drop_after=No
         padding[1, 0] = -xmin
         xmin = 0
     if ymin < 0:
-        padding[2, 0] = -ymin
+        padding[0, 0] = -ymin
         ymin = 0
     if xmax > xmax_img:
         padding[1, 1] = xmax - xmax_img
         xmax = xmax_img
     if ymax > ymax_img:
-        padding[2, 1] = ymax - ymax_img
+        padding[0, 1] = ymax - ymax_img
         ymax = ymax_img  
 
     img_cutout = img_hdu.data[ymin:ymax, xmin:xmax]
@@ -415,7 +415,7 @@ def fits_cut(input_files, coordinates, cutout_size, correct_wcs=False, drop_afte
         for fle in input_files:
             cutout = cutout_hdu_dict[fle]
             if cutout.header.get("EMPTY"):
-                warning.warn("Cutout of {} contains to data and will not be written.".format(fle))
+                warnings.warn("Cutout of {} contains to data and will not be written.".format(fle))
                 continue
 
             cutout_hdus.append(cutout)

astrocut/tests/test_fits_cut.py

@@ -0,0 +1,193 @@
+import numpy as np
+
+from astropy.io import fits
+from astropy import wcs
+from astropy.coordinates import SkyCoord
+import astropy.units as u
+
+from .utils_for_test import create_test_imgs
+from .. import cutouts
+
+
+def test_get_cutout_limits():
+
+    test_img_wcs_kwds = fits.Header(cards=[('NAXIS', 2, 'number of array dimensions'),
+                                       ('NAXIS1', 20, ''),
+                                       ('NAXIS2', 30, ''),
+                                       ('CTYPE1', 'RA---TAN', 'Right ascension, gnomonic projection'),
+                                       ('CTYPE2', 'DEC--TAN', 'Declination, gnomonic projection'),
+                                       ('CRVAL1', 100, '[deg] Coordinate value at reference point'),
+                                       ('CRVAL2', 20, '[deg] Coordinate value at reference point'),
+                                       ('CRPIX1', 10, 'Pixel coordinate of reference point'),
+                                       ('CRPIX2', 15, 'Pixel coordinate of reference point'),
+                                       ('CDELT1', 1.0, '[deg] Coordinate increment at reference point'),
+                                       ('CDELT2', 1.0, '[deg] Coordinate increment at reference point'),
+                                       ('WCSAXES', 2, 'Number of coordinate axes'),
+                                       ('PC1_1', 1, 'Coordinate transformation matrix element'),
+                                       ('PC2_2', 1, 'Coordinate transformation matrix element'),
+                                       ('CUNIT1', 'deg', 'Units of coordinate increment and value'),
+                                       ('CUNIT2', 'deg', 'Units of coordinate increment and value')])
+    test_img_wcs = wcs.WCS(test_img_wcs_kwds)
+
+    center_coord = SkyCoord("100 20", unit='deg')
+    cutout_size = [10,10]
+
+    lims = cutouts._get_cutout_limits(test_img_wcs, center_coord, cutout_size)
+    assert (lims[0,1] - lims[0,0]) == (lims[1,1] - lims[1,0])
+    assert (lims == np.array([[ 4, 14], [ 9, 19]])).all()
+
+    cutout_size = [10,5]
+    lims = cutouts._get_cutout_limits(test_img_wcs, center_coord, cutout_size)
+    assert (lims[0,1] - lims[0,0]) == 10
+    assert (lims[1,1] - lims[1,0]) == 5
+
+    cutout_size = [.1,.1]*u.deg
+    lims = cutouts._get_cutout_limits(test_img_wcs, center_coord, cutout_size)
+    assert (lims[0,1] - lims[0,0]) == (lims[1,1] - lims[1,0])
+    assert (lims[0,1] - lims[0,0]) == 1
+
+    cutout_size = [4,5]*u.deg
+    lims = cutouts._get_cutout_limits(test_img_wcs, center_coord, cutout_size)
+    assert (lims[0,1] - lims[0,0]) == 4
+    assert (lims[1,1] - lims[1,0]) == 5
+
+    center_coord = SkyCoord("90 20", unit='deg')
+    cutout_size = [4,5]*u.deg
+    lims = cutouts._get_cutout_limits(test_img_wcs, center_coord, cutout_size)
+    assert lims[0,0] < 0
+
+    center_coord = SkyCoord("100 5", unit='deg')
+    cutout_size = [4,5]*u.pixel
+    lims = cutouts._get_cutout_limits(test_img_wcs, center_coord, cutout_size)
+    assert lims[1,0] < 0
+
+
+def test_get_cutout_wcs():
+    test_img_wcs_kwds = fits.Header(cards=[('NAXIS', 2, 'number of array dimensions'),
+                                       ('NAXIS1', 20, ''),
+                                       ('NAXIS2', 30, ''),
+                                       ('CTYPE1', 'RA---TAN', 'Right ascension, gnomonic projection'),
+                                       ('CTYPE2', 'DEC--TAN', 'Declination, gnomonic projection'),
+                                       ('CRVAL1', 100, '[deg] Coordinate value at reference point'),
+                                       ('CRVAL2', 20, '[deg] Coordinate value at reference point'),
+                                       ('CRPIX1', 10, 'Pixel coordinate of reference point'),
+                                       ('CRPIX2', 15, 'Pixel coordinate of reference point'),
+                                       ('CDELT1', 1.0, '[deg] Coordinate increment at reference point'),
+                                       ('CDELT2', 1.0, '[deg] Coordinate increment at reference point'),
+                                       ('WCSAXES', 2, 'Number of coordinate axes'),
+                                       ('PC1_1', 1, 'Coordinate transformation matrix element'),
+                                       ('PC2_2', 1, 'Coordinate transformation matrix element'),
+                                       ('CUNIT1', 'deg', 'Units of coordinate increment and value'),
+                                       ('CUNIT2', 'deg', 'Units of coordinate increment and value')])
+    test_img_wcs = wcs.WCS(test_img_wcs_kwds)
+
+    center_coord = SkyCoord("100 20", unit='deg')
+    cutout_size = [4,5]*u.deg
+    lims = cutouts._get_cutout_limits(test_img_wcs, center_coord, cutout_size)
+    cutout_wcs = cutouts._get_cutout_wcs(test_img_wcs, lims)
+    assert (cutout_wcs.wcs.crval == [100,20]).all()
+    assert (cutout_wcs.wcs.crpix == [3,4]).all()
+
+    center_coord = SkyCoord("100 5", unit='deg')
+    cutout_size = [4,5]*u.deg
+    lims = cutouts._get_cutout_limits(test_img_wcs, center_coord, cutout_size)
+    cutout_wcs = cutouts._get_cutout_wcs(test_img_wcs, lims)
+    assert (cutout_wcs.wcs.crval == [100,20]).all()
+    assert (cutout_wcs.wcs.crpix == [3,19]).all()
+
+    center_coord = SkyCoord("110 20", unit='deg')
+    cutout_size = [10,10]*u.deg
+    lims = cutouts._get_cutout_limits(test_img_wcs, center_coord, cutout_size)
+    cutout_wcs = cutouts._get_cutout_wcs(test_img_wcs, lims)
+    assert (cutout_wcs.wcs.crval == [100,20]).all()
+    assert (cutout_wcs.wcs.crpix == [-3,6]).all() 
+
+
+def test_fits_cut(tmpdir):
+
+    test_images = create_test_imgs(50, 6)
+
+    # Single file
+    center_coord = SkyCoord("150.1163213 2.200973097", unit='deg')
+    cutout_size = 10
+    cutout_file = cutouts.fits_cut(test_images, center_coord, cutout_size, single_outfile=True)
+    assert isinstance(cutout_file, str)
+
+    cutout_hdulist = fits.open(cutout_file)
+    assert len(cutout_hdulist) == len(test_images) + 1 # num imgs + primary header
+
+    cut1 = cutout_hdulist[1].data
+    assert cut1.shape == (cutout_size,cutout_size)
+    assert cutout_hdulist[1].data.shape == cutout_hdulist[2].data.shape
+    assert cutout_hdulist[2].data.shape == cutout_hdulist[3].data.shape
+    assert cutout_hdulist[3].data.shape == cutout_hdulist[4].data.shape
+    assert cutout_hdulist[4].data.shape == cutout_hdulist[5].data.shape
+    assert cutout_hdulist[5].data.shape == cutout_hdulist[6].data.shape
+
+    cut_wcs = wcs.WCS(cutout_hdulist[1].header)
+    sra, sdec = cut_wcs.all_pix2world(cutout_size/2,cutout_size/2,0)
+    assert round(float(sra), 4) == round(center_coord.ra.deg, 4)
+    assert round(float(sdec), 4) == round(center_coord.dec.deg, 4)
+
+    cutout_hdulist.close()
+
+    # Multiple files
+    cutout_files = cutouts.fits_cut(test_images, center_coord, cutout_size,
+                                     drop_after="Dummy1", single_outfile=False)
+
+    assert isinstance(cutout_files, list)
+    assert len(cutout_files) == len(test_images)
+
+    cutout_hdulist = fits.open(cutout_files[0])
+    assert len(cutout_hdulist) == 1
+
+    cut1 = cutout_hdulist[0].data
+    assert cut1.shape == (cutout_size,cutout_size)
+
+    cut_wcs = wcs.WCS(cutout_hdulist[0].header)
+    sra,sdec = cut_wcs.all_pix2world(cutout_size/2,cutout_size/2,0)
+    assert round(float(sra), 4) == round(center_coord.ra.deg, 4)
+    assert round(float(sdec), 4) == round(center_coord.dec.deg, 4)
+
+    cutout_hdulist.close()
+
+    # Do an off the edge test
+    center_coord = SkyCoord("150.1163213 2.2005731", unit='deg')
+    cutout_file = cutouts.fits_cut(test_images, center_coord, cutout_size, single_outfile=True)
+    assert isinstance(cutout_file, str)
+
+    cutout_hdulist = fits.open(cutout_file)
+    assert len(cutout_hdulist) == len(test_images) + 1 # num imgs + primary header
+
+    cut1 = cutout_hdulist[1].data
+    assert cut1.shape == (cutout_size,cutout_size)
+    assert np.isnan(cut1[:cutout_size//2,:]).all()
+
+    cutout_hdulist.close()
+
+    # Test when cutout is in some images not others
+
+    # Putting zeros into 2 images
+    for i in range(2):
+        hdu = fits.open(test_images[i], mode="update")
+        hdu[0].data[:20,:] = 0
+        hdu.flush()
+        hdu.close()
+
+
+    center_coord = SkyCoord("150.1163213 2.2007", unit='deg')
+    cutout_file = cutouts.fits_cut(test_images, center_coord, cutout_size, single_outfile=True)
+
+    cutout_hdulist = fits.open(cutout_file)
+    assert len(cutout_hdulist) == len(test_images) + 1 # num imgs + primary header
+    assert (cutout_hdulist[1].data == 0).all()
+    assert (cutout_hdulist[2].data == 0).all()
+    assert ~(cutout_hdulist[3].data == 0).any()
+    assert ~(cutout_hdulist[4].data == 0).any()
+    assert ~(cutout_hdulist[5].data == 0).any()
+    assert ~(cutout_hdulist[6].data == 0).any()
+
+
+    cutout_files = cutouts.fits_cut(test_images, center_coord, cutout_size, single_outfile=False)
+    assert isinstance(cutout_files, list)
+    assert len(cutout_files) == len(test_images) - 2

astrocut/tests/utils_for_test.py

@@ -1,11 +1,10 @@
 import numpy as np
 
-from astropy.io import fits
-
+from astropy.io import fits    
 
 def add_keywords(hdu, extname, time_increment, primary=False):
     """
-    Add a bunch of required keywords (mostly fake values)
+    Add a bunch of required keywords (mostly fake values).
     """
 
     hdu.header['extname'] = 'CAMERA.CCD 1.1 cal'
@@ -62,8 +61,7 @@ def create_test_ffis(img_size, num_images, basename='make_cube-test{:04d}.fits')
     Create test fits files
 
     Write negative values for data array and positive values for error array,
-    with unique values for all the pixels. The header keywords are just
-    populated to make things run but are not tested.
+    with unique values for all the pixels. 
     """
 
     img = np.arange(img_size*img_size, dtype=np.float32).reshape((img_size, img_size))
@@ -75,8 +73,7 @@ def create_test_ffis(img_size, num_images, basename='make_cube-test{:04d}.fits')
 
         primary_hdu = fits.PrimaryHDU()
         add_keywords(primary_hdu, "PRIMARY", i, primary=True)
-
-
+
         hdu = fits.ImageHDU(-img)
         add_keywords(hdu, 'CAMERA.CCD 1.1 cal', i)
 
@@ -89,3 +86,62 @@ def create_test_ffis(img_size, num_images, basename='make_cube-test{:04d}.fits')
         img = img + img_size*img_size
 
     return file_list
+
+
+def add_wcs_nosip_keywords(hdu, img_size):
+    """
+    Adds example wcs keywords without sip distortions to the given header.
+
+    Center coordinate is: 150.1163213, 2.200973097
+    """
+
+    hdu.header.extend([('WCSAXES', 2, 'Number of coordinate axes'),
+                       ('CRPIX1', img_size/2, 'Pixel coordinate of reference point'),
+                       ('CRPIX2', img_size/2, 'Pixel coordinate of reference point'),
+                       ('PC1_1', -1.666667e-05, 'Coordinate transformation matrix element'),
+                       ('PC2_2', 1.666667e-05, 'Coordinate transformation matrix element'),
+                       ('CDELT1', 1.0, '[deg] Coordinate increment at reference point'),
+                       ('CDELT2', 1.0, '[deg] Coordinate increment at reference point'),
+                       ('CUNIT1', 'deg', 'Units of coordinate increment and value'),
+                       ('CUNIT2', 'deg', 'Units of coordinate increment and value'),
+                       ('CTYPE1', 'RA---TAN', 'Right ascension, gnomonic projection'),
+                       ('CTYPE2', 'DEC--TAN', 'Declination, gnomonic projection'),
+                       ('CRVAL1', 150.1163213, '[deg] Coordinate value at reference point'),
+                       ('CRVAL2', 2.200973097, '[deg] Coordinate value at reference point')])
+
+def add_dummy_keywords(hdu):
+    """
+    Adding a number of dummy keywords, basically so the drop_after argument to fits_cut can be tested.
+    """
+    hdu.header.extend([('Dummy1', "Dummy1", 'Dummy1'),
+                       ('Dummy2', 2, 'Dummy2'),
+                       ('Dummy3', "", 'Dummy3')],strip=False)
+
+
+def create_test_imgs(img_size, num_images, dummy_keywords=True, basename='img_{:04d}.fits'):
+    """
+    Create test fits image files, single extension.
+
+    Write unique values for all the pixels. 
+    The header keywords are populated with a simple WCS for testing.
+    """
+
+    img = np.arange(img_size*img_size, dtype=np.float32).reshape((img_size, img_size))
+    file_list = []
+
+    for i in range(num_images):
+
+        file_list.append(basename.format(i))
+
+        primary_hdu = fits.PrimaryHDU(data=img)
+        add_wcs_nosip_keywords(primary_hdu, img_size)
+
+        if dummy_keywords:
+            add_dummy_keywords(primary_hdu)
+
+        hdulist = fits.HDUList([primary_hdu])
+        hdulist.writeto(file_list[-1], overwrite=True, checksum=True)
+
+        img = img + img_size*img_size
+
+    return file_list

docs/astrocut/index.rst

@@ -152,8 +152,7 @@ that all input image are aligned and have the same pixel scale, no checking is d
                 >>> center_coord = SkyCoord("150.0945 2.38681", unit='deg')
                 >>> cutout_size = [200,300]
                 >>> 
-                >>> cutout_file = astrocut.fits_cut(candels_fles, center_coord, cutout_size,
-                ...                                 drop_after="", single_outfile=True)
+                >>> cutout_file = fits_cut(input_files, center_coord, cutout_size, drop_after="", single_outfile=True)
                 >>> print(cutout_file)
                 cutout_150.094500_2.386810_200x300_astrocut.fits