Add convenience function to estimate local backgrounds

[larrybradley]

This would be for aperture photometry using an annulus-type aperture. The function would perform the background estimation, subtraction from the aperture photometry, and error propagation.

[Vb2341]

Since this is an important topic for HST users/the instrument teams, I have went ahead and created some functions that are basically what are described here. I discovered today, however, that you've also recently written the code (for the AAS231 workshop example notebooks) to some of the local background estimation using functions that aren't just the mean of an annulus.

I've went slightly further and implemented some of IRAF/DAOPHOT's background/error computation here using your existing API tools. The outputs here also turn the measurements into magnitudes/mag errors for convenience.

One thing I think is worth noting is that if the pixels inside the annulus are to be sigma clipped, then it seems using astropy.stats.sigma_clipped_stats is quite slow compared to scipy.stats.sigmaclip, see astropy/astropy#7086. In my implementation, I simply extract the pixels of interest before computing the statistics, which appears to be quite a bit faster. That might be a bit more complicated with masked arrays, which sigma_clipped_stats would already handle, so I'm not sure which is the better way forward. I guess the other advantage of using sigmaclip and computing the stats yourself is that it's quite easy to get the standard deviation, which is necessary to account for error in the sky value.

Hope this can be helpful in some way!

[yoachim]

Ah, I was looking to see if there was a way to get medians in an aperture. Would be nice to have these folded back into photutils.

[larrybradley]

@yoachim You can easily get the median of the pixel values in an aperture by converting the aperture to an ApertureMask object and then applying the mask to the data:

from photutils import CircularAperture

aper = CircularAperture((50, 50), r=10)
mask = aper.to_mask(method='center')    # list of ApertureMask

data = np.random.random((100, 100))
aper_data = mask[0].multiply(data)
aper_data_1d = aper_data[mask[0].data > 0]
np.median(aper_data_1d)

You need to use method='center' in the to_mask() method to extract full pixels. The other methods produce partial-pixel masks, which then require the use of weighted statistics (which you can do).

A tutorial notebook is here, which includes an example of sigma-clipped medians in annulus apertures for local background subtraction: https://github.com/astropy/astropy-workshop/blob/master/09-Photutils/photutils_local_backgrounds.ipynb

[yoachim]

Great example notebook, thanks!

[larrybradley]

Local background are now available in the ApertureStats class: https://photutils.readthedocs.io/en/stable/api/photutils.aperture.ApertureStats.html#photutils.aperture.ApertureStats