From f13c877e7719939844c97b1ce1b4581d9b724151 Mon Sep 17 00:00:00 2001 From: Francesco D'Eugenio Date: Tue, 26 Jun 2018 16:25:21 +0200 Subject: [PATCH 1/4] (Re)added functionality to process datacubes. Includes testing. --- photutils/aperture/core.py | 13 +- photutils/aperture/mask.py | 17 +- .../tests/test_aperture_photometry_3D.py | 931 ++++++++++++++++++ photutils/datasets/make.py | 271 ++++- 4 files changed, 1217 insertions(+), 15 deletions(-) create mode 100644 photutils/aperture/tests/test_aperture_photometry_3D.py diff --git a/photutils/aperture/core.py b/photutils/aperture/core.py index f2c8d06b0..3daf82d6d 100644 --- a/photutils/aperture/core.py +++ b/photutils/aperture/core.py @@ -393,21 +393,22 @@ def do_photometry(self, data, error=None, mask=None, method='exact', aperture_sums = [] aperture_sum_errs = [] + output_shape = (1,) if data.ndim==2 else data.shape[0] for mask in self.to_mask(method=method, subpixels=subpixels): data_cutout = mask.cutout(data) if data_cutout is None: - aperture_sums.append(np.nan) + aperture_sums.append(np.repeat(np.nan, output_shape)) else: - aperture_sums.append(np.sum(data_cutout * mask.data)) + aperture_sums.append(np.sum(data_cutout * mask.data, axis=(-2,-1))) if error is not None: error_cutout = mask.cutout(error) if error_cutout is None: - aperture_sum_errs.append(np.nan) + aperture_sum_errs.append(np.repeat(np.nan, output_shape)) else: - aperture_var = np.sum(error_cutout ** 2 * mask.data) + aperture_var = np.sum(error_cutout ** 2 * mask.data, axis=(-2,-1)) aperture_sum_errs.append(np.sqrt(aperture_var)) # handle Quantity objects and input units @@ -726,8 +727,8 @@ def _prepare_photometry_input(data, error, mask, wcs, unit): pass data = np.asanyarray(data) - if data.ndim != 2: - raise ValueError('data must be a 2D array.') + if (data.ndim != 2) and (data.ndim !=3): + raise ValueError('data must be a 2D or 3D array.') if unit is not None: unit = u.Unit(unit, parse_strict='warn') diff --git a/photutils/aperture/mask.py b/photutils/aperture/mask.py index 9c4b0d976..aee1620bc 100644 --- a/photutils/aperture/mask.py +++ b/photutils/aperture/mask.py @@ -178,8 +178,8 @@ def cutout(self, data, fill_value=0., copy=False): """ data = np.asanyarray(data) - if data.ndim != 2: - raise ValueError('data must be a 2D array.') + if (data.ndim != 2) and (data.ndim != 3): + raise ValueError('data must be a 2D or 3D array.') partial_overlap = False if self.bbox.ixmin < 0 or self.bbox.iymin < 0: @@ -189,20 +189,21 @@ def cutout(self, data, fill_value=0., copy=False): # try this for speed -- the result may still be a partial # overlap, in which case the next block will be triggered if copy: - cutout = np.copy(data[self.bbox.slices]) + cutout = np.copy(data[(...,)+self.bbox.slices]) else: - cutout = data[self.bbox.slices] + cutout = data[(...,)+self.bbox.slices] - if partial_overlap or (cutout.shape != self.shape): - slices_large, slices_small = self._overlap_slices(data.shape) + if partial_overlap or (cutout.shape[-2:] != self.shape): + slices_large, slices_small = self._overlap_slices(data.shape[-2:]) if slices_small is None: return None # no overlap # cutout is a copy - cutout = np.zeros(self.shape, dtype=data.dtype) + output_shape = self.shape if data.ndim==2 else (data.shape[0],)+self.shape + cutout = np.zeros(output_shape, dtype=data.dtype) cutout[:] = fill_value - cutout[slices_small] = data[slices_large] + cutout[(...,)+slices_small] = data[(...,)+slices_large] if isinstance(data, u.Quantity): cutout = u.Quantity(cutout, unit=data.unit) diff --git a/photutils/aperture/tests/test_aperture_photometry_3D.py b/photutils/aperture/tests/test_aperture_photometry_3D.py new file mode 100644 index 000000000..e90ca1028 --- /dev/null +++ b/photutils/aperture/tests/test_aperture_photometry_3D.py @@ -0,0 +1,931 @@ +# Licensed under a 3-clause BSD style license - see LICENSE.rst +""" +The tests in this file test the accuracy of the photometric results. +Here we test directly with aperture objects since we are checking the +algorithms in aperture_photometry, not in the wrappers. +""" + +from __future__ import (absolute_import, division, print_function, + unicode_literals) + +import pytest +import numpy as np +from numpy.testing import (assert_allclose, assert_array_equal, + assert_array_less) + +from astropy.coordinates import SkyCoord +from astropy.io import fits +from astropy.nddata import NDData +from astropy.table import Table +from astropy.tests.helper import remote_data +import astropy.units as u +from astropy.utils import minversion +from astropy.wcs.utils import pixel_to_skycoord + +from ..core import aperture_photometry +from ..circle import (CircularAperture, CircularAnnulus, SkyCircularAperture, + SkyCircularAnnulus) +from ..ellipse import (EllipticalAperture, EllipticalAnnulus, + SkyEllipticalAperture, SkyEllipticalAnnulus) +from ..rectangle import (RectangularAperture, RectangularAnnulus, + SkyRectangularAperture, SkyRectangularAnnulus) +from ...datasets import (get_path, make_4gaussians_image, make_wcs, + make_imagehdu, make_4gaussians_cube, make_cubehdu) + +try: + import matplotlib # noqa + HAS_MATPLOTLIB = True +except ImportError: + HAS_MATPLOTLIB = False + +NUMPY_LT_1_14 = not minversion('numpy', '1.14dev') + +APERTURE_CL = [CircularAperture, + CircularAnnulus, + EllipticalAperture, + EllipticalAnnulus, + RectangularAperture, + RectangularAnnulus] + +TEST_APERTURES = list(zip(APERTURE_CL, ((3.,), (3., 5.), + (3., 5., 1.), (3., 5., 4., 1.), + (5, 8, np.pi / 4), + (8, 12, 8, np.pi / 8)))) + + +# THIS TEST HAS BEEN DONE. +@pytest.mark.parametrize(('aperture_class', 'params'), TEST_APERTURES) +def test_outside_array(aperture_class, params): + data = np.ones((100, 10, 10), dtype=np.float) + aperture = aperture_class((-60, 60), *params) + fluxtable = aperture_photometry(data, aperture) + # aperture is fully outside array: + assert np.all(np.isnan(fluxtable['aperture_sum'])) + + +# THIS TEST HAS BEEN DONE. +@pytest.mark.parametrize(('aperture_class', 'params'), TEST_APERTURES) +def test_inside_array_simple(aperture_class, params): + data = np.ones((100, 40, 40), dtype=np.float) + aperture = aperture_class((20., 20.), *params) + table1 = aperture_photometry(data, aperture, method='center', subpixels=10) + table2 = aperture_photometry(data, aperture, method='subpixel', + subpixels=10) + table3 = aperture_photometry(data, aperture, method='exact', subpixels=10) + true_flux = aperture.area() + + if not isinstance(aperture, (RectangularAperture, RectangularAnnulus)): + assert_allclose(table3['aperture_sum'], true_flux) + assert_allclose(table2['aperture_sum'], table3['aperture_sum'], + atol=0.1) + + assert np.all(table1['aperture_sum'] < table3['aperture_sum']) + + +# Unnecessary. +@pytest.mark.skipif('not HAS_MATPLOTLIB') +@pytest.mark.parametrize(('aperture_class', 'params'), TEST_APERTURES) +def test_aperture_plots(aperture_class, params): + # This test should run without any errors, and there is no return + # value. + # TODO: check the content of the plot + aperture = aperture_class((20., 20.), *params) + aperture.plot() + + +# Unnecessary. +def test_aperture_pixel_positions(): + pos1 = (10, 20) + pos2 = u.Quantity((10, 20), unit=u.pixel) + pos3 = ((10, 20, 30), (10, 20, 30)) + pos3_pairs = ((10, 10), (20, 20), (30, 30)) + + r = 3 + ap1 = CircularAperture(pos1, r) + ap2 = CircularAperture(pos2, r) + ap3 = CircularAperture(pos3, r) + + assert_allclose(np.atleast_2d(pos1), ap1.positions) + assert_allclose(np.atleast_2d(pos2.value), ap2.positions) + assert_allclose(pos3_pairs, ap3.positions) + + +class BaseTestAperturePhotometry(object): + + def test_scalar_error(self): + # Scalar error + error = 1. + if not hasattr(self, 'mask'): + mask = None + true_error = np.sqrt(self.area) + else: + mask = self.mask + # 1 masked pixel + true_error = np.sqrt(self.area - 1) + + table1 = aperture_photometry(self.data, + self.aperture, method='center', + mask=mask, error=error) + table2 = aperture_photometry(self.data, + self.aperture, + method='subpixel', subpixels=12, + mask=mask, error=error) + table3 = aperture_photometry(self.data, + self.aperture, method='exact', + mask=mask, error=error) + + if not isinstance(self.aperture, (RectangularAperture, + RectangularAnnulus)): + [assert_allclose(tab3_ap_sum_Col_i, true_flux_i) + for (tab3_ap_sum_Col_i, true_flux_i) + in zip(np.atleast_2d(table3['aperture_sum']), np.atleast_1d(self.true_flux))] + assert_allclose(table2['aperture_sum'], table3['aperture_sum'], + atol=0.1) + assert np.all(table1['aperture_sum'] < table3['aperture_sum']) + + if not isinstance(self.aperture, (RectangularAperture, + RectangularAnnulus)): + [assert_allclose(tab3_ap_sum_err_Col_i, true_flux_i) + for (tab3_ap_sum_err_Col_i, true_flux_i) + in zip(np.atleast_2d(table3['aperture_sum_err']), np.atleast_1d(true_error))] + assert_allclose(table2['aperture_sum'], table3['aperture_sum'], + atol=0.1) + assert np.all(table1['aperture_sum_err'] < table3['aperture_sum_err']) + + def test_array_error(self): + # Array error + error = np.ones(self.data.shape, dtype=np.float) + if not hasattr(self, 'mask'): + mask = None + true_error = np.sqrt(self.area) + else: + mask = self.mask + # 1 masked pixel + true_error = np.sqrt(self.area - 1) + + table1 = aperture_photometry(self.data, + self.aperture, method='center', + mask=mask, error=error) + table2 = aperture_photometry(self.data, + self.aperture, + method='subpixel', subpixels=12, + mask=mask, error=error) + table3 = aperture_photometry(self.data, + self.aperture, method='exact', + mask=mask, error=error) + + if not isinstance(self.aperture, (RectangularAperture, + RectangularAnnulus)): + [assert_allclose(tab3_ap_sum_Col_i, true_flux_i) + for (tab3_ap_sum_Col_i, true_flux_i) + in zip(np.atleast_2d(table3['aperture_sum']), np.atleast_1d(self.true_flux))] + #assert_allclose(table3['aperture_sum'], self.true_flux) + assert_allclose(table2['aperture_sum'], table3['aperture_sum'], + atol=0.1) + assert np.all(table1['aperture_sum'] < table3['aperture_sum']) + + if not isinstance(self.aperture, (RectangularAperture, + RectangularAnnulus)): + [assert_allclose(tab3_ap_sum_err_Col_i, true_flux_i) + for (tab3_ap_sum_err_Col_i, true_flux_i) + in zip(np.atleast_2d(table3['aperture_sum_err']), np.atleast_1d(true_error))] + #assert_allclose(table3['aperture_sum_err'], true_error) + assert_allclose(table2['aperture_sum_err'], + table3['aperture_sum_err'], atol=0.1) + assert np.all(table1['aperture_sum_err'] < table3['aperture_sum_err']) + + +# THIS TEST HAS BEEN DONE. +class TestCircular(BaseTestAperturePhotometry): + + def setup_class(self): + self.data = np.ones((100, 40, 40), dtype=np.float) + position = (20., 20.) + r = 10. + self.aperture = CircularAperture(position, r) + self.area = np.pi * r * r + self.true_flux = self.area + + +# THIS TEST HAS BEEN DONE. +class TestCircularArray(BaseTestAperturePhotometry): + + def setup_class(self): + self.data = np.ones((100, 40, 40), dtype=np.float) + position = ((20., 20.), (25., 25.)) + r = 10. + self.aperture = CircularAperture(position, r) + self.area = np.pi * r * r + self.area = np.array((self.area, ) * 2) + self.true_flux = self.area + + +# THIS TEST HAS BEEN DONE. +class TestCircularAnnulus(BaseTestAperturePhotometry): + + def setup_class(self): + self.data = np.ones((100, 40, 40), dtype=np.float) + position = (20., 20.) + r_in = 8. + r_out = 10. + self.aperture = CircularAnnulus(position, r_in, r_out) + self.area = np.pi * (r_out * r_out - r_in * r_in) + self.true_flux = self.area + + +# THIS TEST HAS BEEN DONE. +class TestCircularAnnulusArray(BaseTestAperturePhotometry): + + def setup_class(self): + self.data = np.ones((100, 40, 40), dtype=np.float) + position = ((20., 20.), (25., 25.)) + r_in = 8. + r_out = 10. + self.aperture = CircularAnnulus(position, r_in, r_out) + self.area = np.pi * (r_out * r_out - r_in * r_in) + self.area = np.array((self.area, ) * 2) + self.true_flux = self.area + + +# THIS TEST HAS BEEN DONE. +class TestElliptical(BaseTestAperturePhotometry): + + def setup_class(self): + self.data = np.ones((100, 40, 40), dtype=np.float) + position = (20., 20.) + a = 10. + b = 5. + theta = -np.pi / 4. + self.aperture = EllipticalAperture(position, a, b, theta) + self.area = np.pi * a * b + self.true_flux = self.area + + +# THIS TEST HAS BEEN DONE. +class TestEllipticalAnnulus(BaseTestAperturePhotometry): + + def setup_class(self): + self.data = np.ones((100, 40, 40), dtype=np.float) + position = (20., 20.) + a_in = 5. + a_out = 8. + b_out = 5. + theta = -np.pi / 4. + self.aperture = EllipticalAnnulus(position, a_in, a_out, b_out, theta) + self.area = (np.pi * (a_out * b_out) - + np.pi * (a_in * b_out * a_in / a_out)) + self.true_flux = self.area + + +# THIS TEST HAS BEEN DONE. +class TestRectangularAperture(BaseTestAperturePhotometry): + + def setup_class(self): + self.data = np.ones((100, 40, 40), dtype=np.float) + position = (20., 20.) + h = 5. + w = 8. + theta = np.pi / 4. + self.aperture = RectangularAperture(position, w, h, theta) + self.area = h * w + self.true_flux = self.area + + +# THIS TEST HAS BEEN DONE. +class TestRectangularAnnulus(BaseTestAperturePhotometry): + + def setup_class(self): + self.data = np.ones((100, 40, 40), dtype=np.float) + position = (20., 20.) + h_out = 8. + w_in = 8. + w_out = 12. + h_in = w_in * h_out / w_out + theta = np.pi / 8. + self.aperture = RectangularAnnulus(position, w_in, w_out, h_out, theta) + self.area = h_out * w_out - h_in * w_in + self.true_flux = self.area + + +# THIS TEST HAS BEEN DONE. +class TestMaskedSkipCircular(BaseTestAperturePhotometry): + + def setup_class(self): + self.data = np.ones((100, 40, 40), dtype=np.float) + self.mask = np.zeros((100, 40, 40), dtype=bool) + self.mask[:, 20, 20] = True + position = (20., 20.) + r = 10. + self.aperture = CircularAperture(position, r) + self.area = np.pi * r * r + self.true_flux = self.area - 1 + + +class BaseTestDifferentData(object): + + def test_basic_circular_aperture_photometry(self): + aperture = CircularAperture(self.position, self.radius) + table = aperture_photometry(self.data, aperture, + method='exact', unit='adu') + + assert_allclose(table['aperture_sum'].value, self.true_flux) + assert table['aperture_sum'].unit, self.fluxunit + + assert np.all(table['xcenter'].value == + np.transpose(self.position)[0]) + assert np.all(table['ycenter'].value == + np.transpose(self.position)[1]) + + +# THIS TEST HAS BEEN DONE. +class TestInputPrimaryHDU(BaseTestDifferentData): + + def setup_class(self): + data = np.ones((100, 40, 40), dtype=np.float) + self.data = fits.ImageHDU(data=data) + self.data.header['BUNIT'] = 'adu' + self.radius = 3 + self.position = (20, 20) + self.true_flux = np.pi * self.radius * self.radius + self.fluxunit = u.adu + + +# THIS TEST HAS BEEN DONE. +class TestInputHDUList(BaseTestDifferentData): + + def setup_class(self): + data0 = np.ones((100, 40, 40), dtype=np.float) + data1 = np.empty((100, 40, 40), dtype=np.float) + data1.fill(2) + self.data = fits.HDUList([fits.ImageHDU(data=data0), + fits.ImageHDU(data=data1)]) + self.radius = 3 + self.position = (20, 20) + # It should stop at the first extension + self.true_flux = np.pi * self.radius * self.radius + + +# THIS TEST HAS BEEN DONE. +class TestInputHDUDifferentBUNIT(BaseTestDifferentData): + + def setup_class(self): + data = np.ones((100, 40, 40), dtype=np.float) + self.data = fits.ImageHDU(data=data) + self.data.header['BUNIT'] = 'Jy' + self.radius = 3 + self.position = (20, 20) + self.true_flux = np.pi * self.radius * self.radius + self.fluxunit = u.adu + + +# THIS TEST HAS BEEN DONE. +class TestInputNDData(BaseTestDifferentData): + + def setup_class(self): + data = np.ones((100, 40, 40), dtype=np.float) + self.data = NDData(data, unit=u.adu) + self.radius = 3 + self.position = [(20, 20), (30, 30)] + self.true_flux = np.pi * self.radius * self.radius + self.fluxunit = u.adu + + +# Needs updating. +@remote_data +def test_wcs_based_photometry_to_catalogue(): + pathcat = get_path('spitzer_example_catalog.xml', location='remote') + pathhdu = get_path('spitzer_example_image.fits', location='remote') + hdu = fits.open(pathhdu) + scale = hdu[0].header['PIXSCAL1'] + + catalog = Table.read(pathcat) + + pos_skycoord = SkyCoord(catalog['l'], catalog['b'], frame='galactic') + + photometry_skycoord = aperture_photometry( + hdu, SkyCircularAperture(pos_skycoord, 4 * u.arcsec)) + + photometry_skycoord_pix = aperture_photometry( + hdu, SkyCircularAperture(pos_skycoord, 4. / scale * u.pixel)) + + assert_allclose(photometry_skycoord['aperture_sum'], + photometry_skycoord_pix['aperture_sum']) + + # Photometric unit conversion is needed to match the catalogue + factor = (1.2 * u.arcsec) ** 2 / u.pixel + converted_aperture_sum = (photometry_skycoord['aperture_sum'] * + factor).to(u.mJy / u.pixel) + + fluxes_catalog = catalog['f4_5'].filled() + + # There shouldn't be large outliers, but some differences is OK, as + # fluxes_catalog is based on PSF photometry, etc. + assert_allclose(fluxes_catalog, converted_aperture_sum.value, rtol=1e0) + + assert(np.mean(np.fabs(((fluxes_catalog - converted_aperture_sum.value) / + fluxes_catalog))) < 0.1) + + +# THIS TEST HAS BEEN DONE. +def test_wcs_based_photometry(): + data = make_4gaussians_cube() + wcs = make_wcs(data.shape) + hdu = make_cubehdu(data, wcs=wcs) + + # hard wired positions in make_4gaussians_image + pos_orig_pixel = u.Quantity(([160., 25., 150., 90.], + [70., 40., 25., 60.]), unit=u.pixel) + + pos_skycoord = pixel_to_skycoord(pos_orig_pixel[0], pos_orig_pixel[1], wcs) + + pos_skycoord_s = pos_skycoord[2] + + photometry_skycoord_circ = aperture_photometry( + hdu, SkyCircularAperture(pos_skycoord, 3 * u.arcsec)) + photometry_skycoord_circ_2 = aperture_photometry( + hdu, SkyCircularAperture(pos_skycoord, 2 * u.arcsec)) + photometry_skycoord_circ_s = aperture_photometry( + hdu, SkyCircularAperture(pos_skycoord_s, 3 * u.arcsec)) + + assert_allclose(np.atleast_2d(photometry_skycoord_circ['aperture_sum'][2]), + np.atleast_2d(photometry_skycoord_circ_s['aperture_sum'])) + + photometry_skycoord_circ_ann = aperture_photometry( + hdu, SkyCircularAnnulus(pos_skycoord, 2 * u.arcsec, 3 * u.arcsec)) + photometry_skycoord_circ_ann_s = aperture_photometry( + hdu, SkyCircularAnnulus(pos_skycoord_s, 2 * u.arcsec, 3 * u.arcsec)) + + assert_allclose(np.atleast_2d(photometry_skycoord_circ_ann['aperture_sum'][2]), + np.atleast_2d(photometry_skycoord_circ_ann_s['aperture_sum'])) + + assert_allclose(photometry_skycoord_circ_ann['aperture_sum'], + photometry_skycoord_circ['aperture_sum'] - + photometry_skycoord_circ_2['aperture_sum']) + + photometry_skycoord_ell = aperture_photometry( + hdu, SkyEllipticalAperture(pos_skycoord, 3 * u.arcsec, + 3.0001 * u.arcsec, 45 * u.arcsec)) + photometry_skycoord_ell_2 = aperture_photometry( + hdu, SkyEllipticalAperture(pos_skycoord, 2 * u.arcsec, + 2.0001 * u.arcsec, 45 * u.arcsec)) + photometry_skycoord_ell_s = aperture_photometry( + hdu, SkyEllipticalAperture(pos_skycoord_s, 3 * u.arcsec, + 3.0001 * u.arcsec, 45 * u.arcsec)) + photometry_skycoord_ell_ann = aperture_photometry( + hdu, SkyEllipticalAnnulus(pos_skycoord, 2 * u.arcsec, 3 * u.arcsec, + 3.0001 * u.arcsec, 45 * u.arcsec)) + photometry_skycoord_ell_ann_s = aperture_photometry( + hdu, SkyEllipticalAnnulus(pos_skycoord_s, 2 * u.arcsec, 3 * u.arcsec, + 3.0001 * u.arcsec, 45 * u.arcsec)) + + assert_allclose(np.atleast_2d(photometry_skycoord_ell['aperture_sum'][2]), + np.atleast_2d(photometry_skycoord_ell_s['aperture_sum'])) + + assert_allclose(np.atleast_2d(photometry_skycoord_ell_ann['aperture_sum'][2]), + np.atleast_2d(photometry_skycoord_ell_ann_s['aperture_sum'])) + + assert_allclose(photometry_skycoord_ell['aperture_sum'], + photometry_skycoord_circ['aperture_sum'], rtol=5e-3) + + assert_allclose(photometry_skycoord_ell_ann['aperture_sum'], + photometry_skycoord_ell['aperture_sum'] - + photometry_skycoord_ell_2['aperture_sum'], rtol=1e-4) + + photometry_skycoord_rec = aperture_photometry( + hdu, SkyRectangularAperture(pos_skycoord, + 6 * u.arcsec, 6 * u.arcsec, + 0 * u.arcsec), + method='subpixel', subpixels=20) + photometry_skycoord_rec_4 = aperture_photometry( + hdu, SkyRectangularAperture(pos_skycoord, + 4 * u.arcsec, 4 * u.arcsec, + 0 * u.arcsec), + method='subpixel', subpixels=20) + photometry_skycoord_rec_s = aperture_photometry( + hdu, SkyRectangularAperture(pos_skycoord_s, + 6 * u.arcsec, 6 * u.arcsec, + 0 * u.arcsec), + method='subpixel', subpixels=20) + photometry_skycoord_rec_ann = aperture_photometry( + hdu, SkyRectangularAnnulus(pos_skycoord, 4 * u.arcsec, 6 * u.arcsec, + 6 * u.arcsec, 0 * u.arcsec), + method='subpixel', subpixels=20) + photometry_skycoord_rec_ann_s = aperture_photometry( + hdu, SkyRectangularAnnulus(pos_skycoord_s, 4 * u.arcsec, 6 * u.arcsec, + 6 * u.arcsec, 0 * u.arcsec), + method='subpixel', subpixels=20) + + assert_allclose(np.atleast_2d(photometry_skycoord_rec['aperture_sum'][2]), + np.atleast_2d(photometry_skycoord_rec_s['aperture_sum'])) + + assert np.all(photometry_skycoord_rec['aperture_sum'] > + photometry_skycoord_circ['aperture_sum']) + + assert_allclose(np.atleast_2d(photometry_skycoord_rec_ann['aperture_sum'][2]), + np.atleast_2d(photometry_skycoord_rec_ann_s['aperture_sum'])) + + assert_allclose(photometry_skycoord_rec_ann['aperture_sum'], + photometry_skycoord_rec['aperture_sum'] - + photometry_skycoord_rec_4['aperture_sum'], rtol=1e-4) + + +# THIS TEST HAS BEEN DONE. +def test_basic_circular_aperture_photometry_unit(): + data1 = np.ones((100, 40, 40), dtype=np.float) + data2 = u.Quantity(data1, unit=u.adu) + + radius = 3 + position = (20, 20) + true_flux = np.pi * radius * radius + unit = u.adu + + table1 = aperture_photometry(data1, CircularAperture(position, radius), + unit=unit) + table2 = aperture_photometry(data2, CircularAperture(position, radius)) + + assert_allclose(table1['aperture_sum'].value, true_flux) + assert_allclose(table2['aperture_sum'].value, true_flux) + assert table1['aperture_sum'].unit == unit + assert table2['aperture_sum'].unit == data2.unit == unit + + +# THIS TEST HAS BEEN DONE. +def test_aperture_photometry_with_error_units(): + """Test aperture_photometry when error has units (see #176).""" + + data1 = np.ones((100, 40, 40), dtype=np.float) + data2 = u.Quantity(data1, unit=u.adu) + error = u.Quantity(data1, unit=u.adu) + radius = 3 + true_flux = np.pi * radius * radius + unit = u.adu + position = (20, 20) + table1 = aperture_photometry(data2, CircularAperture(position, radius), + error=error) + assert_allclose(table1['aperture_sum'].value, true_flux) + assert_allclose(table1['aperture_sum_err'].value, np.sqrt(true_flux)) + assert table1['aperture_sum'].unit == unit + assert table1['aperture_sum_err'].unit == unit + + +# THIS TEST HAS BEEN DONE. +def test_aperture_photometry_inputs_with_mask(): + """ + Test that aperture_photometry does not modify the input + data or error array when a mask is input. + """ + + data = np.ones((100, 5, 5)) + aperture = CircularAperture((2, 2), 2.) + mask = np.zeros_like(data, dtype=bool) + data[:, 2, 2] = 100. # bad pixel + mask[:, 2, 2] = True + error = np.sqrt(data) + data_in = data.copy() + error_in = error.copy() + t1 = aperture_photometry(data, aperture, error=error, mask=mask) + assert_array_equal(data, data_in) + assert_array_equal(error, error_in) + assert_allclose(t1['aperture_sum'][0], 11.5663706144) + t2 = aperture_photometry(data, aperture) + assert_allclose(t2['aperture_sum'][0], 111.566370614) + + +TEST_ELLIPSE_EXACT_APERTURES = [(3.469906, 3.923861394, 3.), + (0.3834415188257778, 0.3834415188257778, 0.3)] + + +# THIS TEST HAS BEEN DONE. +# Unnecessary. +@pytest.mark.parametrize('x,y,r', TEST_ELLIPSE_EXACT_APERTURES) +def test_ellipse_exact_grid(x, y, r): + """ + Test elliptical exact aperture photometry on a grid of pixel positions. + + This is a regression test for the bug discovered in this issue: + https://github.com/astropy/photutils/issues/198 + """ + + data = np.ones((100, 10, 10)) + + aperture = EllipticalAperture((x, y), r, r, 0.) + t = aperture_photometry(data, aperture, method='exact') + actual = t['aperture_sum'][0] / (np.pi * r ** 2) + assert_allclose(actual, 1) + + +# THIS TEST HAS BEEN DONE. +@pytest.mark.parametrize('value', [np.nan, np.inf]) +def test_nan_inf_mask(value): + """Test that nans and infs are properly masked [267].""" + + data = np.ones((100, 9, 9)) + mask = np.zeros_like(data, dtype=bool) + data[:, 4, 4] = value + mask[:, 4, 4] = True + radius = 2. + aper = CircularAperture((4, 4), radius) + tbl = aperture_photometry(data, aper, mask=mask) + desired = (np.pi * radius**2) - 1 + assert_allclose(tbl['aperture_sum'], desired) + + +# THIS TEST HAS BEEN DONE. +def test_aperture_partial_overlap(): + data = np.ones((100, 20, 20)) + error = np.ones((100, 20, 20)) + xypos = [(10, 10), (0, 0), (0, 19), (19, 0), (19, 19)] + r = 5. + aper = CircularAperture(xypos, r=r) + tbl = aperture_photometry(data, aper, error=error) + assert_allclose(tbl['aperture_sum'][0], np.pi * r ** 2) + assert_array_less(tbl['aperture_sum'][1:], np.pi * r ** 2) + + unit = u.MJy / u.sr + tbl = aperture_photometry(data * unit, aper, error=error * unit) + assert_allclose(tbl['aperture_sum'][0].value, np.pi * r ** 2) + assert_array_less(tbl['aperture_sum'][1:].value, np.pi * r ** 2) + assert_array_less(tbl['aperture_sum_err'][1:].value, np.pi * r ** 2) + assert tbl['aperture_sum'].unit == unit + assert tbl['aperture_sum_err'].unit == unit + + +# THIS TEST HAS BEEN DONE. +# Unnecessary. +def test_pixel_aperture_repr(): + aper = CircularAperture((10, 20), r=3.0) + a_repr = '' + a_str = 'Aperture: CircularAperture\npositions: [[10, 20]]\nr: 3.0' + assert repr(aper) == a_repr + assert str(aper) == a_str + + aper = CircularAnnulus((10, 20), r_in=3.0, r_out=5.0) + a_repr = '' + a_str = ('Aperture: CircularAnnulus\npositions: [[10, 20]]\nr_in: 3.0\n' + 'r_out: 5.0') + assert repr(aper) == a_repr + assert str(aper) == a_str + + aper = EllipticalAperture((10, 20), a=5.0, b=3.0, theta=15.0) + a_repr = '' + a_str = ('Aperture: EllipticalAperture\npositions: [[10, 20]]\n' + 'a: 5.0\nb: 3.0\ntheta: 15.0') + assert repr(aper) == a_repr + assert str(aper) == a_str + + aper = EllipticalAnnulus((10, 20), a_in=4.0, a_out=8.0, b_out=4.0, + theta=15.0) + a_repr = ('') + a_str = ('Aperture: EllipticalAnnulus\npositions: [[10, 20]]\na_in: ' + '4.0\na_out: 8.0\nb_out: 4.0\ntheta: 15.0') + assert repr(aper) == a_repr + assert str(aper) == a_str + + aper = RectangularAperture((10, 20), w=5.0, h=3.0, theta=15.0) + a_repr = '' + a_str = ('Aperture: RectangularAperture\npositions: [[10, 20]]\n' + 'w: 5.0\nh: 3.0\ntheta: 15.0') + assert repr(aper) == a_repr + assert str(aper) == a_str + + aper = RectangularAnnulus((10, 20), w_in=4.0, w_out=8.0, h_out=4.0, + theta=15.0) + a_repr = ('') + a_str = ('Aperture: RectangularAnnulus\npositions: [[10, 20]]\n' + 'w_in: 4.0\nw_out: 8.0\nh_out: 4.0\ntheta: 15.0') + assert repr(aper) == a_repr + assert str(aper) == a_str + + +# THIS TEST HAS BEEN DONE. +# Unnecessary. +def test_sky_aperture_repr(): + s = SkyCoord([1, 2], [3, 4], unit='deg') + + aper = SkyCircularAperture(s, r=3*u.pix) + if NUMPY_LT_1_14: + a_repr = (', r=3.0 pix)>') + a_str = ('Aperture: SkyCircularAperture\npositions: \n' + 'r: 3.0 pix') + else: + a_repr = (', r=3.0 pix)>') + a_str = ('Aperture: SkyCircularAperture\npositions: \n' + 'r: 3.0 pix') + + assert repr(aper) == a_repr + assert str(aper) == a_str + + aper = SkyCircularAnnulus(s, r_in=3.*u.pix, r_out=5*u.pix) + if NUMPY_LT_1_14: + a_repr = (', r_in=3.0 pix, ' + 'r_out=5.0 pix)>') + a_str = ('Aperture: SkyCircularAnnulus\npositions: \n' + 'r_in: 3.0 pix\nr_out: 5.0 pix') + else: + a_repr = (', r_in=3.0 pix, r_out=5.0 pix)>') + a_str = ('Aperture: SkyCircularAnnulus\npositions: \n' + 'r_in: 3.0 pix\nr_out: 5.0 pix') + + assert repr(aper) == a_repr + assert str(aper) == a_str + + aper = SkyEllipticalAperture(s, a=3*u.pix, b=5*u.pix, theta=15*u.deg) + if NUMPY_LT_1_14: + a_repr = (', a=3.0 pix, b=5.0 pix,' + ' theta=15.0 deg)>') + a_str = ('Aperture: SkyEllipticalAperture\npositions: \n' + 'a: 3.0 pix\nb: 5.0 pix\ntheta: 15.0 deg') + else: + a_repr = (', a=3.0 pix, b=5.0 pix,' + ' theta=15.0 deg)>') + a_str = ('Aperture: SkyEllipticalAperture\npositions: \n' + 'a: 3.0 pix\nb: 5.0 pix\ntheta: 15.0 deg') + + assert repr(aper) == a_repr + assert str(aper) == a_str + + aper = SkyEllipticalAnnulus(s, a_in=3*u.pix, a_out=5*u.pix, b_out=3*u.pix, + theta=15*u.deg) + if NUMPY_LT_1_14: + a_repr = (', a_in=3.0 pix, ' + 'a_out=5.0 pix, b_out=3.0 pix, theta=15.0 deg)>') + a_str = ('Aperture: SkyEllipticalAnnulus\npositions: \n' + 'a_in: 3.0 pix\na_out: 5.0 pix\nb_out: 3.0 pix\n' + 'theta: 15.0 deg') + else: + a_repr = (', a_in=3.0 pix, ' + 'a_out=5.0 pix, b_out=3.0 pix, theta=15.0 deg)>') + a_str = ('Aperture: SkyEllipticalAnnulus\npositions: \n' + 'a_in: 3.0 pix\na_out: 5.0 pix\nb_out: 3.0 pix\n' + 'theta: 15.0 deg') + + assert repr(aper) == a_repr + assert str(aper) == a_str + + aper = SkyRectangularAperture(s, w=3*u.pix, h=5*u.pix, theta=15*u.deg) + if NUMPY_LT_1_14: + a_repr = (', w=3.0 pix, h=5.0 pix' + ', theta=15.0 deg)>') + a_str = ('Aperture: SkyRectangularAperture\npositions: \n' + 'w: 3.0 pix\nh: 5.0 pix\ntheta: 15.0 deg') + else: + a_repr = (', w=3.0 pix, h=5.0 pix' + ', theta=15.0 deg)>') + a_str = ('Aperture: SkyRectangularAperture\npositions: \n' + 'w: 3.0 pix\nh: 5.0 pix\ntheta: 15.0 deg') + + assert repr(aper) == a_repr + assert str(aper) == a_str + + aper = SkyRectangularAnnulus(s, w_in=3*u.pix, w_out=3.4*u.pix, + h_out=5*u.pix, theta=15*u.deg) + if NUMPY_LT_1_14: + a_repr = (', w_in=3.0 pix, ' + 'w_out=3.4 pix, h_out=5.0 pix, theta=15.0 deg)>') + a_str = ('Aperture: SkyRectangularAnnulus\npositions: \n' + 'w_in: 3.0 pix\nw_out: 3.4 pix\nh_out: 5.0 pix\n' + 'theta: 15.0 deg') + else: + a_repr = (', w_in=3.0 pix, ' + 'w_out=3.4 pix, h_out=5.0 pix, theta=15.0 deg)>') + a_str = ('Aperture: SkyRectangularAnnulus\npositions: \n' + 'w_in: 3.0 pix\nw_out: 3.4 pix\nh_out: 5.0 pix\n' + 'theta: 15.0 deg') + + assert repr(aper) == a_repr + assert str(aper) == a_str + + +# THIS TEST HAS BEEN DONE. +# Unnecessary. +def test_rectangular_bbox(): + # odd sizes + width = 7 + height = 3 + a = RectangularAperture((50, 50), w=width, h=height, theta=0) + assert a.bounding_boxes[0].shape == (height, width) + + a = RectangularAperture((50.5, 50.5), w=width, h=height, theta=0) + assert a.bounding_boxes[0].shape == (height + 1, width + 1) + + a = RectangularAperture((50, 50), w=width, h=height, theta=90.*np.pi/180.) + assert a.bounding_boxes[0].shape == (width, height) + + # even sizes + width = 8 + height = 4 + a = RectangularAperture((50, 50), w=width, h=height, theta=0) + assert a.bounding_boxes[0].shape == (height + 1, width + 1) + + a = RectangularAperture((50.5, 50.5), w=width, h=height, theta=0) + assert a.bounding_boxes[0].shape == (height, width) + + a = RectangularAperture((50.5, 50.5), w=width, h=height, + theta=90.*np.pi/180.) + assert a.bounding_boxes[0].shape == (width, height) + + +# THIS TEST HAS BEEN DONE. +# Unnecessary. +def test_elliptical_bbox(): + # integer axes + a = 7 + b = 3 + ap = EllipticalAperture((50, 50), a=a, b=b, theta=0) + assert ap.bounding_boxes[0].shape == (2*b + 1, 2*a + 1) + + ap = EllipticalAperture((50.5, 50.5), a=a, b=b, theta=0) + assert ap.bounding_boxes[0].shape == (2*b, 2*a) + + ap = EllipticalAperture((50, 50), a=a, b=b, theta=90.*np.pi/180.) + assert ap.bounding_boxes[0].shape == (2*a + 1, 2*b + 1) + + # fractional axes + a = 7.5 + b = 4.5 + ap = EllipticalAperture((50, 50), a=a, b=b, theta=0) + assert ap.bounding_boxes[0].shape == (2*b, 2*a) + + ap = EllipticalAperture((50.5, 50.5), a=a, b=b, theta=0) + assert ap.bounding_boxes[0].shape == (2*b + 1, 2*a + 1) + + ap = EllipticalAperture((50, 50), a=a, b=b, theta=90.*np.pi/180.) + assert ap.bounding_boxes[0].shape == (2*a, 2*b) + + +# THIS TEST HAS BEEN DONE. +def test_to_sky_pixel(): + data = make_4gaussians_cube() + wcs = make_wcs(data.shape) + + ap = CircularAperture(((12.3, 15.7), (48.19, 98.14)), r=3.14) + ap2 = ap.to_sky(wcs).to_pixel(wcs) + assert_allclose(ap.positions, ap2.positions) + assert_allclose(ap.r, ap2.r) + + ap = CircularAnnulus(((12.3, 15.7), (48.19, 98.14)), r_in=3.14, + r_out=5.32) + ap2 = ap.to_sky(wcs).to_pixel(wcs) + assert_allclose(ap.positions, ap2.positions) + assert_allclose(ap.r_in, ap2.r_in) + assert_allclose(ap.r_out, ap2.r_out) + + ap = EllipticalAperture(((12.3, 15.7), (48.19, 98.14)), a=3.14, b=5.32, + theta=103.*np.pi/180.) + ap2 = ap.to_sky(wcs).to_pixel(wcs) + assert_allclose(ap.positions, ap2.positions) + assert_allclose(ap.a, ap2.a) + assert_allclose(ap.b, ap2.b) + assert_allclose(ap.theta, ap2.theta) + + ap = EllipticalAnnulus(((12.3, 15.7), (48.19, 98.14)), a_in=3.14, + a_out=15.32, b_out=4.89, theta=103.*np.pi/180.) + ap2 = ap.to_sky(wcs).to_pixel(wcs) + assert_allclose(ap.positions, ap2.positions) + assert_allclose(ap.a_in, ap2.a_in) + assert_allclose(ap.a_out, ap2.a_out) + assert_allclose(ap.b_out, ap2.b_out) + assert_allclose(ap.theta, ap2.theta) + + ap = RectangularAperture(((12.3, 15.7), (48.19, 98.14)), w=3.14, h=5.32, + theta=103.*np.pi/180.) + ap2 = ap.to_sky(wcs).to_pixel(wcs) + assert_allclose(ap.positions, ap2.positions) + assert_allclose(ap.w, ap2.w) + assert_allclose(ap.h, ap2.h) + assert_allclose(ap.theta, ap2.theta) + + ap = RectangularAnnulus(((12.3, 15.7), (48.19, 98.14)), w_in=3.14, + w_out=15.32, h_out=4.89, theta=103.*np.pi/180.) + ap2 = ap.to_sky(wcs).to_pixel(wcs) + assert_allclose(ap.positions, ap2.positions) + assert_allclose(ap.w_in, ap2.w_in) + assert_allclose(ap.w_out, ap2.w_out) + assert_allclose(ap.h_out, ap2.h_out) + assert_allclose(ap.theta, ap2.theta) diff --git a/photutils/datasets/make.py b/photutils/datasets/make.py index a75c271b9..337ddc7fb 100644 --- a/photutils/datasets/make.py +++ b/photutils/datasets/make.py @@ -21,7 +21,9 @@ 'make_random_models_table', 'make_random_gaussians_table', 'make_model_sources_image', 'make_gaussian_sources_image', 'make_4gaussians_image', 'make_100gaussians_image', - 'make_wcs', 'make_imagehdu'] + 'make_4gaussians_cube', 'make_100gaussians_cube', + 'make_wcs', 'make_imagehdu', 'make_cubehdu', + ] def apply_poisson_noise(data, random_state=None): @@ -556,6 +558,7 @@ def make_gaussian_sources_image(shape, source_table, oversample=1): oversample=oversample) + def make_4gaussians_image(noise=True): """ Make an example image containing four 2D Gaussians plus a constant @@ -670,7 +673,168 @@ def make_100gaussians_image(noise=True): return data + +def make_4gaussians_cube(noise=True): + """ + Make an example datacube containing four 2D Gaussians plus a constant + background. + + The background has a mean of 5. + + If ``noise`` is `True`, then Gaussian noise with a mean of 0 and a + standard deviation of 5 is added to the output image. + + Parameters + ---------- + noise : bool, optional + Whether to include noise in the output image (default is + `True`). + + Returns + ------- + datacube : 3D `~numpy.ndarray` + Image containing four 2D Gaussian sources. + + See Also + -------- + make_100gaussians_cube + + Examples + -------- + .. plot:: + :include-source: + + from photutils import datasets + image = datasets.make_4gaussians_image() + plt.imshow(image, origin='lower', interpolation='nearest') + """ + + table = Table() + table['amplitude'] = [50, 70, 150, 210] + table['x_mean'] = [160, 25, 150, 90] + table['y_mean'] = [70, 40, 25, 60] + table['x_stddev'] = [15.2, 5.1, 3., 8.1] + table['y_stddev'] = [2.6, 2.5, 3., 4.7] + table['theta'] = np.array([145., 20., 0., 60.]) * np.pi / 180. + + shape = (1000, 100, 200) + data = np.empty(shape) + data[:, :, :] = make_gaussian_sources_image(shape[-2:], table)[None, :, :] + 5. + + if noise: + data += make_noise_image(shape[-2:], type='gaussian', mean=0., + stddev=5., random_state=12345)[None, :, :] + + return data + + +def make_100gaussians_cube(noise=True): + """ + Make an example datacube containing 100 2D Gaussians plus a constant + background. + + The background has a mean of 5. + + If ``noise`` is `True`, then Gaussian noise with a mean of 0 and a + standard deviation of 2 is added to the output image. + + Parameters + ---------- + noise : bool, optional + Whether to include noise in the output image (default is + `True`). + + Returns + ------- + datacube : 3D `~numpy.ndarray` + Datacube containing 100 2D Gaussian sources. + + See Also + -------- + make_4gaussians_cube + + Examples + -------- + .. plot:: + :include-source: + + from photutils import datasets + image = datasets.make_100gaussians_image() + plt.imshow(image, origin='lower', interpolation='nearest') + """ + + n_sources = 100 + flux_range = [500, 1000] + xmean_range = [0, 500] + ymean_range = [0, 300] + xstddev_range = [1, 5] + ystddev_range = [1, 5] + params = OrderedDict([('flux', flux_range), + ('x_mean', xmean_range), + ('y_mean', ymean_range), + ('x_stddev', xstddev_range), + ('y_stddev', ystddev_range), + ('theta', [0, 2*np.pi])]) + + sources = make_random_gaussians_table(n_sources, params, + random_state=12345) + + shape = (1000, 300, 500) + data = np.empty(shape) + data[:, :, :] = make_gaussian_sources_image(shape[-2:], table)[None, :, :] + 5. + + if noise: + data += make_noise_image(shape[-2:], type='gaussian', mean=0., + stddev=5., random_state=12345)[None, :, :] + + return data + + def make_wcs(shape, galactic=False): + """ + Create a simple celestial WCS object in either the ICRS or Galactic + coordinate frame. A wrapper to both `_make_wcs_2D` and `_make_wcs_3D`. + + Parameters + ---------- + shape : 2-tuple of int + The shape of the 2D array to be used with the output + `~astropy.wcs.WCS` object. + + galactic : bool, optional + If `True`, then the output WCS will be in the Galactic + coordinate frame. If `False` (default), then the output WCS + will be in the ICRS coordinate frame. + + Returns + ------- + wcs : `~astropy.wcs.WCS` object + The world coordinate system (WCS) transformation. + + See Also + -------- + make_imagehdu, make_cubehdu + + Examples + -------- + >>> from photutils.datasets import make_wcs + >>> shape = (100, 100) + >>> wcs = make_wcs(shape) + >>> print(wcs.wcs.crpix) # doctest: +FLOAT_CMP + [50. 50.] + >>> print(wcs.wcs.crval) # doctest: +FLOAT_CMP + [197.8925 -1.36555556] + """ + + if len(shape) == 2: + return _make_wcs_2D(shape, galactic=galactic) + elif len(shape) == 3: + return _make_wcs_3D(shape, galactic=galactic) + else: + raise ValueError('`shape` must correspond to a 2D or 3D array') + + +def _make_wcs_2D(shape, galactic=False): """ Create a simple celestial WCS object in either the ICRS or Galactic coordinate frame. @@ -725,6 +889,64 @@ def make_wcs(shape, galactic=False): return wcs +def _make_wcs_3D(shape, galactic=False): + """ + Create a simple celestial WCS object in either the ICRS or Galactic + coordinate frame. + + Parameters + ---------- + shape : 3-tuple of int + The shape of the 3D array to be used with the output + `~astropy.wcs.WCS` object. + + galactic : bool, optional + If `True`, then the output WCS will be in the Galactic + coordinate frame. If `False` (default), then the output WCS + will be in the ICRS coordinate frame. + + Returns + ------- + wcs : `~astropy.wcs.WCS` object + The world coordinate system (WCS) transformation. + + See Also + -------- + make_imagehdu + + Examples + -------- + >>> from photutils.datasets import make_wcs + >>> shape = (1000, 100, 100) + >>> wcs = make_wcs(shape) + >>> print(wcs.wcs.crpix) # doctest: +FLOAT_CMP + [50. 50. 500.] + >>> print(wcs.wcs.crval) # doctest: +FLOAT_CMP + [197.8925 -1.36555556 4861.] + """ + + wcs = WCS(naxis=3) + rho = np.pi / 3. + scale = 0.1 / 3600. + delta_wave = 0.5 + wcs._naxis1 = shape[2] # nx + wcs._naxis2 = shape[1] # ny + wcs._naxis3 = shape[0] # nz + wcs.wcs.crpix = [shape[2] / 2, shape[1] / 2, shape[0] / 2] # 1-indexed (x, y, z) + wcs.wcs.crval = [197.8925, -1.36555556, 4861.] + wcs.wcs.cunit = ['deg', 'deg', 'nm'] + wcs.wcs.cd = [[-scale * np.cos(rho), scale * np.sin(rho), 0.], + [ scale * np.sin(rho), scale * np.cos(rho), 0.], + [ 0., 0., delta_wave]] + if not galactic: + wcs.wcs.radesys = 'ICRS' + wcs.wcs.ctype = ['RA---TAN', 'DEC--TAN', 'WAVE'] + else: + wcs.wcs.ctype = ['GLON-CAR', 'GLAT-CAR', 'WAVE'] + + return wcs + + def make_imagehdu(data, wcs=None): """ Create a FITS `~astropy.io.fits.ImageHDU` containing the input 2D @@ -769,3 +991,50 @@ def make_imagehdu(data, wcs=None): header = None return fits.ImageHDU(data, header=header) + + + +def make_cubehdu(data, wcs=None): + """ + Create a FITS `~astropy.io.fits.ImageHDU` containing the input 3D + datacube. + + Parameters + ---------- + data : 3D array-like + The input 3D data. + + wcs : `~astropy.wcs.WCS`, optional + The world coordinate system (WCS) transformation to include in + the output FITS header. + + Returns + ------- + image_hdu : `~astropy.io.fits.ImageHDU` + The FITS `~astropy.io.fits.ImageHDU`. + + See Also + -------- + make_wcs + + Examples + -------- + >>> from photutils.datasets import make_cubehdu, make_wcs + >>> shape = (1000, 100, 100) + >>> data = np.ones(shape) + >>> wcs = make_wcs(shape) + >>> hdu = make_cubehdu(data, wcs=wcs) + >>> print(hdu.data.shape) + (1000, 100, 100) + """ + + data = np.asanyarray(data) + if data.ndim != 3: + raise ValueError('data must be a 3D array') + + if wcs is not None: + header = wcs.to_header() + else: + header = None + + return fits.ImageHDU(data, header=header) From ef539f587c0edac8c6da5accbf45e31168ec7c2f Mon Sep 17 00:00:00 2001 From: Francesco D'Eugenio Date: Tue, 3 Jul 2018 09:40:41 +0200 Subject: [PATCH 2/4] (Re)added functionality to process datacubes. Includes testing. Rebased on most recent master --- photutils/aperture/core.py | 13 +- photutils/aperture/mask.py | 17 +- .../tests/test_aperture_photometry_3D.py | 931 ++++++++++++++++++ photutils/datasets/make.py | 271 ++++- 4 files changed, 1217 insertions(+), 15 deletions(-) create mode 100644 photutils/aperture/tests/test_aperture_photometry_3D.py diff --git a/photutils/aperture/core.py b/photutils/aperture/core.py index f2c8d06b0..3daf82d6d 100644 --- a/photutils/aperture/core.py +++ b/photutils/aperture/core.py @@ -393,21 +393,22 @@ def do_photometry(self, data, error=None, mask=None, method='exact', aperture_sums = [] aperture_sum_errs = [] + output_shape = (1,) if data.ndim==2 else data.shape[0] for mask in self.to_mask(method=method, subpixels=subpixels): data_cutout = mask.cutout(data) if data_cutout is None: - aperture_sums.append(np.nan) + aperture_sums.append(np.repeat(np.nan, output_shape)) else: - aperture_sums.append(np.sum(data_cutout * mask.data)) + aperture_sums.append(np.sum(data_cutout * mask.data, axis=(-2,-1))) if error is not None: error_cutout = mask.cutout(error) if error_cutout is None: - aperture_sum_errs.append(np.nan) + aperture_sum_errs.append(np.repeat(np.nan, output_shape)) else: - aperture_var = np.sum(error_cutout ** 2 * mask.data) + aperture_var = np.sum(error_cutout ** 2 * mask.data, axis=(-2,-1)) aperture_sum_errs.append(np.sqrt(aperture_var)) # handle Quantity objects and input units @@ -726,8 +727,8 @@ def _prepare_photometry_input(data, error, mask, wcs, unit): pass data = np.asanyarray(data) - if data.ndim != 2: - raise ValueError('data must be a 2D array.') + if (data.ndim != 2) and (data.ndim !=3): + raise ValueError('data must be a 2D or 3D array.') if unit is not None: unit = u.Unit(unit, parse_strict='warn') diff --git a/photutils/aperture/mask.py b/photutils/aperture/mask.py index 9c4b0d976..aee1620bc 100644 --- a/photutils/aperture/mask.py +++ b/photutils/aperture/mask.py @@ -178,8 +178,8 @@ def cutout(self, data, fill_value=0., copy=False): """ data = np.asanyarray(data) - if data.ndim != 2: - raise ValueError('data must be a 2D array.') + if (data.ndim != 2) and (data.ndim != 3): + raise ValueError('data must be a 2D or 3D array.') partial_overlap = False if self.bbox.ixmin < 0 or self.bbox.iymin < 0: @@ -189,20 +189,21 @@ def cutout(self, data, fill_value=0., copy=False): # try this for speed -- the result may still be a partial # overlap, in which case the next block will be triggered if copy: - cutout = np.copy(data[self.bbox.slices]) + cutout = np.copy(data[(...,)+self.bbox.slices]) else: - cutout = data[self.bbox.slices] + cutout = data[(...,)+self.bbox.slices] - if partial_overlap or (cutout.shape != self.shape): - slices_large, slices_small = self._overlap_slices(data.shape) + if partial_overlap or (cutout.shape[-2:] != self.shape): + slices_large, slices_small = self._overlap_slices(data.shape[-2:]) if slices_small is None: return None # no overlap # cutout is a copy - cutout = np.zeros(self.shape, dtype=data.dtype) + output_shape = self.shape if data.ndim==2 else (data.shape[0],)+self.shape + cutout = np.zeros(output_shape, dtype=data.dtype) cutout[:] = fill_value - cutout[slices_small] = data[slices_large] + cutout[(...,)+slices_small] = data[(...,)+slices_large] if isinstance(data, u.Quantity): cutout = u.Quantity(cutout, unit=data.unit) diff --git a/photutils/aperture/tests/test_aperture_photometry_3D.py b/photutils/aperture/tests/test_aperture_photometry_3D.py new file mode 100644 index 000000000..e90ca1028 --- /dev/null +++ b/photutils/aperture/tests/test_aperture_photometry_3D.py @@ -0,0 +1,931 @@ +# Licensed under a 3-clause BSD style license - see LICENSE.rst +""" +The tests in this file test the accuracy of the photometric results. +Here we test directly with aperture objects since we are checking the +algorithms in aperture_photometry, not in the wrappers. +""" + +from __future__ import (absolute_import, division, print_function, + unicode_literals) + +import pytest +import numpy as np +from numpy.testing import (assert_allclose, assert_array_equal, + assert_array_less) + +from astropy.coordinates import SkyCoord +from astropy.io import fits +from astropy.nddata import NDData +from astropy.table import Table +from astropy.tests.helper import remote_data +import astropy.units as u +from astropy.utils import minversion +from astropy.wcs.utils import pixel_to_skycoord + +from ..core import aperture_photometry +from ..circle import (CircularAperture, CircularAnnulus, SkyCircularAperture, + SkyCircularAnnulus) +from ..ellipse import (EllipticalAperture, EllipticalAnnulus, + SkyEllipticalAperture, SkyEllipticalAnnulus) +from ..rectangle import (RectangularAperture, RectangularAnnulus, + SkyRectangularAperture, SkyRectangularAnnulus) +from ...datasets import (get_path, make_4gaussians_image, make_wcs, + make_imagehdu, make_4gaussians_cube, make_cubehdu) + +try: + import matplotlib # noqa + HAS_MATPLOTLIB = True +except ImportError: + HAS_MATPLOTLIB = False + +NUMPY_LT_1_14 = not minversion('numpy', '1.14dev') + +APERTURE_CL = [CircularAperture, + CircularAnnulus, + EllipticalAperture, + EllipticalAnnulus, + RectangularAperture, + RectangularAnnulus] + +TEST_APERTURES = list(zip(APERTURE_CL, ((3.,), (3., 5.), + (3., 5., 1.), (3., 5., 4., 1.), + (5, 8, np.pi / 4), + (8, 12, 8, np.pi / 8)))) + + +# THIS TEST HAS BEEN DONE. +@pytest.mark.parametrize(('aperture_class', 'params'), TEST_APERTURES) +def test_outside_array(aperture_class, params): + data = np.ones((100, 10, 10), dtype=np.float) + aperture = aperture_class((-60, 60), *params) + fluxtable = aperture_photometry(data, aperture) + # aperture is fully outside array: + assert np.all(np.isnan(fluxtable['aperture_sum'])) + + +# THIS TEST HAS BEEN DONE. +@pytest.mark.parametrize(('aperture_class', 'params'), TEST_APERTURES) +def test_inside_array_simple(aperture_class, params): + data = np.ones((100, 40, 40), dtype=np.float) + aperture = aperture_class((20., 20.), *params) + table1 = aperture_photometry(data, aperture, method='center', subpixels=10) + table2 = aperture_photometry(data, aperture, method='subpixel', + subpixels=10) + table3 = aperture_photometry(data, aperture, method='exact', subpixels=10) + true_flux = aperture.area() + + if not isinstance(aperture, (RectangularAperture, RectangularAnnulus)): + assert_allclose(table3['aperture_sum'], true_flux) + assert_allclose(table2['aperture_sum'], table3['aperture_sum'], + atol=0.1) + + assert np.all(table1['aperture_sum'] < table3['aperture_sum']) + + +# Unnecessary. +@pytest.mark.skipif('not HAS_MATPLOTLIB') +@pytest.mark.parametrize(('aperture_class', 'params'), TEST_APERTURES) +def test_aperture_plots(aperture_class, params): + # This test should run without any errors, and there is no return + # value. + # TODO: check the content of the plot + aperture = aperture_class((20., 20.), *params) + aperture.plot() + + +# Unnecessary. +def test_aperture_pixel_positions(): + pos1 = (10, 20) + pos2 = u.Quantity((10, 20), unit=u.pixel) + pos3 = ((10, 20, 30), (10, 20, 30)) + pos3_pairs = ((10, 10), (20, 20), (30, 30)) + + r = 3 + ap1 = CircularAperture(pos1, r) + ap2 = CircularAperture(pos2, r) + ap3 = CircularAperture(pos3, r) + + assert_allclose(np.atleast_2d(pos1), ap1.positions) + assert_allclose(np.atleast_2d(pos2.value), ap2.positions) + assert_allclose(pos3_pairs, ap3.positions) + + +class BaseTestAperturePhotometry(object): + + def test_scalar_error(self): + # Scalar error + error = 1. + if not hasattr(self, 'mask'): + mask = None + true_error = np.sqrt(self.area) + else: + mask = self.mask + # 1 masked pixel + true_error = np.sqrt(self.area - 1) + + table1 = aperture_photometry(self.data, + self.aperture, method='center', + mask=mask, error=error) + table2 = aperture_photometry(self.data, + self.aperture, + method='subpixel', subpixels=12, + mask=mask, error=error) + table3 = aperture_photometry(self.data, + self.aperture, method='exact', + mask=mask, error=error) + + if not isinstance(self.aperture, (RectangularAperture, + RectangularAnnulus)): + [assert_allclose(tab3_ap_sum_Col_i, true_flux_i) + for (tab3_ap_sum_Col_i, true_flux_i) + in zip(np.atleast_2d(table3['aperture_sum']), np.atleast_1d(self.true_flux))] + assert_allclose(table2['aperture_sum'], table3['aperture_sum'], + atol=0.1) + assert np.all(table1['aperture_sum'] < table3['aperture_sum']) + + if not isinstance(self.aperture, (RectangularAperture, + RectangularAnnulus)): + [assert_allclose(tab3_ap_sum_err_Col_i, true_flux_i) + for (tab3_ap_sum_err_Col_i, true_flux_i) + in zip(np.atleast_2d(table3['aperture_sum_err']), np.atleast_1d(true_error))] + assert_allclose(table2['aperture_sum'], table3['aperture_sum'], + atol=0.1) + assert np.all(table1['aperture_sum_err'] < table3['aperture_sum_err']) + + def test_array_error(self): + # Array error + error = np.ones(self.data.shape, dtype=np.float) + if not hasattr(self, 'mask'): + mask = None + true_error = np.sqrt(self.area) + else: + mask = self.mask + # 1 masked pixel + true_error = np.sqrt(self.area - 1) + + table1 = aperture_photometry(self.data, + self.aperture, method='center', + mask=mask, error=error) + table2 = aperture_photometry(self.data, + self.aperture, + method='subpixel', subpixels=12, + mask=mask, error=error) + table3 = aperture_photometry(self.data, + self.aperture, method='exact', + mask=mask, error=error) + + if not isinstance(self.aperture, (RectangularAperture, + RectangularAnnulus)): + [assert_allclose(tab3_ap_sum_Col_i, true_flux_i) + for (tab3_ap_sum_Col_i, true_flux_i) + in zip(np.atleast_2d(table3['aperture_sum']), np.atleast_1d(self.true_flux))] + #assert_allclose(table3['aperture_sum'], self.true_flux) + assert_allclose(table2['aperture_sum'], table3['aperture_sum'], + atol=0.1) + assert np.all(table1['aperture_sum'] < table3['aperture_sum']) + + if not isinstance(self.aperture, (RectangularAperture, + RectangularAnnulus)): + [assert_allclose(tab3_ap_sum_err_Col_i, true_flux_i) + for (tab3_ap_sum_err_Col_i, true_flux_i) + in zip(np.atleast_2d(table3['aperture_sum_err']), np.atleast_1d(true_error))] + #assert_allclose(table3['aperture_sum_err'], true_error) + assert_allclose(table2['aperture_sum_err'], + table3['aperture_sum_err'], atol=0.1) + assert np.all(table1['aperture_sum_err'] < table3['aperture_sum_err']) + + +# THIS TEST HAS BEEN DONE. +class TestCircular(BaseTestAperturePhotometry): + + def setup_class(self): + self.data = np.ones((100, 40, 40), dtype=np.float) + position = (20., 20.) + r = 10. + self.aperture = CircularAperture(position, r) + self.area = np.pi * r * r + self.true_flux = self.area + + +# THIS TEST HAS BEEN DONE. +class TestCircularArray(BaseTestAperturePhotometry): + + def setup_class(self): + self.data = np.ones((100, 40, 40), dtype=np.float) + position = ((20., 20.), (25., 25.)) + r = 10. + self.aperture = CircularAperture(position, r) + self.area = np.pi * r * r + self.area = np.array((self.area, ) * 2) + self.true_flux = self.area + + +# THIS TEST HAS BEEN DONE. +class TestCircularAnnulus(BaseTestAperturePhotometry): + + def setup_class(self): + self.data = np.ones((100, 40, 40), dtype=np.float) + position = (20., 20.) + r_in = 8. + r_out = 10. + self.aperture = CircularAnnulus(position, r_in, r_out) + self.area = np.pi * (r_out * r_out - r_in * r_in) + self.true_flux = self.area + + +# THIS TEST HAS BEEN DONE. +class TestCircularAnnulusArray(BaseTestAperturePhotometry): + + def setup_class(self): + self.data = np.ones((100, 40, 40), dtype=np.float) + position = ((20., 20.), (25., 25.)) + r_in = 8. + r_out = 10. + self.aperture = CircularAnnulus(position, r_in, r_out) + self.area = np.pi * (r_out * r_out - r_in * r_in) + self.area = np.array((self.area, ) * 2) + self.true_flux = self.area + + +# THIS TEST HAS BEEN DONE. +class TestElliptical(BaseTestAperturePhotometry): + + def setup_class(self): + self.data = np.ones((100, 40, 40), dtype=np.float) + position = (20., 20.) + a = 10. + b = 5. + theta = -np.pi / 4. + self.aperture = EllipticalAperture(position, a, b, theta) + self.area = np.pi * a * b + self.true_flux = self.area + + +# THIS TEST HAS BEEN DONE. +class TestEllipticalAnnulus(BaseTestAperturePhotometry): + + def setup_class(self): + self.data = np.ones((100, 40, 40), dtype=np.float) + position = (20., 20.) + a_in = 5. + a_out = 8. + b_out = 5. + theta = -np.pi / 4. + self.aperture = EllipticalAnnulus(position, a_in, a_out, b_out, theta) + self.area = (np.pi * (a_out * b_out) - + np.pi * (a_in * b_out * a_in / a_out)) + self.true_flux = self.area + + +# THIS TEST HAS BEEN DONE. +class TestRectangularAperture(BaseTestAperturePhotometry): + + def setup_class(self): + self.data = np.ones((100, 40, 40), dtype=np.float) + position = (20., 20.) + h = 5. + w = 8. + theta = np.pi / 4. + self.aperture = RectangularAperture(position, w, h, theta) + self.area = h * w + self.true_flux = self.area + + +# THIS TEST HAS BEEN DONE. +class TestRectangularAnnulus(BaseTestAperturePhotometry): + + def setup_class(self): + self.data = np.ones((100, 40, 40), dtype=np.float) + position = (20., 20.) + h_out = 8. + w_in = 8. + w_out = 12. + h_in = w_in * h_out / w_out + theta = np.pi / 8. + self.aperture = RectangularAnnulus(position, w_in, w_out, h_out, theta) + self.area = h_out * w_out - h_in * w_in + self.true_flux = self.area + + +# THIS TEST HAS BEEN DONE. +class TestMaskedSkipCircular(BaseTestAperturePhotometry): + + def setup_class(self): + self.data = np.ones((100, 40, 40), dtype=np.float) + self.mask = np.zeros((100, 40, 40), dtype=bool) + self.mask[:, 20, 20] = True + position = (20., 20.) + r = 10. + self.aperture = CircularAperture(position, r) + self.area = np.pi * r * r + self.true_flux = self.area - 1 + + +class BaseTestDifferentData(object): + + def test_basic_circular_aperture_photometry(self): + aperture = CircularAperture(self.position, self.radius) + table = aperture_photometry(self.data, aperture, + method='exact', unit='adu') + + assert_allclose(table['aperture_sum'].value, self.true_flux) + assert table['aperture_sum'].unit, self.fluxunit + + assert np.all(table['xcenter'].value == + np.transpose(self.position)[0]) + assert np.all(table['ycenter'].value == + np.transpose(self.position)[1]) + + +# THIS TEST HAS BEEN DONE. +class TestInputPrimaryHDU(BaseTestDifferentData): + + def setup_class(self): + data = np.ones((100, 40, 40), dtype=np.float) + self.data = fits.ImageHDU(data=data) + self.data.header['BUNIT'] = 'adu' + self.radius = 3 + self.position = (20, 20) + self.true_flux = np.pi * self.radius * self.radius + self.fluxunit = u.adu + + +# THIS TEST HAS BEEN DONE. +class TestInputHDUList(BaseTestDifferentData): + + def setup_class(self): + data0 = np.ones((100, 40, 40), dtype=np.float) + data1 = np.empty((100, 40, 40), dtype=np.float) + data1.fill(2) + self.data = fits.HDUList([fits.ImageHDU(data=data0), + fits.ImageHDU(data=data1)]) + self.radius = 3 + self.position = (20, 20) + # It should stop at the first extension + self.true_flux = np.pi * self.radius * self.radius + + +# THIS TEST HAS BEEN DONE. +class TestInputHDUDifferentBUNIT(BaseTestDifferentData): + + def setup_class(self): + data = np.ones((100, 40, 40), dtype=np.float) + self.data = fits.ImageHDU(data=data) + self.data.header['BUNIT'] = 'Jy' + self.radius = 3 + self.position = (20, 20) + self.true_flux = np.pi * self.radius * self.radius + self.fluxunit = u.adu + + +# THIS TEST HAS BEEN DONE. +class TestInputNDData(BaseTestDifferentData): + + def setup_class(self): + data = np.ones((100, 40, 40), dtype=np.float) + self.data = NDData(data, unit=u.adu) + self.radius = 3 + self.position = [(20, 20), (30, 30)] + self.true_flux = np.pi * self.radius * self.radius + self.fluxunit = u.adu + + +# Needs updating. +@remote_data +def test_wcs_based_photometry_to_catalogue(): + pathcat = get_path('spitzer_example_catalog.xml', location='remote') + pathhdu = get_path('spitzer_example_image.fits', location='remote') + hdu = fits.open(pathhdu) + scale = hdu[0].header['PIXSCAL1'] + + catalog = Table.read(pathcat) + + pos_skycoord = SkyCoord(catalog['l'], catalog['b'], frame='galactic') + + photometry_skycoord = aperture_photometry( + hdu, SkyCircularAperture(pos_skycoord, 4 * u.arcsec)) + + photometry_skycoord_pix = aperture_photometry( + hdu, SkyCircularAperture(pos_skycoord, 4. / scale * u.pixel)) + + assert_allclose(photometry_skycoord['aperture_sum'], + photometry_skycoord_pix['aperture_sum']) + + # Photometric unit conversion is needed to match the catalogue + factor = (1.2 * u.arcsec) ** 2 / u.pixel + converted_aperture_sum = (photometry_skycoord['aperture_sum'] * + factor).to(u.mJy / u.pixel) + + fluxes_catalog = catalog['f4_5'].filled() + + # There shouldn't be large outliers, but some differences is OK, as + # fluxes_catalog is based on PSF photometry, etc. + assert_allclose(fluxes_catalog, converted_aperture_sum.value, rtol=1e0) + + assert(np.mean(np.fabs(((fluxes_catalog - converted_aperture_sum.value) / + fluxes_catalog))) < 0.1) + + +# THIS TEST HAS BEEN DONE. +def test_wcs_based_photometry(): + data = make_4gaussians_cube() + wcs = make_wcs(data.shape) + hdu = make_cubehdu(data, wcs=wcs) + + # hard wired positions in make_4gaussians_image + pos_orig_pixel = u.Quantity(([160., 25., 150., 90.], + [70., 40., 25., 60.]), unit=u.pixel) + + pos_skycoord = pixel_to_skycoord(pos_orig_pixel[0], pos_orig_pixel[1], wcs) + + pos_skycoord_s = pos_skycoord[2] + + photometry_skycoord_circ = aperture_photometry( + hdu, SkyCircularAperture(pos_skycoord, 3 * u.arcsec)) + photometry_skycoord_circ_2 = aperture_photometry( + hdu, SkyCircularAperture(pos_skycoord, 2 * u.arcsec)) + photometry_skycoord_circ_s = aperture_photometry( + hdu, SkyCircularAperture(pos_skycoord_s, 3 * u.arcsec)) + + assert_allclose(np.atleast_2d(photometry_skycoord_circ['aperture_sum'][2]), + np.atleast_2d(photometry_skycoord_circ_s['aperture_sum'])) + + photometry_skycoord_circ_ann = aperture_photometry( + hdu, SkyCircularAnnulus(pos_skycoord, 2 * u.arcsec, 3 * u.arcsec)) + photometry_skycoord_circ_ann_s = aperture_photometry( + hdu, SkyCircularAnnulus(pos_skycoord_s, 2 * u.arcsec, 3 * u.arcsec)) + + assert_allclose(np.atleast_2d(photometry_skycoord_circ_ann['aperture_sum'][2]), + np.atleast_2d(photometry_skycoord_circ_ann_s['aperture_sum'])) + + assert_allclose(photometry_skycoord_circ_ann['aperture_sum'], + photometry_skycoord_circ['aperture_sum'] - + photometry_skycoord_circ_2['aperture_sum']) + + photometry_skycoord_ell = aperture_photometry( + hdu, SkyEllipticalAperture(pos_skycoord, 3 * u.arcsec, + 3.0001 * u.arcsec, 45 * u.arcsec)) + photometry_skycoord_ell_2 = aperture_photometry( + hdu, SkyEllipticalAperture(pos_skycoord, 2 * u.arcsec, + 2.0001 * u.arcsec, 45 * u.arcsec)) + photometry_skycoord_ell_s = aperture_photometry( + hdu, SkyEllipticalAperture(pos_skycoord_s, 3 * u.arcsec, + 3.0001 * u.arcsec, 45 * u.arcsec)) + photometry_skycoord_ell_ann = aperture_photometry( + hdu, SkyEllipticalAnnulus(pos_skycoord, 2 * u.arcsec, 3 * u.arcsec, + 3.0001 * u.arcsec, 45 * u.arcsec)) + photometry_skycoord_ell_ann_s = aperture_photometry( + hdu, SkyEllipticalAnnulus(pos_skycoord_s, 2 * u.arcsec, 3 * u.arcsec, + 3.0001 * u.arcsec, 45 * u.arcsec)) + + assert_allclose(np.atleast_2d(photometry_skycoord_ell['aperture_sum'][2]), + np.atleast_2d(photometry_skycoord_ell_s['aperture_sum'])) + + assert_allclose(np.atleast_2d(photometry_skycoord_ell_ann['aperture_sum'][2]), + np.atleast_2d(photometry_skycoord_ell_ann_s['aperture_sum'])) + + assert_allclose(photometry_skycoord_ell['aperture_sum'], + photometry_skycoord_circ['aperture_sum'], rtol=5e-3) + + assert_allclose(photometry_skycoord_ell_ann['aperture_sum'], + photometry_skycoord_ell['aperture_sum'] - + photometry_skycoord_ell_2['aperture_sum'], rtol=1e-4) + + photometry_skycoord_rec = aperture_photometry( + hdu, SkyRectangularAperture(pos_skycoord, + 6 * u.arcsec, 6 * u.arcsec, + 0 * u.arcsec), + method='subpixel', subpixels=20) + photometry_skycoord_rec_4 = aperture_photometry( + hdu, SkyRectangularAperture(pos_skycoord, + 4 * u.arcsec, 4 * u.arcsec, + 0 * u.arcsec), + method='subpixel', subpixels=20) + photometry_skycoord_rec_s = aperture_photometry( + hdu, SkyRectangularAperture(pos_skycoord_s, + 6 * u.arcsec, 6 * u.arcsec, + 0 * u.arcsec), + method='subpixel', subpixels=20) + photometry_skycoord_rec_ann = aperture_photometry( + hdu, SkyRectangularAnnulus(pos_skycoord, 4 * u.arcsec, 6 * u.arcsec, + 6 * u.arcsec, 0 * u.arcsec), + method='subpixel', subpixels=20) + photometry_skycoord_rec_ann_s = aperture_photometry( + hdu, SkyRectangularAnnulus(pos_skycoord_s, 4 * u.arcsec, 6 * u.arcsec, + 6 * u.arcsec, 0 * u.arcsec), + method='subpixel', subpixels=20) + + assert_allclose(np.atleast_2d(photometry_skycoord_rec['aperture_sum'][2]), + np.atleast_2d(photometry_skycoord_rec_s['aperture_sum'])) + + assert np.all(photometry_skycoord_rec['aperture_sum'] > + photometry_skycoord_circ['aperture_sum']) + + assert_allclose(np.atleast_2d(photometry_skycoord_rec_ann['aperture_sum'][2]), + np.atleast_2d(photometry_skycoord_rec_ann_s['aperture_sum'])) + + assert_allclose(photometry_skycoord_rec_ann['aperture_sum'], + photometry_skycoord_rec['aperture_sum'] - + photometry_skycoord_rec_4['aperture_sum'], rtol=1e-4) + + +# THIS TEST HAS BEEN DONE. +def test_basic_circular_aperture_photometry_unit(): + data1 = np.ones((100, 40, 40), dtype=np.float) + data2 = u.Quantity(data1, unit=u.adu) + + radius = 3 + position = (20, 20) + true_flux = np.pi * radius * radius + unit = u.adu + + table1 = aperture_photometry(data1, CircularAperture(position, radius), + unit=unit) + table2 = aperture_photometry(data2, CircularAperture(position, radius)) + + assert_allclose(table1['aperture_sum'].value, true_flux) + assert_allclose(table2['aperture_sum'].value, true_flux) + assert table1['aperture_sum'].unit == unit + assert table2['aperture_sum'].unit == data2.unit == unit + + +# THIS TEST HAS BEEN DONE. +def test_aperture_photometry_with_error_units(): + """Test aperture_photometry when error has units (see #176).""" + + data1 = np.ones((100, 40, 40), dtype=np.float) + data2 = u.Quantity(data1, unit=u.adu) + error = u.Quantity(data1, unit=u.adu) + radius = 3 + true_flux = np.pi * radius * radius + unit = u.adu + position = (20, 20) + table1 = aperture_photometry(data2, CircularAperture(position, radius), + error=error) + assert_allclose(table1['aperture_sum'].value, true_flux) + assert_allclose(table1['aperture_sum_err'].value, np.sqrt(true_flux)) + assert table1['aperture_sum'].unit == unit + assert table1['aperture_sum_err'].unit == unit + + +# THIS TEST HAS BEEN DONE. +def test_aperture_photometry_inputs_with_mask(): + """ + Test that aperture_photometry does not modify the input + data or error array when a mask is input. + """ + + data = np.ones((100, 5, 5)) + aperture = CircularAperture((2, 2), 2.) + mask = np.zeros_like(data, dtype=bool) + data[:, 2, 2] = 100. # bad pixel + mask[:, 2, 2] = True + error = np.sqrt(data) + data_in = data.copy() + error_in = error.copy() + t1 = aperture_photometry(data, aperture, error=error, mask=mask) + assert_array_equal(data, data_in) + assert_array_equal(error, error_in) + assert_allclose(t1['aperture_sum'][0], 11.5663706144) + t2 = aperture_photometry(data, aperture) + assert_allclose(t2['aperture_sum'][0], 111.566370614) + + +TEST_ELLIPSE_EXACT_APERTURES = [(3.469906, 3.923861394, 3.), + (0.3834415188257778, 0.3834415188257778, 0.3)] + + +# THIS TEST HAS BEEN DONE. +# Unnecessary. +@pytest.mark.parametrize('x,y,r', TEST_ELLIPSE_EXACT_APERTURES) +def test_ellipse_exact_grid(x, y, r): + """ + Test elliptical exact aperture photometry on a grid of pixel positions. + + This is a regression test for the bug discovered in this issue: + https://github.com/astropy/photutils/issues/198 + """ + + data = np.ones((100, 10, 10)) + + aperture = EllipticalAperture((x, y), r, r, 0.) + t = aperture_photometry(data, aperture, method='exact') + actual = t['aperture_sum'][0] / (np.pi * r ** 2) + assert_allclose(actual, 1) + + +# THIS TEST HAS BEEN DONE. +@pytest.mark.parametrize('value', [np.nan, np.inf]) +def test_nan_inf_mask(value): + """Test that nans and infs are properly masked [267].""" + + data = np.ones((100, 9, 9)) + mask = np.zeros_like(data, dtype=bool) + data[:, 4, 4] = value + mask[:, 4, 4] = True + radius = 2. + aper = CircularAperture((4, 4), radius) + tbl = aperture_photometry(data, aper, mask=mask) + desired = (np.pi * radius**2) - 1 + assert_allclose(tbl['aperture_sum'], desired) + + +# THIS TEST HAS BEEN DONE. +def test_aperture_partial_overlap(): + data = np.ones((100, 20, 20)) + error = np.ones((100, 20, 20)) + xypos = [(10, 10), (0, 0), (0, 19), (19, 0), (19, 19)] + r = 5. + aper = CircularAperture(xypos, r=r) + tbl = aperture_photometry(data, aper, error=error) + assert_allclose(tbl['aperture_sum'][0], np.pi * r ** 2) + assert_array_less(tbl['aperture_sum'][1:], np.pi * r ** 2) + + unit = u.MJy / u.sr + tbl = aperture_photometry(data * unit, aper, error=error * unit) + assert_allclose(tbl['aperture_sum'][0].value, np.pi * r ** 2) + assert_array_less(tbl['aperture_sum'][1:].value, np.pi * r ** 2) + assert_array_less(tbl['aperture_sum_err'][1:].value, np.pi * r ** 2) + assert tbl['aperture_sum'].unit == unit + assert tbl['aperture_sum_err'].unit == unit + + +# THIS TEST HAS BEEN DONE. +# Unnecessary. +def test_pixel_aperture_repr(): + aper = CircularAperture((10, 20), r=3.0) + a_repr = '' + a_str = 'Aperture: CircularAperture\npositions: [[10, 20]]\nr: 3.0' + assert repr(aper) == a_repr + assert str(aper) == a_str + + aper = CircularAnnulus((10, 20), r_in=3.0, r_out=5.0) + a_repr = '' + a_str = ('Aperture: CircularAnnulus\npositions: [[10, 20]]\nr_in: 3.0\n' + 'r_out: 5.0') + assert repr(aper) == a_repr + assert str(aper) == a_str + + aper = EllipticalAperture((10, 20), a=5.0, b=3.0, theta=15.0) + a_repr = '' + a_str = ('Aperture: EllipticalAperture\npositions: [[10, 20]]\n' + 'a: 5.0\nb: 3.0\ntheta: 15.0') + assert repr(aper) == a_repr + assert str(aper) == a_str + + aper = EllipticalAnnulus((10, 20), a_in=4.0, a_out=8.0, b_out=4.0, + theta=15.0) + a_repr = ('') + a_str = ('Aperture: EllipticalAnnulus\npositions: [[10, 20]]\na_in: ' + '4.0\na_out: 8.0\nb_out: 4.0\ntheta: 15.0') + assert repr(aper) == a_repr + assert str(aper) == a_str + + aper = RectangularAperture((10, 20), w=5.0, h=3.0, theta=15.0) + a_repr = '' + a_str = ('Aperture: RectangularAperture\npositions: [[10, 20]]\n' + 'w: 5.0\nh: 3.0\ntheta: 15.0') + assert repr(aper) == a_repr + assert str(aper) == a_str + + aper = RectangularAnnulus((10, 20), w_in=4.0, w_out=8.0, h_out=4.0, + theta=15.0) + a_repr = ('') + a_str = ('Aperture: RectangularAnnulus\npositions: [[10, 20]]\n' + 'w_in: 4.0\nw_out: 8.0\nh_out: 4.0\ntheta: 15.0') + assert repr(aper) == a_repr + assert str(aper) == a_str + + +# THIS TEST HAS BEEN DONE. +# Unnecessary. +def test_sky_aperture_repr(): + s = SkyCoord([1, 2], [3, 4], unit='deg') + + aper = SkyCircularAperture(s, r=3*u.pix) + if NUMPY_LT_1_14: + a_repr = (', r=3.0 pix)>') + a_str = ('Aperture: SkyCircularAperture\npositions: \n' + 'r: 3.0 pix') + else: + a_repr = (', r=3.0 pix)>') + a_str = ('Aperture: SkyCircularAperture\npositions: \n' + 'r: 3.0 pix') + + assert repr(aper) == a_repr + assert str(aper) == a_str + + aper = SkyCircularAnnulus(s, r_in=3.*u.pix, r_out=5*u.pix) + if NUMPY_LT_1_14: + a_repr = (', r_in=3.0 pix, ' + 'r_out=5.0 pix)>') + a_str = ('Aperture: SkyCircularAnnulus\npositions: \n' + 'r_in: 3.0 pix\nr_out: 5.0 pix') + else: + a_repr = (', r_in=3.0 pix, r_out=5.0 pix)>') + a_str = ('Aperture: SkyCircularAnnulus\npositions: \n' + 'r_in: 3.0 pix\nr_out: 5.0 pix') + + assert repr(aper) == a_repr + assert str(aper) == a_str + + aper = SkyEllipticalAperture(s, a=3*u.pix, b=5*u.pix, theta=15*u.deg) + if NUMPY_LT_1_14: + a_repr = (', a=3.0 pix, b=5.0 pix,' + ' theta=15.0 deg)>') + a_str = ('Aperture: SkyEllipticalAperture\npositions: \n' + 'a: 3.0 pix\nb: 5.0 pix\ntheta: 15.0 deg') + else: + a_repr = (', a=3.0 pix, b=5.0 pix,' + ' theta=15.0 deg)>') + a_str = ('Aperture: SkyEllipticalAperture\npositions: \n' + 'a: 3.0 pix\nb: 5.0 pix\ntheta: 15.0 deg') + + assert repr(aper) == a_repr + assert str(aper) == a_str + + aper = SkyEllipticalAnnulus(s, a_in=3*u.pix, a_out=5*u.pix, b_out=3*u.pix, + theta=15*u.deg) + if NUMPY_LT_1_14: + a_repr = (', a_in=3.0 pix, ' + 'a_out=5.0 pix, b_out=3.0 pix, theta=15.0 deg)>') + a_str = ('Aperture: SkyEllipticalAnnulus\npositions: \n' + 'a_in: 3.0 pix\na_out: 5.0 pix\nb_out: 3.0 pix\n' + 'theta: 15.0 deg') + else: + a_repr = (', a_in=3.0 pix, ' + 'a_out=5.0 pix, b_out=3.0 pix, theta=15.0 deg)>') + a_str = ('Aperture: SkyEllipticalAnnulus\npositions: \n' + 'a_in: 3.0 pix\na_out: 5.0 pix\nb_out: 3.0 pix\n' + 'theta: 15.0 deg') + + assert repr(aper) == a_repr + assert str(aper) == a_str + + aper = SkyRectangularAperture(s, w=3*u.pix, h=5*u.pix, theta=15*u.deg) + if NUMPY_LT_1_14: + a_repr = (', w=3.0 pix, h=5.0 pix' + ', theta=15.0 deg)>') + a_str = ('Aperture: SkyRectangularAperture\npositions: \n' + 'w: 3.0 pix\nh: 5.0 pix\ntheta: 15.0 deg') + else: + a_repr = (', w=3.0 pix, h=5.0 pix' + ', theta=15.0 deg)>') + a_str = ('Aperture: SkyRectangularAperture\npositions: \n' + 'w: 3.0 pix\nh: 5.0 pix\ntheta: 15.0 deg') + + assert repr(aper) == a_repr + assert str(aper) == a_str + + aper = SkyRectangularAnnulus(s, w_in=3*u.pix, w_out=3.4*u.pix, + h_out=5*u.pix, theta=15*u.deg) + if NUMPY_LT_1_14: + a_repr = (', w_in=3.0 pix, ' + 'w_out=3.4 pix, h_out=5.0 pix, theta=15.0 deg)>') + a_str = ('Aperture: SkyRectangularAnnulus\npositions: \n' + 'w_in: 3.0 pix\nw_out: 3.4 pix\nh_out: 5.0 pix\n' + 'theta: 15.0 deg') + else: + a_repr = (', w_in=3.0 pix, ' + 'w_out=3.4 pix, h_out=5.0 pix, theta=15.0 deg)>') + a_str = ('Aperture: SkyRectangularAnnulus\npositions: \n' + 'w_in: 3.0 pix\nw_out: 3.4 pix\nh_out: 5.0 pix\n' + 'theta: 15.0 deg') + + assert repr(aper) == a_repr + assert str(aper) == a_str + + +# THIS TEST HAS BEEN DONE. +# Unnecessary. +def test_rectangular_bbox(): + # odd sizes + width = 7 + height = 3 + a = RectangularAperture((50, 50), w=width, h=height, theta=0) + assert a.bounding_boxes[0].shape == (height, width) + + a = RectangularAperture((50.5, 50.5), w=width, h=height, theta=0) + assert a.bounding_boxes[0].shape == (height + 1, width + 1) + + a = RectangularAperture((50, 50), w=width, h=height, theta=90.*np.pi/180.) + assert a.bounding_boxes[0].shape == (width, height) + + # even sizes + width = 8 + height = 4 + a = RectangularAperture((50, 50), w=width, h=height, theta=0) + assert a.bounding_boxes[0].shape == (height + 1, width + 1) + + a = RectangularAperture((50.5, 50.5), w=width, h=height, theta=0) + assert a.bounding_boxes[0].shape == (height, width) + + a = RectangularAperture((50.5, 50.5), w=width, h=height, + theta=90.*np.pi/180.) + assert a.bounding_boxes[0].shape == (width, height) + + +# THIS TEST HAS BEEN DONE. +# Unnecessary. +def test_elliptical_bbox(): + # integer axes + a = 7 + b = 3 + ap = EllipticalAperture((50, 50), a=a, b=b, theta=0) + assert ap.bounding_boxes[0].shape == (2*b + 1, 2*a + 1) + + ap = EllipticalAperture((50.5, 50.5), a=a, b=b, theta=0) + assert ap.bounding_boxes[0].shape == (2*b, 2*a) + + ap = EllipticalAperture((50, 50), a=a, b=b, theta=90.*np.pi/180.) + assert ap.bounding_boxes[0].shape == (2*a + 1, 2*b + 1) + + # fractional axes + a = 7.5 + b = 4.5 + ap = EllipticalAperture((50, 50), a=a, b=b, theta=0) + assert ap.bounding_boxes[0].shape == (2*b, 2*a) + + ap = EllipticalAperture((50.5, 50.5), a=a, b=b, theta=0) + assert ap.bounding_boxes[0].shape == (2*b + 1, 2*a + 1) + + ap = EllipticalAperture((50, 50), a=a, b=b, theta=90.*np.pi/180.) + assert ap.bounding_boxes[0].shape == (2*a, 2*b) + + +# THIS TEST HAS BEEN DONE. +def test_to_sky_pixel(): + data = make_4gaussians_cube() + wcs = make_wcs(data.shape) + + ap = CircularAperture(((12.3, 15.7), (48.19, 98.14)), r=3.14) + ap2 = ap.to_sky(wcs).to_pixel(wcs) + assert_allclose(ap.positions, ap2.positions) + assert_allclose(ap.r, ap2.r) + + ap = CircularAnnulus(((12.3, 15.7), (48.19, 98.14)), r_in=3.14, + r_out=5.32) + ap2 = ap.to_sky(wcs).to_pixel(wcs) + assert_allclose(ap.positions, ap2.positions) + assert_allclose(ap.r_in, ap2.r_in) + assert_allclose(ap.r_out, ap2.r_out) + + ap = EllipticalAperture(((12.3, 15.7), (48.19, 98.14)), a=3.14, b=5.32, + theta=103.*np.pi/180.) + ap2 = ap.to_sky(wcs).to_pixel(wcs) + assert_allclose(ap.positions, ap2.positions) + assert_allclose(ap.a, ap2.a) + assert_allclose(ap.b, ap2.b) + assert_allclose(ap.theta, ap2.theta) + + ap = EllipticalAnnulus(((12.3, 15.7), (48.19, 98.14)), a_in=3.14, + a_out=15.32, b_out=4.89, theta=103.*np.pi/180.) + ap2 = ap.to_sky(wcs).to_pixel(wcs) + assert_allclose(ap.positions, ap2.positions) + assert_allclose(ap.a_in, ap2.a_in) + assert_allclose(ap.a_out, ap2.a_out) + assert_allclose(ap.b_out, ap2.b_out) + assert_allclose(ap.theta, ap2.theta) + + ap = RectangularAperture(((12.3, 15.7), (48.19, 98.14)), w=3.14, h=5.32, + theta=103.*np.pi/180.) + ap2 = ap.to_sky(wcs).to_pixel(wcs) + assert_allclose(ap.positions, ap2.positions) + assert_allclose(ap.w, ap2.w) + assert_allclose(ap.h, ap2.h) + assert_allclose(ap.theta, ap2.theta) + + ap = RectangularAnnulus(((12.3, 15.7), (48.19, 98.14)), w_in=3.14, + w_out=15.32, h_out=4.89, theta=103.*np.pi/180.) + ap2 = ap.to_sky(wcs).to_pixel(wcs) + assert_allclose(ap.positions, ap2.positions) + assert_allclose(ap.w_in, ap2.w_in) + assert_allclose(ap.w_out, ap2.w_out) + assert_allclose(ap.h_out, ap2.h_out) + assert_allclose(ap.theta, ap2.theta) diff --git a/photutils/datasets/make.py b/photutils/datasets/make.py index a75c271b9..337ddc7fb 100644 --- a/photutils/datasets/make.py +++ b/photutils/datasets/make.py @@ -21,7 +21,9 @@ 'make_random_models_table', 'make_random_gaussians_table', 'make_model_sources_image', 'make_gaussian_sources_image', 'make_4gaussians_image', 'make_100gaussians_image', - 'make_wcs', 'make_imagehdu'] + 'make_4gaussians_cube', 'make_100gaussians_cube', + 'make_wcs', 'make_imagehdu', 'make_cubehdu', + ] def apply_poisson_noise(data, random_state=None): @@ -556,6 +558,7 @@ def make_gaussian_sources_image(shape, source_table, oversample=1): oversample=oversample) + def make_4gaussians_image(noise=True): """ Make an example image containing four 2D Gaussians plus a constant @@ -670,7 +673,168 @@ def make_100gaussians_image(noise=True): return data + +def make_4gaussians_cube(noise=True): + """ + Make an example datacube containing four 2D Gaussians plus a constant + background. + + The background has a mean of 5. + + If ``noise`` is `True`, then Gaussian noise with a mean of 0 and a + standard deviation of 5 is added to the output image. + + Parameters + ---------- + noise : bool, optional + Whether to include noise in the output image (default is + `True`). + + Returns + ------- + datacube : 3D `~numpy.ndarray` + Image containing four 2D Gaussian sources. + + See Also + -------- + make_100gaussians_cube + + Examples + -------- + .. plot:: + :include-source: + + from photutils import datasets + image = datasets.make_4gaussians_image() + plt.imshow(image, origin='lower', interpolation='nearest') + """ + + table = Table() + table['amplitude'] = [50, 70, 150, 210] + table['x_mean'] = [160, 25, 150, 90] + table['y_mean'] = [70, 40, 25, 60] + table['x_stddev'] = [15.2, 5.1, 3., 8.1] + table['y_stddev'] = [2.6, 2.5, 3., 4.7] + table['theta'] = np.array([145., 20., 0., 60.]) * np.pi / 180. + + shape = (1000, 100, 200) + data = np.empty(shape) + data[:, :, :] = make_gaussian_sources_image(shape[-2:], table)[None, :, :] + 5. + + if noise: + data += make_noise_image(shape[-2:], type='gaussian', mean=0., + stddev=5., random_state=12345)[None, :, :] + + return data + + +def make_100gaussians_cube(noise=True): + """ + Make an example datacube containing 100 2D Gaussians plus a constant + background. + + The background has a mean of 5. + + If ``noise`` is `True`, then Gaussian noise with a mean of 0 and a + standard deviation of 2 is added to the output image. + + Parameters + ---------- + noise : bool, optional + Whether to include noise in the output image (default is + `True`). + + Returns + ------- + datacube : 3D `~numpy.ndarray` + Datacube containing 100 2D Gaussian sources. + + See Also + -------- + make_4gaussians_cube + + Examples + -------- + .. plot:: + :include-source: + + from photutils import datasets + image = datasets.make_100gaussians_image() + plt.imshow(image, origin='lower', interpolation='nearest') + """ + + n_sources = 100 + flux_range = [500, 1000] + xmean_range = [0, 500] + ymean_range = [0, 300] + xstddev_range = [1, 5] + ystddev_range = [1, 5] + params = OrderedDict([('flux', flux_range), + ('x_mean', xmean_range), + ('y_mean', ymean_range), + ('x_stddev', xstddev_range), + ('y_stddev', ystddev_range), + ('theta', [0, 2*np.pi])]) + + sources = make_random_gaussians_table(n_sources, params, + random_state=12345) + + shape = (1000, 300, 500) + data = np.empty(shape) + data[:, :, :] = make_gaussian_sources_image(shape[-2:], table)[None, :, :] + 5. + + if noise: + data += make_noise_image(shape[-2:], type='gaussian', mean=0., + stddev=5., random_state=12345)[None, :, :] + + return data + + def make_wcs(shape, galactic=False): + """ + Create a simple celestial WCS object in either the ICRS or Galactic + coordinate frame. A wrapper to both `_make_wcs_2D` and `_make_wcs_3D`. + + Parameters + ---------- + shape : 2-tuple of int + The shape of the 2D array to be used with the output + `~astropy.wcs.WCS` object. + + galactic : bool, optional + If `True`, then the output WCS will be in the Galactic + coordinate frame. If `False` (default), then the output WCS + will be in the ICRS coordinate frame. + + Returns + ------- + wcs : `~astropy.wcs.WCS` object + The world coordinate system (WCS) transformation. + + See Also + -------- + make_imagehdu, make_cubehdu + + Examples + -------- + >>> from photutils.datasets import make_wcs + >>> shape = (100, 100) + >>> wcs = make_wcs(shape) + >>> print(wcs.wcs.crpix) # doctest: +FLOAT_CMP + [50. 50.] + >>> print(wcs.wcs.crval) # doctest: +FLOAT_CMP + [197.8925 -1.36555556] + """ + + if len(shape) == 2: + return _make_wcs_2D(shape, galactic=galactic) + elif len(shape) == 3: + return _make_wcs_3D(shape, galactic=galactic) + else: + raise ValueError('`shape` must correspond to a 2D or 3D array') + + +def _make_wcs_2D(shape, galactic=False): """ Create a simple celestial WCS object in either the ICRS or Galactic coordinate frame. @@ -725,6 +889,64 @@ def make_wcs(shape, galactic=False): return wcs +def _make_wcs_3D(shape, galactic=False): + """ + Create a simple celestial WCS object in either the ICRS or Galactic + coordinate frame. + + Parameters + ---------- + shape : 3-tuple of int + The shape of the 3D array to be used with the output + `~astropy.wcs.WCS` object. + + galactic : bool, optional + If `True`, then the output WCS will be in the Galactic + coordinate frame. If `False` (default), then the output WCS + will be in the ICRS coordinate frame. + + Returns + ------- + wcs : `~astropy.wcs.WCS` object + The world coordinate system (WCS) transformation. + + See Also + -------- + make_imagehdu + + Examples + -------- + >>> from photutils.datasets import make_wcs + >>> shape = (1000, 100, 100) + >>> wcs = make_wcs(shape) + >>> print(wcs.wcs.crpix) # doctest: +FLOAT_CMP + [50. 50. 500.] + >>> print(wcs.wcs.crval) # doctest: +FLOAT_CMP + [197.8925 -1.36555556 4861.] + """ + + wcs = WCS(naxis=3) + rho = np.pi / 3. + scale = 0.1 / 3600. + delta_wave = 0.5 + wcs._naxis1 = shape[2] # nx + wcs._naxis2 = shape[1] # ny + wcs._naxis3 = shape[0] # nz + wcs.wcs.crpix = [shape[2] / 2, shape[1] / 2, shape[0] / 2] # 1-indexed (x, y, z) + wcs.wcs.crval = [197.8925, -1.36555556, 4861.] + wcs.wcs.cunit = ['deg', 'deg', 'nm'] + wcs.wcs.cd = [[-scale * np.cos(rho), scale * np.sin(rho), 0.], + [ scale * np.sin(rho), scale * np.cos(rho), 0.], + [ 0., 0., delta_wave]] + if not galactic: + wcs.wcs.radesys = 'ICRS' + wcs.wcs.ctype = ['RA---TAN', 'DEC--TAN', 'WAVE'] + else: + wcs.wcs.ctype = ['GLON-CAR', 'GLAT-CAR', 'WAVE'] + + return wcs + + def make_imagehdu(data, wcs=None): """ Create a FITS `~astropy.io.fits.ImageHDU` containing the input 2D @@ -769,3 +991,50 @@ def make_imagehdu(data, wcs=None): header = None return fits.ImageHDU(data, header=header) + + + +def make_cubehdu(data, wcs=None): + """ + Create a FITS `~astropy.io.fits.ImageHDU` containing the input 3D + datacube. + + Parameters + ---------- + data : 3D array-like + The input 3D data. + + wcs : `~astropy.wcs.WCS`, optional + The world coordinate system (WCS) transformation to include in + the output FITS header. + + Returns + ------- + image_hdu : `~astropy.io.fits.ImageHDU` + The FITS `~astropy.io.fits.ImageHDU`. + + See Also + -------- + make_wcs + + Examples + -------- + >>> from photutils.datasets import make_cubehdu, make_wcs + >>> shape = (1000, 100, 100) + >>> data = np.ones(shape) + >>> wcs = make_wcs(shape) + >>> hdu = make_cubehdu(data, wcs=wcs) + >>> print(hdu.data.shape) + (1000, 100, 100) + """ + + data = np.asanyarray(data) + if data.ndim != 3: + raise ValueError('data must be a 3D array') + + if wcs is not None: + header = wcs.to_header() + else: + header = None + + return fits.ImageHDU(data, header=header) From 7b915971f68c4415b827441c447bb809d4c9e34d Mon Sep 17 00:00:00 2001 From: Francesco D'Eugenio Date: Tue, 3 Jul 2018 12:50:34 +0200 Subject: [PATCH 3/4] Fixed python2.7 compatibility issue --- photutils/aperture/mask.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/photutils/aperture/mask.py b/photutils/aperture/mask.py index aee1620bc..10575eea3 100644 --- a/photutils/aperture/mask.py +++ b/photutils/aperture/mask.py @@ -189,9 +189,9 @@ def cutout(self, data, fill_value=0., copy=False): # try this for speed -- the result may still be a partial # overlap, in which case the next block will be triggered if copy: - cutout = np.copy(data[(...,)+self.bbox.slices]) + cutout = np.copy(data[(Ellipsis,)+self.bbox.slices]) else: - cutout = data[(...,)+self.bbox.slices] + cutout = data[(Ellipsis,)+self.bbox.slices] if partial_overlap or (cutout.shape[-2:] != self.shape): slices_large, slices_small = self._overlap_slices(data.shape[-2:]) @@ -203,7 +203,7 @@ def cutout(self, data, fill_value=0., copy=False): output_shape = self.shape if data.ndim==2 else (data.shape[0],)+self.shape cutout = np.zeros(output_shape, dtype=data.dtype) cutout[:] = fill_value - cutout[(...,)+slices_small] = data[(...,)+slices_large] + cutout[(Ellipsis,)+slices_small] = data[(Ellipsis,)+slices_large] if isinstance(data, u.Quantity): cutout = u.Quantity(cutout, unit=data.unit) From 91153dc9562863c6967817cdcb4f1e1ff37b90e5 Mon Sep 17 00:00:00 2001 From: Francesco D'Eugenio Date: Fri, 11 Jan 2019 11:26:05 +0100 Subject: [PATCH 4/4] Updated to increase test coverage. --- photutils/datasets/make.py | 7 +-- photutils/datasets/tests/test_make.py | 73 ++++++++++++++++++++++++++- 2 files changed, 76 insertions(+), 4 deletions(-) diff --git a/photutils/datasets/make.py b/photutils/datasets/make.py index f4f8e32f5..08555ff10 100644 --- a/photutils/datasets/make.py +++ b/photutils/datasets/make.py @@ -155,7 +155,7 @@ def make_noise_image(shape, type='gaussian', mean=None, stddev=None, image = prng.poisson(lam=mean, size=shape) else: raise ValueError('Invalid type: {0}. Use one of ' - '{"gaussian", "poisson"}.'.format(type)) + '{{"gaussian", "poisson"}}.'.format(type)) return image @@ -863,11 +863,12 @@ def make_100gaussians_cube(noise=True): shape = (1000, 300, 500) data = np.empty(shape) - data[:, :, :] = make_gaussian_sources_image(shape[-2:], table)[None, :, :] + 5. + data[:, :, :] = make_gaussian_sources_image( + shape[-2:], sources)[None, :, :] + 5. if noise: data += make_noise_image(shape[-2:], type='gaussian', mean=0., - stddev=5., random_state=12345)[None, :, :] + stddev=2., random_state=12345)[None, :, :] return data diff --git a/photutils/datasets/tests/test_make.py b/photutils/datasets/tests/test_make.py index 55539c2d4..ee6f6479b 100644 --- a/photutils/datasets/tests/test_make.py +++ b/photutils/datasets/tests/test_make.py @@ -11,8 +11,10 @@ from .. import (make_noise_image, apply_poisson_noise, make_gaussian_sources_image, make_random_gaussians_table, make_4gaussians_image, make_100gaussians_image, + make_4gaussians_cube, make_100gaussians_cube, make_random_models_table, make_model_sources_image, - make_wcs, make_gaussian_prf_sources_image) + make_imagehdu, make_cubehdu, make_wcs, + make_gaussian_prf_sources_image) try: import scipy # noqa @@ -50,6 +52,13 @@ def test_make_noise_image_poisson(): assert_allclose(image.mean(), 1., atol=1.) +def test_make_noise_image_unktype(): + """Test if ValueError raises if type is not valid.""" + + with pytest.raises(ValueError): + shape = (100, 100) + make_noise_image(shape, 'invalid_type', mean=0., stddev=1.) + def test_make_noise_image_nomean(): """Test if ValueError raises if mean is not input.""" @@ -156,6 +165,22 @@ def test_make_100gaussians_image(): assert_allclose(image.sum(), data_sum, rtol=1.e-6) +def test_make_4gaussians_cube(): + shape = (1000, 100, 200) + data_sum = 176219180.59091496 + cube = make_4gaussians_cube() + assert cube.shape == shape + assert_allclose(cube.sum(), data_sum, rtol=1.e-6) + + +def test_make_100gaussians_cube(): + shape = (1000, 300, 500) + data_sum = 826182245.01251709 + cube = make_100gaussians_cube() + assert cube.shape == shape + assert_allclose(cube.sum(), data_sum, rtol=1.e-6) + + def test_make_random_models_table(): model = Moffat2D(amplitude=1) param_ranges = {'x_0': (0, 300), 'y_0': (0, 500), @@ -169,7 +194,9 @@ def test_make_random_models_table(): def test_make_wcs(): + shape = (100, 200) + wcs = make_wcs(shape) if astropy_version < '3.1': @@ -183,3 +210,47 @@ def test_make_wcs(): wcs = make_wcs(shape, galactic=True) assert wcs.wcs.ctype[0] == 'GLON-CAR' assert wcs.wcs.ctype[1] == 'GLAT-CAR' + + +def test_make_wcs_3D(): + + shape = (1000, 100, 200) + + wcs = make_wcs(shape) + + if astropy_version < '3.1': + assert wcs._naxis1 == shape[2] + assert wcs._naxis2 == shape[1] + else: + assert wcs.pixel_shape == shape[1:][::-1] + + assert wcs.wcs.radesys == 'ICRS' + + wcs = make_wcs(shape, galactic=True) + assert wcs.wcs.ctype[0] == 'GLON-CAR' + assert wcs.wcs.ctype[1] == 'GLAT-CAR' + + +def test_make_wcs_shape(): + """Test if ValueError raises if shape is not 2D or 3D.""" + + with pytest.raises(ValueError): + shape = (100, 200, 1, 1) + wcs = make_wcs(shape) + + +def test_make_imagehdu_wrongshape(): + """Test if ValueError raises if shape is not 3D.""" + + with pytest.raises(ValueError): + data = np.empty(10) + image = make_imagehdu(data) + + +def test_make_cubehdu_wrongshape(): + """Test if ValueError raises if shape is not 3D.""" + + with pytest.raises(ValueError): + shape = (5, 2) + data = np.empty(10).reshape(*shape) + cube = make_cubehdu(data)