Curve of Growth: Support different aperture shapes?

[pllim]

As of Jdaviz v4.1, we have this way of calculating curve of growth in our Aperture Photometry plugin. In replacing that with what is in photutils/profiles/curve_of_growth.py as of photutils v2.1, we lost the ability to compute the growth using selected aperture shape (particularly ellipse and rectangle). However, adding this feature to existing class might break photutils API, so maybe a new class or function in photutils/profiles/curve_of_growth.py instead if such a feature is useful for photutils proper?

Here is the code from Jdaviz (needs rewriting to be more general):

def _curve_of_growth(data, centroid, aperture, final_sum, wcs=None, background=0,
                     n_datapoints=10, pixarea_fac=None, display_unit=None, equivalencies=[]):
    """Calculate curve of growth for aperture photometry.

    Parameters
    ----------
    data : ndarray or `~astropy.units.Quantity`
        Data for the calculation.

    centroid : tuple of int
        ``ApertureStats`` centroid or desired center in ``(x, y)``.

    aperture : obj
        ``photutils`` aperture to use, except its center will be
        changed to the given ``centroid``. This is because the aperture
        might be hand-drawn and a more accurate centroid has been
        recalculated separately.

    final_sum : float or `~astropy.units.Quantity`
        Aperture sum that is already calculated in the
        main plugin above.

    wcs : obj or `None`
        Supported WCS objects or `None`.

    background : float or `~astropy.units.Quantity`
        Background to subtract, if any. Unit must match ``data``.

    n_datapoints : int
        Number of data points in the curve.

    pixarea_fac : float or `None`
        For ``flux_unit/sr`` to ``flux_unit`` conversion.

    display_unit : str or None
        (For cubeviz only to deal with display unit conversion). Desired unit
        for output. If unit is a surface brightness, a Flux unit will be
        returned if pixarea_fac is provided.

    Returns
    -------
    x_arr : ndarray
        Data for X-axis of the curve.

    sum_arr : ndarray or `~astropy.units.Quantity`
        Data for Y-axis of the curve.

    x_label, y_label : str
        X- and Y-axis labels, respectively.

    Raises
    ------
    TypeError
        Unsupported aperture.

    """
    n_datapoints += 1  # n + 1

    # determined desired unit for output sum array and y label
    # cubeviz only to handle unit conversion display unit changes
    if display_unit is not None:
        sum_unit = u.Unit(display_unit)
    else:
        if isinstance(data, u.Quantity):
            sum_unit = data.unit
        else:
            sum_unit = None
    if sum_unit and pixarea_fac is not None:
        # multiply data unit by its solid angle to convert sum in sb to sum in flux
        sum_unit *= check_if_unit_is_per_solid_angle(sum_unit, return_unit=True)

    if hasattr(aperture, 'to_pixel'):
        aperture = aperture.to_pixel(wcs)

    if isinstance(aperture, CircularAperture):
        x_label = 'Radius (pix)'
        x_arr = np.linspace(0, aperture.r, num=n_datapoints)[1:]
        aper_list = [CircularAperture(centroid, cur_r) for cur_r in x_arr[:-1]]
    elif isinstance(aperture, EllipticalAperture):
        x_label = 'Semimajor axis (pix)'
        x_arr = np.linspace(0, aperture.a, num=n_datapoints)[1:]
        a_arr = x_arr[:-1]
        b_arr = aperture.b * a_arr / aperture.a
        aper_list = [EllipticalAperture(centroid, cur_a, cur_b, theta=aperture.theta)
                     for (cur_a, cur_b) in zip(a_arr, b_arr)]
    elif isinstance(aperture, RectangularAperture):
        x_label = 'Width (pix)'
        x_arr = np.linspace(0, aperture.w, num=n_datapoints)[1:]
        w_arr = x_arr[:-1]
        h_arr = aperture.h * w_arr / aperture.w
        aper_list = [RectangularAperture(centroid, cur_w, cur_h, theta=aperture.theta)
                     for (cur_w, cur_h) in zip(w_arr, h_arr)]
    else:
        raise TypeError(f'Unsupported aperture: {aperture}')

    sum_arr = [ApertureStats(data, cur_aper, wcs=wcs, local_bkg=background).sum
               for cur_aper in aper_list]
    if isinstance(sum_arr[0], u.Quantity):
        sum_arr = u.Quantity(sum_arr)
    else:
        sum_arr = np.array(sum_arr)
    if pixarea_fac is not None:
        sum_arr = sum_arr * pixarea_fac
    if isinstance(final_sum, u.Quantity):
        final_sum = flux_conversion_general(final_sum.value, final_sum.unit,
                                            sum_arr.unit, equivalencies)
    sum_arr = np.append(sum_arr, final_sum)

    if sum_unit is None:
        y_label = 'Value'
    else:
        y_label = sum_unit.to_string()
        # bqplot does not like Quantity
        sum_arr = flux_conversion_general(sum_arr.value, sum_arr.unit, sum_unit,
                                          equivalencies, with_unit=False)

    return x_arr, sum_arr, x_label, y_label

🐱