Merge branch 'develop' of https://github.com/roman-corgi/corgidrp int…

…o improve_err

corgidrp/fluxcal.py

@@ -8,6 +8,7 @@
 from astropy.coordinates import SkyCoord
 import corgidrp
 from photutils.aperture import CircularAperture
+from photutils.background import LocalBackground
 from photutils.psf import fit_2dgaussian
 from scipy import integrate
 import urllib
@@ -244,11 +245,11 @@ def calculate_band_irradiance(filter_curve, calspec_flux, filter_wavelength):
 
     return irrad
 
-def aper_phot(image, encircled_radius, frac_enc_energy = 1., method = 'exact', subpixels = 5):
+def aper_phot(image, encircled_radius, frac_enc_energy = 1., method = 'exact', subpixels = 5, background_sub = False, r_in = 5 , r_out = 10):
     """
     returns the flux in photo-electrons of a point source at the target Ra/Dec position
     and using a circular aperture by applying aperture_photometry of photutils.
-    We assume that background subtraction is already done.
+    Background subtraction can be done optionally using a user defined circular annulus.
     
     Args:
         image (corgidrp.data.Image): combined source exposure image
@@ -258,29 +259,42 @@ def aper_phot(image, encircled_radius, frac_enc_energy = 1., method = 'exact', s
         default is 'exact'. For detailed description see https://photutils.readthedocs.io/en/stable/api/photutils.aperture.CircularAnnulus.html
         subpixels (int): For the 'subpixel' method, resample pixels by this factor in each dimension. That is, each pixel is divided 
                          into subpixels**2 subpixels. This keyword is ignored unless method='subpixel', default is 5
+        background_sub (boolean): background can be determine from a circular annulus (default: False).
+        r_in (float): inner radius of circular annulus in pixels, (default: 5)
+        r_out (float): outer radius of circular annulus in pixels, (default: 10)
     
     Returns:
-        float: integrated flux of the point source in unit photo-electrons and corresponding error
+        list: integrated flux of the point source in unit photo-electrons and corresponding error and optional local background value
     """
     #calculate the x and y pixel positions using the RA/Dec target position and applying WCS conversion
     if frac_enc_energy <=0 or frac_enc_energy > 1:
         raise ValueError("frac_enc_energy {0} should be within 0 < fee <= 1".format(str(frac_enc_energy)))
     ra = image.pri_hdr['RA']
     dec = image.pri_hdr['DEC']
+    dat = image.data.copy()
 
     target_skycoord = SkyCoord(ra = ra, dec = dec, unit='deg')
     w = wcs.WCS(image.ext_hdr)
     pix = wcs.utils.skycoord_to_pixel(target_skycoord, w, origin = 1)
+    if background_sub:
+        #This is essentially the median in a circular annulus 
+        bkg = LocalBackground(r_in, r_out)
+        back = bkg(dat, pix[0], pix[1], mask = image.dq.astype(bool))
+        dat -= back
+
     aper = CircularAperture(pix, encircled_radius)
     aperture_sums, aperture_sums_errs = \
-        aper.do_photometry(image.data, error = image.err[0], mask = image.dq.astype(bool), method = method, subpixels = subpixels)
-    return aperture_sums[0]/frac_enc_energy, aperture_sums_errs[0]/frac_enc_energy
+        aper.do_photometry(dat, error = image.err[0], mask = image.dq.astype(bool), method = method, subpixels = subpixels)
+    if background_sub:
+        return [aperture_sums[0]/frac_enc_energy, aperture_sums_errs[0]/frac_enc_energy, back]
+    else:
+        return [aperture_sums[0]/frac_enc_energy, aperture_sums_errs[0]/frac_enc_energy]
 
-def phot_by_gauss2d_fit(image, fwhm, fit_shape = None):
+def phot_by_gauss2d_fit (image, fwhm, fit_shape = None, background_sub = False, r_in = 5 , r_out = 10):
     """
     returns the flux in photo-electrons of a point source at the target Ra/Dec position
     and using a circular aperture by applying aperture_photometry of photutils
-    We assume that background subtraction is already done.
+    Background subtraction can be done optionally using a user defined circular annulus.
     
     Args:
         image (corgidrp.data.Image): combined source exposure image
@@ -289,34 +303,48 @@ def phot_by_gauss2d_fit(image, fwhm, fit_shape = None):
             The shape of the fitting region. If a scalar, then it is assumed
             to be a square. If `None`, then the shape of the input ``data``.
             It must be an odd value and should be much bigger than fwhm.
+            background_sub (boolean): background can be determine from a circular annulus (default: False).
+            r_in (float): inner radius of circular annulus in pixels, (default: 5)
+            r_out (float): outer radius of circular annulus in pixels, (default: 10)
     
     Returns:
-        float: integrated flux of the Gaussian2d fit of the point source in unit photo-electrons and corresponding error
+        list: integrated flux of the Gaussian2d fit of the point source in unit photo-electrons and corresponding error and optional local background value
     """
     ra = image.pri_hdr['RA']
     dec = image.pri_hdr['DEC']
+    dat = image.data.copy()
 
     target_skycoord = SkyCoord(ra = ra, dec = dec, unit='deg')
     w = wcs.WCS(image.ext_hdr)
     pix = wcs.utils.skycoord_to_pixel(target_skycoord, w, origin = 1)
+    if background_sub:
+        #This is essentially the median in a circular annulus 
+        bkg = LocalBackground(r_in, r_out)
+        back = bkg(dat, pix[0], pix[1], mask = image.dq.astype(bool))
+        dat -= back
+
     # fit_2dgaussian: error weighting raises exception if error is zero
     err = image.err[0]
     err[err == 0] = np.finfo(np.float32).eps
 
     if fit_shape == None:
-        fit_shape = np.shape(image.data)[0] -1
+        fit_shape = np.shape(dat)[0] -1
 
-    psf_phot = fit_2dgaussian(image.data, xypos = pix, fwhm = fwhm, fit_shape = fit_shape, mask = image.dq.astype(bool), error = err)
+    psf_phot = fit_2dgaussian(dat, xypos = pix, fwhm = fwhm, fit_shape = fit_shape, mask = image.dq.astype(bool), error = err)
     flux = psf_phot.results['flux_fit'][0]
     flux_err = psf_phot.results['flux_err'][0]
-    return flux, flux_err
+    if background_sub:
+        return [flux, flux_err, back]
+    else:
+        return [flux, flux_err]
 
-def calibrate_fluxcal_aper(dataset_or_image, encircled_radius, frac_enc_energy = 1., method = 'exact', subpixels = 5):
+def calibrate_fluxcal_aper(dataset_or_image, encircled_radius, frac_enc_energy = 1., method = 'exact', subpixels = 5, background_sub = False, r_in = 5 , r_out = 10):
     """
     fills the FluxcalFactors calibration product values for one filter band,
     calculates the flux calibration factors by aperture photometry.
     The band flux values are divided by the found photoelectrons.
     Propagates also errors to flux calibration factor calfile.
+    Background subtraction can be done optionally using a user defined circular annulus.
     
     Args:
         dataset_or_image (corgidrp.data.Dataset or corgidrp.data.Image): dataset with one image or image as combined source exposure
@@ -326,6 +354,9 @@ def calibrate_fluxcal_aper(dataset_or_image, encircled_radius, frac_enc_energy =
         default is 'exact'. For detailed description see https://photutils.readthedocs.io/en/stable/api/photutils.aperture.CircularAnnulus.html
         subpixels (int): For the 'subpixel' method, resample pixels by this factor in each dimension. That is, each pixel is divided 
                          into subpixels**2 subpixels. This keyword is ignored unless method='subpixel', default is 5
+        background_sub (boolean): background can be determine from a circular annulus (default: False).
+        r_in (float): inner radius of circular annulus in pixels, (default: 5)
+        r_out (float): outer radius of circular annulus in pixels, (default: 10)
     
     Returns:
         corgidrp.data.FluxcalFactor: FluxcalFactor calibration object with the value and error of the corresponding filter
@@ -348,22 +379,27 @@ def calibrate_fluxcal_aper(dataset_or_image, encircled_radius, frac_enc_energy =
     # calculate the flux from the user given CALSPEC file binned on the wavelength grid of the filter
     flux_ref = read_cal_spec(calspec_filepath, wave)
     flux = calculate_band_flux(filter_trans, flux_ref, wave)
-
-    ap_sum, ap_sum_err = aper_phot(image, encircled_radius, frac_enc_energy, method = method, subpixels = subpixels)
-
+    aper_phot_result = aper_phot(image, encircled_radius, frac_enc_energy, method = method, subpixels = subpixels, background_sub = background_sub, r_in = r_in, r_out = r_out)
+    ap_sum = aper_phot_result[0]
+    ap_sum_err = aper_phot_result[1]
+    exthdr = image.ext_hdr
+    if background_sub:
+        exthdr['LOCBACK'] = aper_phot_result[2]
     fluxcal_fac = flux/ap_sum
     fluxcal_fac_err = flux/ap_sum**2 * ap_sum_err
 
-    fluxcal = corgidrp.data.FluxcalFactor(np.array([[fluxcal_fac]]), err = np.array([[[fluxcal_fac_err]]]), pri_hdr = image.pri_hdr, ext_hdr = image.ext_hdr, input_dataset = dataset)
+    exthdr['HISTORY'] = "Flux calibration factor was determined by aperture photometry"
+    fluxcal = corgidrp.data.FluxcalFactor(np.array([[fluxcal_fac]]), err = np.array([[[fluxcal_fac_err]]]), pri_hdr = image.pri_hdr, ext_hdr = exthdr, input_dataset = dataset)
 
     return fluxcal
 
-def calibrate_fluxcal_gauss2d(dataset_or_image, fwhm, fit_shape = None):
+def calibrate_fluxcal_gauss2d(dataset_or_image, fwhm, fit_shape = None, background_sub = False, r_in = 5 , r_out = 10):
     """
     fills the FluxcalFactors calibration product values for one filter band,
     calculates the flux calibration factors by fitting a 2D Gaussian.
     The band flux values are divided by the found photoelectrons.
     Propagates also errors to flux calibration factor calfile.
+    Background subtraction can be done optionally using a user defined circular annulus.
     
     Args:
         dataset_or_image (corgidrp.data.Dataset or corgidrp.data.Image): dataset with one image or image as combined source exposure
@@ -372,7 +408,10 @@ def calibrate_fluxcal_gauss2d(dataset_or_image, fwhm, fit_shape = None):
             The shape of the fitting region. If a scalar, then it is assumed
             to be a square. If `None`, then the shape of the input ``data``.
             It must be an odd value and should be much bigger than fwhm.
-    
+        background_sub (boolean): background can be determine from a circular annulus (default: False).
+        r_in (float): inner radius of circular annulus in pixels, (default: 5)
+        r_out (float): outer radius of circular annulus in pixels, (default: 10)
+        
     Returns:
         corgidrp.data.FluxcalFactor: FluxcalFactor calibration object with the value and error of the corresponding filter
     """
@@ -394,11 +433,15 @@ def calibrate_fluxcal_gauss2d(dataset_or_image, fwhm, fit_shape = None):
     flux_ref = read_cal_spec(calspec_filepath, wave)
     flux = calculate_band_flux(filter_trans, flux_ref, wave)
 
-    flux_sum, flux_sum_err = phot_by_gauss2d_fit(image, fwhm, fit_shape = fit_shape)
+    flux_sum = phot_by_gauss2d_fit(image, fwhm, fit_shape = fit_shape, background_sub = background_sub, r_in = r_in, r_out = r_out)
 
-    fluxcal_fac = flux/flux_sum
-    fluxcal_fac_err = flux/flux_sum**2 * flux_sum_err
+    fluxcal_fac = flux/flux_sum[0]
+    fluxcal_fac_err = flux/flux_sum[0]**2 * flux_sum[1]
 
-    fluxcal = corgidrp.data.FluxcalFactor(np.array([[fluxcal_fac]]), err = np.array([[[fluxcal_fac_err]]]), pri_hdr = image.pri_hdr, ext_hdr = image.ext_hdr, input_dataset = dataset)
+    exthdr = image.ext_hdr
+    if background_sub:
+        exthdr['LOCBACK'] = flux_sum[2]
+    exthdr['HISTORY'] = "Flux calibration factor was determined by a Gaussian 2D fit photometry"
+    fluxcal = corgidrp.data.FluxcalFactor(np.array([[fluxcal_fac]]), err = np.array([[[fluxcal_fac_err]]]), pri_hdr = image.pri_hdr, ext_hdr = exthdr, input_dataset = dataset)
 
     return fluxcal

corgidrp/mocks.py

@@ -1860,10 +1860,10 @@ def add_defect(sch_imgs, prob, ori, temp):
                     hdul.writeto(filename, overwrite = True)
 
 
-def create_flux_image(star_flux, fwhm, cal_factor, filedir=None, color_cor = 1., platescale=21.8, add_gauss_noise=True, noise_scale=1., file_save=False):
+def create_flux_image(star_flux, fwhm, cal_factor, filedir=None, color_cor = 1., platescale=21.8, add_gauss_noise=True, noise_scale=1., background = 0., file_save=False):
     """
     Create simulated data for absolute flux calibration. This is a point source in the image center with a 2D-Gaussian PSF
-    and Gaussian noise
+    and Gaussian noise and a background.
 
     Args:
         star_flux (float): flux of point source in erg/(s*cm^2*AA)
@@ -1874,6 +1874,7 @@ def create_flux_image(star_flux, fwhm, cal_factor, filedir=None, color_cor = 1.,
         platescale (float): The plate scale of the created image data (default: 21.8 [mas/pixel])
         add_gauss_noise (boolean): Argument to determine if Gaussian noise should be added to the data (default: True)
         noise_scale (float): spread of the Gaussian noise
+        background (float): optional additive background value
         file_save (boolean): save the simulated Image or not (default: False)
 
     Returns:
@@ -1943,6 +1944,8 @@ def create_flux_image(star_flux, fwhm, cal_factor, filedir=None, color_cor = 1.,
     # inject the star into the image
     sim_data[ymin:ymax + 1, xmin:xmax + 1] += psf
 
+    #add a background
+    sim_data += background
     if add_gauss_noise:
         # add Gaussian random noise
         noise_rng = np.random.default_rng(10)

tests/test_fluxcal.py

@@ -59,6 +59,7 @@ def test_flux_calc():
 def test_colorcor():
     """
     test that the pivot reference wavelengths is close to the center of the bandpass
+    and the color correction is calculated correctly
     """
 
     lambda_piv = fluxcal.calculate_pivot_lambda(transmission, wave)
@@ -164,13 +165,16 @@ def test_abs_fluxcal():
     sigma = fwhm/(2.*np.sqrt(2*np.log(2)))
     radius = 3.* sigma
 
+    #Test the aperture photometry
     #The error of one pixel is 1, so the error of the aperture sum should be: 
     error_sum = np.sqrt(np.pi * radius * radius)
-    flux_el_ap, flux_err_ap = fluxcal.aper_phot(flux_image, radius, 0.997)
+    [flux_el_ap, flux_err_ap] = fluxcal.aper_phot(flux_image, radius, 0.997)
     #200 is the input count of photo electrons
     assert flux_el_ap == pytest.approx(200, rel = 0.05)
     assert flux_err_ap == pytest.approx(error_sum, rel = 0.05)
 
+
+    #Test the 2D Gaussian fit photometry
     #The error of one pixel is 1, so the error of a circular aperture with radius of 2 sigma should be about:
     error_gauss = np.sqrt(np.pi * 4 * sigma * sigma)
     flux_el_gauss, flux_err_gauss = fluxcal.phot_by_gauss2d_fit(flux_image, fwhm, fit_shape = 41)
@@ -190,12 +194,12 @@ def test_abs_fluxcal():
     err_fluxcal_ap = band_flux/flux_el_ap**2*flux_err_ap
     assert fluxcal_factor.fluxcal_err == pytest.approx(err_fluxcal_ap)
     assert fluxcal_factor.filename == 'sim_fluxcal_FluxcalFactor_3C_ND475.fits'
-    fluxcal_factor = fluxcal.calibrate_fluxcal_gauss2d(dataset, fwhm)
-    assert fluxcal_factor.filter == '3C'
-    assert fluxcal_factor.fluxcal_fac == pytest.approx(cal_factor,rel = 0.05)
+    fluxcal_factor_gauss = fluxcal.calibrate_fluxcal_gauss2d(dataset, fwhm)
+    assert fluxcal_factor_gauss.filter == '3C'
+    assert fluxcal_factor_gauss.fluxcal_fac == pytest.approx(cal_factor,rel = 0.05)
     #divisive error propagation of the 2D Gaussian fit phot error
     err_fluxcal_gauss = band_flux/flux_el_gauss**2*flux_err_gauss
-    assert fluxcal_factor.fluxcal_err == pytest.approx(err_fluxcal_gauss)
+    assert fluxcal_factor_gauss.fluxcal_err == pytest.approx(err_fluxcal_gauss)
 
     #test the flux conversion computation.
     old_ind = corgidrp.track_individual_errors
@@ -222,8 +226,34 @@ def test_abs_fluxcal():
     assert output_dataset[0].err[0,512, 512] == pytest.approx(err_comb)
     assert output_dataset[0].err[1,512, 512] == pytest.approx(err_layer2)
 
+    #Test the optional background subtraction#
+    flux_image_back = create_flux_image(band_flux, fwhm, cal_factor, background = 3, filedir=datadir, file_save=True)
+
+    [flux_back, flux_err_back, back] = fluxcal.aper_phot(flux_image_back, radius, background_sub = True)
+    #calculated median background should be close to the input
+    assert back == pytest.approx(3, abs = 0.03)
+    #the found values should be close to the ones without background subtraction
+    assert flux_back == pytest.approx(flux_el_ap, abs = 1)
+    assert flux_err_back == pytest.approx(flux_err_ap, abs = 0.03)
+
+    [flux_back, flux_err_back, back] = fluxcal.phot_by_gauss2d_fit(flux_image_back, fwhm, background_sub = True)
+    #calculated median background should be close to the input
+    assert back == pytest.approx(3, abs = 0.03)
+    #the found values should be close to the ones without background subtraction
+    assert flux_back == pytest.approx(flux_el_gauss, abs = 1)
+    assert flux_err_back == pytest.approx(flux_err_gauss, abs = 0.03)
+
+    #Also test again the generation of the cal file now with a background subtraction
+    fluxcal_factor_back = fluxcal.calibrate_fluxcal_aper(flux_image_back, radius, 0.997, background_sub = True)
+    assert fluxcal_factor_back.fluxcal_fac == pytest.approx(fluxcal_factor.fluxcal_fac)
+    assert fluxcal_factor_back.ext_hdr["LOCBACK"] == back
+    fluxcal_factor_back_gauss = fluxcal.calibrate_fluxcal_gauss2d(flux_image_back, fwhm, background_sub = True)
+    assert fluxcal_factor_back_gauss.fluxcal_fac == pytest.approx(fluxcal_factor_gauss.fluxcal_fac)
+    assert fluxcal_factor_back_gauss.ext_hdr["LOCBACK"] == back
+
     corgidrp.track_individual_errors = old_ind
 
+
 if __name__ == '__main__':
     test_get_filter_name()
     test_flux_calc()