Use centroid for aperture photometry

iop4lib/db/aperphotresult.py

@@ -39,6 +39,8 @@ class AperPhotResult(models.Model):
     flux_counts_err = models.FloatField(null=True, blank=True)
 
     ## extra fields
+    x_px = models.FloatField(null=True, help_text="used pixel position of the source in the image, x coordinate.")
+    y_px = models.FloatField(null=True, help_text="used pixel position of the source in the image, y coordinate.")
     fwhm = models.FloatField(null=True, blank=True)
     photopolresults = models.ManyToManyField("PhotoPolResult", related_name='aperphotresults', help_text="The PhotoPolResult(s) this AperPhotResult has been used for.")
     modified = models.DateTimeField(auto_now=True)

iop4lib/instruments/instrument.py

@@ -17,6 +17,7 @@
 import datetime
 import glob
 import astrometry
+from photutils.centroids import centroid_com, centroid_sources
 
 # iop4lib imports
 from iop4lib.enums import *
@@ -518,35 +519,25 @@ def compute_aperture_photometry(cls, redf, aperpix, r_in, r_out):
         from photutils.utils import calc_total_error
         from astropy.stats import SigmaClip
 
-        if redf.mdata.shape[0] == 1024: # andor cameras
-            bkg_box_size = 128
-        elif redf.mdata.shape[0] == 2048: # andor cameras
-            bkg_box_size = 256
-        elif redf.mdata.shape[0] == 800: # cafos
-            bkg_box_size = 100
-        elif redf.mdata.shape[0] == 900: # dipol polarimetry
-            bkg_box_size = 90
-        elif redf.mdata.shape[0] == 896: # dipol polarimetry (also)
-            bkg_box_size = 112
-        elif redf.mdata.shape[0] == 4144: # dipol photometry
-            bkg_box_size = 518
-        else:
-            logger.warning(f"Image size {redf.mdata.shape[0]} not expected.")
-            bkg_box_size = redf.mdata.shape[0]//10
+        bkg_box_size = redf.mdata.shape[0]//10
 
         bkg = get_bkg(redf.mdata, filter_size=1, box_size=bkg_box_size)
         img = redf.mdata
 
-        if np.sum(redf.mdata <= 0.0) >= 1:
-            logger.debug(f"{redf}: {np.sum(img <= 0.0):.0f} px < 0  ({math.sqrt(np.sum(img <= 0.0)):.0f} px2) in image.")
-
         error = calc_total_error(img, bkg.background_rms, cls.gain_e_adu)
 
         for astrosource in redf.sources_in_field.all():
             for pairs, wcs in (('O', redf.wcs1), ('E', redf.wcs2)) if redf.has_pairs else (('O',redf.wcs),):
 
-                ap = CircularAperture(astrosource.coord.to_pixel(wcs), r=aperpix)
-                annulus = CircularAnnulus(astrosource.coord.to_pixel(wcs), r_in=r_in, r_out=r_out)
+                wcs_px_pos = astrosource.coord.to_pixel(wcs)
+
+                logger.debug(f"{redf}: computing aperture photometry for {astrosource}")
+
+                centroid_px_pos = centroid_sources(img, xpos=wcs_px_pos[0], ypos=wcs_px_pos[1], box_size=11, centroid_func=centroid_com)
+                centroid_px_pos = (centroid_px_pos[0][0], centroid_px_pos[1][0])
+
+                ap = CircularAperture(centroid_px_pos, r=aperpix)
+                annulus = CircularAnnulus(centroid_px_pos, r_in=r_in, r_out=r_out)
 
                 annulus_stats = ApertureStats(redf.mdata, annulus, error=error, sigma_clip=SigmaClip(sigma=5.0, maxiters=10))
                 ap_stats = ApertureStats(redf.mdata, ap, error=error)
@@ -561,6 +552,7 @@ def compute_aperture_photometry(cls, redf, aperpix, r_in, r_out):
                                       astrosource=astrosource, 
                                       aperpix=aperpix, 
                                       r_in=r_in, r_out=r_out,
+                                      x_px=centroid_px_pos[0], y_px=centroid_px_pos[1],
                                       pairs=pairs, 
                                       bkg_flux_counts=bkg_flux_counts, bkg_flux_counts_err=bkg_flux_counts_err,
                                       flux_counts=flux_counts, flux_counts_err=flux_counts_err)

iop4lib/utils/__init__.py

@@ -307,7 +307,7 @@ def fit_gaussian(px_start, redf=None, data=None, sigma_start=7, r_max=None, r_se
 def get_angle_from_history(redf: 'ReducedFit' = None, 
                            target_src: 'AstroSource' = None, 
                            calibrated_fits: Sequence['ReducedFit'] = None, 
-                           n_max_fits=20) -> (float, float):
+                           n_max_fits=20) -> Union[float, float]:
     """ Compute the average rotation angle from a list of already calibrated fits.
 
     To compute it, it checks the WCS in the header of each calibrated fit.