Dev

opticam/align.py

@@ -103,13 +103,11 @@ def align_batch(
     background_rmss = {}
 
     for file in batch:
-        data = np.asarray(
-            get_data(
-                file,
-                flat_corrector=flat_corrector,
-                rebin_factor=rebin_factor,
-                return_error=False,
-                remove_cosmic_rays=remove_cosmic_rays),
+        data = get_data(
+            file,
+            flat_corrector=flat_corrector,
+            rebin_factor=rebin_factor,
+            remove_cosmic_rays=remove_cosmic_rays,
             )
 
         # calculate and subtract background

opticam/background/global_background.py

@@ -10,7 +10,10 @@ class BaseBackground(ABC):
     Base class for OPTICAM background estimators.
     """
 
-    def __init__(self, box_size: int | Tuple[int, int]):
+    def __init__(
+        self,
+        box_size: int | Tuple[int, int],
+        ):
         """
         Initialize a background estimator.
 
@@ -28,7 +31,10 @@ def __init__(self, box_size: int | Tuple[int, int]):
         self.box_size = box_size
 
     @abstractmethod
-    def __call__(self, image: NDArray) -> Background2D:
+    def __call__(
+        self,
+        image: NDArray,
+        ) -> Background2D:
         """
         Compute the 2D background for an image.
 

opticam/background/local_background.py

@@ -39,7 +39,6 @@ def __init__(
     def __call__(
         self,
         data: NDArray,
-        error: NDArray,
         position: NDArray,
         semimajor_axis: float,
         semiminor_axis: float,
@@ -51,8 +50,6 @@ def __call__(
         ----------
         data : NDArray
             The image data.
-        error : NDArray
-            The error in the image data.
         semimajor_axis : float
             The (unscaled) semi-major axis of the aperture.
         semiminor_axis : float
@@ -79,7 +76,6 @@ class DefaultLocalBackground(BaseLocalBackground):
     def __call__(
         self,
         data: NDArray,
-        error: NDArray,
         position: NDArray,
         semimajor_axis: float,
         semiminor_axis: float | None = None,
@@ -125,7 +121,6 @@ def __call__(
         stats = ApertureStats(
             data,
             annulus,
-            error=error,
             sigma_clip=self.sigma_clip,
             )
 

opticam/fitting/models.py

@@ -0,0 +1,83 @@
+import numpy as np
+from numpy.typing import NDArray
+
+
+def gaussian(
+    x: NDArray,
+    amp: float,
+    mu: float,
+    std: float,
+    ) -> NDArray:
+    """
+    Gaussian function.
+
+    Parameters
+    ----------
+    x : NDArray
+        The function inputs.
+    amplitude : float
+        The height of the function peak.
+    mu : float
+        The mean of the Gaussian.
+    std : float
+        The standard deviation of the Gaussian.
+
+    Returns
+    -------
+    NDArray
+        The Gaussian evaluated at the inputs.
+    """
+
+    return amp * np.exp(-.5 * ((x - mu) / std)**2)
+
+
+def power_law(
+    x: NDArray,
+    norm: float,
+    exp: float,
+    ) -> NDArray:
+    """
+    Power law function.
+
+    Parameters
+    ----------
+    x : NDArray
+        The function inputs.
+    norm : float
+        The normalisation.
+    exp : float
+        The exponent.
+
+    Returns
+    -------
+    NDArray
+        The power law evaluated at the inputs.
+    """
+
+    return norm*x**exp
+
+
+def straight_line(
+    x: NDArray,
+    m: float,
+    c: float,
+    ) -> NDArray:
+    """
+    Straight line function.
+
+    Parameters
+    ----------
+    x : NDArray
+        The function inputs.
+    m : float
+        The gradient.
+    c : float
+        The intercept.
+
+    Returns
+    -------
+    NDArray
+        The straight line evaluated at the inputs.
+    """
+
+    return m*x + c

opticam/fitting/routines.py

@@ -0,0 +1,96 @@
+from typing import Dict
+
+import numpy as np
+from numpy.typing import NDArray
+from scipy.optimize import curve_fit
+
+from opticam.fitting.models import power_law, straight_line
+
+def fit_rms_vs_flux(
+    data: Dict,
+    ) -> Dict[str, Dict[str, NDArray]]:
+    """
+    Iteratively fit a straight line (in log space) to the RMS vs flux plots for each catalog. This can be used to
+    identify variable sources and good comparison sources.
+
+    Parameters
+    ----------
+    data : Dict
+        The RMS vs flux data.
+
+    Returns
+    -------
+    Dict[str, Dict[str, NDArray]]
+        The power law fits for each filter `{filter: {'flux': NDArray, 'rms': NDArray}}`.
+    """
+
+    pl_fits = {}
+
+    for fltr in data.keys():
+        rms, flux = [], []
+        for values in data[fltr].values():
+            rms.append(values['rms'])
+            flux.append(values['flux'])
+
+        order = np.argsort(flux)
+        x = np.array(flux)[order]
+        y = np.array(rms)[order]
+
+        try:
+            converged = False
+            prev, prev_err = None, None
+            while not converged:
+                log_x = np.log10(flux)
+                log_y = np.log10(rms)
+
+                popt, pcov = curve_fit(
+                    straight_line,
+                    log_x,
+                    log_y,
+                    )
+                perr = np.sqrt(np.diag(pcov))
+
+                if prev is not None and prev_err is not None:
+                    # assume fit has converged is params are within 20%
+                    # do not use perr since it may be very large
+                    converged = np.allclose(popt, prev, rtol=0.2, atol=0)
+
+                # remove largest outliers
+                model = straight_line(log_x, *popt)
+                r = log_y - model
+                i = np.argmax(r)
+
+                rms.pop(i)
+                flux.pop(i)
+
+                prev = popt
+                prev_err = perr
+        except:
+            popt, pcov = curve_fit(
+                    straight_line,
+                    np.log10(x),
+                    np.log10(y),
+                    )
+
+        # get model prediction band using Monte Carlo method
+        N = 1000
+        y_models = np.zeros((N, x.size))
+        for i in range(N):
+            rng = np.random.default_rng(i)
+            a, b = rng.multivariate_normal(popt, pcov)
+            y_models[i] += power_law(x, 10**a, b)
+        y_model_err = 5 * np.std(y_models, axis=0)  # 5 sigma error
+
+        y_model = power_law(
+            x,
+            10**popt[1],
+            popt[0],
+            )
+
+        pl_fits[fltr] = {
+            'flux': x,
+            'rms': y_model,
+            'err': y_model_err,
+        }
+
+    return pl_fits