Coverage for python / lsst / source / injection / utils / test_utils.py: 25%

1# This file is part of source_injection. 

2# 

3# Developed for the LSST Data Management System. 

4# This product includes software developed by the LSST Project 

5# (https://www.lsst.org). 

6# See the COPYRIGHT file at the top-level directory of this distribution 

7# for details of code ownership. 

8# 

9# This program is free software: you can redistribute it and/or modify 

10# it under the terms of the GNU General Public License as published by 

11# the Free Software Foundation, either version 3 of the License, or 

12# (at your option) any later version. 

13# 

14# This program is distributed in the hope that it will be useful, 

15# but WITHOUT ANY WARRANTY; without even the implied warranty of 

16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 

17# GNU General Public License for more details. 

18# 

19# You should have received a copy of the GNU General Public License 

20# along with this program. If not, see <https://www.gnu.org/licenses/>. 

21 

22 

23import numpy as np 

24 

25from lsst.afw.geom import makeCdMatrix, makeSkyWcs 

26from lsst.afw.image import makePhotoCalibFromCalibZeroPoint 

27from lsst.geom import Box2I, Extent2I, Point2D, Point2I, SpherePoint, degrees 

28from lsst.ip.isr.isrTask import IsrTask 

29from lsst.meas.algorithms.testUtils import plantSources 

30from lsst.pipe.base import Pipeline 

31from lsst.pipe.base.pipelineIR import ContractIR, LabeledSubset 

32from lsst.pipe.tasks.calibrate import CalibrateTask 

33from lsst.pipe.tasks.characterizeImage import CharacterizeImageTask 

34from lsst.source.injection import generate_injection_catalog 

35 

36 

37def make_test_exposure(): 

38 """Make a test exposure with a PSF attached and stars placed randomly. 

39 

40 This function generates a noisy synthetic image with Gaussian PSFs injected 

41 into the frame. 

42 The exposure is returned with a WCS, PhotoCalib and PSF attached. 

43 

44 Returns 

45 ------- 

46 exposure : `lsst.afw.image.Exposure` 

47 Exposure with calibs attached and stars placed randomly. 

48 """ 

49 # Inspired by meas_algorithms test_dynamicDetection.py. 

50 xy0 = Point2I(12345, 67890) # xy0 for image 

51 dims = Extent2I(2345, 2345) # Dimensions of image 

52 bbox = Box2I(xy0, dims) # Bounding box of image 

53 sigma = 3.21 # PSF sigma 

54 buffer = 4.0 # Buffer for star centers around edge 

55 n_sigma = 5.0 # Number of PSF sigmas for kernel 

56 sky = 12345.6 # Initial sky level 

57 num_stars = 100 # Number of stars 

58 noise = np.sqrt(sky) * np.pi * sigma**2 # Poisson noise per PSF 

59 faint = 1.0 * noise # Faintest level for star fluxes 

60 bright = 100.0 * noise # Brightest level for star fluxes 

61 star_bbox = Box2I(bbox) # Area on image in which we can put star centers 

62 star_bbox.grow(-int(buffer * sigma)) # Shrink star_bbox 

63 pixel_scale = 1.0e-5 * degrees # Pixel scale (1E-5 deg = 0.036 arcsec) 

64 

65 # Make an exposure with a PSF attached; place stars randomly. 

66 rng = np.random.default_rng(12345) 

67 stars = [ 

68 (xx, yy, ff, sigma) 

69 for xx, yy, ff in zip( 

70 rng.uniform(star_bbox.getMinX(), star_bbox.getMaxX(), num_stars), 

71 rng.uniform(star_bbox.getMinY(), star_bbox.getMaxY(), num_stars), 

72 np.linspace(faint, bright, num_stars), 

73 ) 

74 ] 

75 exposure = plantSources(bbox, 2 * int(n_sigma * sigma) + 1, sky, stars, True) 

76 

77 # Set WCS and PhotoCalib. 

78 exposure.setWcs( 

79 makeSkyWcs( 

80 crpix=Point2D(0, 0), 

81 crval=SpherePoint(0, 0, degrees), 

82 cdMatrix=makeCdMatrix(scale=pixel_scale), 

83 ) 

84 ) 

85 exposure.setPhotoCalib(makePhotoCalibFromCalibZeroPoint(1e10, 1e8)) 

86 return exposure 

87 

88 

89def make_test_injection_catalog(wcs, bbox): 

90 """Make a test source injection catalog. 

91 

92 This function generates a test source injection catalog consisting of 30 

93 star-like sources of varying magnitude. 

94 

95 Parameters 

96 ---------- 

97 wcs : `lsst.afw.geom.SkyWcs` 

98 WCS associated with the exposure. 

99 bbox : `lsst.geom.Box2I` 

100 Bounding box of the exposure. 

101 

102 Returns 

103 ------- 

104 injection_catalog : `astropy.table.Table` 

105 Source injection catalog. 

106 """ 

107 radec0 = wcs.pixelToSky(bbox.getBeginX(), bbox.getBeginY()) 

108 radec1 = wcs.pixelToSky(bbox.getEndX(), bbox.getEndY()) 

109 ra_lim = sorted([radec0.getRa().asDegrees(), radec1.getRa().asDegrees()]) 

110 dec_lim = sorted([radec0.getDec().asDegrees(), radec1.getDec().asDegrees()]) 

111 injection_catalog = generate_injection_catalog( 

112 ra_lim=ra_lim, 

113 dec_lim=dec_lim, 

114 wcs=wcs, 

115 number=10, 

116 source_type="DeltaFunction", 

117 mag=[10.0, 15.0, 20.0], 

118 ) 

119 return injection_catalog 

120 

121 

122def make_test_reference_pipeline(): 

123 """Make a test reference pipeline containing initial single-frame tasks.""" 

124 reference_pipeline = Pipeline("reference_pipeline") 

125 reference_pipeline.addTask(IsrTask, "isr") 

126 reference_pipeline.addTask(CharacterizeImageTask, "characterizeImage") 

127 reference_pipeline.addTask(CalibrateTask, "calibrate") 

128 reference_pipeline._pipelineIR.labeled_subsets["test_subset"] = LabeledSubset("test_subset", set(), None) 

129 reference_pipeline.addLabelToSubset("test_subset", "isr") 

130 reference_pipeline.addLabelToSubset("test_subset", "characterizeImage") 

131 reference_pipeline.addLabelToSubset("test_subset", "calibrate") 

132 isr_out = "isr.connections.ConnectionsClass(config=isr).outputExposure.name" 

133 char_inp = "characterizeImage.connections.ConnectionsClass(config=characterizeImage).exposure.name" 

134 char_out = "characterizeImage.connections.ConnectionsClass(config=characterizeImage).characterized.name" 

135 calib_inp = "calibrate.connections.ConnectionsClass(config=calibrate).exposure.name" 

136 # When injecting into a postISRCCD, contract1 should be violated. 

137 # The make_injection_pipeline utility should be robust to this. 

138 contract1 = ContractIR( 

139 contract=f"{isr_out} == {char_inp}", 

140 msg="isr output == characterizeImage input", 

141 ) 

142 contract2 = ContractIR( 

143 contract=f"{char_out} == {calib_inp}", 

144 msg="characterizeImage output == calibrate input", 

145 ) 

146 reference_pipeline._pipelineIR.contracts = [contract1, contract2] 

147 return reference_pipeline