Skip to content

Commit 5ca5aae

Browse files
PFLegetPFLeget
committed
Compute TEX metric.
1 parent f28766c commit 5ca5aae

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):
@@ -410,10 +460,13 @@ def update_psf_stats(
410460
psfE1 = (psfXX - psfYY)/(psfXX + psfYY)
411461
psfE2 = 2*psfXY/(psfXX + psfYY)
412462

413-
psfStarDeltaE1Median = np.median(starE1 - psfE1)
414-
psfStarDeltaE1Scatter = sigmaMad(starE1 - psfE1, scale='normal')
415-
psfStarDeltaE2Median = np.median(starE2 - psfE2)
416-
psfStarDeltaE2Scatter = sigmaMad(starE2 - psfE2, scale='normal')
463+
e1Residuals = starE1 - psfE1
464+
e2Residuals = starE2 - psfE2
465+
466+
psfStarDeltaE1Median = np.median(e1Residuals)
467+
psfStarDeltaE1Scatter = sigmaMad(e1Residuals, scale='normal')
468+
psfStarDeltaE2Median = np.median(e2Residuals)
469+
psfStarDeltaE2Scatter = sigmaMad(e2Residuals, scale='normal')
417470

418471
psfStarDeltaSizeMedian = np.median(starSize - psfSize)
419472
psfStarDeltaSizeScatter = sigmaMad(starSize - psfSize, scale='normal')
@@ -436,6 +489,178 @@ def update_psf_stats(
436489
)
437490
summary.maxDistToNearestPsf = float(maxDistToNearestPsf)
438491

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

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

Lines changed: 12 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1560,6 +1560,18 @@ def run(self, visitSummaries):
15601560
visitEntry["obsStart"] = visitEntry["expMidpt"] - 0.5 * np.timedelta64(int(expTime * 1E9), "ns")
15611561
expTime_days = expTime / (60*60*24)
15621562
visitEntry["obsStartMJD"] = visitEntry["expMidptMJD"] - 0.5 * expTime_days
1563+
visitEntry["psfTE1e1"] = visitRow["psfTE1e1"]
1564+
visitEntry["psfTE1e2"] = visitRow["psfTE1e2"]
1565+
visitEntry["psfTE1ex"] = visitRow["psfTE1ex"]
1566+
visitEntry["psfTE2e1"] = visitRow["psfTE2e1"]
1567+
visitEntry["psfTE2e2"] = visitRow["psfTE2e2"]
1568+
visitEntry["psfTE2ex"] = visitRow["psfTE2ex"]
1569+
visitEntry["psfTE3e1"] = visitRow["psfTE3e1"]
1570+
visitEntry["psfTE3e2"] = visitRow["psfTE3e2"]
1571+
visitEntry["psfTE3ex"] = visitRow["psfTE3ex"]
1572+
visitEntry["psfTE4e1"] = visitRow["psfTE4e1"]
1573+
visitEntry["psfTE4e2"] = visitRow["psfTE4e2"]
1574+
visitEntry["psfTE4ex"] = visitRow["psfTE4ex"]
15631575
visitEntries.append(visitEntry)
15641576

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

0 commit comments

Comments
 (0)