/
multiBand.py
1447 lines (1179 loc) · 64.5 KB
/
multiBand.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
#!/usr/bin/env python
#
# LSST Data Management System
# Copyright 2008-2015 AURA/LSST.
#
# This product includes software developed by the
# LSST Project (http://www.lsst.org/).
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the LSST License Statement and
# the GNU General Public License along with this program. If not,
# see <https://www.lsstcorp.org/LegalNotices/>.
#
from __future__ import absolute_import, division, print_function
from builtins import zip
from builtins import range
import numpy
from lsst.coadd.utils.coaddDataIdContainer import ExistingCoaddDataIdContainer
from lsst.pipe.base import CmdLineTask, Struct, TaskRunner, ArgumentParser, ButlerInitializedTaskRunner
from lsst.pex.config import Config, Field, ListField, ConfigurableField, RangeField, ConfigField
from lsst.meas.algorithms import SourceDetectionTask
from lsst.meas.base import SingleFrameMeasurementTask, ApplyApCorrTask, CatalogCalculationTask
from lsst.meas.deblender import SourceDeblendTask
from lsst.pipe.tasks.coaddBase import getSkyInfo, scaleVariance
from lsst.meas.astrom import DirectMatchTask, denormalizeMatches
from lsst.pipe.tasks.fakes import BaseFakeSourcesTask
from lsst.pipe.tasks.setPrimaryFlags import SetPrimaryFlagsTask
from lsst.pipe.tasks.propagateVisitFlags import PropagateVisitFlagsTask
import lsst.afw.image as afwImage
import lsst.afw.table as afwTable
import lsst.afw.math as afwMath
import lsst.afw.geom as afwGeom
import lsst.afw.detection as afwDetect
from lsst.daf.base import PropertyList
"""
New dataset types:
* deepCoadd_det: detections from what used to be processCoadd (tract, patch, filter)
* deepCoadd_mergeDet: merged detections (tract, patch)
* deepCoadd_meas: measurements of merged detections (tract, patch, filter)
* deepCoadd_ref: reference sources (tract, patch)
All of these have associated *_schema catalogs that require no data ID and hold no records.
In addition, we have a schema-only dataset, which saves the schema for the PeakRecords in
the mergeDet, meas, and ref dataset Footprints:
* deepCoadd_peak_schema
"""
def _makeGetSchemaCatalogs(datasetSuffix):
"""Construct a getSchemaCatalogs instance method
These are identical for most of the classes here, so we'll consolidate
the code.
datasetSuffix: Suffix of dataset name, e.g., "src" for "deepCoadd_src"
"""
def getSchemaCatalogs(self):
"""Return a dict of empty catalogs for each catalog dataset produced by this task."""
src = afwTable.SourceCatalog(self.schema)
if hasattr(self, "algMetadata"):
src.getTable().setMetadata(self.algMetadata)
return {self.config.coaddName + "Coadd_" + datasetSuffix: src}
return getSchemaCatalogs
def _makeMakeIdFactory(datasetName):
"""Construct a makeIdFactory instance method
These are identical for all the classes here, so this consolidates
the code.
datasetName: Dataset name without the coadd name prefix, e.g., "CoaddId" for "deepCoaddId"
"""
def makeIdFactory(self, dataRef):
"""Return an IdFactory for setting the detection identifiers
The actual parameters used in the IdFactory are provided by
the butler (through the provided data reference.
"""
expBits = dataRef.get(self.config.coaddName + datasetName + "_bits")
expId = int(dataRef.get(self.config.coaddName + datasetName))
return afwTable.IdFactory.makeSource(expId, 64 - expBits)
return makeIdFactory
def getShortFilterName(name):
"""Given a longer, camera-specific filter name (e.g. "HSC-I") return its shorthand name ("i").
"""
# I'm not sure if this is the way this is supposed to be implemented, but it seems to work,
# and its the only way I could get it to work.
return afwImage.Filter(name).getFilterProperty().getName()
##############################################################################################################
class DetectCoaddSourcesConfig(Config):
"""!
\anchor DetectCoaddSourcesConfig_
\brief Configuration parameters for the DetectCoaddSourcesTask
"""
doScaleVariance = Field(dtype=bool, default=True, doc="Scale variance plane using empirical noise?")
detection = ConfigurableField(target=SourceDetectionTask, doc="Source detection")
coaddName = Field(dtype=str, default="deep", doc="Name of coadd")
mask = ListField(dtype=str, default=["DETECTED", "BAD", "SAT", "NO_DATA", "INTRP"],
doc="Mask planes for pixels to ignore when scaling variance")
doInsertFakes = Field(dtype=bool, default=False,
doc="Run fake sources injection task")
insertFakes = ConfigurableField(target=BaseFakeSourcesTask,
doc="Injection of fake sources for testing "
"purposes (must be retargeted)")
def setDefaults(self):
Config.setDefaults(self)
self.detection.thresholdType = "pixel_stdev"
self.detection.isotropicGrow = True
# Coadds are made from background-subtracted CCDs, so background subtraction should be very basic
self.detection.background.useApprox = False
self.detection.background.binSize = 4096
self.detection.background.undersampleStyle = 'REDUCE_INTERP_ORDER'
## \addtogroup LSST_task_documentation
## \{
## \page DetectCoaddSourcesTask
## \ref DetectCoaddSourcesTask_ "DetectCoaddSourcesTask"
## \copybrief DetectCoaddSourcesTask
## \}
class DetectCoaddSourcesTask(CmdLineTask):
"""!
\anchor DetectCoaddSourcesTask_
\brief Detect sources on a coadd
\section pipe_tasks_multiBand_Contents Contents
- \ref pipe_tasks_multiBand_DetectCoaddSourcesTask_Purpose
- \ref pipe_tasks_multiBand_DetectCoaddSourcesTask_Initialize
- \ref pipe_tasks_multiBand_DetectCoaddSourcesTask_Run
- \ref pipe_tasks_multiBand_DetectCoaddSourcesTask_Config
- \ref pipe_tasks_multiBand_DetectCoaddSourcesTask_Debug
- \ref pipe_tasks_multiband_DetectCoaddSourcesTask_Example
\section pipe_tasks_multiBand_DetectCoaddSourcesTask_Purpose Description
Command-line task that detects sources on a coadd of exposures obtained with a single filter.
Coadding individual visits requires each exposure to be warped. This introduces covariance in the noise
properties across pixels. Before detection, we correct the coadd variance by scaling the variance plane
in the coadd to match the observed variance. This is an approximate approach -- strictly, we should
propagate the full covariance matrix -- but it is simple and works well in practice.
After scaling the variance plane, we detect sources and generate footprints by delegating to the \ref
SourceDetectionTask_ "detection" subtask.
\par Inputs:
deepCoadd{tract,patch,filter}: ExposureF
\par Outputs:
deepCoadd_det{tract,patch,filter}: SourceCatalog (only parent Footprints)
\n deepCoadd_calexp{tract,patch,filter}: Variance scaled, background-subtracted input
exposure (ExposureF)
\n deepCoadd_calexp_background{tract,patch,filter}: BackgroundList
\par Data Unit:
tract, patch, filter
DetectCoaddSourcesTask delegates most of its work to the \ref SourceDetectionTask_ "detection" subtask.
You can retarget this subtask if you wish.
\section pipe_tasks_multiBand_DetectCoaddSourcesTask_Initialize Task initialization
\copydoc \_\_init\_\_
\section pipe_tasks_multiBand_DetectCoaddSourcesTask_Run Invoking the Task
\copydoc run
\section pipe_tasks_multiBand_DetectCoaddSourcesTask_Config Configuration parameters
See \ref DetectCoaddSourcesConfig_ "DetectSourcesConfig"
\section pipe_tasks_multiBand_DetectCoaddSourcesTask_Debug Debug variables
The \link lsst.pipe.base.cmdLineTask.CmdLineTask command line task\endlink interface supports a
flag \c -d to import \b debug.py from your \c PYTHONPATH; see \ref baseDebug for more about \b debug.py
files.
DetectCoaddSourcesTask has no debug variables of its own because it relegates all the work to
\ref SourceDetectionTask_ "SourceDetectionTask"; see the documetation for
\ref SourceDetectionTask_ "SourceDetectionTask" for further information.
\section pipe_tasks_multiband_DetectCoaddSourcesTask_Example A complete example of using DetectCoaddSourcesTask
DetectCoaddSourcesTask is meant to be run after assembling a coadded image in a given band. The purpose of
the task is to update the background, detect all sources in a single band and generate a set of parent
footprints. Subsequent tasks in the multi-band processing procedure will merge sources across bands and,
eventually, perform forced photometry. Command-line usage of DetectCoaddSourcesTask expects a data
reference to the coadd to be processed. A list of the available optional arguments can be obtained by
calling detectCoaddSources.py with the `--help` command line argument:
\code
detectCoaddSources.py --help
\endcode
To demonstrate usage of the DetectCoaddSourcesTask in the larger context of multi-band processing, we
will process HSC data in the [ci_hsc](https://github.com/lsst/ci_hsc) package. Assuming one has followed
steps 1 - 4 at \ref pipeTasks_multiBand, one may detect all the sources in each coadd as follows:
\code
detectCoaddSources.py $CI_HSC_DIR/DATA --id patch=5,4 tract=0 filter=HSC-I
\endcode
that will process the HSC-I band data. The results are written to
`$CI_HSC_DIR/DATA/deepCoadd-results/HSC-I`.
It is also necessary to run:
\code
detectCoaddSources.py $CI_HSC_DIR/DATA --id patch=5,4 tract=0 filter=HSC-R
\endcode
to generate the sources catalogs for the HSC-R band required by the next step in the multi-band
processing procedure: \ref MergeDetectionsTask_ "MergeDetectionsTask".
"""
_DefaultName = "detectCoaddSources"
ConfigClass = DetectCoaddSourcesConfig
getSchemaCatalogs = _makeGetSchemaCatalogs("det")
makeIdFactory = _makeMakeIdFactory("CoaddId")
@classmethod
def _makeArgumentParser(cls):
parser = ArgumentParser(name=cls._DefaultName)
parser.add_id_argument("--id", "deepCoadd", help="data ID, e.g. --id tract=12345 patch=1,2 filter=r",
ContainerClass=ExistingCoaddDataIdContainer)
return parser
def __init__(self, schema=None, **kwargs):
"""!
\brief Initialize the task. Create the \ref SourceDetectionTask_ "detection" subtask.
Keyword arguments (in addition to those forwarded to CmdLineTask.__init__):
\param[in] schema: initial schema for the output catalog, modified-in place to include all
fields set by this task. If None, the source minimal schema will be used.
\param[in] **kwargs: keyword arguments to be passed to lsst.pipe.base.task.Task.__init__
"""
CmdLineTask.__init__(self, **kwargs)
if schema is None:
schema = afwTable.SourceTable.makeMinimalSchema()
if self.config.doInsertFakes:
self.makeSubtask("insertFakes")
self.schema = schema
self.makeSubtask("detection", schema=self.schema)
def run(self, patchRef):
"""!
\brief Run detection on a coadd.
Invokes \ref runDetection and then uses \ref write to output the
results.
\param[in] patchRef: data reference for patch
"""
exposure = patchRef.get(self.config.coaddName + "Coadd", immediate=True)
results = self.runDetection(exposure, self.makeIdFactory(patchRef))
self.write(exposure, results, patchRef)
return results
def runDetection(self, exposure, idFactory):
"""!
\brief Run detection on an exposure.
First scale the variance plane to match the observed variance
using \ref scaleVariance. Then invoke the \ref SourceDetectionTask_ "detection" subtask to
detect sources.
\param[in] exposure: Exposure on which to detect
\param[in] idFactory: IdFactory to set source identifiers
\return a pipe.base.Struct with fields
- sources: catalog of detections
- backgrounds: list of backgrounds
"""
if self.config.doScaleVariance:
scaleVariance(exposure.getMaskedImage(), self.config.mask, log=self.log)
backgrounds = afwMath.BackgroundList()
if self.config.doInsertFakes:
self.insertFakes.run(exposure, background=backgrounds)
table = afwTable.SourceTable.make(self.schema, idFactory)
detections = self.detection.makeSourceCatalog(table, exposure)
sources = detections.sources
fpSets = detections.fpSets
if fpSets.background:
backgrounds.append(fpSets.background)
return Struct(sources=sources, backgrounds=backgrounds)
def write(self, exposure, results, patchRef):
"""!
\brief Write out results from runDetection.
\param[in] exposure: Exposure to write out
\param[in] results: Struct returned from runDetection
\param[in] patchRef: data reference for patch
"""
coaddName = self.config.coaddName + "Coadd"
patchRef.put(results.backgrounds, coaddName + "_calexp_background")
patchRef.put(results.sources, coaddName + "_det")
patchRef.put(exposure, coaddName + "_calexp")
##############################################################################################################
class MergeSourcesRunner(TaskRunner):
"""!
\anchor MergeSourcesRunner_
\brief Task runner for the \ref MergeSourcesTask_ "MergeSourcesTask". Required because the run method
requires a list of dataRefs rather than a single dataRef.
"""
def makeTask(self, parsedCmd=None, args=None):
"""!
\brief Provide a butler to the Task constructor.
\param[in] parsedCmd the parsed command
\param[in] args tuple of a list of data references and kwargs (un-used)
\throws RuntimeError if both parsedCmd & args are None
"""
if parsedCmd is not None:
butler = parsedCmd.butler
elif args is not None:
dataRefList, kwargs = args
butler = dataRefList[0].getButler()
else:
raise RuntimeError("Neither parsedCmd or args specified")
return self.TaskClass(config=self.config, log=self.log, butler=butler)
@staticmethod
def getTargetList(parsedCmd, **kwargs):
"""!
\brief Provide a list of patch references for each patch.
The patch references within the list will have different filters.
\param[in] parsedCmd the parsed command
\param **kwargs key word arguments (unused)
\throws RuntimeError if multiple references are provided for the same combination of tract, patch and
filter
"""
refList = {} # Will index this as refList[tract][patch][filter] = ref
for ref in parsedCmd.id.refList:
tract = ref.dataId["tract"]
patch = ref.dataId["patch"]
filter = ref.dataId["filter"]
if not tract in refList:
refList[tract] = {}
if not patch in refList[tract]:
refList[tract][patch] = {}
if filter in refList[tract][patch]:
raise RuntimeError("Multiple versions of %s" % (ref.dataId,))
refList[tract][patch][filter] = ref
return [(list(p.values()), kwargs) for t in refList.values() for p in t.values()]
class MergeSourcesConfig(Config):
"""!
\anchor MergeSourcesConfig_
\brief Configuration for merging sources.
"""
priorityList = ListField(dtype=str, default=[],
doc="Priority-ordered list of bands for the merge.")
coaddName = Field(dtype=str, default="deep", doc="Name of coadd")
def validate(self):
Config.validate(self)
if len(self.priorityList) == 0:
raise RuntimeError("No priority list provided")
class MergeSourcesTask(CmdLineTask):
"""!
\anchor MergeSourcesTask_
\brief A base class for merging source catalogs.
Merging detections (MergeDetectionsTask) and merging measurements (MergeMeasurementsTask) are
so similar that it makes sense to re-use the code, in the form of this abstract base class.
NB: Do not use this class directly. Instead use one of the child classes that inherit from
MergeSourcesTask such as \ref MergeDetectionsTask_ "MergeDetectionsTask" or \ref MergeMeasurementsTask_
"MergeMeasurementsTask"
Sub-classes should set the following class variables:
* `_DefaultName`: name of Task
* `inputDataset`: name of dataset to read
* `outputDataset`: name of dataset to write
* `getSchemaCatalogs` to the result of `_makeGetSchemaCatalogs(outputDataset)`
In addition, sub-classes must implement the mergeCatalogs method.
"""
_DefaultName = None
ConfigClass = MergeSourcesConfig
RunnerClass = MergeSourcesRunner
inputDataset = None
outputDataset = None
getSchemaCatalogs = None
@classmethod
def _makeArgumentParser(cls):
"""!
\brief Create a suitable ArgumentParser.
We will use the ArgumentParser to get a provide a list of data
references for patches; the RunnerClass will sort them into lists
of data references for the same patch
"""
parser = ArgumentParser(name=cls._DefaultName)
parser.add_id_argument("--id", "deepCoadd_" + cls.inputDataset,
ContainerClass=ExistingCoaddDataIdContainer,
help="data ID, e.g. --id tract=12345 patch=1,2 filter=g^r^i")
return parser
def getInputSchema(self, butler=None, schema=None):
"""!
\brief Obtain the input schema either directly or froma butler reference.
\param[in] butler butler reference to obtain the input schema from
\param[in] schema the input schema
"""
if schema is None:
assert butler is not None, "Neither butler nor schema specified"
schema = butler.get(self.config.coaddName + "Coadd_" + self.inputDataset + "_schema",
immediate=True).schema
return schema
def __init__(self, butler=None, schema=None, **kwargs):
"""!
\brief Initialize the task.
Keyword arguments (in addition to those forwarded to CmdLineTask.__init__):
\param[in] schema the schema of the detection catalogs used as input to this one
\param[in] butler a butler used to read the input schema from disk, if schema is None
Derived classes should use the getInputSchema() method to handle the additional
arguments and retreive the actual input schema.
"""
CmdLineTask.__init__(self, **kwargs)
def run(self, patchRefList):
"""!
\brief Merge coadd sources from multiple bands. Calls \ref mergeCatalogs which must be defined in
subclasses that inherit from MergeSourcesTask.
\param[in] patchRefList list of data references for each filter
"""
catalogs = dict(self.readCatalog(patchRef) for patchRef in patchRefList)
mergedCatalog = self.mergeCatalogs(catalogs, patchRefList[0])
self.write(patchRefList[0], mergedCatalog)
def readCatalog(self, patchRef):
"""!
\brief Read input catalog.
We read the input dataset provided by the 'inputDataset'
class variable.
\param[in] patchRef data reference for patch
\return tuple consisting of the filter name and the catalog
"""
filterName = patchRef.dataId["filter"]
catalog = patchRef.get(self.config.coaddName + "Coadd_" + self.inputDataset, immediate=True)
self.log.info("Read %d sources for filter %s: %s" % (len(catalog), filterName, patchRef.dataId))
return filterName, catalog
def mergeCatalogs(self, catalogs, patchRef):
"""!
\brief Merge multiple catalogs. This function must be defined in all subclasses that inherit from
MergeSourcesTask.
\param[in] catalogs dict mapping filter name to source catalog
\return merged catalog
"""
raise NotImplementedError()
def write(self, patchRef, catalog):
"""!
\brief Write the output.
\param[in] patchRef data reference for patch
\param[in] catalog catalog
We write as the dataset provided by the 'outputDataset'
class variable.
"""
patchRef.put(catalog, self.config.coaddName + "Coadd_" + self.outputDataset)
# since the filter isn't actually part of the data ID for the dataset we're saving,
# it's confusing to see it in the log message, even if the butler simply ignores it.
mergeDataId = patchRef.dataId.copy()
del mergeDataId["filter"]
self.log.info("Wrote merged catalog: %s" % (mergeDataId,))
def writeMetadata(self, dataRefList):
"""!
\brief No metadata to write, and not sure how to write it for a list of dataRefs.
"""
pass
class CullPeaksConfig(Config):
"""!
\anchor CullPeaksConfig_
\brief Configuration for culling garbage peaks after merging footprints.
Peaks may also be culled after detection or during deblending; this configuration object
only deals with culling after merging Footprints.
These cuts are based on three quantities:
- nBands: the number of bands in which the peak was detected
- peakRank: the position of the peak within its family, sorted from brightest to faintest.
- peakRankNormalized: the peak rank divided by the total number of peaks in the family.
The formula that identifie peaks to cull is:
nBands < nBandsSufficient
AND (rank >= rankSufficient)
AND (rank >= rankConsider OR rank >= rankNormalizedConsider)
To disable peak culling, simply set nBandsSafe=1.
"""
nBandsSufficient = RangeField(dtype=int, default=2, min=1,
doc="Always keep peaks detected in this many bands")
rankSufficient = RangeField(dtype=int, default=20, min=1,
doc="Always keep this many peaks in each family")
rankConsidered = RangeField(dtype=int, default=30, min=1,
doc=("Keep peaks with less than this rank that also match the "
"rankNormalizedConsidered condition."))
rankNormalizedConsidered = RangeField(dtype=float, default=0.7, min=0.0,
doc=("Keep peaks with less than this normalized rank that"
" also match the rankConsidered condition."))
class MergeDetectionsConfig(MergeSourcesConfig):
"""!
\anchor MergeDetectionsConfig_
\brief Configuration parameters for the MergeDetectionsTask.
"""
minNewPeak = Field(dtype=float, default=1,
doc="Minimum distance from closest peak to create a new one (in arcsec).")
maxSamePeak = Field(dtype=float, default=0.3,
doc="When adding new catalogs to the merge, all peaks less than this distance "
" (in arcsec) to an existing peak will be flagged as detected in that catalog.")
cullPeaks = ConfigField(dtype=CullPeaksConfig, doc="Configuration for how to cull peaks.")
skyFilterName = Field(dtype=str, default="sky",
doc="Name of `filter' used to label sky objects (e.g. flag merge_peak_sky is set)\n"
"(N.b. should be in MergeMeasurementsConfig.pseudoFilterList)")
skySourceRadius = Field(dtype=float, default=8,
doc="Radius, in pixels, of sky objects")
skyGrowDetectedFootprints = Field(dtype=int, default=0,
doc="Number of pixels to grow the detected footprint mask "
"when adding sky objects")
nSkySourcesPerPatch = Field(dtype=int, default=100,
doc="Try to add this many sky objects to the mergeDet list, which will\n"
"then be measured along with the detected objects in sourceMeasurementTask")
nTrialSkySourcesPerPatch = Field(dtype=int, default=None, optional=True,
doc="Maximum number of trial sky object positions\n"
"(default: nSkySourcesPerPatch*nTrialSkySourcesPerPatchMultiplier)")
nTrialSkySourcesPerPatchMultiplier = Field(dtype=int, default=5,
doc="Set nTrialSkySourcesPerPatch to\n"
" nSkySourcesPerPatch*nTrialSkySourcesPerPatchMultiplier\n"
"if nTrialSkySourcesPerPatch is None")
## \addtogroup LSST_task_documentation
## \{
## \page MergeDetectionsTask
## \ref MergeDetectionsTask_ "MergeDetectionsTask"
## \copybrief MergeDetectionsTask
## \}
class MergeDetectionsTask(MergeSourcesTask):
"""!
\anchor MergeDetectionsTask_
\brief Merge coadd detections from multiple bands.
\section pipe_tasks_multiBand_Contents Contents
- \ref pipe_tasks_multiBand_MergeDetectionsTask_Purpose
- \ref pipe_tasks_multiBand_MergeDetectionsTask_Init
- \ref pipe_tasks_multiBand_MergeDetectionsTask_Run
- \ref pipe_tasks_multiBand_MergeDetectionsTask_Config
- \ref pipe_tasks_multiBand_MergeDetectionsTask_Debug
- \ref pipe_tasks_multiband_MergeDetectionsTask_Example
\section pipe_tasks_multiBand_MergeDetectionsTask_Purpose Description
Command-line task that merges sources detected in coadds of exposures obtained with different filters.
To perform photometry consistently across coadds in multiple filter bands, we create a master catalog of
sources from all bands by merging the sources (peaks & footprints) detected in each coadd, while keeping
track of which band each source originates in.
The catalog merge is performed by \ref getMergedSourceCatalog. Spurious peaks detected around bright
objects are culled as described in \ref CullPeaksConfig_.
\par Inputs:
deepCoadd_det{tract,patch,filter}: SourceCatalog (only parent Footprints)
\par Outputs:
deepCoadd_mergeDet{tract,patch}: SourceCatalog (only parent Footprints)
\par Data Unit:
tract, patch
MergeDetectionsTask subclasses \ref MergeSourcesTask_ "MergeSourcesTask".
\section pipe_tasks_multiBand_MergeDetectionsTask_Init Task initialisation
\copydoc \_\_init\_\_
\section pipe_tasks_multiBand_MergeDetectionsTask_Run Invoking the Task
\copydoc run
\section pipe_tasks_multiBand_MergeDetectionsTask_Config Configuration parameters
See \ref MergeDetectionsConfig_
\section pipe_tasks_multiBand_MergeDetectionsTask_Debug Debug variables
The \link lsst.pipe.base.cmdLineTask.CmdLineTask command line task\endlink interface supports a flag \c -d
to import \b debug.py from your \c PYTHONPATH; see \ref baseDebug for more about \b debug.py files.
MergeDetectionsTask has no debug variables.
\section pipe_tasks_multiband_MergeDetectionsTask_Example A complete example of using MergeDetectionsTask
MergeDetectionsTask is meant to be run after detecting sources in coadds generated for the chosen subset
of the available bands.
The purpose of the task is to merge sources (peaks & footprints) detected in the coadds generated from the
chosen subset of filters.
Subsequent tasks in the multi-band processing procedure will deblend the generated master list of sources
and, eventually, perform forced photometry.
Command-line usage of MergeDetectionsTask expects data references for all the coadds to be processed.
A list of the available optional arguments can be obtained by calling mergeCoaddDetections.py with the
`--help` command line argument:
\code
mergeCoaddDetections.py --help
\endcode
To demonstrate usage of the DetectCoaddSourcesTask in the larger context of multi-band processing, we
will process HSC data in the [ci_hsc](https://github.com/lsst/ci_hsc) package. Assuming one has finished
step 5 at \ref pipeTasks_multiBand, one may merge the catalogs of sources from each coadd as follows:
\code
mergeCoaddDetections.py $CI_HSC_DIR/DATA --id patch=5,4 tract=0 filter=HSC-I^HSC-R
\endcode
This will merge the HSC-I & -R band parent source catalogs and write the results to
`$CI_HSC_DIR/DATA/deepCoadd-results/merged/0/5,4/mergeDet-0-5,4.fits`.
The next step in the multi-band processing procedure is
\ref MeasureMergedCoaddSourcesTask_ "MeasureMergedCoaddSourcesTask"
"""
ConfigClass = MergeDetectionsConfig
_DefaultName = "mergeCoaddDetections"
inputDataset = "det"
outputDataset = "mergeDet"
makeIdFactory = _makeMakeIdFactory("MergedCoaddId")
def __init__(self, butler=None, schema=None, **kwargs):
"""!
\brief Initialize the merge detections task.
A \ref FootprintMergeList_ "FootprintMergeList" will be used to
merge the source catalogs.
Additional keyword arguments (forwarded to MergeSourcesTask.__init__):
\param[in] schema the schema of the detection catalogs used as input to this one
\param[in] butler a butler used to read the input schema from disk, if schema is None
\param[in] **kwargs keyword arguments to be passed to MergeSourcesTask.__init__
The task will set its own self.schema attribute to the schema of the output merged catalog.
"""
MergeSourcesTask.__init__(self, butler=butler, schema=schema, **kwargs)
self.schema = self.getInputSchema(butler=butler, schema=schema)
filterNames = [getShortFilterName(name) for name in self.config.priorityList]
if self.config.nSkySourcesPerPatch > 0:
filterNames += [self.config.skyFilterName]
self.merged = afwDetect.FootprintMergeList(self.schema, filterNames)
def mergeCatalogs(self, catalogs, patchRef):
"""!
\brief Merge multiple catalogs.
After ordering the catalogs and filters in priority order,
\ref getMergedSourceCatalog of the \ref FootprintMergeList_ "FootprintMergeList" created by
\ref \_\_init\_\_ is used to perform the actual merging. Finally, \ref cullPeaks is used to remove
garbage peaks detected around bright objects.
\param[in] catalogs
\param[in] patchRef
\param[out] mergedList
"""
# Convert distance to tract coordinate
skyInfo = getSkyInfo(coaddName=self.config.coaddName, patchRef=patchRef)
tractWcs = skyInfo.wcs
peakDistance = self.config.minNewPeak / tractWcs.pixelScale().asArcseconds()
samePeakDistance = self.config.maxSamePeak / tractWcs.pixelScale().asArcseconds()
# Put catalogs, filters in priority order
orderedCatalogs = [catalogs[band] for band in self.config.priorityList if band in catalogs.keys()]
orderedBands = [getShortFilterName(band) for band in self.config.priorityList
if band in catalogs.keys()]
mergedList = self.merged.getMergedSourceCatalog(orderedCatalogs, orderedBands, peakDistance,
self.schema, self.makeIdFactory(patchRef),
samePeakDistance)
#
# Add extra sources that correspond to blank sky
#
skySourceFootprints = self.getSkySourceFootprints(
mergedList, skyInfo, self.config.skyGrowDetectedFootprints)
if skySourceFootprints:
key = mergedList.schema.find("merge_footprint_%s" % self.config.skyFilterName).key
for foot in skySourceFootprints:
s = mergedList.addNew()
s.setFootprint(foot)
s.set(key, True)
self.log.info("Added %d sky sources (%.0f%% of requested)",
len(skySourceFootprints),
100*len(skySourceFootprints)/float(self.config.nSkySourcesPerPatch))
# Sort Peaks from brightest to faintest
for record in mergedList:
record.getFootprint().sortPeaks()
self.log.info("Merged to %d sources" % len(mergedList))
# Attempt to remove garbage peaks
self.cullPeaks(mergedList)
return mergedList
def cullPeaks(self, catalog):
"""!
\brief Attempt to remove garbage peaks (mostly on the outskirts of large blends).
\param[in] catalog Source catalog
"""
keys = [item.key for item in self.merged.getPeakSchema().extract("merge.peak.*").values()]
totalPeaks = 0
culledPeaks = 0
for parentSource in catalog:
# Make a list copy so we can clear the attached PeakCatalog and append the ones we're keeping
# to it (which is easier than deleting as we iterate).
keptPeaks = parentSource.getFootprint().getPeaks()
oldPeaks = list(keptPeaks)
keptPeaks.clear()
familySize = len(oldPeaks)
totalPeaks += familySize
for rank, peak in enumerate(oldPeaks):
if ((rank < self.config.cullPeaks.rankSufficient) or
(self.config.cullPeaks.nBandsSufficient > 1 and
sum([peak.get(k) for k in keys]) >= self.config.cullPeaks.nBandsSufficient) or
(rank < self.config.cullPeaks.rankConsidered and
rank < self.config.cullPeaks.rankNormalizedConsidered * familySize)):
keptPeaks.append(peak)
else:
culledPeaks += 1
self.log.info("Culled %d of %d peaks" % (culledPeaks, totalPeaks))
def getSchemaCatalogs(self):
"""!
Return a dict of empty catalogs for each catalog dataset produced by this task.
\param[out] dictionary of empty catalogs
"""
mergeDet = afwTable.SourceCatalog(self.schema)
peak = afwDetect.PeakCatalog(self.merged.getPeakSchema())
return {self.config.coaddName + "Coadd_mergeDet": mergeDet,
self.config.coaddName + "Coadd_peak": peak}
def getSkySourceFootprints(self, mergedList, skyInfo, growDetectedFootprints=0):
"""!
\brief Return a list of Footprints of sky objects which don't overlap with anything in mergedList
\param mergedList The merged Footprints from all the input bands
\param skyInfo A description of the patch
\param growDetectedFootprints The number of pixels to grow the detected footprints
"""
if self.config.nSkySourcesPerPatch <= 0:
return []
skySourceRadius = self.config.skySourceRadius
nSkySourcesPerPatch = self.config.nSkySourcesPerPatch
nTrialSkySourcesPerPatch = self.config.nTrialSkySourcesPerPatch
if nTrialSkySourcesPerPatch is None:
nTrialSkySourcesPerPatch = self.config.nTrialSkySourcesPerPatchMultiplier*nSkySourcesPerPatch
#
# We are going to find circular Footprints that don't intersect any pre-existing Footprints,
# and the easiest way to do this is to generate a Mask containing all the detected pixels (as
# merged by this task).
#
patchBBox = skyInfo.patchInfo.getOuterBBox()
mask = afwImage.Mask(patchBBox)
detectedMask = mask.getPlaneBitMask("DETECTED")
for s in mergedList:
foot = s.getFootprint()
if growDetectedFootprints > 0:
foot.dilate(growDetectedFootprints)
foot.spans.setMask(mask, detectedMask)
xmin, ymin = patchBBox.getMin()
xmax, ymax = patchBBox.getMax()
# Avoid objects partially off the image
xmin += skySourceRadius + 1
ymin += skySourceRadius + 1
xmax -= skySourceRadius + 1
ymax -= skySourceRadius + 1
skySourceFootprints = []
maskToSpanSet = afwGeom.SpanSet.fromMask(mask)
for i in range(nTrialSkySourcesPerPatch):
if len(skySourceFootprints) == nSkySourcesPerPatch:
break
x = int(numpy.random.uniform(xmin, xmax))
y = int(numpy.random.uniform(ymin, ymax))
spans = afwGeom.SpanSet.fromShape(int(skySourceRadius),
offset=(x, y))
foot = afwDetect.Footprint(spans, patchBBox)
foot.setPeakSchema(self.merged.getPeakSchema())
if not foot.spans.overlaps(maskToSpanSet):
foot.addPeak(x, y, 0)
foot.getPeaks()[0].set("merge_peak_%s" % self.config.skyFilterName, True)
skySourceFootprints.append(foot)
return skySourceFootprints
class MeasureMergedCoaddSourcesConfig(Config):
"""!
\anchor MeasureMergedCoaddSourcesConfig_
\brief Configuration parameters for the MeasureMergedCoaddSourcesTask
"""
doDeblend = Field(dtype=bool, default=True, doc="Deblend sources?")
deblend = ConfigurableField(target=SourceDeblendTask, doc="Deblend sources")
measurement = ConfigurableField(target=SingleFrameMeasurementTask, doc="Source measurement")
setPrimaryFlags = ConfigurableField(target=SetPrimaryFlagsTask, doc="Set flags for primary tract/patch")
doPropagateFlags = Field(
dtype=bool, default=True,
doc="Whether to match sources to CCD catalogs to propagate flags (to e.g. identify PSF stars)"
)
propagateFlags = ConfigurableField(target=PropagateVisitFlagsTask, doc="Propagate visit flags to coadd")
doMatchSources = Field(dtype=bool, default=True, doc="Match sources to reference catalog?")
match = ConfigurableField(target=DirectMatchTask, doc="Matching to reference catalog")
doWriteMatchesDenormalized = Field(
dtype=bool,
default=False,
doc=("Write reference matches in denormalized format? "
"This format uses more disk space, but is more convenient to read."),
)
coaddName = Field(dtype=str, default="deep", doc="Name of coadd")
checkUnitsParseStrict = Field(
doc="Strictness of Astropy unit compatibility check, can be 'raise', 'warn' or 'silent'",
dtype=str,
default="raise",
)
doApCorr = Field(
dtype=bool,
default=True,
doc="Apply aperture corrections"
)
applyApCorr = ConfigurableField(
target=ApplyApCorrTask,
doc="Subtask to apply aperture corrections"
)
doRunCatalogCalculation = Field(
dtype=bool,
default=True,
doc='Run catalogCalculation task'
)
catalogCalculation = ConfigurableField(
target=CatalogCalculationTask,
doc="Subtask to run catalogCalculation plugins on catalog"
)
def setDefaults(self):
Config.setDefaults(self)
self.deblend.propagateAllPeaks = True
self.measurement.plugins.names |= ['base_InputCount']
# The following line must be set if clipped pixel flags are to be added to the output table
# The clipped mask plane is added by running SafeClipAssembleCoaddTask
self.measurement.plugins['base_PixelFlags'].masksFpAnywhere = ['CLIPPED']
## \addtogroup LSST_task_documentation
## \{
## \page MeasureMergedCoaddSourcesTask
## \ref MeasureMergedCoaddSourcesTask_ "MeasureMergedCoaddSourcesTask"
## \copybrief MeasureMergedCoaddSourcesTask
## \}
class MeasureMergedCoaddSourcesTask(CmdLineTask):
"""!
\anchor MeasureMergedCoaddSourcesTask_
\brief Deblend sources from master catalog in each coadd seperately and measure.
\section pipe_tasks_multiBand_Contents Contents
- \ref pipe_tasks_multiBand_MeasureMergedCoaddSourcesTask_Purpose
- \ref pipe_tasks_multiBand_MeasureMergedCoaddSourcesTask_Initialize
- \ref pipe_tasks_multiBand_MeasureMergedCoaddSourcesTask_Run
- \ref pipe_tasks_multiBand_MeasureMergedCoaddSourcesTask_Config
- \ref pipe_tasks_multiBand_MeasureMergedCoaddSourcesTask_Debug
- \ref pipe_tasks_multiband_MeasureMergedCoaddSourcesTask_Example
\section pipe_tasks_multiBand_MeasureMergedCoaddSourcesTask_Purpose Description
Command-line task that uses peaks and footprints from a master catalog to perform deblending and
measurement in each coadd.
Given a master input catalog of sources (peaks and footprints), deblend and measure each source on the
coadd. Repeating this procedure with the same master catalog across multiple coadds will generate a
consistent set of child sources.
The deblender retains all peaks and deblends any missing peaks (dropouts in that band) as PSFs. Source
properties are measured and the \c is-primary flag (indicating sources with no children) is set. Visit
flags are propagated to the coadd sources.
Optionally, we can match the coadd sources to an external reference catalog.
\par Inputs:
deepCoadd_mergeDet{tract,patch}: SourceCatalog
\n deepCoadd_calexp{tract,patch,filter}: ExposureF
\par Outputs:
deepCoadd_meas{tract,patch,filter}: SourceCatalog
\par Data Unit:
tract, patch, filter
MeasureMergedCoaddSourcesTask delegates most of its work to a set of sub-tasks:
<DL>
<DT> \ref SourceDeblendTask_ "deblend"
<DD> Deblend all the sources from the master catalog.</DD>
<DT> \ref SingleFrameMeasurementTask_ "measurement"
<DD> Measure source properties of deblended sources.</DD>
<DT> \ref SetPrimaryFlagsTask_ "setPrimaryFlags"
<DD> Set flag 'is-primary' as well as related flags on sources. 'is-primary' is set for sources that are
not at the edge of the field and that have either not been deblended or are the children of deblended
sources</DD>
<DT> \ref PropagateVisitFlagsTask_ "propagateFlags"
<DD> Propagate flags set in individual visits to the coadd.</DD>
<DT> \ref DirectMatchTask_ "match"
<DD> Match input sources to a reference catalog (optional).
</DD>
</DL>
These subtasks may be retargeted as required.
\section pipe_tasks_multiBand_MeasureMergedCoaddSourcesTask_Initialize Task initialization
\copydoc \_\_init\_\_
\section pipe_tasks_multiBand_MeasureMergedCoaddSourcesTask_Run Invoking the Task
\copydoc run
\section pipe_tasks_multiBand_MeasureMergedCoaddSourcesTask_Config Configuration parameters
See \ref MeasureMergedCoaddSourcesConfig_
\section pipe_tasks_multiBand_MeasureMergedCoaddSourcesTask_Debug Debug variables
The \link lsst.pipe.base.cmdLineTask.CmdLineTask command line task\endlink interface supports a
flag \c -d to import \b debug.py from your \c PYTHONPATH; see \ref baseDebug for more about \b debug.py
files.
MeasureMergedCoaddSourcesTask has no debug variables of its own because it delegates all the work to
the various sub-tasks. See the documetation for individual sub-tasks for more information.
\section pipe_tasks_multiband_MeasureMergedCoaddSourcesTask_Example A complete example of using MeasureMergedCoaddSourcesTask