Source code for lsst.sims.maf.metrics.crowdingMetric

import numpy as np
from scipy.interpolate import interp1d
from lsst.sims.maf.metrics import BaseMetric
import healpy as hp

# Modifying from Knut Olson's fork at:
# https://github.com/knutago/sims_maf_contrib/blob/master/tutorials/CrowdingMetric.ipynb

__all__ = ['CrowdingM5Metric', 'CrowdingMagUncertMetric', 'NstarsMetric']


def _compCrowdError(magVector, lumFunc, seeing, singleMag=None):
    """
    Compute the photometric crowding error given the luminosity function and best seeing.

    Parameters
    ----------
    magVector : np.array
        Stellar magnitudes.
    lumFunc : np.array
        Stellar luminosity function.
    seeing : float
        The best seeing conditions. Assuming forced-photometry can use the best seeing conditions
        to help with confusion errors.
    singleMag : float (None)
        If singleMag is None, the crowding error is calculated for each mag in magVector. If
        singleMag is a float, the crowding error is interpolated to that single value.

    Returns
    -------
    np.array
        Magnitude uncertainties.

    Equation from Olsen, Blum, & Rigaut 2003, AJ, 126, 452
    """
    lumAreaArcsec = 3600.0 ** 2
    lumVector = 10 ** (-0.4 * magVector)
    coeff = np.sqrt(np.pi / lumAreaArcsec) * seeing / 2.
    myInt = (np.add.accumulate((lumVector ** 2 * lumFunc)[::-1]))[::-1]
    temp = np.sqrt(myInt) / lumVector
    if singleMag is not None:
        interp = interp1d(magVector, temp)
        temp = interp(singleMag)
    crowdError = coeff * temp
    return crowdError


[docs]class CrowdingM5Metric(BaseMetric): """Return the magnitude at which the photometric error exceeds crowding_error threshold. """ def __init__(self, crowding_error=0.1, filtername='r', seeingCol='seeingFwhmGeom', metricName=None, maps=['StellarDensityMap'], **kwargs): """ Parameters ---------- crowding_error : float, opt The magnitude uncertainty from crowding in magnitudes. Default 0.1 mags. filtername: str, opt The bandpass in which to calculate the crowding limit. Default r. seeingCol : str, opt The name of the seeing column. m5Col : str, opt The name of the m5 depth column. maps : list of str, opt Names of maps required for the metric. Returns ------- float The magnitude of a star which has a photometric error of `crowding_error` """ cols = [seeingCol] units = 'mag' self.crowding_error = crowding_error self.filtername = filtername self.seeingCol = seeingCol if metricName is None: metricName = 'Crowding to Precision %.2f' % (crowding_error) super().__init__(col=cols, maps=maps, units=units, metricName=metricName, **kwargs)
[docs] def run(self, dataSlice, slicePoint=None): # Set magVector to the same length as starLumFunc (lower edge of mag bins) magVector = slicePoint[f'starMapBins_{self.filtername}'][1:] # Pull up density of stars at this point in the sky lumFunc = slicePoint[f'starLumFunc_{self.filtername}'] # Calculate the crowding error using the best seeing value (in any filter?) crowdError = _compCrowdError(magVector, lumFunc, seeing=min(dataSlice[self.seeingCol])) # Locate at which point crowding error is greater than user-defined limit aboveCrowd = np.where(crowdError >= self.crowding_error)[0] if np.size(aboveCrowd) == 0: result = max(magVector) else: crowdMag = magVector[max(aboveCrowd[0]-1, 0)] result = crowdMag return result
[docs]class NstarsMetric(BaseMetric): """Return the number of stars visible above some uncertainty limit, taking image depth and crowding into account. """ def __init__(self, crowding_error=0.1, filtername='r', seeingCol='seeingFwhmGeom', m5Col='fiveSigmaDepth', metricName=None, maps=['StellarDensityMap'], ignore_crowding=False, **kwargs): """ Parameters ---------- crowding_error : float, opt The magnitude uncertainty from crowding in magnitudes. Default 0.1 mags. filtername: str, opt The bandpass in which to calculate the crowding limit. Default r. seeingCol : str, opt The name of the seeing column. m5Col : str, opt The name of the m5 depth column. maps : list of str, opt Names of maps required for the metric. ignore_crowding : bool (False) Ignore the cowding limit. Returns ------- float The number of stars above the error limit """ cols = [seeingCol, m5Col] units = 'N stars' self.crowding_error = crowding_error self.m5Col = m5Col self.filtername = filtername self.seeingCol = seeingCol self.ignore_crowding = ignore_crowding if metricName is None: metricName = 'N stars to Precision %.2f' % (crowding_error) super().__init__(col=cols, maps=maps, units=units, metricName=metricName, **kwargs)
[docs] def run(self, dataSlice, slicePoint=None): pix_area = hp.nside2pixarea(slicePoint['nside'], degrees=True) # Set magVector to the same length as starLumFunc (lower edge of mag bins) magVector = slicePoint[f'starMapBins_{self.filtername}'][1:] # Pull up density of stars at this point in the sky lumFunc = slicePoint[f'starLumFunc_{self.filtername}'] # Calculate the crowding error using the best seeing value (in any filter?) crowdError = _compCrowdError(magVector, lumFunc, seeing=min(dataSlice[self.seeingCol])) # Locate at which point crowding error is greater than user-defined limit aboveCrowd = np.where(crowdError >= self.crowding_error)[0] if np.size(aboveCrowd) == 0: crowdMag = max(magVector) else: crowdMag = magVector[max(aboveCrowd[0]-1, 0)] # Compute the coadded depth, and the mag where that depth hits the error specified coadded_depth = 1.25 * np.log10(np.sum(10.**(.8*dataSlice[self.m5Col]))) mag_limit = -2.5*np.log10(1./(self.crowding_error*(1.09*5)))+coadded_depth # Use the shallower depth, crowding or coadded if self.ignore_crowding: min_mag = mag_limit else: min_mag = np.min([crowdMag, mag_limit]) # Interpolate to the number of stars result = np.interp(min_mag, slicePoint[f'starMapBins_{self.filtername}'][1:], slicePoint[f'starLumFunc_{self.filtername}']) * pix_area return result
[docs]class CrowdingMagUncertMetric(BaseMetric): """ Given a stellar magnitude, calculate the mean uncertainty on the magnitude from crowding. """ def __init__(self, rmag=20., seeingCol='seeingFwhmGeom', units='mag', metricName=None, filtername='r', maps=['StellarDensityMap'], **kwargs): """ Parameters ---------- rmag : float The magnitude of the star to consider. Returns ------- float The uncertainty in magnitudes caused by crowding for a star of rmag. """ self.filtername = filtername self.seeingCol = seeingCol self.rmag = rmag if metricName is None: metricName = 'CrowdingError at %.2f' % (rmag) super().__init__(col=[seeingCol], maps=maps, units=units, metricName=metricName, **kwargs)
[docs] def run(self, dataSlice, slicePoint=None): magVector = slicePoint[f'starMapBins_{self.filtername}'][1:] lumFunc = slicePoint[f'starLumFunc_{self.filtername}'] # Magnitude uncertainty given crowding dmagCrowd = _compCrowdError(magVector, lumFunc, dataSlice[self.seeingCol], singleMag=self.rmag) result = np.mean(dmagCrowd) return result