Coverage for metadetect / lsst / tests / test_lsst_skysub.py: 19%

1import sys 

2import os 

3import numpy as np 

4import pytest 

5import tqdm 

6import logging 

7import metadetect.lsst.skysub as lsst_skysub 

8import metadetect.lsst.measure as lsst_measure 

9from lsst.utils import getPackageDir 

10 

11try: 

12 getPackageDir('descwl_shear_sims') 

13 skip_tests_on_simulations = False 

14except LookupError: 

15 skip_tests_on_simulations = True 

16 

17logging.basicConfig( 

18 stream=sys.stdout, 

19 level=logging.WARN, 

20) 

21 

22 

23def make_lsst_sim(rng, gal_type, sky_n_sigma, star_density=0): 

24 import descwl_shear_sims 

25 

26 coadd_dim = 251 

27 

28 # the EDGE region is 5 pixels wide but, give a bit more space because the 

29 # sky sub seems to fail with a lsst.pex.exceptions.wrappers.LengthError, 

30 # presumably due to an object near the edge 

31 

32 buff = 20 

33 

34 if gal_type == 'fixed': 

35 galaxy_catalog = descwl_shear_sims.galaxies.FixedGalaxyCatalog( 

36 rng=rng, 

37 coadd_dim=coadd_dim, 

38 buff=buff, 

39 layout='random', 

40 mag=22, 

41 hlr=0.5, 

42 ) 

43 elif gal_type == 'wldeblend': 

44 galaxy_catalog = descwl_shear_sims.galaxies.WLDeblendGalaxyCatalog( 

45 rng=rng, 

46 coadd_dim=coadd_dim, 

47 buff=buff, 

48 ) 

49 else: 

50 raise ValueError('bad gal type: %s' % gal_type) 

51 

52 if star_density > 0: 

53 stars = descwl_shear_sims.stars.StarCatalog( 

54 rng=rng, 

55 coadd_dim=coadd_dim, 

56 buff=buff, 

57 density=star_density, 

58 ) 

59 else: 

60 stars = None 

61 

62 psf = descwl_shear_sims.psfs.make_fixed_psf(psf_type='gauss') 

63 

64 sim_data = descwl_shear_sims.make_sim( 

65 rng=rng, 

66 galaxy_catalog=galaxy_catalog, 

67 star_catalog=stars, 

68 coadd_dim=coadd_dim, 

69 g1=0.02, 

70 g2=0.00, 

71 psf=psf, 

72 se_dim=coadd_dim, 

73 sky_n_sigma=sky_n_sigma, 

74 ) 

75 return sim_data 

76 

77 

78def make_exp(dims): 

79 import lsst.afw.image as afw_image 

80 exp = afw_image.ExposureF(width=dims[1], height=dims[0]) 

81 return exp 

82 

83 

84def show_mask(exp): 

85 import lsst.afw.display as afw_display 

86 display = afw_display.getDisplay(backend='ds9') 

87 display.mtv(exp.mask) 

88 input('hit a key') 

89 

90 

91def show_image(exp): 

92 import lsst.afw.display as afw_display 

93 display = afw_display.getDisplay(backend='ds9') 

94 display.mtv(exp) 

95 display.scale('log', 'minmax') 

96 input('hit a key') 

97 

98 

99def check_skysub(meanvals, errvals, image_noise, true_sky): 

100 

101 meansky = meanvals.mean() 

102 stdsky = meanvals.std() 

103 errsky = stdsky/np.sqrt(meanvals.size) 

104 errmean = errvals.mean() 

105 

106 tol = image_noise / 10 

107 

108 # we want the uncertainty on the mean sky to be small enough for a 

109 # meaningful test 

110 assert errsky < tol / 3 

111 

112 # our sims have no sky in them 

113 

114 print('image_noise:', image_noise) 

115 print('tol:', tol) 

116 print('true sky value is', true_sky) 

117 print('mean of all trials: %g +/- %g' % (meansky, errsky)) 

118 print('sky error from trials:', stdsky) 

119 print('mean predicted error:', errmean) 

120 

121 assert np.abs(meansky - true_sky) < tol 

122 

123 

124def test_skysub_smoke(): 

125 seed = 5 

126 rng = np.random.RandomState(seed) 

127 

128 dims = [1000, 1000] 

129 noise = 1.0 

130 skyval = 0.0 

131 

132 exp = make_exp(dims) 

133 exp.image.array[:, :] = rng.normal(scale=noise, size=dims, loc=skyval) 

134 exp.variance.array[:, :] = noise**2 

135 

136 lsst_skysub.determine_and_subtract_sky(exp) 

137 

138 meta = exp.getMetadata() 

139 assert 'BGMEAN' in meta 

140 assert 'BGVAR' in meta 

141 

142 

143def test_skysub_pure_noise(): 

144 """ 

145 check the mean sky over all trials is within a fraction of the noise level 

146 """ 

147 seed = 5 

148 rng = np.random.RandomState(seed) 

149 

150 dims = [1000, 1000] 

151 noise = 1.0 

152 skyval = 0.0 

153 

154 ntrial = 100 

155 meanvals = np.zeros(ntrial) 

156 errvals = np.zeros(ntrial) 

157 

158 for itrial in range(ntrial): 

159 exp = make_exp(dims) 

160 exp.image.array[:, :] = rng.normal(scale=noise, size=dims, loc=skyval) 

161 exp.variance.array[:, :] = noise**2 

162 

163 lsst_skysub.determine_and_subtract_sky(exp) 

164 

165 meta = exp.getMetadata() 

166 

167 meanvals[itrial] = meta['BGMEAN'] 

168 errvals[itrial] = np.sqrt(meta['BGVAR']) 

169 

170 check_skysub(meanvals, errvals, noise, true_sky=0) 

171 

172 

173@pytest.mark.skipif( 

174 skip_tests_on_simulations, 

175 reason='descwl_shear_sims not available' 

176) 

177def test_skysub_sim_smoke(): 

178 seed = 812 

179 rng = np.random.RandomState(seed) 

180 sim = make_lsst_sim(rng, gal_type='fixed', sky_n_sigma=3) 

181 

182 exp = sim['band_data']['i'][0] 

183 lsst_skysub.determine_and_subtract_sky(exp) 

184 

185 meta = exp.getMetadata() 

186 assert 'BGMEAN' in meta 

187 assert 'BGVAR' in meta 

188 

189 

190@pytest.mark.skipif( 

191 skip_tests_on_simulations, 

192 reason='descwl_shear_sims not available' 

193) 

194@pytest.mark.parametrize('sky_n_sigma', [-0.5, -2.0, -100.0]) 

195def test_skysub_sim_fixed_gal(sky_n_sigma): 

196 """ 

197 check the measured mean sky over all trials is within 1/10 of the noise 

198 level 

199 """ 

200 seed = 184 

201 rng = np.random.RandomState(seed) 

202 

203 ntrial = 20 

204 meanvals = np.zeros(ntrial) 

205 errvals = np.zeros(ntrial) 

206 

207 for itrial in range(ntrial): 

208 sim = make_lsst_sim(rng, gal_type='fixed', sky_n_sigma=sky_n_sigma) 

209 

210 exp = sim['band_data']['i'][0] 

211 

212 noise = np.sqrt(np.median(exp.variance.array)) 

213 true_sky = sky_n_sigma * noise 

214 

215 if False: 

216 show_image(exp) 

217 

218 lsst_skysub.iterate_detection_and_skysub( 

219 exposure=exp, thresh=5, 

220 ) 

221 meta = exp.getMetadata() 

222 if 'BGMEAN' not in meta: 

223 raise RuntimeError('sky sub failed') 

224 

225 sky_meas = meta['BGMEAN'] 

226 

227 meanvals[itrial] = sky_meas 

228 errvals[itrial] = np.sqrt(meta['BGVAR']) 

229 

230 image_noise = np.median(exp.variance.array) 

231 check_skysub(meanvals, errvals, image_noise, true_sky=true_sky) 

232 

233 

234@pytest.mark.skipif( 

235 skip_tests_on_simulations, 

236 reason='descwl_shear_sims not available' 

237) 

238@pytest.mark.skipif( 

239 "CATSIM_DIR" not in os.environ, 

240 reason='simulation input data is not present', 

241) 

242@pytest.mark.parametrize('star_density', [20.0]) 

243@pytest.mark.parametrize('sky_n_sigma', [-0.5, -2.0, -100.0]) 

244def test_skysub_sim_wldeblend_gal(star_density, sky_n_sigma): 

245 """ 

246 check the measured mean sky over all trials is within 1/10 of the noise 

247 level 

248 

249 putting a nominal test at stellar density of 20 but really need to do a 

250 full shear recover test to explore this better 

251 """ 

252 seed = 312 

253 rng = np.random.RandomState(seed) 

254 

255 ntrial = 20 

256 meanvals = np.zeros(ntrial) 

257 errvals = np.zeros(ntrial) 

258 

259 for itrial in tqdm.trange(ntrial): 

260 sim = make_lsst_sim( 

261 rng, gal_type='wldeblend', star_density=star_density, 

262 sky_n_sigma=sky_n_sigma, 

263 ) 

264 

265 exp = sim['band_data']['i'][0] 

266 

267 noise = np.sqrt(np.median(exp.variance.array)) 

268 true_sky = sky_n_sigma * noise 

269 

270 if False: 

271 show_image(exp) 

272 

273 if True: 

274 lsst_skysub.iterate_detection_and_skysub( 

275 exposure=exp, thresh=5, 

276 ) 

277 meta = exp.getMetadata() 

278 sky_meas = meta['BGMEAN'] 

279 else: 

280 # this one is for debugging; we do the iterations ourselves so we 

281 # can display the result 

282 _, _ = lsst_measure.detect_and_deblend( 

283 exposure=exp, thresh=5, 

284 ) 

285 if False: 

286 show_mask(exp) 

287 

288 lsst_skysub.determine_and_subtract_sky(exp) 

289 sky_meas = exp.getMetadata()['BGMEAN'] 

290 

291 _, _ = lsst_measure.detect_and_deblend( 

292 exposure=exp, thresh=5, 

293 ) 

294 if False: 

295 show_mask(exp) 

296 

297 lsst_skysub.determine_and_subtract_sky(exp) 

298 

299 meta = exp.getMetadata() 

300 sky_meas += meta['BGMEAN'] 

301 

302 meanvals[itrial] = sky_meas 

303 errvals[itrial] = np.sqrt(meta['BGVAR']) 

304 

305 image_noise = np.median(exp.variance.array) 

306 

307 check_skysub(meanvals, errvals, image_noise, true_sky=true_sky)