Skip to content

Commit 7edf784

Browse files
PFLegetPFLeget
committed
Compute TEX metric.
1 parent 511e8cf commit 7edf784

File tree

2 files changed

+241
-4
lines changed

2 files changed

+241
-4
lines changed

‎python/lsst/pipe/tasks/computeExposureSummaryStats.py‎

Lines changed: 229 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -35,6 +35,7 @@
3535
import lsst.afw.image as afwImage
3636
import lsst.geom as geom
3737
from lsst.meas.algorithms import ScienceSourceSelectorTask
38+
from lsst.meas.algorithms.computeExPsf import ComputeExPsfTask, ComputeExPsfConfig
3839
from lsst.utils.timer import timeMethod
3940
import lsst.ip.isr as ipIsr
4041

@@ -137,6 +138,22 @@ class ComputeExposureSummaryStatsConfig(pexConfig.Config):
137138
doc="Signal-to-noise ratio for computing the magnitude limit depth.",
138139
default=5.0
139140
)
141+
psfTE1TreecorrConfig = pexConfig.ConfigField(
142+
dtype=ComputeExPsfConfig,
143+
doc="Treecorr config for computing scalar value of TE1.",
144+
)
145+
psfTE2TreecorrConfig = pexConfig.ConfigField(
146+
dtype=ComputeExPsfConfig,
147+
doc="Treecorr config for computing scalar value of TE2.",
148+
)
149+
psfTE3TreecorrConfig = pexConfig.ConfigField(
150+
dtype=ComputeExPsfConfig,
151+
doc="Treecorr config for computing scalar value of TE3.",
152+
)
153+
psfTE4TreecorrConfig = pexConfig.ConfigField(
154+
dtype=ComputeExPsfConfig,
155+
doc="Treecorr config for computing scalar value of TE4.",
156+
)
140157

141158
def setDefaults(self):
142159
super().setDefaults()
@@ -159,6 +176,22 @@ def setDefaults(self):
159176
self.starSelector.signalToNoise.fluxField = "slot_PsfFlux_instFlux"
160177
self.starSelector.signalToNoise.errField = "slot_PsfFlux_instFluxErr"
161178

179+
min_theta = [1e-6, 5.0, 1e-6, 5.0]
180+
max_theta = [1.0, 180.0, 5.0, 12.0]
181+
TExTreecorrConfig = [
182+
self.psfTE1TreecorrConfig,
183+
self.psfTE2TreecorrConfig,
184+
self.psfTE3TreecorrConfig,
185+
self.psfTE4TreecorrConfig,
186+
]
187+
188+
for texc, mint, maxt in zip(TExTreecorrConfig, min_theta, max_theta):
189+
texc.treecorr.min_sep = mint / 60.0
190+
texc.treecorr.max_sep = maxt / 60.0
191+
texc.treecorr.nbins = 1
192+
texc.treecorr.bin_type = "Linear"
193+
texc.treecorr.sep_units = "degree"
194+
162195

163196
class ComputeExposureSummaryStatsTask(pipeBase.Task):
164197
"""Task to compute exposure summary statistics.
@@ -200,6 +233,22 @@ class ComputeExposureSummaryStatsTask(pipeBase.Task):
200233
- psfTraceRadiusDelta
201234
- psfApFluxDelta
202235
236+
These quantities are computed as part of:
237+
https://rubinobs.atlassian.net/browse/DM-40780
238+
239+
- psfTE1e1
240+
- psfTE1e2
241+
- psfTE1ex
242+
- psfTE2e1
243+
- psfTE2e2
244+
- psfTE2ex
245+
- psfTE3e1
246+
- psfTE3e2
247+
- psfTE3ex
248+
- psfTE4e1
249+
- psfTE4e2
250+
- psfTE4ex
251+
203252
This quantity is computed based on the aperture correction map, the
204253
psfSigma, and the image mask to assess the robustness of the aperture
205254
corrections across a given detector:
@@ -219,6 +268,7 @@ def __init__(self, **kwargs):
219268
super().__init__(**kwargs)
220269

221270
self.makeSubtask("starSelector")
271+
self.isTEXComputationDone = False
222272

223273
@timeMethod
224274
def run(self, exposure, sources, background):
@@ -408,10 +458,13 @@ def update_psf_stats(
408458
psfE1 = (psfXX - psfYY)/(psfXX + psfYY)
409459
psfE2 = 2*psfXY/(psfXX + psfYY)
410460

411-
psfStarDeltaE1Median = np.median(starE1 - psfE1)
412-
psfStarDeltaE1Scatter = sigmaMad(starE1 - psfE1, scale='normal')
413-
psfStarDeltaE2Median = np.median(starE2 - psfE2)
414-
psfStarDeltaE2Scatter = sigmaMad(starE2 - psfE2, scale='normal')
461+
e1Residuals = starE1 - psfE1
462+
e2Residuals = starE2 - psfE2
463+
464+
psfStarDeltaE1Median = np.median(e1Residuals)
465+
psfStarDeltaE1Scatter = sigmaMad(e1Residuals, scale='normal')
466+
psfStarDeltaE2Median = np.median(e2Residuals)
467+
psfStarDeltaE2Scatter = sigmaMad(e2Residuals, scale='normal')
415468

416469
psfStarDeltaSizeMedian = np.median(starSize - psfSize)
417470
psfStarDeltaSizeScatter = sigmaMad(starSize - psfSize, scale='normal')
@@ -434,6 +487,178 @@ def update_psf_stats(
434487
)
435488
summary.maxDistToNearestPsf = float(maxDistToNearestPsf)
436489

490+
def comp_psf_TEX_visit_level(self, summary, sources, sources_is_astropy=False):
491+
"""Compute all summary-statistic fields at visit level for TEx metric of PSF.
492+
493+
Parameters
494+
----------
495+
summary : `lsst.afw.image.ExposureSummaryStats`
496+
Summary object to update in-place.
497+
sources : `lsst.afw.table.SourceCatalog` or `astropy.table.Table`
498+
Catalog for quantities that are computed from source table columns.
499+
If `None`, these quantities will be reset (generally to NaN).
500+
The type of this table must correspond to the
501+
``sources_is_astropy`` argument.
502+
sources_is_astropy : `bool`, optional
503+
Whether ``sources`` is an `astropy.table.Table` instance instead
504+
of an `lsst.afw.table.Catalog` instance. Default is `False` (the
505+
latter).
506+
"""
507+
508+
if self.isTEXComputationDone:
509+
510+
summary.psfTE1e1 = self.psfTE1e1
511+
summary.psfTE1e2 = self.psfTE1e2
512+
summary.psfTE1ex = self.psfTE1ex
513+
summary.psfTE2e1 = self.psfTE2e1
514+
summary.psfTE2e2 = self.psfTE2e2
515+
summary.psfTE2ex = self.psfTE2ex
516+
summary.psfTE3e1 = self.psfTE3e1
517+
summary.psfTE3e2 = self.psfTE3e2
518+
summary.psfTE3ex = self.psfTE3ex
519+
summary.psfTE4e1 = self.psfTE4e1
520+
summary.psfTE4e2 = self.psfTE4e2
521+
summary.psfTE4ex = self.psfTE4ex
522+
523+
else:
524+
525+
self.isTEXComputationDone = True
526+
527+
nan = float("nan")
528+
summary.psfTE1e1, self.psfTE1e1 = nan, nan
529+
summary.psfTE1e2, self.psfTE1e2 = nan, nan
530+
summary.psfTE1ex, self.psfTE1ex = nan, nan
531+
summary.psfTE2e1, self.psfTE2e1 = nan, nan
532+
summary.psfTE2e2, self.psfTE2e2 = nan, nan
533+
summary.psfTE2ex, self.psfTE2ex = nan, nan
534+
summary.psfTE3e1, self.psfTE3e1 = nan, nan
535+
summary.psfTE3e2, self.psfTE3e2 = nan, nan
536+
summary.psfTE3ex, self.psfTE3ex = nan, nan
537+
summary.psfTE4e1, self.psfTE4e1 = nan, nan
538+
summary.psfTE4e2, self.psfTE4e2 = nan, nan
539+
summary.psfTE4ex, self.psfTE4ex = nan, nan
540+
541+
psf_mask = self.starSelector.run(sources).selected
542+
543+
nPsfStarsUsedInStats = psf_mask.sum()
544+
545+
if nPsfStarsUsedInStats == 0:
546+
# No stars to measure statistics, so we must return the defaults
547+
# of 0 stars and NaN values.
548+
return
549+
550+
if sources_is_astropy:
551+
psf_cat = sources[psf_mask]
552+
else:
553+
psf_cat = sources[psf_mask].copy(deep=True)
554+
555+
starXX = psf_cat[self.config.starShape + '_xx']
556+
starYY = psf_cat[self.config.starShape + '_yy']
557+
starXY = psf_cat[self.config.starShape + '_xy']
558+
psfXX = psf_cat[self.config.psfShape + '_xx']
559+
psfYY = psf_cat[self.config.psfShape + '_yy']
560+
psfXY = psf_cat[self.config.psfShape + '_xy']
561+
562+
starE1 = (starXX - starYY)/(starXX + starYY)
563+
starE2 = 2*starXY/(starXX + starYY)
564+
565+
psfE1 = (psfXX - psfYY)/(psfXX + psfYY)
566+
psfE2 = 2*psfXY/(psfXX + psfYY)
567+
568+
e1Residuals = starE1 - psfE1
569+
e2Residuals = starE2 - psfE2
570+
571+
# Comp TEx
572+
ra = psf_cat["coord_ra"].to(units.deg)
573+
dec = psf_cat["coord_dec"].to(units.deg)
574+
575+
TExTreecorrConfig = {
576+
"TE1": self.config.psfTE1TreecorrConfig,
577+
"TE2": self.config.psfTE2TreecorrConfig,
578+
}
579+
580+
gatherE12Stat = {
581+
"TE1": {"E1": np.nan, "E2": np.nan, "Ex": np.nan, },
582+
"TE2": {"E1": np.nan, "E2": np.nan, "Ex": np.nan, },
583+
}
584+
585+
isNotNan = np.array([True] * len(ra))
586+
isNotNan &= np.isfinite(ra)
587+
isNotNan &= np.isfinite(dec)
588+
isNotNan &= np.isfinite(e1Residuals)
589+
isNotNan &= np.isfinite(e2Residuals)
590+
591+
if np.sum(isNotNan) >= 2:
592+
# TE1 and TE2 computation, over visit.
593+
for TEX in ["TE1", "TE2"]:
594+
595+
task = ComputeExPsfTask(TExTreecorrConfig[TEX])
596+
output = task.run(
597+
e1Residuals[isNotNan], e2Residuals[isNotNan],
598+
ra[isNotNan], dec[isNotNan],
599+
units="degree",
600+
)
601+
602+
gatherE12Stat[TEX]["E1"] = output.metric_E1
603+
gatherE12Stat[TEX]["E2"] = output.metric_E2
604+
gatherE12Stat[TEX]["Ex"] = output.metric_Ex
605+
606+
# TE3 and TE4 loop over detector and then median on visit.
607+
608+
TExTreecorrConfig = {
609+
"TE3": self.config.psfTE3TreecorrConfig,
610+
"TE4": self.config.psfTE4TreecorrConfig,
611+
}
612+
613+
gatherE34Stat = {
614+
"TE3": {"E1": [], "E2": [], "Ex": [], },
615+
"TE4": {"E1": [], "E2": [], "Ex": [], },
616+
}
617+
# calibrateImage run at detector level,
618+
# need to wait second run of PSF to run this.
619+
if "detector" in psf_cat.colnames:
620+
detectorIds = list(set(psf_cat["detector"]))
621+
for TEX in ["TE3", "TE4"]:
622+
for ccdId in detectorIds:
623+
isccdId = (ccdId == psf_cat["detector"])
624+
mask = (isccdId & isNotNan)
625+
if np.sum(mask) >= 2:
626+
task = ComputeExPsfTask(TExTreecorrConfig[TEX])
627+
output = task.run(
628+
e1Residuals[mask], e2Residuals[mask],
629+
ra[mask], dec[mask],
630+
units="degree",
631+
)
632+
gatherE34Stat[TEX]["E1"].append(output.metric_E1)
633+
gatherE34Stat[TEX]["E2"].append(output.metric_E2)
634+
gatherE34Stat[TEX]["Ex"].append(output.metric_Ex)
635+
636+
summary.psfTE1e1 = gatherE12Stat["TE1"]["E1"]
637+
summary.psfTE1e2 = gatherE12Stat["TE1"]["E2"]
638+
summary.psfTE1ex = gatherE12Stat["TE1"]["Ex"]
639+
summary.psfTE2e1 = gatherE12Stat["TE2"]["E1"]
640+
summary.psfTE2e2 = gatherE12Stat["TE2"]["E2"]
641+
summary.psfTE2ex = gatherE12Stat["TE2"]["Ex"]
642+
summary.psfTE3e1 = np.nanmedian(gatherE34Stat["TE3"]["E1"])
643+
summary.psfTE3e2 = np.nanmedian(gatherE34Stat["TE3"]["E2"])
644+
summary.psfTE3ex = np.nanmedian(gatherE34Stat["TE3"]["Ex"])
645+
summary.psfTE4e1 = np.nanmedian(gatherE34Stat["TE4"]["E1"])
646+
summary.psfTE4e2 = np.nanmedian(gatherE34Stat["TE4"]["E2"])
647+
summary.psfTE4ex = np.nanmedian(gatherE34Stat["TE4"]["Ex"])
648+
649+
self.psfTE1e1 = summary.psfTE1e1
650+
self.psfTE1e2 = summary.psfTE1e2
651+
self.psfTE1ex = summary.psfTE1ex
652+
self.psfTE2e1 = summary.psfTE2e1
653+
self.psfTE2e2 = summary.psfTE2e2
654+
self.psfTE2ex = summary.psfTE2ex
655+
self.psfTE3e1 = summary.psfTE3e1
656+
self.psfTE3e2 = summary.psfTE3e2
657+
self.psfTE3ex = summary.psfTE3ex
658+
self.psfTE4e1 = summary.psfTE4e1
659+
self.psfTE4e2 = summary.psfTE4e2
660+
self.psfTE4ex = summary.psfTE4ex
661+
437662
def update_wcs_stats(self, summary, wcs, bbox, visitInfo):
438663
"""Compute all summary-statistic fields that depend on the WCS model.
439664

‎python/lsst/pipe/tasks/postprocess.py‎

Lines changed: 12 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1568,6 +1568,18 @@ def run(self, visitSummaries):
15681568
visitEntry["obsStart"] = visitEntry["expMidpt"] - 0.5 * np.timedelta64(int(expTime * 1E9), "ns")
15691569
expTime_days = expTime / (60*60*24)
15701570
visitEntry["obsStartMJD"] = visitEntry["expMidptMJD"] - 0.5 * expTime_days
1571+
visitEntry["psfTE1e1"] = visitRow["psfTE1e1"]
1572+
visitEntry["psfTE1e2"] = visitRow["psfTE1e2"]
1573+
visitEntry["psfTE1ex"] = visitRow["psfTE1ex"]
1574+
visitEntry["psfTE2e1"] = visitRow["psfTE2e1"]
1575+
visitEntry["psfTE2e2"] = visitRow["psfTE2e2"]
1576+
visitEntry["psfTE2ex"] = visitRow["psfTE2ex"]
1577+
visitEntry["psfTE3e1"] = visitRow["psfTE3e1"]
1578+
visitEntry["psfTE3e2"] = visitRow["psfTE3e2"]
1579+
visitEntry["psfTE3ex"] = visitRow["psfTE3ex"]
1580+
visitEntry["psfTE4e1"] = visitRow["psfTE4e1"]
1581+
visitEntry["psfTE4e2"] = visitRow["psfTE4e2"]
1582+
visitEntry["psfTE4ex"] = visitRow["psfTE4ex"]
15711583
visitEntries.append(visitEntry)
15721584

15731585
# TODO: DM-30623, Add programId, exposureType, cameraTemp,

0 commit comments

Comments
 (0)