fix bug: aperture estimation, position correction (#174)

Author: juanep97

iop4lib/instruments/instrument.py

@@ -520,10 +520,17 @@ def compute_aperture_photometry(cls, redf, aperpix, r_in, r_out):
         from photutils.utils import calc_total_error
         from astropy.stats import SigmaClip
 
-        bkg_box_size = redf.mdata.shape[0]//10
+        img = redf.mdata
+
+        # centroid_sources and fit_fwhm need bkg-substracted data
 
+        # opt 1) estimate bkg using more complex algs
+        bkg_box_size = redf.mdata.shape[0]//10
         bkg = get_bkg(redf.mdata, filter_size=1, box_size=bkg_box_size)
-        img = redf.mdata
+        median = bkg.background_median
+
+        # opt 2) estimate bkg using simple sigma clip
+        # mean, median, std = sigma_clipped_stats(img, sigma=3.0)
 
         error = calc_total_error(img, bkg.background_rms, cls.gain_e_adu)
 
@@ -544,17 +551,31 @@ def compute_aperture_photometry(cls, redf, aperpix, r_in, r_out):
                     logger.warning(f"{redf}: ({pairs}) image of {astrosource.name} is too close to the border, skipping aperture photometry.")
                     continue
 
+                # # correct position using centroid
+                # # choose a box size that is somewhat larger than the aperture
+                # # in case of pairs, cap box size so that it is somewhat smaller than the distance between pairs
+
+                # box_size = math.ceil(1.6 * aperpix)//2 * 2 + 1
+
+                # if redf.has_pairs:
+                #     box_size_max = math.ceil(np.linalg.norm(Instrument.by_name(redf.instrument).disp_sign_mean))//2 * 2 - 1
+                #     box_size = min(box_size, box_size_max)
+
                 # correct position using centroid
-                # choose a box size that is somewhat larger than the aperture
-                # in case of pairs, cap box size that is somewhat smaller than the distance between pairs
+                # choose a box size around 12.0 arcsec (0.4'' is excellent seeing)
+                # in case of pairs, cap box size so that it is somewhat smaller than the distance between pairs
 
-                box_size = math.ceil(1.6 * aperpix)//2 * 2 + 1
+                arcsec_px = redf.pixscale.to(u.arcsec/u.pix).value
+
+                box_size = math.ceil( ( 12 / arcsec_px ) ) // 2 * 2 + 1
+
+                logger.debug(f"{box_size=}")
 
                 if redf.has_pairs:
-                    _box_size_max = math.ceil(np.linalg.norm(Instrument.by_name(redf.instrument).disp_sign_mean))//2 * 2 - 1
-                    box_size = min(box_size, _box_size_max)
+                    box_size_max = math.ceil(np.linalg.norm(Instrument.by_name(redf.instrument).disp_sign_mean))//2 * 2 - 1
+                    box_size = min(box_size, box_size_max)
 
-                centroid_px_pos = centroid_sources(img, xpos=wcs_px_pos[0], ypos=wcs_px_pos[1], box_size=box_size, centroid_func=centroid_2dg)
+                centroid_px_pos = centroid_sources(img-median, xpos=wcs_px_pos[0], ypos=wcs_px_pos[1], box_size=box_size, centroid_func=centroid_2dg)
                 centroid_px_pos = (centroid_px_pos[0][0], centroid_px_pos[1][0])
 
                 # check that the centroid position is within the borders of the image
@@ -580,14 +601,16 @@ def compute_aperture_photometry(cls, redf, aperpix, r_in, r_out):
                 flux_counts = ap_stats.sum - annulus_stats.mean*ap_stats.sum_aper_area.value # TODO: check if i should use mean!
                 flux_counts_err = ap_stats.sum_err
 
-                apres = AperPhotResult.create(reducedfit=redf, 
-                                      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)
+                apres = AperPhotResult.create(
+                    reducedfit=redf, 
+                    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,
+                )
 
                 apres_L.append(apres)
 
@@ -709,8 +732,10 @@ def estimate_common_apertures(cls, reducedfits, reductionmethod=None, fit_boxsiz
         It fits the target source profile in the fields and returns some multiples of the fwhm which are used as the aperture and as the inner and outer radius of the annulus for local bkg estimation).
         """
 
-        from iop4lib.utils import fit_gaussian
-        from iop4lib.db import AstroSource
+        from photutils.psf import fit_fwhm
+        from iop4lib.utils.sourcedetection import get_bkg
+
+        # select the target source from the list of files (usually, there is only one)
 
         astrosource_S = set.union(*[set(redf.sources_in_field.all()) for redf in reducedfits])
         target_L = [astrosource for astrosource in astrosource_S if not astrosource.is_calibrator]
@@ -723,30 +748,78 @@ def estimate_common_apertures(cls, reducedfits, reductionmethod=None, fit_boxsiz
             target = astrosource_S.pop()
         else:
             return np.nan, np.nan, np.nan, {'mean_fwhm':np.nan, 'sigma':np.nan}
+
+        logger.debug(f"{target=}")
+
+        # Iterate over each file, get the fwhm of the source
 
         fwhm_L = list()
 
         for redf in reducedfits:
+
             try:
-                gaussian = fit_gaussian(px_start=redf.wcs.world_to_pixel(target.coord), redf=redf)
-                fwhm = (2*np.sqrt(2*math.log(2))) * np.sqrt(gaussian[0].x_stddev.value**2+gaussian[0].y_stddev.value**2)
-                if not (fwhm_min < fwhm < fwhm_max):
-                    logger.warning(f"ReducedFit {redf.id} {target.name}: fwhm = {fwhm} px, skipping this reduced fit")
-                    continue
-                logger.debug(f"ReducedFit {redf.id} [{target.name}] Gaussian FWHM: {fwhm:.1f} px")
-                fwhm_L.append(fwhm)
+
+                img = redf.mdata
+
+                # centroid_sources and fit_fwhm need bkg-substracted data
+
+                # opt 1) estimate bkg using more complex algs
+                bkg_box_size = redf.mdata.shape[0]//10
+                bkg = get_bkg(redf.mdata, filter_size=1, box_size=bkg_box_size)
+                median = bkg.background_median
+
+                # opt 2) estimate bkg using simple sigma clip
+                # mean, median, std = sigma_clipped_stats(img, sigma=3.0)
+
+                for pairs, wcs in (('O', redf.wcs1), ('E', redf.wcs2)) if redf.has_pairs else (('O',redf.wcs),):
+
+                    wcs_px_pos = target.coord.to_pixel(wcs)
+
+                    # correct position using centroid
+                    # choose a box size around 12.0 arcsec (0.4'' is excellent seeing)
+                    # in case of pairs, cap box size so that it is somewhat smaller than the distance between pairs
+
+                    arcsec_px = redf.pixscale.to(u.arcsec/u.pix).value
+
+                    box_size = math.ceil( ( 12 / arcsec_px ) ) // 2 * 2 + 1
+
+                    logger.debug(f"{box_size=}")
+
+                    if redf.has_pairs:
+                        box_size_max = math.ceil(np.linalg.norm(Instrument.by_name(redf.instrument).disp_sign_mean))//2 * 2 - 1
+                        box_size = min(box_size, box_size_max)
+
+                    centroid_px_pos = centroid_sources(img-median, xpos=wcs_px_pos[0], ypos=wcs_px_pos[1], box_size=box_size, centroid_func=centroid_2dg)
+                    centroid_px_pos = (centroid_px_pos[0][0], centroid_px_pos[1][0])
+
+                    # log the difference between the WCS and the centroid
+                    wcs_diff = np.sqrt((centroid_px_pos[0] - wcs_px_pos[0])**2 + (centroid_px_pos[1] - wcs_px_pos[1])**2)
+                    
+                    logger.debug(f"ReducedFit {redf.id}: {target.name} {pairs}: WCS centroid distance = {wcs_diff:.1f} px")
+
+                    fwhm = fit_fwhm(img-median, xypos=centroid_px_pos, fit_shape=7)[0]
+
+                    if not (fwhm_min < fwhm < fwhm_max):
+                        logger.warning(f"ReducedFit {redf.id} {target.name}: fwhm = {fwhm} px, skipping this reduced fit")
+                        continue
+
+                    logger.debug(f"ReducedFit {redf.id} [{target.name}] Gaussian FWHM: {fwhm} px")
+
+                    fwhm_L.append(fwhm)
+
             except Exception as e:
                 logger.warning(f"ReducedFit {redf.id} {target.name}: error in gaussian fit, skipping this reduced fit ({e})")
                 continue
 
         if len(fwhm_L) > 0:
             mean_fwhm = np.mean(fwhm_L)
         else:
-            logger.error(f"Could not find an appropriate aperture for Reduced Fits {[redf.id for redf in reducedfits]}, using standard fwhm of 3.5px")
+            logger.error(f"Could not find an appropriate aperture for Reduced Fits {[redf.id for redf in reducedfits]}, using standard fwhm of {fwhm_default} px")
             mean_fwhm = fwhm_default
 
         sigma = mean_fwhm / (2*np.sqrt(2*math.log(2)))
 
+        
         return 3.0*sigma, 5.0*sigma, 9.0*sigma, {'mean_fwhm':mean_fwhm, 'sigma':sigma}