From 8b2dab72351eb2a0b8ec7bd9e62d2ef212f7d0fa Mon Sep 17 00:00:00 2001 From: Madelyn Madsen Date: Wed, 22 Jul 2020 13:11:06 -0500 Subject: [PATCH 1/3] Debug aij photometry --- stellarphot/photometry.py | 5 +++-- stellarphot/source_detection.py | 7 ++++--- 2 files changed, 7 insertions(+), 5 deletions(-) diff --git a/stellarphot/photometry.py b/stellarphot/photometry.py index 578eb94d..23955a32 100644 --- a/stellarphot/photometry.py +++ b/stellarphot/photometry.py @@ -517,12 +517,12 @@ def calculate_noise(gain=1.0, read_noise=0.0, dark_current_per_sec=0.0, exposure=0, include_digitization=False): - if annulus_area == 0: + if annulus_area.all() == 0: area_ratio = aperture_area else: area_ratio = aperture_area * (1 + aperture_area / annulus_area) - poisson_source = gain * flux + poisson_source = gain * flux / u.adu sky = area_ratio * gain * sky_per_pix dark = area_ratio * dark_current_per_sec * exposure @@ -533,4 +533,5 @@ def calculate_noise(gain=1.0, read_noise=0.0, dark_current_per_sec=0.0, if include_digitization: digitization = area_ratio * (gain * 0.289) ** 2 + return np.sqrt(poisson_source + sky + dark + rn_error + digitization) diff --git a/stellarphot/source_detection.py b/stellarphot/source_detection.py index 667f6608..f4954d31 100644 --- a/stellarphot/source_detection.py +++ b/stellarphot/source_detection.py @@ -5,17 +5,18 @@ from photutils import fit_2dgaussian from astropy.nddata.utils import Cutout2D from astropy.stats import gaussian_sigma_to_fwhm +from astropy import units as u __all__ = ['source_detection', 'compute_fwhm'] def compute_fwhm(ccd, sources, fwhm_estimate=5, - x_column='xcentroid', y_column='ycentroid'): + x_column='xcenter', y_column='ycenter'): fwhm_x = [] fwhm_y = [] for source in sources: - x = source[x_column] - y = source[y_column] + x = source[x_column] / u.pixel + y = source[y_column] / u.pixel cutout = Cutout2D(ccd, (x, y), 5 * fwhm_estimate) fit = fit_2dgaussian(cutout.data) fwhm_x.append(gaussian_sigma_to_fwhm * fit.x_stddev) From 4c880779fe24a4a340ea2fc2993ab48d6cf2d31d Mon Sep 17 00:00:00 2001 From: Madelyn Madsen Date: Wed, 22 Jul 2020 13:16:05 -0500 Subject: [PATCH 2/3] Remove blank line --- stellarphot/photometry.py | 1 - 1 file changed, 1 deletion(-) diff --git a/stellarphot/photometry.py b/stellarphot/photometry.py index 23955a32..07e40e94 100644 --- a/stellarphot/photometry.py +++ b/stellarphot/photometry.py @@ -533,5 +533,4 @@ def calculate_noise(gain=1.0, read_noise=0.0, dark_current_per_sec=0.0, if include_digitization: digitization = area_ratio * (gain * 0.289) ** 2 - return np.sqrt(poisson_source + sky + dark + rn_error + digitization) From 90023c88b8be0961d29939fe4b05c2c7bc2af1bc Mon Sep 17 00:00:00 2001 From: Matt Craig Date: Wed, 22 Jul 2020 15:15:11 -0500 Subject: [PATCH 3/3] Make unit workarounds more robust This is really a temporary workaround until we get something sensible in place for units. --- stellarphot/photometry.py | 15 +++++++++++++-- stellarphot/source_detection.py | 17 +++++++++++++--- stellarphot/tests/test_detection.py | 30 ++++++++++++++++++++++------- 3 files changed, 50 insertions(+), 12 deletions(-) diff --git a/stellarphot/photometry.py b/stellarphot/photometry.py index 07e40e94..49f73d56 100644 --- a/stellarphot/photometry.py +++ b/stellarphot/photometry.py @@ -517,12 +517,23 @@ def calculate_noise(gain=1.0, read_noise=0.0, dark_current_per_sec=0.0, exposure=0, include_digitization=False): - if annulus_area.all() == 0: + try: + no_annulus = (annulus_area == 0).all() + except AttributeError: + no_annulus = annulus_area == 0 + + if no_annulus: area_ratio = aperture_area else: area_ratio = aperture_area * (1 + aperture_area / annulus_area) - poisson_source = gain * flux / u.adu + # Units in here are an absolute mess + poisson_source = gain * flux + + try: + poisson_source = poisson_source.value + except AttributeError: + pass sky = area_ratio * gain * sky_per_pix dark = area_ratio * dark_current_per_sec * exposure diff --git a/stellarphot/source_detection.py b/stellarphot/source_detection.py index f4954d31..3884da21 100644 --- a/stellarphot/source_detection.py +++ b/stellarphot/source_detection.py @@ -15,8 +15,17 @@ def compute_fwhm(ccd, sources, fwhm_estimate=5, fwhm_x = [] fwhm_y = [] for source in sources: - x = source[x_column] / u.pixel - y = source[y_column] / u.pixel + x = source[x_column] + y = source[y_column] + + # Cutout2D needs no units on the center position, so remove unit + # if it is present. + try: + x = x.value + y = y.value + except AttributeError: + pass + cutout = Cutout2D(ccd, (x, y), 5 * fwhm_estimate) fit = fit_2dgaussian(cutout.data) fwhm_x.append(gaussian_sigma_to_fwhm * fit.x_stddev) @@ -63,7 +72,9 @@ def source_detection(ccd, fwhm=8, sigma=3.0, iters=5, mean, median, std = sigma_clipped_stats(ccd, sigma=sigma, maxiters=iters) daofind = DAOStarFinder(fwhm=fwhm, threshold=threshold * std) sources = daofind(ccd - median) + print(sources) if find_fwhm: - x, y = compute_fwhm(ccd, sources, fwhm_estimate=fwhm) + x, y = compute_fwhm(ccd, sources, fwhm_estimate=fwhm, + x_column='xcentroid', y_column='ycentroid') sources['FWHM'] = (x + y) / 2 return sources diff --git a/stellarphot/tests/test_detection.py b/stellarphot/tests/test_detection.py index 5a8f13b7..68d7a33f 100644 --- a/stellarphot/tests/test_detection.py +++ b/stellarphot/tests/test_detection.py @@ -4,8 +4,11 @@ from astropy.utils.data import get_pkg_data_filename from astropy.table import Table from astropy.stats import gaussian_sigma_to_fwhm +from astropy import units as u + from photutils.datasets import make_gaussian_sources_image, make_noise_image -from stellarphot.source_detection import source_detection + +from stellarphot.source_detection import source_detection, compute_fwhm from stellarphot.photometry import photutils_stellar_photometry @@ -31,6 +34,23 @@ def image(self): return self._stars + self._noise +@pytest.mark.parametrize('units', [u.pixel, None]) +def test_compute_fwhm(units): + fake_image = FakeImage() + sources = fake_image.sources + if units is not None: + # It turns out having a unit on a column is not the same as + # things in the column having units. The construct below ensures + # that the source table values have units. + # Do not try: sources['x_mean'] = sources['x_mean'] * units + # Turns out individual values do NOT have units in that case. + sources['x_mean'] = [v * units for v in sources['x_mean']] + sources['y_mean'] = [v * units for v in sources['y_mean']] + + fwhm_x, fwhm_y = compute_fwhm(fake_image.image, sources, + x_column='x_mean', y_column='y_mean') + + def test_detect_source_number_location(): """ Make sure we detect the sources in the input table.... @@ -40,12 +60,10 @@ def test_detect_source_number_location(): found_sources = source_detection(fake_image.image, fwhm=2 * sources['x_stddev'].mean(), threshold=10) - print(found_sources.colnames) # Sort by flux so we can reliably match them sources.sort('amplitude') found_sources.sort('flux') - print(found_sources) # Do we have the right number of sources? assert len(sources) == len(found_sources) @@ -74,8 +92,7 @@ def test_aperture_photometry_no_outlier_rejection(): phot.sort('aperture_sum') sources.sort('amplitude') found_sources.sort('flux') - print(found_sources['xcentroid']) - print(phot['background_per_pixel', 'aperture_sum_err', 'net_flux', 'aperture_sum']) + for inp, out in zip(sources, phot): stdev = inp['x_stddev'] expected_flux = (inp['amplitude'] * 2 * np.pi * @@ -131,8 +148,7 @@ def test_aperture_photometry_with_outlier_rejection(reject): phot.sort('aperture_sum') sources.sort('amplitude') found_sources.sort('flux') - print(found_sources['xcentroid']) - print(phot['background_per_pixel', 'aperture_sum_err', 'net_flux', 'aperture_sum']) + for inp, out in zip(sources, phot): stdev = inp['x_stddev'] expected_flux = (inp['amplitude'] * 2 * np.pi *