Coverage for python / lsst / obs / subaru / arcturus.py: 0%

1"""Process the scattered light data from Arcturus""" 

2 

3import numpy as np 

4import multiprocessing 

5try: 

6 pyplot 

7except NameError: 

8 import matplotlib.pyplot as pyplot 

9 pyplot.interactive(1) 

10 

11import lsst.afw.display as afwDisplay 

12import lsst.analysis.utils as utils 

13 

14position = { # position of Arcturus on 16x16 binned "showCamera" images 

15 906142: (8, 8), 

16 906144: (8, 116), 

17 906146: (6, 228), 

18 906148: (5, 341), 

19 906150: (5, 455), 

20 906152: (8, 566), 

21 906154: (9, 682), 

22 906156: (8, 798), 

23 906158: (8, 917), 

24 906160: (8, 1039), 

25 906162: (8, 1038 + 108), 

26 906164: (8, 1038 + 2*108), 

27 906186: (8, 1038 + 3*108), 

28 906188: (8, 1038 + 4*108), 

29 906190: (8, 1038 + 5*108), 

30} 

31 

32 

33def showFrames(mos, frame0=1, R=23, subtractSky=True): 

34 visits = sorted(position.keys()) 

35 

36 frame = frame0 - 1 

37 for v in visits: 

38 frame += 1 

39 

40 xc, yc = position[v] 

41 xc -= mos[v].getX0() 

42 yc -= mos[v].getY0() 

43 

44 im = mos[v] 

45 if subtractSky: 

46 im = im.clone() 

47 im[2121:2230, 590:830] = np.nan # QE for this chip is bad 

48 

49 ima = im.getArray() 

50 im[:] -= np.percentile(ima, 25) 

51 

52 disp = afwDisplay.Display(frame=frame) 

53 disp.mtv(im, title=v) 

54 disp.dot("o", xc, yc, size=R, 

55 ctype=afwDisplay.GREEN if yc < mos[v].getHeight() else afwDisplay.RED) 

56 

57 

58class MedianFilterImageWorker(object): 

59 

60 def __init__(self, mos, medianN): 

61 self.mos = mos 

62 self.medianN = medianN 

63 

64 def __call__(self, v): 

65 print("Median filtering visit %d" % v) 

66 return v, utils.medianFilterImage(self.mos[-v], self.medianN) 

67 

68 

69def prepareFrames(mos, frame0=1, R=100, medianN=23, onlyVisits=[], nJob=1, force=False): 

70 """Prepare frames to have their radial profiles measured. 

71Subtract the first quartile 

72Set a radius R circle about the centre of Arcturus to NaN 

73Set a chip with a bad QE value to NaN 

74Median filter with a medianN x medianN filter 

75 

76The results are set as mos[-visit] 

77 

78If onlyVisits is specified, only process those chips [n.b. frame0 is still 

79obeyed]; otherwise process every visit in the positions dict 

80 """ 

81 visits = sorted(position.keys()) 

82 

83 visits0 = visits[:] # needed to keep frame numbers aligned 

84 if not force: 

85 visits = [v for v in visits if -v not in mos] # don't process ones we've already seen 

86 if onlyVisits: 

87 visits = [v for v in visits if v in onlyVisits] 

88 

89 width, height = mos[visits[0]].getDimensions() 

90 X, Y = np.meshgrid(np.arange(width), np.arange(height)) 

91 

92 for v in visits: 

93 im = mos[v].clone() 

94 mos[-v] = im 

95 im[2121:2230, 590:830] = np.nan # QE for this chip is bad 

96 

97 ima = im.getArray() 

98 im[:] -= np.percentile(ima, 25) 

99 

100 xc, yc = position[v] 

101 xc -= im.getX0() 

102 yc -= im.getY0() 

103 

104 ima[np.hypot(X - xc, Y - yc) < R] = np.nan 

105 # 

106 # multiprocessing 

107 # 

108 worker = MedianFilterImageWorker(mos, medianN) 

109 if nJob > 1: 

110 pool = multiprocessing.Pool(max([len(visits), nJob])) 

111 

112 # We use map_async(...).get(9999) instead of map(...) to workaround 

113 # a python but in handling ^C in subprocesses 

114 # (http://bugs.python.org/issue8296) 

115 for v, im in pool.map_async(worker, visits).get(9999): 

116 mos[-v] = im 

117 

118 pool.close() 

119 pool.join() 

120 else: 

121 for v in visits: 

122 mos[-v] = worker(v) 

123 # 

124 # Display 

125 # 

126 frame = frame0 - 1 

127 for v in visits0: 

128 frame += 1 

129 

130 if v not in visits: 

131 continue 

132 

133 im = mos[-v] 

134 

135 disp = afwDisplay.Display(frame=frame) 

136 if True: 

137 disp.mtv(im, title=v) 

138 else: 

139 disp.erase() 

140 

141 xc, yc = position[v] 

142 xc -= im.getX0() 

143 yc -= im.getY0() 

144 

145 disp.dot("o", xc, yc, size=R, ctype=afwDisplay.GREEN if yc < im.getHeight() else afwDisplay.RED) 

146 

147 

148def OLDprepareFrames(mos, frame0=1, R=23, medianN=23, onlyVisits=[], force=False): 

149 """Prepare frames to have their radial profiles measured. 

150 

151 #. Subtract the first quartile 

152 #. Set a radius R circle about the centre of Arcturus to NaN 

153 #. Set a chip with a bad QE value to NaN 

154 #. Median filter with a medianN x medianN filter 

155 

156 The result is set as mos[-visit] 

157 

158 If onlyVisits is specified, only process those chips [n.b. frame0 is still 

159 obeyed] 

160 """ 

161 

162 visits = sorted(position.keys()) 

163 

164 width, height = mos[visits[0]].getDimensions() 

165 

166 X, Y = np.meshgrid(np.arange(width), np.arange(height)) 

167 

168 frame = frame0 - 1 

169 for v in visits: 

170 frame += 1 

171 if onlyVisits and v not in onlyVisits: 

172 continue 

173 

174 print("Processing %d" % v) 

175 

176 im = mos[v].clone() 

177 im[2121:2230, 590:830] = np.nan # QE for this chip is bad 

178 

179 ima = im.getArray() 

180 im[:] -= np.percentile(ima, 25) 

181 

182 xc, yc = position[v] 

183 xc -= im.getX0() 

184 yc -= im.getY0() 

185 

186 ima[np.hypot(X - xc, Y - yc) < R] = np.nan 

187 

188 if force or -v not in mos: 

189 im = utils.medianFilterImage(im, medianN) 

190 mos[-v] = im 

191 else: 

192 im = mos[-v] 

193 

194 disp = afwDisplay.Display(frame=frame) 

195 if True: 

196 disp.mtv(im, title=v) 

197 else: 

198 disp.erase() 

199 disp.dot("o", xc, yc, size=R, ctype=afwDisplay.GREEN if yc < im.getHeight() else afwDisplay.RED) 

200 

201 

202def annularAverage(im, nBin=100, rmin=None, rmax=None, median=False): 

203 """Return image im's radial profile (more accurately the tuple 

204 (rbar, prof), binned into nBin bins. 

205 

206 r is measured from (0, 0) in the image, accounting properly for XY0 

207 

208 If rmax is provided it's the maximum value of r to consider 

209 """ 

210 

211 width, height = im.getDimensions() 

212 

213 if not rmax: 

214 rmax = 0.5*min(width, height) 

215 

216 bins = np.linspace(0, rmax, nBin) 

217 binWidth = bins[1] - bins[0] 

218 

219 rbar = np.empty_like(bins) 

220 prof = np.empty_like(bins) 

221 

222 X, Y = np.meshgrid(np.arange(width), np.arange(height)) 

223 R = np.hypot(X + im.getX0(), Y + im.getY0()).flatten() 

224 imProf = im.getArray().flatten() 

225 

226 if rmin: 

227 imProf = imProf[R > rmin] 

228 R = R[R > rmin] 

229 

230 for i in range(len(bins)): 

231 inBin = np.abs(R - bins[i]) < 0.5*binWidth 

232 rbar[i] = np.mean(R[inBin]) 

233 if np.sum(inBin): 

234 ivals = imProf[inBin] 

235 if median: 

236 val = np.median(ivals) 

237 else: 

238 val = np.mean(ivals[np.isfinite(ivals)]) 

239 else: 

240 val = 0 

241 

242 prof[i] = val 

243 

244 return rbar, prof 

245 

246 

247def plotAnnularAverages(mos, visits, bin=16, useSmoothed=True, weightByR=True, showSum=True, **kwargs): 

248 pyplot.clf() 

249 # 

250 # Find the spacing between exposures of Arcturus 

251 # 

252 allVisits = sorted(position.keys()) 

253 spacing = np.hypot(*position[allVisits[1]]) - np.hypot(*position[allVisits[0]]) 

254 

255 sumProf = None 

256 for v in visits: 

257 rbar, prof = annularAverage(mos[-v if useSmoothed else v], **kwargs) 

258 rbar *= bin 

259 prof *= spacing*bin # scale up as if we had exposures of Arcturus with 1 pixel steps 

260 # 

261 # Allow for the extra phase space further away from the boresight 

262 if weightByR: 

263 r = bin*(np.hypot(*position[v]) + 0.3*spacing) # really a geometric mean? 

264 prof *= 2*np.pi*r 

265 

266 print("%d %07.2g" % (v, np.sum(prof))) 

267 

268 if showSum: 

269 if sumProf is None: 

270 sumProf = np.zeros_like(prof) 

271 

272 sumProf += prof 

273 

274 pyplot.plot(rbar, prof, label=str(v)) 

275 

276 if sumProf is not None: 

277 pyplot.plot(rbar, sumProf, label="sum") 

278 

279 legend = pyplot.legend(loc="best", ncol=2, 

280 borderpad=0.1, labelspacing=0, handletextpad=0, handlelength=2, borderaxespad=0) 

281 

282 if legend: 

283 legend.draggable(True) 

284 

285 pyplot.show()