Source code for rubin_sim.maf.metrics.incremental_template_metric

import numpy as np
from rubin_scheduler.utils import SURVEY_START_MJD, Site, SysEngVals, calc_season

from .base_metric import BaseMetric

__all__ = ["TemplateTime"]


[docs] class TemplateTime(BaseMetric): """Find the time at which we expect to hit incremental template availability. Note that there are some complications to real template generation that make this an approximation and not an exact answer -- one aspect is that templates are generated in `patches` and not per pixel. However, it may be possible to generate parts of these patches at about 5 arcsecond scales, which implies running with a healpix slicer at nside=512 or 1024. Parameters ---------- n_images_in_template : `int`, opt Number of qualified visits required for template generation. Default 3. seeing_threshold : `float`, opt The threshold for counting "good seeing" images, in arcseconds. This is at zenith, at 500 nm - it will be scaled to the slicer's best possible value at the relevant declination and in the template's bandpass. seeing_percentile : `float`, opt When evaluating the local seeing distribution, if there `n_images_in_template` better than `seeing_percentile`, go ahead and make the template. This allows template creation in locations with poor seeing compared to Maximum percentile seeing to allow in the qualified images (0 - 100). The maximum of seeing_threshold (scaled) or the seeing at seeing_percentil will be used as the cutoff for template creation. m5_threshold : `float` or None, opt The threshold for counting an image as "good m5". This is also scaled to the expected best minimum airmass/seeing given the slicepoints' declination. m5_percentile : `float`, opt Maximum percentile m5 to allow in the qualified images (0 - 100). The minimum of m5_threshold (scaled) or m5 resulting from m5_percentile will be used as the cutoff for template creation. seeing_col : `str`, opt Name of the seeing column to use. m5_col : `str`, opt Name of the five sigma depth columns. night_col : `str`, opt Name of the column describing the night of the visit. mjd_col : `str`, opt Name of column describing time of the visit band_col : `str`, opt Name of column describing band mjd_start : `float`, opt The starting time of the survey, to use to count "nights before template creation", etc. site_latiude_rad : `float`, opt The latitude of the site, used to determine minimum best airmass. """ def __init__( self, n_images_in_template=3, seeing_threshold=0.8, m5_threshold=None, seeing_percentile=50, m5_percentile=50, seeing_col="seeingFwhmEff", m5_col="fiveSigmaDepth", night_col="night", mjd_col="observationStartMJD", band_col="band", mjd_start=None, filter_col=None, site_latitude_rad=Site("LSST").latitude_rad, **kwargs, ): self.n_images_in_template = n_images_in_template self.seeing_threshold = seeing_threshold if m5_threshold is None: # Update these - currently just m5 defaults for v1.9 m5_threshold = {"u": 23.7, "g": 24.9, "r": 24.5, "i": 24.1, "z": 23.5, "y": 22.5} self.m5_threshold = m5_threshold self.seeing_percentile = seeing_percentile self.m5_percentile = m5_percentile self.seeing_col = seeing_col self.m5_col = m5_col self.night_col = night_col self.mjd_col = mjd_col self.band_col = band_col if filter_col is not None: self.band_col = filter_col if "metric_name" in kwargs: self.metric_name = kwargs["metric_name"] del kwargs["metric_name"] else: self.metric_name = "TemplateTime" if mjd_start is None: self.mjd_start = SURVEY_START_MJD else: self.mjd_start = mjd_start self.site_latitude_rad = site_latitude_rad sev = SysEngVals() self.eff_wavelens = dict([(b, sev.eff_wavelengths[b]) for b in "ugrizy"]) self.extinction_coeffs = { "u": -0.46, "g": -0.21, "r": -0.12, "i": -0.07, "z": -0.06, "y": -0.10, } super().__init__( col=[self.seeing_col, self.m5_col, self.night_col, self.mjd_col, self.band_col], metric_name=self.metric_name, units="days", **kwargs, )
[docs] def run(self, data_slice, slice_point): result = {} # Bail if not enough visits at all if len(data_slice) < self.n_images_in_template: return self.badval # Check that the visits are sorted in time data_slice.sort(order=self.mjd_col) min_z_possible = np.abs(slice_point["dec"] - self.site_latitude_rad) min_airmass_possible = 1.0 / np.cos(min_z_possible) best_seeing_possible = self.seeing_threshold * np.power(min_airmass_possible, 0.6) current_band = np.unique(data_slice[self.band_col])[0] best_seeing_possible = best_seeing_possible * np.power(500 / self.eff_wavelens[current_band], 0.3) percentile_seeing = np.percentile(data_slice[self.seeing_col], self.seeing_percentile) bench_seeing = np.max([best_seeing_possible, percentile_seeing]) try: m5_threshold = self.m5_threshold[current_band] except TypeError: m5_threshold = self.m5_threshold best_m5 = ( m5_threshold - self.extinction_coeffs[current_band] * (1 - min_airmass_possible) + 2.5 * np.log10(self.seeing_threshold / best_seeing_possible) ) percentile_m5 = np.percentile(data_slice[self.m5_col], 100 - self.m5_percentile) bench_m5 = np.min([best_m5, percentile_m5]) seeing_ok = np.where(data_slice[self.seeing_col] <= bench_seeing, True, False) m5_ok = np.where(data_slice[self.m5_col] >= bench_m5, True, False) ok_template_input = np.where(seeing_ok & m5_ok)[0] if ( len(ok_template_input) < self.n_images_in_template ): # If seeing_ok and/or m5_ok are "false", returned as bad value return self.badval idx_template_created = ok_template_input[ self.n_images_in_template - 1 ] # Last image needed for template idx_template_inputs = ok_template_input[: self.n_images_in_template] # Images included in template night_template_created = data_slice[self.night_col][idx_template_created] frac_into_season_template_created = calc_season( np.degrees(slice_point["ra"]), data_slice[self.mjd_col][idx_template_created], mjd_start=self.mjd_start, ) where_template = ( data_slice[self.night_col] > night_template_created ) # of later images where we have a template n_images_with_template = np.sum(where_template) template_m5 = 1.25 * np.log10(np.sum(10.0 ** (0.8 * data_slice[self.m5_col][idx_template_inputs]))) # derive variance-weighted PSF width # 1. normalize the weights template_f0 = np.sqrt(np.sum(25 * 10 ** (0.8 * data_slice[self.m5_col][idx_template_inputs]))) # 2. compute per-input noise image_noise = template_f0 / 5 * 10 ** (-0.4 * data_slice[self.m5_col]) # 3. used variance-weighted sum of squares to derive template seeing template_seeing = np.sqrt( np.sum( (data_slice[self.seeing_col][idx_template_inputs] / image_noise[idx_template_inputs]) ** 2.0 ) ) delta_m5 = -2.5 * np.log10(np.sqrt(1.0 + 10 ** (-0.8 * (template_m5 - data_slice[self.m5_col])))) diff_m5s = data_slice[self.m5_col] + delta_m5 n_alerts_per_diffim = 1e4 * 10 ** (0.6 * (diff_m5s - 24.7)) result["bench_seeing"] = bench_seeing result["bench_m5"] = bench_m5 result["night_template_created"] = night_template_created result["frac_into_season_template_created"] = frac_into_season_template_created result["n_images_until_template"] = idx_template_created + 1 result["n_images_with_template"] = n_images_with_template result["template_m5"] = template_m5 result["template_seeing"] = template_seeing result["total_alerts"] = np.sum(n_alerts_per_diffim[where_template]) result["template_input_m5s"] = data_slice[self.m5_col][idx_template_inputs] result["template_input_seeing"] = data_slice[self.seeing_col][idx_template_inputs] result["fraction_better_template_seeing"] = np.sum( (data_slice[self.seeing_col][where_template] > template_seeing) ) / np.sum(where_template) result["diffim_lc"] = { "mjd": data_slice[self.mjd_col][where_template], "night": data_slice[self.night_col][where_template], "band": data_slice[self.band_col][where_template], "diff_m5": diff_m5s[where_template], "science_m5": data_slice[self.m5_col][where_template], "template_m5": template_m5 * np.ones(np.sum(where_template)), "science_seeing": data_slice[self.seeing_col][where_template], "template_seeing": template_seeing * np.ones(np.sum(where_template)), "n_alerts": n_alerts_per_diffim[where_template], } return result
def reduce_night_template_created(self, metric_val): # returns night of template creation return metric_val["night_template_created"] def reduce_n_images_until_template( self, metric_val ): # returns number of images needed to complete template return metric_val["n_images_until_template"] def reduce_n_images_with_template(self, metric_val): # returns number of images with a template return metric_val["n_images_with_template"] def reduce_template_m5(self, metric_val): # calculated coadded m5 of resulting template return metric_val["template_m5"] def reduce_template_seeing(self, metric_val): # calculated seeing of resulting template return metric_val["template_seeing"] def reduce_fraction_better_template_seeing(self, metric_val): # calculated seeing of resulting template return metric_val["fraction_better_template_seeing"] def reduce_total_alerts(self, metric_val): return metric_val["total_alerts"] def reduce_frac_into_season(self, metric_val): return metric_val["frac_into_season_template_created"]