diff --git a/docs/reduction.rst b/docs/reduction.rst
index 46c5336..2787c8a 100644
--- a/docs/reduction.rst
+++ b/docs/reduction.rst
@@ -357,6 +357,55 @@ it shows you how to load in the rather human-unreadable |hiper| log
files (huge numbers of columns and rows). While observing |plog| is
helpful for plotting the sky level as twilight approaches.
+PSF Photometry
+==============
+In crowded fields, aperture photometry can be difficult. The effect
+of nearby stars on the sky background estimates can be minimised by
+masking the nearby stars (as shown in |fig-setaper|), but this does
+not solve the problem of target apertures containing flux from more
+than one object. For crowded fields where this is an issue PSF
+photometry is often the best solution.
+
+|reduce| can carry out PSF photometry, using `photutils `_
+as the underlying engine. This is an optional install, so unless you
+installed the pipeline with the `psf` extra you will need to install
+photutils separately. To use PSF photometry, all one needs to do is
+set the `use_psf` parameter in the reduce file to `yes`.
+
+.. Note::
+ Because PSF photometry uses photutils as the underlying engine,
+ it is not possible to carry out PSF photometry with apertures
+ that lie in different windows of the CCD. This is because we
+ need to send photutils a contiguous image to work with,
+ and this is not possible if the apertures lie in different windows.
+
+In PSF photometry, :ref:`aperture positioning `
+is carried out in exactly the same way as for aperture photometry
+to find the positions of the stars in the frame. The PSF model
+is then fit to the reference stars to find the shape of the PSF.
+The PSF model is then fit to all stars defined in the aperture file
+with the shape held fixed to find the fluxes. When running |reduce|
+you can display the residual image after the PSF model has been
+subtracted from the data to check the quality of the fit.
+
+Optionally you can let the positions of the stars vary from the locations
+found by the aperture positioning step, by setting the `positions` parameter
+in the `psf_photom` section of the reduction file to `variable`. This
+repositioning does not respect :ref:`linked apertures `
+and usually gives poorer results, so it is not recommended.
+
+Creating aperture files for the crowded fields where PSF photometry is needed
+can be tricky. Overlapping sources can make using |setaper| difficult or
+impossible. In this case, it is best to use the |psfaper| script.
+This script allows you to zoom into the region where you want to
+carry out PSF photometry and add a small number of well-isolated bright
+reference stars. The script then fits a PSF model to the reference stars
+to constrain the shape of the PSF and carries out multiple iterations
+of a FIND-FIT-SUBTRACT loop to find all the stars within a region of interest.
+The aperture file created with |psfaper| can (and should) be hand-edited with
+|setaper| to add any extra features such as masking and linking.
+
+
Customisation
=============
diff --git a/hipercam/psf_reduction.py b/hipercam/psf_reduction.py
index 4840870..b554bb4 100644
--- a/hipercam/psf_reduction.py
+++ b/hipercam/psf_reduction.py
@@ -4,174 +4,386 @@
import numpy as np
from astropy import units as u
from astropy.modeling import Fittable2DModel, Parameter
-from astropy.modeling.fitting import LevMarLSQFitter
-from astropy.stats import SigmaClip, gaussian_fwhm_to_sigma, gaussian_sigma_to_fwhm
+from astropy.modeling.utils import ellipse_extent
+from astropy.stats import SigmaClip
from astropy.table import Table
-from photutils.background import MMMBackground
-from photutils.psf import DAOGroup, BasicPSFPhotometry
-from scipy.special import erf
+from photutils.background import LocalBackground, MedianBackground
+from photutils.psf import GaussianPRF, PSFPhotometry
+from photutils.psf.functional_models import FLOAT_EPSILON, GAUSSIAN_FWHM_TO_SIGMA
import hipercam as hcam
-from hipercam.reduction import moveApers
-
# Stuff below here are helper routines that are not exported
-class IntegratedGaussianPRF2(Fittable2DModel):
- r"""
- Circular Gaussian model integrated over pixels.
- Because it is integrated, this model is considered a PRF, *not* a
- PSF (see :ref:`psf-terminology` for more about the terminology used
- here.)
- This model is a Gaussian *integrated* over an area of
- ``1`` (in units of the model input coordinates, e.g., 1
- pixel). This is in contrast to the apparently similar
- `astropy.modeling.functional_models.Gaussian2D`, which is the value
- of a 2D Gaussian *at* the input coordinates, with no integration.
- So this model is equivalent to assuming the PSF is Gaussian at a
- *sub-pixel* level.
+class MoffatPSF(Fittable2DModel):
+ r"""
+ Standard Moffat model, but with parameter names that work with photutils
- Based on the same named model in photutils, but with different
- std deviations in x and y.
+ Different FWHM parameters apply in x and y, the fwhm property gives
+ the average of both.
Parameters
----------
- sigma_x : float
- Width of the Gaussian PSF in x.
- sigma_y : float
- Width of the Gaussian PSF in y.
flux : float, optional
- Total integrated flux over the entire PSF
+ Total integrated flux over the entire PSF.
+
x_0 : float, optional
- Position of the peak in x direction.
+ Position of the peak along the x axis.
+
y_0 : float, optional
- Position of the peak in y direction.
+ Position of the peak along the y axis.
+
+ x_fwhm : float, optional
+ Full-Width at Half-Maximum of profile in x
+
+ y_fwhm : float, optional
+ Full-Width at Half-Maximum of profile in y
+
+ beta : float, optional
+ Beta parameter of Moffat profile
+
+ theta : float, optional
+ The counterclockwise rotation angle either as a float (in
+ degrees) or a `~astropy.units.Quantity` angle (optional).
+
+ bbox_factor : float, optional
+ The multiple of the FWHM used to define the bounding box limits.
Notes
-----
This model is evaluated according to the following formula:
.. math::
+ f(x, y) = F \frac{beta - 1}{\pi \alpha^2}
+ \left[1 + a(x-x_0)**2 + +b(x-x_0)(y-y_0) + c(y-y_0)**2 \right]**(-beta)
- f(x, y) =
- \frac{F}{4}
- \left[
- {\rm erf} \left(\frac{x - x_0 + 0.5}
- {\sqrt{2} \sigma_x} \right) -
- {\rm erf} \left(\frac{x - x_0 - 0.5}
- {\sqrt{2} \sigma_x} \right)
- \right]
- \left[
- {\rm erf} \left(\frac{y - y_0 + 0.5}
- {\sqrt{2} \sigma_y} \right) -
- {\rm erf} \left(\frac{y - y_0 - 0.5}
- {\sqrt{2} \sigma_y} \right)
- \right]
-
- where ``erf`` denotes the error function and ``F`` the total
- integrated flux.
- """
-
- flux = Parameter(default=1)
- x_0 = Parameter(default=0)
- y_0 = Parameter(default=0)
- sigma_x = Parameter(default=1, fixed=False)
- sigma_y = Parameter(default=1, fixed=False)
-
- _erf = None
-
- @property
- def bounding_box(self):
- halfwidth = 2 * (self.sigma_x + self.sigma_y)
- return (
- (int(self.y_0 - halfwidth), int(self.y_0 + halfwidth)),
- (int(self.x_0 - halfwidth), int(self.x_0 + halfwidth)),
- )
+ where :math:`F` is the total integrated flux, :math:`(x_{0},
+ y_{0})` is the position of the peak. Note that :math:`beta` must be
+ greater than 1.
- @property
- def fwhm(self):
- return u.Quantity(0.5 * (self.sigma_x + self.sigma_y)) * gaussian_sigma_to_fwhm
-
- def evaluate(self, x, y, flux, x_0, y_0, sigma_x, sigma_y):
- """Model function Gaussian PSF model."""
- return (
- flux
- / 4
- * (
- (
- erf((x - x_0 + 0.5) / (np.sqrt(2) * sigma_x))
- - erf((x - x_0 - 0.5) / (np.sqrt(2) * sigma_x))
- )
- * (
- erf((y - y_0 + 0.5) / (np.sqrt(2) * sigma_y))
- - erf((y - y_0 - 0.5) / (np.sqrt(2) * sigma_y))
- )
- )
- )
+ The parameters :math:`a`, :math:`b` and :math:`c` are given by:
+ .. math::
+ a = \cos(\theta)^2 / \alpha_x^2 + \sin(\theta)^2 / \alpha_y^2
+ b = \sin(2\theta) / \alpha_x^2 - \sin(2\theta) / \alpha_y^2
+ c = \sin(\theta)^2 / \alpha_x^2 + \cos(\theta)^2 / \alpha_y^2
-class MoffatPSF(Fittable2DModel):
- """
- Standard Moffat model, but with parameter names that work with photutils
+ where :math:`\alpha_x` and :math:`\alpha_y` are chosen to given the correct
+ FWHM in the x and y directions respectively. :math:`\theta` is the anti-clockwise
+ rotation of the PSF. :math:`alpha` is the geometric mean of :math:`\alpha_x` and
+ :math:`\alpha_y`. The FWHM of the Moffat profile is given by:
- Different FWHM parameters apply in x and y, the fwhm property gives
- the average of both.
+ .. math::
- Parameters
- ----------
- fwhm_x : float
- Full-Width at Half-Maximum of profile in x
- fwhm_y : float
- Full-Width at Half-Maximum of profile in y
- beta : float
- Beta parameter of Moffat profile
- flux : float (default 1)
- Total integrated flux over the entire PSF
- x_0 : float (default 0)
- Position of peak in x direction
- y_0 : float (default 0)
- Position of peak in y direction
+ \rm{FWHM} = 2 \alpha \sqrt{2^{1 / beta} - 1}
- Notes
- -----
- This model is evaluated according to the following formula:
+ The model is normalized such that, for :math:`beta > 1`:
- .. math::
- f(x, y) = F [1 + ((x-x_o)/\alpha_x)**2 + ((y-y_o)/\alpha_y)**2]**(-beta)
+ .. math::
- where ``F`` is the total integrated flux, and ``alpha_i`` is chosen
- so the profile has the appropriate FWHM.
+ \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} f(x, y)
+ \,dx \,dy = F
"""
flux = Parameter(default=1)
x_0 = Parameter(default=0)
y_0 = Parameter(default=0)
- beta = Parameter(default=2.5, fixed=False)
- fwhm_x = Parameter(default=12, fixed=False)
- fwhm_y = Parameter(default=12, fixed=False)
-
- fit_deriv = None
+ x_fwhm = Parameter(default=12, bounds=(FLOAT_EPSILON, None), fixed=False)
+ y_fwhm = Parameter(default=12, bounds=(FLOAT_EPSILON, None), fixed=False)
+ beta = Parameter(default=2.5, bounds=(1.0, None), fixed=False)
+ theta = Parameter(default=0, bounds=(-90, 90), fixed=False)
+
+ def __init__(
+ self,
+ *,
+ flux=flux.default,
+ x_0=x_0.default,
+ y_0=y_0.default,
+ x_fwhm=x_fwhm.default,
+ y_fwhm=y_fwhm.default,
+ beta=beta.default,
+ theta=theta.default,
+ bbox_factor=15.5,
+ **kwargs,
+ ):
+ super().__init__(
+ flux=flux,
+ x_0=x_0,
+ y_0=y_0,
+ x_fwhm=x_fwhm,
+ y_fwhm=y_fwhm,
+ beta=beta,
+ theta=theta,
+ **kwargs,
+ )
+ self.bbox_factor = bbox_factor
@property
def fwhm(self):
- return u.Quantity(0.5 * (self.fwhm_x + self.fwhm_y))
+ return u.Quantity(0.5 * (self.x_fwhm + self.y_fwhm))
+
+ def _calc_bounding_box(self, factor=5.5):
+ """
+ Calculate a bounding box defining the limits of the model.
+
+ The limits are adjusted for rotation.
+
+ Parameters
+ ----------
+ factor : float, optional
+ The multiple of the x and y standard deviations (sigma) used
+ to define the limits.
+
+ Returns
+ -------
+ bbox : tuple
+ A bounding box defining the ((y_min, y_max), (x_min, x_max))
+ limits of the model.
+ """
+ a = factor * self.x_fwhm * GAUSSIAN_FWHM_TO_SIGMA
+ b = factor * self.y_fwhm * GAUSSIAN_FWHM_TO_SIGMA
+ theta = self.theta
+ if not isinstance(theta, u.Quantity):
+ theta = np.deg2rad(theta)
+
+ dx, dy = ellipse_extent(a, b, theta)
+ return ((self.y_0 - dy, self.y_0 + dy), (self.x_0 - dx, self.x_0 + dx))
@property
def bounding_box(self):
- halfwidth = 3 * self.fwhm
- return (
- (int(self.y_0 - halfwidth), int(self.y_0 + halfwidth)),
- (int(self.x_0 - halfwidth), int(self.x_0 + halfwidth)),
+ """
+ The bounding box of the model.
+
+ Examples
+ --------
+ >>> from hipercam.psf_reduction import MoffatPSF
+ >>> model = MoffatPSF(x_0=0, y_0=0, x_fwhm=5, y_fwhm=3, theta=20)
+ >>> model.bounding_box # doctest: +FLOAT_CMP
+ ModelBoundingBox(
+ intervals={
+ x: Interval(lower=-11.232523683597986, upper=11.232523683597986)
+ y: Interval(lower=-7.701096862895421, upper=7.701096862895421)
+ }
+ model=MoffatPSF(inputs=('x', 'y'))
+ order='C'
+ )
+ >>> model.bbox_factor = 7
+ >>> model.bounding_box # doctest: +FLOAT_CMP
+ ModelBoundingBox(
+ intervals={
+ x: Interval(lower=-14.295939233670163, upper=14.295939233670163)
+ y: Interval(lower=-9.801396007321445, upper=9.801396007321445)
+ }
+ model=MoffatPSF(inputs=('x', 'y'))
+ order='C'
+ )
+ """
+ return self._calc_bounding_box(factor=self.bbox_factor)
+
+ def evaluate(self, x, y, flux, x_0, y_0, x_fwhm, y_fwhm, beta, theta):
+ """
+ Calculate the value of the 2D Moffat model at the input
+ coordinates for the given model parameters.
+
+ Parameters
+ ----------
+ x, y : float or array_like
+ The x and y coordinates at which to evaluate the model.
+
+ flux : float
+ Total integrated flux over the entire PSF.
+
+ x_0, y_0 : float
+ Position of the peak along the x and y axes.
+
+ x_fwhm, y_fwhm : float
+ FWHM of the Gaussian along the x and y axes.
+
+ beta : float, optional
+ Beta parameter of Moffat profile
+
+ theta : float
+ The counterclockwise rotation angle either as a float (in
+ degrees) or a `~astropy.units.Quantity` angle (optional).
+
+ Returns
+ -------
+ result : `~numpy.ndarray`
+ The value of the model evaluated at the input coordinates.
+ """
+ # find useful derived properties of theta
+ if not isinstance(theta, u.Quantity):
+ theta = np.deg2rad(theta)
+ cost2 = np.cos(theta) ** 2
+ sint2 = np.sin(theta) ** 2
+ sin2t = np.sin(2.0 * theta)
+
+ # alpha_sq terms for x and y
+ alpha_sq_x = x_fwhm**2 / 4 / (2 ** (1 / beta) - 1)
+ alpha_sq_y = y_fwhm**2 / 4 / (2 ** (1 / beta) - 1)
+
+ # terms a, b, c
+ a = cost2 / alpha_sq_x + sint2 / alpha_sq_y
+ b = sin2t / alpha_sq_x - sin2t / alpha_sq_y
+ c = sint2 / alpha_sq_x + cost2 / alpha_sq_y
+ x_term = a * (x - x_0) ** 2
+ xy_term = b * (x - x_0) * (y - y_0)
+ y_term = c * (y - y_0) ** 2
+ profile = (1 + x_term + y_term + xy_term) ** (-beta)
+ normalisation = (beta - 1) / np.pi / np.sqrt(alpha_sq_x * alpha_sq_y)
+ return normalisation * flux * profile
+
+ @staticmethod
+ def fit_deriv(x, y, flux, x_0, y_0, x_fwhm, y_fwhm, beta, theta):
+ """
+ Calculate the partial derivatives of the 2D Moffat function
+ with respect to the parameters.
+
+ Parameters
+ ----------
+ x, y : float or array_like
+ The x and y coordinates at which to evaluate the model.
+
+ flux : float
+ Total integrated flux over the entire PSF.
+
+ x_0, y_0 : float
+ Position of the peak along the x and y axes.
+
+ x_fwhm, y_fwhm : float
+ FWHM of the Gaussian along the x and y axes.
+
+ beta : float, optional
+ Beta parameter of Moffat profile
+
+ theta : float
+ The counterclockwise rotation angle either as a float (in
+ degrees) or a `~astropy.units.Quantity` angle (optional).
+
+ Returns
+ -------
+ result : list of `~numpy.ndarray`
+ The list of partial derivatives with respect to each
+ parameter.
+ """
+ if not isinstance(theta, u.Quantity):
+ theta = np.deg2rad(theta)
+ else:
+ theta = theta.to_value(u.rad)
+
+ alpha_sq_x = x_fwhm**2 / 4 / (2 ** (1 / beta) - 1)
+ alpha_sq_y = y_fwhm**2 / 4 / (2 ** (1 / beta) - 1)
+ alpha_x = np.sqrt(alpha_sq_x)
+ alpha_y = np.sqrt(alpha_sq_y)
+
+ x1 = beta - 1
+ x2 = 1 / np.pi
+ x3 = 1 / alpha_x
+ x4 = 1 / alpha_y
+ x5 = x - x_0
+ x6 = x5**2
+ x7 = alpha_x ** (-2)
+ x8 = np.cos(theta)
+ x9 = x8**2
+ x10 = alpha_y ** (-2)
+ x11 = np.sin(theta)
+ x12 = x11**2
+ x13 = x10 * x12 + x7 * x9
+ x14 = y - y_0
+ x15 = x14**2
+ x16 = x10 * x9 + x12 * x7
+ x17 = 2 * theta
+ x18 = np.sin(x17)
+ x19 = -x10 * x18 + x18 * x7
+ x20 = x14 * x19
+ x21 = x13 * x6 + x15 * x16 + x20 * x5 + 1
+ x22 = x21 ** (-beta)
+ x23 = x1 * x2 * x22 * x3 * x4
+ x24 = 1 / x21
+ x25 = 2 / alpha_x**3
+ x26 = x14 * x5
+ x27 = alpha_y ** (-3)
+ x28 = 2 * x27
+ x29 = 2 * x11 * x8
+ x30 = x10 * x29 - x29 * x7
+ x31 = np.cos(x17)
+
+ dg_dflux = x23
+ dg_dx_0 = -flux * beta * x23 * x24 * (x13 * (-2 * x + 2 * x_0) - x20)
+ dg_dy_0 = -flux * beta * x23 * x24 * (x16 * (-2 * y + 2 * y_0) - x19 * x5)
+ dg_dalpha_x = (
+ -flux
+ * beta
+ * x23
+ * x24
+ * (-x12 * x15 * x25 - x18 * x25 * x26 - x25 * x6 * x9)
+ - flux * x1 * x2 * x22 * x4 * x7
+ )
+ dg_dalpha_y = (
+ -flux
+ * beta
+ * x23
+ * x24
+ * (-x12 * x28 * x6 + 2 * x14 * x18 * x27 * x5 - x15 * x28 * x9)
+ - flux * x1 * x10 * x2 * x22 * x3
+ )
+ partial_dg_dbeta = flux * x2 * x22 * x3 * x4 - flux * x23 * np.log(x21)
+ dg_dtheta = (
+ (
+ -flux
+ * beta
+ * x23
+ * x24
+ * (-x15 * x30 + x26 * (-2 * x10 * x31 + 2 * x31 * x7) + x30 * x6)
+ )
+ * np.pi
+ / 180
)
- def evaluate(self, x, y, flux, x_0, y_0, beta, fwhm_x, fwhm_y):
- alpha_sq_x = fwhm_x**2 / 4 / (2 ** (1 / beta) - 1)
- alpha_sq_y = fwhm_y**2 / 4 / (2 ** (1 / beta) - 1)
- x_term = (x - x_0) ** 2 / alpha_sq_x
- y_term = (y - y_0) ** 2 / alpha_sq_y
- prof = flux / (1 + x_term + y_term) ** beta
- return (beta - 1) * prof / np.pi / np.sqrt(alpha_sq_x * alpha_sq_y)
+ # jacobians
+ dg_dx_fwhm = alpha_x * dg_dalpha_x / x_fwhm
+ dg_dy_fwhm = alpha_y * dg_dalpha_y / y_fwhm
+ dalpha_x_dbeta = (
+ 2 ** (1 / beta) * np.sqrt(x_fwhm**2 / (2 ** (1 / beta) - 1)) * np.log(2)
+ )
+ dalpha_x_dbeta /= 4 * beta**2 * (2 ** (1 / beta) - 1)
+ dalpha_y_dbeta = (
+ 2 ** (1 / beta) * np.sqrt(y_fwhm**2 / (2 ** (1 / beta) - 1)) * np.log(2)
+ )
+ dalpha_y_dbeta /= 4 * beta**2 * (2 ** (1 / beta) - 1)
+ dg_dbeta = (
+ partial_dg_dbeta
+ + dg_dalpha_x * dalpha_x_dbeta
+ + dg_dalpha_y * dalpha_y_dbeta
+ )
+
+ return dg_dflux, dg_dx_0, dg_dy_0, dg_dx_fwhm, dg_dy_fwhm, dg_dbeta, dg_dtheta
+
+
+def create_psf_model(photom_results, method, fixed_positions=False):
+ if method == "moffat":
+ psf_model = MoffatPSF(flux=1)
+ for param in ["x_fwhm", "y_fwhm", "theta", "beta"]:
+ colname = param + "_fit"
+ if colname in photom_results.colnames:
+ setattr(psf_model, param, np.median(photom_results[colname]))
+ psf_model.beta.fixed = True
+ else:
+ psf_model = GaussianPRF(flux=1)
+ for param in ["x_fwhm", "y_fwhm", "theta"]:
+ colname = param + "_fit"
+ if colname in photom_results.colnames:
+ setattr(psf_model, param, np.median(photom_results[colname]))
+
+ if fixed_positions:
+ psf_model.x_0.fixed = True
+ psf_model.y_0.fixed = True
+
+ psf_model.x_fwhm.fixed = True
+ psf_model.y_fwhm.fixed = True
+ psf_model.theta.fixed = True
+ return psf_model
def extractFluxPSF(cnam, ccd, bccd, rccd, read, gain, ccdwin, rfile, store):
@@ -313,6 +525,43 @@ def extractFluxPSF(cnam, ccd, bccd, rccd, read, gain, ccdwin, rfile, store):
"cmax": 0,
}
return results
+
+ # we need at least one reference aperture to proceed
+ nref = sum(1 for ap in ccdaper.values() if ap.ref)
+ if nref == 0:
+ print(
+ (
+ " *** WARNING: CCD {:s}: no reference apertures"
+ " defined within window; no extraction possible"
+ ).format(cnam)
+ )
+
+ # set flag to indicate no extraction
+ flag = hcam.NO_EXTRACTION
+
+ # return empty results
+ for apnam, aper in ccdaper.items():
+ info = store[apnam]
+ results[apnam] = {
+ "x": aper.x,
+ "xe": info["xe"],
+ "y": aper.y,
+ "ye": info["ye"],
+ "fwhm": info["fwhm"],
+ "fwhme": info["fwhme"],
+ "beta": info["beta"],
+ "betae": info["betae"],
+ "counts": 0.0,
+ "countse": -1,
+ "sky": 0.0,
+ "skye": 0.0,
+ "nsky": 0,
+ "nrej": 0,
+ "flag": flag,
+ "cmax": 0,
+ }
+ return results
+
wnam = wnames.pop()
# PSF params are in binned pixels, so find binning
@@ -323,21 +572,24 @@ def extractFluxPSF(cnam, ccd, bccd, rccd, read, gain, ccdwin, rfile, store):
if method == "moffat":
# fix beta if initial aperture tweak used Gaussian profiles
beta = mbeta if mbeta > 0.0 else 2.5
- psf_model = MoffatPSF(beta=beta, fwhm_x=mfwhm / bin_fac, fwhm_y=mfwhm / bin_fac)
- # TODO: adjust psf_model based on reference stars
- psf_model.fwhm_x.fixed = True
- psf_model.fwhm_y.fixed = True
- psf_model.beta.fixed = True
- extra_output_cols = ["fwhm_x", "fwhm_y"]
+ psf_model = MoffatPSF(beta=beta, x_fwhm=mfwhm / bin_fac, y_fwhm=mfwhm / bin_fac)
+ # let the shape vary when we constrain psf_model based on reference stars
+ psf_model.x_fwhm.fixed = False
+ psf_model.y_fwhm.fixed = False
+ psf_model.beta.fixed = False
+ psf_model.theta.fixed = False
+ psf_model.theta.bounds = (-90, 90)
else:
- psf_model = IntegratedGaussianPRF2(
- sigma_x=mfwhm * gaussian_fwhm_to_sigma / bin_fac,
- sigma_y=mfwhm * gaussian_fwhm_to_sigma / bin_fac,
+ psf_model = GaussianPRF(
+ x_fwhm=mfwhm / bin_fac,
+ y_fwhm=mfwhm / bin_fac,
+ theta=0.0,
)
- # TODO: adjust psf model based on reference stars
- psf_model.sigma_x.fixed = True
- psf_model.sigma_y.fixed = True
- extra_output_cols = ["sigma_x", "sigma_y"]
+ # let the shape vary when we constrain psf_model based on reference stars
+ psf_model.x_fwhm.fixed = False
+ psf_model.y_fwhm.fixed = False
+ psf_model.theta.fixed = False
+ psf_model.theta.bounds = (-90, 90)
# force photometry only at aperture positions
# this means PSF shape and positions are fixed, we are only fitting flux
@@ -345,30 +597,15 @@ def extractFluxPSF(cnam, ccd, bccd, rccd, read, gain, ccdwin, rfile, store):
psf_model.x_0.fixed = True
psf_model.y_0.fixed = True
- # create instances for PSF photometry
- gfac = float(rfile["psf_photom"]["gfac"])
- sclip = float(rfile["sky"]["thresh"])
- daogroup = DAOGroup(gfac * mfwhm / bin_fac)
- mmm_bkg = MMMBackground(sigma_clip=SigmaClip(sclip))
- fitter = LevMarLSQFitter()
- fitshape_box_size = int(2 * int(float(rfile["psf_photom"]["fit_half_width"])) + 1)
- fitshape = (fitshape_box_size, fitshape_box_size)
-
- photometry_task = BasicPSFPhotometry(
- group_maker=daogroup,
- bkg_estimator=mmm_bkg,
- psf_model=psf_model,
- fitter=fitter,
- fitshape=fitshape,
- extra_output_cols=extra_output_cols,
- )
-
# initialise flag
flag = hcam.ALL_OK
# extract Windows relevant for these apertures
wdata = ccd[wnam]
wraw = rccd[wnam]
+ wbias = bccd[wnam]
+ wread = read[wnam]
+ wgain = gain[wnam]
# extract sub-windows that include all of the apertures, plus a little
# extra around the edges.
@@ -378,18 +615,74 @@ def extractFluxPSF(cnam, ccd, bccd, rccd, read, gain, ccdwin, rfile, store):
y2 = max([ap.y + ap.rsky2 + wdata.ybin for ap in ccdaper.values()])
# extract sub-Windows
- swdata = wdata.window(x1, x2, y1, y2)
- swraw = wraw.window(x1, x2, y1, y2)
+ box_limits = (x1, x2, y1, y2)
+ swdata = wdata.window(*box_limits)
+ swraw = wraw.window(*box_limits)
+ swbias = wbias.window(*box_limits)
+ swread = wread.window(*box_limits)
+ swgain = wgain.window(*box_limits)
- # compute pixel positions of apertures in windows
+ # This is where the pure-debiassed data are used
+ sigma = np.sqrt(swread.data**2 + np.maximum(0, swbias.data) / swgain.data)
+
+ # compute pixel positions of apertures in windows (twice, once for reference aps only)
xpos, ypos = zip(
*((swdata.x_pixel(ap.x), swdata.y_pixel(ap.y)) for ap in ccdaper.values())
)
+ xpos = np.array(xpos)
+ ypos = np.array(ypos)
+
positions = Table(names=["x_0", "y_0"], data=(xpos, ypos))
+ # and positions of PSF stars in the sub-window
+ ref_idx = [i for i, ap in enumerate(ccdaper.values()) if ap.ref]
+ psf_positions = Table(names=["x_0", "y_0"], data=(xpos[ref_idx], ypos[ref_idx]))
+ # assign all PSF stars to the same group, so the same PSF params are used for all
+ psf_positions["id"] = 1
+
+ # create instances for PSF photometry
+ sclip = float(rfile["sky"]["thresh"])
+ bkg = LocalBackground(
+ rfile["extraction"][cnam][7], # inner sky
+ rfile["extraction"][cnam][10], # outer sky
+ bkg_estimator=MedianBackground(sigma_clip=SigmaClip(sclip)),
+ )
+ fitshape_box_size = int(2 * int(float(rfile["psf_photom"]["fit_half_width"])) + 1)
+ fit_shape = (fitshape_box_size, fitshape_box_size)
+
+ # aperture radius is only used to guess initial flux
+ aperture_radius = 0.5 * (
+ rfile["extraction"][cnam][3] + rfile["extraction"][cnam][4]
+ )
+ photometry_task = PSFPhotometry(
+ psf_model=psf_model,
+ aperture_radius=aperture_radius,
+ fit_shape=fit_shape,
+ localbkg_estimator=bkg,
+ xy_bounds=rfile["apertures"]["fit_max_shift"],
+ )
# do the PSF photometry
- photom_results = photometry_task(swdata.data, init_guesses=positions)
- slevel = mmm_bkg(swdata.data)
+ photom_results = photometry_task(
+ swdata.data, error=sigma, init_params=psf_positions
+ )
+
+ # now create the PSF model from fitting the reference stars
+ psf_model = create_psf_model(
+ photom_results, method, rfile["psf_photom"]["positions"] == "fixed"
+ )
+
+ # re-do the PSF photometry with the fixed PSF shape and all stars
+ photometry_task = PSFPhotometry(
+ psf_model=psf_model,
+ aperture_radius=aperture_radius,
+ fit_shape=fit_shape,
+ localbkg_estimator=bkg,
+ xy_bounds=rfile["apertures"]["fit_max_shift"],
+ )
+ photom_results = photometry_task(swdata.data, error=sigma, init_params=positions)
+
+ # store PSF task and settings for display by external routines
+ store[f"psf_model_{cnam}"] = (photometry_task, wnam, box_limits)
# unpack the results and check apertures
for apnam, aper in ccdaper.items():
@@ -409,7 +702,10 @@ def extractFluxPSF(cnam, ccd, bccd, rccd, read, gain, ccdwin, rfile, store):
"ambiguous lookup for this aperture in PSF photometry"
)
else:
+ # USE PSF PHOTOMETRY FLAGS TO SET EXTRACTION FLAG
result_row = result_row[0]
+ if result_row['flags'] > 0:
+ flag |= hcam.NO_EXTRACTION
# compute X, Y arrays over the sub-window relative to the centre
# of the aperture and the distance squared from the centre (Rsq)
@@ -445,11 +741,20 @@ def extractFluxPSF(cnam, ccd, bccd, rccd, read, gain, ccdwin, rfile, store):
counts = result_row["flux_fit"]
try:
- countse = result_row["flux_unc"]
+ countse = result_row["flux_err"]
except KeyError:
raise hcam.HipercamError(
"unable to find errors on solution, model does not depend on params"
)
+ slevel = result_row["local_bkg"]
+
+ # check flags from PSF photom
+ if result_row["flags"] == 2:
+ flag = hcam.TARGET_AT_EDGE
+ elif result_row["flags"] > 8:
+ # fit failed
+ flag = hcam.NO_EXTRACTION
+
info = store[apnam]
results[apnam] = {
diff --git a/hipercam/reduction.py b/hipercam/reduction.py
index 5d98f51..fca9542 100644
--- a/hipercam/reduction.py
+++ b/hipercam/reduction.py
@@ -4,34 +4,34 @@
from collections import OrderedDict
import sys
import warnings
-import numpy as np
-import hipercam as hcam
-from hipercam import utils, fitting
-from hipercam.ccd import crop_calib_frame_to
+import numpy as np
from trm.pgplot import (
- pgpap,
- pgsubp,
- pgsci,
- pgsch,
+ pgbbuf,
+ pgbox,
+ pgdraw,
+ pgebuf,
pgenv,
- pglab,
- pgqvp,
- pgvstd,
pgeras,
pgerry,
- pgpt,
+ pglab,
+ pgmove,
pgpanl,
- pgbbuf,
- pgebuf,
+ pgpap,
+ pgpt,
+ pgpt1,
+ pgqvp,
+ pgsch,
+ pgsci,
+ pgsubp,
pgsvp,
pgswin,
- pgbox,
- pgmove,
- pgdraw,
- pgpt1,
+ pgvstd,
)
+import hipercam as hcam
+from hipercam import fitting, utils
+from hipercam.ccd import crop_calib_frame_to
NaN = float("NaN")
@@ -86,8 +86,7 @@ def read(cls, filename, readcal=True):
else:
raise hcam.HipercamError(
- "found entry line before"
- " any section = \n{:s}".format(line)
+ "found entry line before any section = \n{:s}".format(line)
)
# process it. this is a matter of checking entries and
@@ -624,11 +623,11 @@ def crop(self, mccd):
self.gain = crop_calib_frame_to(mccd, self.gain, "gain")
-def setup_plots(rfile, ccds, nx, plot_lims, implot=True, lplot=True):
+def setup_plots(rfile, ccds, nx, plot_lims, implot=True, lplot=True, psfplot=True):
"""
Perform initial setup of image and lightcurve plots for reduction outputs
"""
- imdev = lcdev = spanel = tpanel = xpanel = ypanel = lpanel = None
+ imdev = lcdev = psfdev = spanel = tpanel = xpanel = ypanel = lpanel = None
if implot:
xlo, xhi, ylo, yhi = plot_lims
# image plot
@@ -655,6 +654,31 @@ def setup_plots(rfile, ccds, nx, plot_lims, implot=True, lplot=True):
pgsch(hcam.pgp.Params["axis.label.ch"])
pglab("X", "Y", "CCD {:s}".format(cnam))
+ if psfplot:
+ xlo, xhi, ylo, yhi = plot_lims
+ # plot of image minus PSF model
+ psfdev = hcam.pgp.Device(rfile["general"]["psfdevice"])
+ iwidth = rfile["general"]["iwidth"]
+ iheight = rfile["general"]["iheight"]
+
+ if iwidth > 0 and iheight > 0:
+ pgpap(iwidth, iheight / iwidth)
+ # set up panels and axes
+ nccd = len(ccds)
+ ny = nccd // nx if nccd % nx == 0 else nccd // nx + 1
+
+ # slice up viewport
+ pgsubp(nx, ny)
+
+ # plot axes, labels, titles. Happens once only
+ for cnam in ccds:
+ pgsci(hcam.pgp.Params["axis.ci"])
+ pgsch(hcam.pgp.Params["axis.number.ch"])
+ pgenv(xlo, xhi, ylo, yhi, 1, 0)
+ pgsci(hcam.pgp.Params["axis.label.ci"])
+ pgsch(hcam.pgp.Params["axis.label.ch"])
+ pglab("X", "Y", "CCD {:s}".format(cnam))
+
if lplot:
# open the light curve plot
lcdev = hcam.pgp.Device(rfile["general"]["ldevice"])
@@ -813,7 +837,7 @@ def setup_plots(rfile, ccds, nx, plot_lims, implot=True, lplot=True):
rfile["light"]["y2"],
)
- return imdev, lcdev, spanel, tpanel, xpanel, ypanel, lpanel
+ return imdev, lcdev, psfdev, spanel, tpanel, xpanel, ypanel, lpanel
def setup_plot_buffers(rfile):
@@ -855,10 +879,13 @@ def decorate(func):
def update_plots(
results,
rfile,
+ store,
implot,
lplot,
+ psfplot,
imdev,
lcdev,
+ psfdev,
pccd,
ccds,
msub,
@@ -902,12 +929,18 @@ def update_plots(
lplot : bool
whether to plot light curves
+ psfplot : bool
+ whether to plot PSF residuals
+
imdev : hipercam.pgp.Device
represents the image plot device (if any)
lcdev : hipercam.pgp.Device
represents the light curve plot device (if any)
+ psfdev : hipercam.pgp.Device
+ represents the psf residuals plot device (if any)
+
pccd : MCCD [only matters if imdev]
current debiassed, flat-fielded frame. In the case of frame
groups, this will be the last one added.
@@ -1001,6 +1034,67 @@ def update_plots(
# end of CCD display loop
print(message)
+ if psfplot:
+ # plot the PSF residuals
+ psfdev.select()
+
+ # display the CCDs chosen
+ message = "; "
+ for nc, cnam in enumerate(ccds):
+ ccd = pccd[cnam]
+
+ if ccd.is_data():
+ # this should be data as opposed to a blank frame
+ # between data frames that occur with nskip > 0
+
+ # set to the correct panel and then plot CCD
+ ix = (nc % nx) + 1
+ iy = nc // nx + 1
+ pgpanl(ix, iy)
+ psf_model, box_wnam, box_limits = store[cnam].get(
+ f"psf_model_{cnam}", None
+ )
+
+ if psf_model is not None:
+ data_minus_model = ccd.copy()
+
+ # find the indices of the stamp used for PSF photometry
+ window = data_minus_model[box_wnam]
+ x1, x2, y1, y2 = box_limits
+ # code copied from window class
+ winh = super(hcam.Window, window).window(x1, x2, y1, y2)
+ ix1 = (winh.llx - window.llx) // window.xbin
+ iy1 = (winh.lly - window.lly) // window.ybin
+ iy2 = iy1 + winh.ny
+ ix2 = ix1 + winh.nx
+
+ # extract a view of the data array
+ data = window.data[iy1:iy2, ix1:ix2]
+ resid = psf_model.make_residual_image(data)
+ data_minus_model[box_wnam].data[iy1:iy2, ix1:ix2] = resid
+ vmin, vmax = ccd.percentile((plo, phi), xlo, xhi, ylo, yhi)
+ hcam.pgp.pCcd(
+ data_minus_model,
+ iset="d",
+ dlo=vmin,
+ dhi=vmax,
+ xlo=xlo,
+ xhi=xhi,
+ ylo=ylo,
+ yhi=yhi,
+ )
+
+ # plot apertures
+ if cnam in rfile.aper and len(rfile.aper[cnam]) > 0:
+ # find the result for this ccd name
+ res = next(res for (this_cnam, res) in results if this_cnam == cnam)
+ # we always use the last result in the parallel group
+ last_res = res[-1]
+ ccd_aper = last_res[2]
+ hcam.pgp.pCcdAper(ccd_aper)
+
+ # end of PSF display loop
+
if lplot:
# plot the light curve
@@ -2535,7 +2629,7 @@ def extractFlux(cnam, ccd, bccd, rccd, read, gain, ccdwin, rfile, store):
if extype == "optimal":
# optimal extraction. Need the profile
warnings.warn(
- "Transmission plot is not reliable" " with optimal extraction"
+ "Transmission plot is not reliable with optimal extraction"
)
mbeta = store["mbeta"]
@@ -3662,9 +3756,7 @@ def write_header(self):
# of the HiPERCAM reduction software, and was generated using the following
# command-line inputs to 'reduce':
#
-""".format(
- hipercam_version=self.hipercam_version
- )
+""".format(hipercam_version=self.hipercam_version)
)
# second, list the command-line inputs to the logfile
diff --git a/hipercam/scripts/genred.py b/hipercam/scripts/genred.py
index 222e880..32998ec 100644
--- a/hipercam/scripts/genred.py
+++ b/hipercam/scripts/genred.py
@@ -424,6 +424,7 @@ def genred(args=None):
gsec["iwidth"] = 0.0
gsec["iheight"] = 0.0
gsec["toffset"] = 0
+ gsec["psfdevice"] = "3/xs"
gsec["skipbadt"] = "yes"
gsec["ncpu"] = 1
gsec["ngroup"] = 1
@@ -461,7 +462,6 @@ def genred(args=None):
psfsec = tvals["psf_photom"] = {}
psfsec["use_psf"] = "no"
psfsec["psf_model"] = "moffat"
- psfsec["gfac"] = 3.0
psfsec["fit_half_width"] = 15.0
psfsec["positions"] = "fixed"
@@ -867,6 +867,7 @@ def genred(args=None):
iwidth = {gsec["iwidth"]} # image curve plot width, inches, 0 to let program choose
iheight = {gsec["iheight"]} # image curve plot height, inches
+psfdevice = {gsec["psfdevice"]} # PGPLOT plot device for PSF residual plots [if psfplot True]
toffset = {gsec["toffset"]} # integer offset to subtract from the MJD
# skip points with bad times in plots. HiPERCAM has a problem in not
@@ -1045,26 +1046,21 @@ def genred(args=None):
[extraction]
{extraction}
-# The next lines are specific to the PSF photometry option. 'gfac' is
-# used to label the sources according to groups, such that stars
-# closer than 'gfac' times the FWHM are labelled in the same
-# group. Each group has a PSF model fit independently. The reason
-# behind the construction of groups is to reduce the dimensionality of
-# the fitting procedure. Usually you want closely seperated stars to
-# be fit simultaneously, but too large a value will mean fitting a
-# model with many free parameters, which can fail to converge. The
-# size of the box over which data is collected for fitting is set by
+# The next lines are specific to the PSF photometry option.
+# In PSF photometry, the PSF model is fit to the reference stars and
+# then the shape of the PSF is held fixed and only the scaling of the PSF
+# [and optionally the position] is fit for other sources.
+# The size of the box over which data is collected for fitting is set by
# 'fit_half_width'. Finally, 'positions' determines whether the star's
# positions should be considered variable in the PSF fitting. If this
# is set to fixed, the positions are held at the locations found in
# the aperture repositioning step, otherwise the positions are refined
# during PSF fitting. This step can fail for PSF photometry of faint
-# sources.
+# sources, generally fixed positions yields better results.
[psf_photom]
use_psf = {psfsec["use_psf"]} # 'yes' or 'no'
psf_model = {psfsec["psf_model"]} # 'gaussian' or 'moffat'
-gfac = {psfsec["gfac"]} # multiple of the FWHM to use in grouping objects
fit_half_width = {psfsec["fit_half_width"]} # size of window used to collect the data to do the fitting
positions = {psfsec["positions"]} # 'fixed' or 'variable'
diff --git a/hipercam/scripts/psfaper.py b/hipercam/scripts/psfaper.py
index 022d891..f6c774b 100644
--- a/hipercam/scripts/psfaper.py
+++ b/hipercam/scripts/psfaper.py
@@ -1,25 +1,30 @@
-import sys
import os
+import sys
import warnings
from copy import deepcopy
-import numpy as np
import matplotlib as mpl
-
-from photutils.psf import DAOPhotPSFPhotometry, IntegratedGaussianPRF
+import numpy as np
+from astropy.stats import SigmaClip
+from astropy.table import Table
from photutils.background import MADStdBackgroundRMS, MMMBackground
-from astropy.stats import gaussian_fwhm_to_sigma, SigmaClip
-
+from photutils.detection import DAOStarFinder
+from photutils.psf import (
+ GaussianPRF,
+ IterativePSFPhotometry,
+ PSFPhotometry,
+ SourceGrouper,
+)
from trm import cline
from trm.cline import Cline
import hipercam as hcam
-from hipercam.psf_reduction import MoffatPSF
+from hipercam.psf_reduction import MoffatPSF, create_psf_model
# re-configure the cursors: backend specific.
# aim to get rid of irritating 'hand' icon in
# favour of something pointier.
-
+
backend = mpl.get_backend()
if backend == "Qt4Agg" or "Qt5Agg":
@@ -216,7 +221,6 @@ def psfaper(args=None):
# get input section
with Cline("HIPERCAM_ENV", ".hipercam", command, args) as cl:
-
# register parameters
cl.register("mccd", Cline.LOCAL, Cline.PROMPT)
cl.register("aper", Cline.LOCAL, Cline.PROMPT)
@@ -263,9 +267,7 @@ def psfaper(args=None):
else:
# create empty container
mccdaper = hcam.MccdAper()
- print(
- "No file called {:s} exists; " "will create from scratch".format(aper)
- )
+ print("No file called {:s} exists; will create from scratch".format(aper))
# define the panel grid
nxdef = cl.get_default("nx", 3)
@@ -342,7 +344,7 @@ def psfaper(args=None):
iset = cl.get_value(
"iset",
- "set intensity a(utomatically)," " d(irectly) or with p(ercentiles)?",
+ "set intensity a(utomatically), d(irectly) or with p(ercentiles)?",
"a",
lvals=["a", "A", "d", "D", "p", "P"],
)
@@ -371,21 +373,21 @@ def psfaper(args=None):
shbox = cl.get_value(
"shbox",
- "half width of box for initial" " location of target [unbinned pixels]",
+ "half width of box for initial location of target [unbinned pixels]",
11.0,
2.0,
)
smooth = cl.get_value(
"smooth",
- "FWHM for smoothing for initial object" " detection [binned pixels]",
+ "FWHM for smoothing for initial object detection [binned pixels]",
6.0,
)
fhbox = cl.get_value(
- "fhbox", "half width of box for profile fit" " [unbinned pixels]", 21.0, 3.0
+ "fhbox", "half width of box for profile fit [unbinned pixels]", 21.0, 3.0
)
- read = cl.get_value("read", "readout noise, RMS ADU", 3.0)
- gain = cl.get_value("gain", "gain, ADU/e-", 1.0)
+ read = cl.get_value("read", "readout noise, RMS ADU", 3.5)
+ gain = cl.get_value("gain", "gain, ADU/e-", 1.1)
rejthresh = cl.get_value(
"rejthresh", "RMS rejection threshold for sky fitting", 4.0
)
@@ -465,14 +467,16 @@ def psfaper(args=None):
if cnam in mccdaper:
# plot any pre-existing apertures, keeping track of
# the plot objects
- pobjs[cnam] = hcam.mpl.pCcdAper(axes, mccdaper[cnam])
+ pobjs[cnam] = hcam.Group(list)
+ for key, value in hcam.mpl.pCcdAper(axes, mccdaper[cnam]):
+ pobjs[cnam][key] = value
else:
# add in an empty CcdApers for any CCD not already present
mccdaper[cnam] = hcam.CcdAper()
# and an empty container for any new plot objects
- pobjs[cnam] = hcam.Group(tuple)
+ pobjs[cnam] = hcam.Group(list)
print(
"""
@@ -530,7 +534,7 @@ def psfaper(args=None):
)
# squeeze space a bit
- plt.subplots_adjust(wspace=0.1, hspace=0.1)
+ plt.tight_layout()
# finally show stuff ....
plt.show()
@@ -572,7 +576,6 @@ def __init__(
pobjs,
apernam,
):
-
# save the inputs, tack on event handlers.
self.fig = fig
self.fig.canvas.mpl_connect("key_press_event", self._keyPressEvent)
@@ -771,7 +774,7 @@ def _delete(self):
print(' deleted aperture "{:s}"'.format(apnam))
else:
- print(" found no aperture near enough " "the cursor position for deletion")
+ print(" found no aperture near enough the cursor position for deletion")
def _find_aper(self):
"""Finds the nearest aperture to the currently selected position,
@@ -795,14 +798,20 @@ def _find_aper(self):
def _find_stars_and_quit(self):
"""
- Runs the PSF photometry on each ROI, and writes the aperture file
+ Runs PSF photometry on the region of interest, finds the stars, adds them as apertures.
+
+ For each CCD the PSF shape parameters are estimated from the reference stars,
+ at which point the PSF shape parameters are held fixed and multiple iterations
+ of the FIND-FIT-SUBTRACT loop are run to find more stars.
+
+ The final apertures are then saved to the aperture file and the program exits.
"""
for cnam, ccdaper in self.mccdaper.items():
xlo, xhi = self.anams[cnam].get_xlim()
ylo, yhi = self.anams[cnam].get_ylim()
if cnam not in self.psf_data:
- warnings.warn("no reference stars for CCD{} - skipping".format(cnam))
+ warnings.warn("no (new) reference stars for CCD{} - skipping".format(cnam))
continue
fwhm = -1
@@ -824,6 +833,7 @@ def _find_stars_and_quit(self):
xlocs, ylocs = daophot(
cnam,
self.mccd[cnam],
+ ccdaper,
xlo,
xhi,
ylo,
@@ -835,6 +845,8 @@ def _find_stars_and_quit(self):
self.gfac,
self.thresh,
self.rejthresh,
+ self.read,
+ self.gain,
)
# let's find which of our positions best matches our reference apertures
@@ -868,12 +880,45 @@ def _find_stars_and_quit(self):
def daophot(
- cnam, ccd, xlo, xhi, ylo, yhi, niters, method, fwhm, beta, gfac, thresh, rejthresh
+ cnam, ccd, ccdaper, xlo, xhi, ylo, yhi, niters, method, fwhm, beta, gfac, thresh, rejthresh, read, gain
):
"""
- Perform iterative PSF photometry and star finding on region of CCD
+ Perform PSF photometry on region of CCD.
+
+ The PSF shape parameters are determined from the reference stars, at which point
+ the PSF shape parameters are held fixed and multiple iterations of the FIND-FIT-SUBTRACT
+ loop are run to find more stars.
+
+ Parameters
+ ----------
+ cnam : str
+ name of CCD
+ ccd : hcam.Ccd
+ CCD to perform photometry on
+ ccdaper : hcam.CcdAper
+ reference apertures to use for determining PSF shape parameters
+ xlo, xhi, ylo, yhi : float
+ limits of region to perform photometry on, in device coordinates
+ niters : int
+ number of iterations of FIND-FIT-SUBTRACT to perform
+ method : str
+ PSF model to use, either 'm' for Moffat, 'g' for Gaussian
+ fwhm : float
+ initial FWHM to use for PSF model, in unbinned pixels
+ beta : float
+ initial beta to use for Moffat PSF model, ignored if method is 'g'
+ gfac : float
+ multiple of FWHM used to group stars for fitting.
+ Stars within gfac*FWHM of each other are fitted simultaneously in the FIND-FIT-SUBTRACT loop.
+ thresh : float
+ threshold for object detection in FIND-FIT-SUBTRACT, in multiples of the background RMS
+ rejthresh : float
+ RMS rejection threshold for sky fitting, in multiples of the background RMS.
+ read : float
+ readout noise, RMS ADU, for assigning uncertainties
+ gain : float
+ gain, ADU/count, for assigning uncertainties
"""
- print(xlo, ylo, xhi, yhi)
# first check we are within a single window
wnam1 = ccd.inside(xlo, ylo, 2)
wnam2 = ccd.inside(xhi, yhi, 2)
@@ -882,11 +927,11 @@ def daophot(
"PSF photometry cannot currently be run across seperate windows"
)
wnam = wnam1
- print(wnam)
+
# background stats from whole windpw
# estimate background RMS
wind = ccd[wnam]
-
+
rms_func = MADStdBackgroundRMS(sigma_clip=SigmaClip(sigma=rejthresh))
bkg_rms = rms_func(wind.data)
bkg_func = MMMBackground(sigma_clip=SigmaClip(sigma=rejthresh))
@@ -895,15 +940,23 @@ def daophot(
# crop window to ROI
wind = ccd[wnam].window(xlo, xhi, ylo, yhi)
+ # uncertainties
+ sigma = np.sqrt(read**2 + np.maximum(0, wind.data) / gain)
# correct FWHM for binning
fwhm /= wind.xbin
if method == "m":
- psf_model = MoffatPSF(fwhm, beta)
+ psf_model = MoffatPSF(beta=beta, x_fwhm=fwhm, y_fwhm=fwhm)
print(" FWHM = {:.1f}, BETA={:.1f}".format(fwhm, beta))
+ psf_model_name = 'moffat'
else:
- psf_model = IntegratedGaussianPRF(sigma=fwhm * gaussian_fwhm_to_sigma)
+ psf_model = GaussianPRF(x_fwhm=fwhm, y_fwhm=fwhm)
+ psf_model.x_fwhm.fixed = False
+ psf_model.y_fwhm.fixed = False
+ psf_model.theta.fixed = False
+ psf_model.theta.bounds = (-90, 90)
print(" FWHM = {:.1f}".format(fwhm))
+ psf_model_name = 'gaussian'
# define region to extract around positions for fits
fitshape = int(5 * fwhm)
@@ -911,30 +964,58 @@ def daophot(
if fitshape % 2 == 0:
fitshape += 1
- photometry_task = DAOPhotPSFPhotometry(
- gfac * fwhm,
- thresh * bkg_rms,
- fwhm,
- psf_model,
- fitshape,
- niters=niters,
- sigma=rejthresh,
+ # Step 1: fit the reference stars to determine the PSF parameters.
+ # get pixel positions of reference apertures
+ xpos, ypos = list(zip(*[(wind.x_pixel(aper.x), wind.y_pixel(aper.y)) for aper in ccdaper.values()]))
+ reference_positions = Table(names=['x_0', 'y_0'], data=(xpos, ypos))
+
+ photometry_task = PSFPhotometry(
+ psf_model=psf_model,
+ aperture_radius=1.7*fwhm,
+ fit_shape=fitshape,
+ )
+ # do the PSF photometry
+ photom_results = photometry_task(
+ wind.data - bkg, error=sigma, init_params=reference_positions
+ )
+ colnames = [col for col in photom_results.colnames if 'init' not in col]
+ print(photom_results[colnames])
+
+ # now create the PSF model from fitting the reference stars
+ psf_model = create_psf_model(
+ photom_results, psf_model_name, fixed_positions=True
+ )
+ print(psf_model)
+ # Step 2: run the FIND-FIT-SUBTRACT loop to find all stars in the region.
+ photometry_task = IterativePSFPhotometry(
+ psf_model=psf_model,
+ fit_shape=fitshape,
+ finder=DAOStarFinder(threshold=thresh * bkg_rms, fwhm=fwhm),
+ grouper=SourceGrouper(gfac * fwhm),
+ maxiters=niters,
+ mode='all',
+ aperture_radius=1.7*fwhm,
+ localbkg_estimator=None,
)
with warnings.catch_warnings():
warnings.simplefilter("ignore")
- results = photometry_task(wind.data - bkg)
+ results = photometry_task(wind.data - bkg, error=sigma)
+ colnames = [col for col in results.colnames if 'fit' in col]
+ print(results[colnames])
# filter out junk fits
+ results = results[results["flags"] == 0]
+
tiny = 1e-30
bad_errs = (
- (results["flux_unc"] < tiny)
- | (results["x_0_unc"] < tiny)
- | (results["y_0_unc"] < tiny)
+ (results["flux_err"] < tiny)
+ | (results["x_err"] < tiny)
+ | (results["y_err"] < tiny)
)
results = results[~bad_errs]
- results.write("table_{}.fits".format(cnam))
+ results.write("table_{}.fits".format(cnam), overwrite=True)
print(" found {} stars".format(len(results)))
xlocs, ylocs = results["x_fit"], results["y_fit"]
diff --git a/hipercam/scripts/reduce.py b/hipercam/scripts/reduce.py
index eed59fc..cfd264c 100644
--- a/hipercam/scripts/reduce.py
+++ b/hipercam/scripts/reduce.py
@@ -4,7 +4,6 @@
from functools import partial
import numpy as np
-
from trm import cline
from trm.cline import Cline
@@ -155,6 +154,10 @@ def reduce(args=None):
implot : bool
flag to indicate you want to plot images.
+ psfplot : bool
+ flag to indicate you want to plot images after subtraction of PSF model.
+ inherits all other parameters from the image plot
+
ccd : string [if implot]
CCD(s) to plot, '0' for all, '1 3' to plot '1' and '3' only, etc.
@@ -225,6 +228,7 @@ def reduce(args=None):
cl.register("tkeep", Cline.GLOBAL, Cline.PROMPT)
cl.register("lplot", Cline.LOCAL, Cline.PROMPT)
cl.register("implot", Cline.LOCAL, Cline.PROMPT)
+ cl.register("psfplot", Cline.LOCAL, Cline.PROMPT)
cl.register("ccd", Cline.LOCAL, Cline.PROMPT)
cl.register("nx", Cline.LOCAL, Cline.PROMPT)
cl.register("msub", Cline.GLOBAL, Cline.PROMPT)
@@ -318,8 +322,14 @@ def reduce(args=None):
lplot = cl.get_value("lplot", "do you want to plot light curves?", True)
implot = cl.get_value("implot", "do you want to plot images?", True)
+ if rfile['psf_photom'].get('use_psf', 'no').lower().strip() != 'no':
+ psfplot = cl.get_value(
+ "psfplot", "do you want to plot residual of PSF model?", True
+ )
+ else:
+ psfplot = False
- if implot:
+ if implot or psfplot:
# define the panel grid. first get the labels and maximum
# dimensions
ccdinf = spooler.get_ccd_pars(source, resource)
@@ -403,8 +413,8 @@ def reduce(args=None):
else:
plot_lims = None
- imdev, lcdev, spanel, tpanel, xpanel, ypanel, lpanel = setup_plots(
- rfile, ccds, nx, plot_lims, implot, lplot
+ imdev, lcdev, psfdev, spanel, tpanel, xpanel, ypanel, lpanel = setup_plots(
+ rfile, ccds, nx, plot_lims, implot, lplot, psfplot
)
# a couple of initialisations
@@ -468,10 +478,13 @@ def reduce(args=None):
update_plots(
results,
rfile,
+ processor.store,
implot,
lplot,
+ psfplot,
imdev,
lcdev,
+ psfdev,
pccd,
ccds,
msub,
@@ -535,10 +548,13 @@ def reduce(args=None):
update_plots(
results,
rfile,
+ processor.store,
implot,
lplot,
+ psfplot,
imdev,
lcdev,
+ psfdev,
pccd,
ccds,
msub,
@@ -711,10 +727,13 @@ def reduce(args=None):
update_plots(
results,
rfile,
+ processor.store,
implot,
lplot,
+ psfplot,
imdev,
lcdev,
+ psfdev,
pccds[-1],
ccds,
msub,
@@ -759,10 +778,13 @@ def reduce(args=None):
update_plots(
results,
rfile,
+ processor.store,
implot,
lplot,
+ psfplot,
imdev,
lcdev,
+ psfdev,
pccd,
ccds,
msub,
diff --git a/pyproject.toml b/pyproject.toml
index 1908dac..8e0e506 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -40,6 +40,7 @@ dependencies = [
[project.optional-dependencies]
dev = ["check-manifest"]
test = ["coverage", "nose"]
+psf = ["photutils>=2.1"]
[project.urls]
Homepage = "https://github.com/HiPERCAM/hipercam"
@@ -94,6 +95,7 @@ ncal = "hipercam.scripts.ncal:ncal"
pbands = "hipercam.scripts.pbands:pbands"
pfolder = "hipercam.scripts.pfolder:pfolder"
plog = "hipercam.scripts.plog:plog"
+psfaper = "hipercam.scripts.psfaper:psfaper"
redanal = "hipercam.scripts.redanal:redanal"
redplt = "hipercam.scripts.redplt:redplt"
reduce = "hipercam.scripts.reduce:reduce"