Hide keyboard shortcuts

Hot-keys on this page

r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

1"""Annotations for sky maps.""" 

2import io 

3import os 

4import tempfile 

5 

6from astropy.io import fits 

7from astropy import table 

8from celery import group 

9from celery.exceptions import Ignore 

10from ligo.skymap.tool import ligo_skymap_flatten 

11from ligo.skymap.tool import ligo_skymap_from_samples 

12from ligo.skymap.tool import ligo_skymap_plot 

13from ligo.skymap.tool import ligo_skymap_plot_volume 

14from matplotlib import pyplot as plt 

15import numpy as np 

16 

17from . import gracedb 

18from . import lvalert 

19from ..import app 

20from ..jinja import env 

21from ..util.cmdline import handling_system_exit 

22from ..util.tempfile import NamedTemporaryFile 

23 

24 

25@app.task(ignore_result=True, shared=False) 

26def annotate_fits(filecontents, versioned_filename, graceid, tags): 

27 """Perform annotations on a sky map. 

28 

29 This function downloads a FITS file and then generates and uploads all 

30 derived images as well as an HTML dump of the FITS header. 

31 """ 

32 filebase = versioned_filename.partition('.fits')[0] 

33 header_msg = ( 

34 'FITS headers for <a href="/api/superevents/{graceid}/files/' 

35 '{versioned_filename}">{versioned_filename}</a>').format( 

36 graceid=graceid, versioned_filename=versioned_filename) 

37 allsky_msg = ( 

38 'Mollweide projection of <a href="/api/superevents/{graceid}/files/' 

39 '{versioned_filename}">{versioned_filename}</a>').format( 

40 graceid=graceid, versioned_filename=versioned_filename) 

41 volume_msg = ( 

42 'Volume rendering of <a href="/api/superevents/{graceid}/files/' 

43 '{versioned_filename}">{versioned_filename}</a>').format( 

44 graceid=graceid, versioned_filename=versioned_filename) 

45 

46 group( 

47 fits_header.s(versioned_filename) | 

48 gracedb.upload.s( 

49 filebase + '.html', graceid, header_msg, tags), 

50 

51 plot_allsky.s() | 

52 gracedb.upload.s( 

53 filebase + '.png', graceid, allsky_msg, tags), 

54 

55 annotate_fits_volume.s( 

56 filebase + '.volume.png', graceid, volume_msg, tags) 

57 ).delay(filecontents) 

58 

59 

60def is_3d_fits_file(filecontents): 

61 """Determine if a FITS file has distance information.""" 

62 with NamedTemporaryFile(content=filecontents) as fitsfile: 

63 return 'DISTNORM' in table.Table.read(fitsfile.name).colnames 

64 

65 

66@app.task(ignore_result=True, shared=False) 

67def annotate_fits_volume(filecontents, *args): 

68 """Perform annotations that are specific to 3D sky maps.""" 

69 if is_3d_fits_file(filecontents): 

70 ( 

71 plot_volume.s(filecontents) 

72 | 

73 gracedb.upload.s(*args) 

74 ).apply_async() 

75 

76 

77@app.task(shared=False) 

78def fits_header(filecontents, filename): 

79 """Dump FITS header to HTML.""" 

80 template = env.get_template('fits_header.jinja2') 

81 with NamedTemporaryFile(content=filecontents) as fitsfile, \ 

82 fits.open(fitsfile.name) as hdus: 

83 return template.render(filename=filename, hdus=hdus) 

84 

85 

86@app.task(shared=False) 

87def plot_allsky(filecontents, ra=None, dec=None): 

88 """Plot a Mollweide projection of a sky map using the command-line tool 

89 :doc:`ligo-skymap-plot <ligo/skymap/tool/ligo_skymap_plot>`. 

90 """ 

91 # Note: plt.style.context added as workaround for 

92 # https://github.com/astropy/astropy/issues/8004. 

93 with NamedTemporaryFile(mode='rb', suffix='.png') as pngfile, \ 

94 NamedTemporaryFile(content=filecontents) as fitsfile, \ 

95 plt.style.context({'text.usetex': False}, after_reset=True), \ 

96 handling_system_exit(): 

97 if ra is not None and dec is not None: 

98 ligo_skymap_plot.main([fitsfile.name, '-o', pngfile.name, 

99 '--annotate', '--radec', str(ra), str(dec)]) 

100 else: 

101 ligo_skymap_plot.main([fitsfile.name, '-o', pngfile.name, 

102 '--annotate', '--contour', '50', '90']) 

103 return pngfile.read() 

104 

105 

106@app.task(priority=1, queue='openmp', shared=False) 

107def plot_volume(filecontents): 

108 """Plot a 3D volume rendering of a sky map using the command-line tool 

109 :doc:`ligo-skymap-plot-volume <ligo/skymap/tool/ligo_skymap_plot_volume>`. 

110 """ 

111 # Note: plt.style.context added as workaround for 

112 # https://github.com/astropy/astropy/issues/8004. 

113 with NamedTemporaryFile(mode='rb', suffix='.png') as pngfile, \ 

114 NamedTemporaryFile(content=filecontents) as fitsfile, \ 

115 plt.style.context({'text.usetex': False}, after_reset=True), \ 

116 handling_system_exit(): 

117 ligo_skymap_plot_volume.main([fitsfile.name, '-o', 

118 pngfile.name, '--annotate']) 

119 return pngfile.read() 

120 

121 

122@app.task(shared=False) 

123def flatten(filecontents, filename): 

124 """Convert a HEALPix FITS file from multi-resolution UNIQ indexing to the 

125 more common IMPLICIT indexing using the command-line tool 

126 :doc:`ligo-skymap-flatten <ligo/skymap/tool/ligo_skymap_flatten>`. 

127 """ 

128 with NamedTemporaryFile(content=filecontents) as infile, \ 

129 tempfile.TemporaryDirectory() as tmpdir, \ 

130 handling_system_exit(): 

131 outfilename = os.path.join(tmpdir, filename) 

132 ligo_skymap_flatten.main([infile.name, outfilename]) 

133 return open(outfilename, 'rb').read() 

134 

135 

136@app.task(shared=False, queue='openmp') 

137def skymap_from_samples(samplefilecontents): 

138 """Generate multi-resolution fits file from samples.""" 

139 with NamedTemporaryFile(content=samplefilecontents) as samplefile, \ 

140 tempfile.TemporaryDirectory() as tmpdir, \ 

141 handling_system_exit(): 

142 ligo_skymap_from_samples.main( 

143 ['-j', '-o', tmpdir, samplefile.name]) 

144 with open(os.path.join(tmpdir, 'skymap.fits'), 'rb') as f: 

145 return f.read() 

146 

147 

148def plot_bayes_factor(logb, 

149 values=(1, 3, 5), 

150 labels=('', 'strong', 'very strong'), 

151 xlim=7, title=None, palette='RdYlBu'): 

152 """Visualize a Bayes factor as a `bullet graph`_. 

153 

154 Make a bar chart of a log Bayes factor as compared to a set of subjective 

155 threshold values. By default, use the thresholds from 

156 Kass & Raftery (1995). 

157 

158 .. _`bullet graph`: https://en.wikipedia.org/wiki/Bullet_graph 

159 

160 Parameters 

161 ---------- 

162 logb : float 

163 The natural logarithm of the Bayes factor. 

164 values : list 

165 A list of floating point values for human-friendly confidence levels. 

166 labels : list 

167 A list of string labels for human-friendly confidence levels. 

168 xlim : float 

169 Limits of plot (`-xlim` to `+xlim`). 

170 title : str 

171 Title for plot. 

172 palette : str 

173 Color palette. 

174 

175 Returns 

176 ------- 

177 fig : Matplotlib figure 

178 ax : Matplotlib axes 

179 

180 Example 

181 ------- 

182 .. plot:: 

183 :alt: Example Bayes factor plot 

184 

185 from gwcelery.tasks.skymaps import plot_bayes_factor 

186 plot_bayes_factor(6.3, title='GWCelery is awesome') 

187 

188 """ 

189 with plt.style.context('seaborn-notebook'): 

190 fig, ax = plt.subplots(figsize=(6, 1.7), tight_layout=True) 

191 ax.set_xlim(-xlim, xlim) 

192 ax.set_ylim(-0.5, 0.5) 

193 ax.set_yticks([]) 

194 ax.set_title(title) 

195 ax.set_ylabel(r'$\ln\,B$', rotation=0, 

196 rotation_mode='anchor', 

197 ha='right', va='center') 

198 

199 # Add human-friendly labels 

200 ticks = (*(-x for x in reversed(values)), 0, *values) 

201 ticklabels = ( 

202 *(f'{s}\nevidence\nagainst'.strip() for s in reversed(labels)), '', 

203 *(f'{s}\nevidence\nfor'.strip() for s in labels)) 

204 ax.set_xticks(ticks) 

205 ax.set_xticklabels(ticklabels) 

206 plt.setp(ax.get_xticklines(), visible=False) 

207 plt.setp(ax.get_xticklabels()[:len(ticks) // 2], ha='right') 

208 plt.setp(ax.get_xticklabels()[len(ticks) // 2:], ha='left') 

209 

210 # Plot colored bands for confidence thresholds 

211 fmt = plt.FuncFormatter(lambda x, _: f'{x:+g}'.replace('+0', '0')) 

212 ax2 = ax.twiny() 

213 ax2.set_xlim(*ax.get_xlim()) 

214 ax2.set_xticks(ticks) 

215 ax2.xaxis.set_major_formatter(fmt) 

216 levels = (-xlim, *ticks, xlim) 

217 colors = plt.get_cmap(palette)(np.arange(1, len(levels)) / len(levels)) 

218 ax.barh(0, np.diff(levels), 1, levels[:-1], 

219 linewidth=plt.rcParams['xtick.major.width'], 

220 color=colors, edgecolor='white') 

221 

222 # Plot bar for log Bayes factor value 

223 ax.barh(0, logb, 0.5, color='black', 

224 linewidth=plt.rcParams['xtick.major.width'], 

225 edgecolor='white') 

226 

227 for ax_ in fig.axes: 

228 ax_.grid(False) 

229 for spine in ax_.spines.values(): 

230 spine.set_visible(False) 

231 

232 return fig, ax 

233 

234 

235@app.task(shared=False) 

236def plot_coherence(filecontents): 

237 """LVAlert handler to plot the coherence Bayes factor. 

238 

239 Parameters 

240 ---------- 

241 contents : str, bytes 

242 The contents of the FITS file. 

243 

244 Returns 

245 ------- 

246 png : bytes 

247 The contents of a PNG file. 

248 

249 Notes 

250 ----- 

251 Under the hood, this just calls :meth:`plot_bayes_factor`. 

252 

253 """ 

254 with NamedTemporaryFile(content=filecontents) as fitsfile: 

255 header = fits.getheader(fitsfile, 1) 

256 try: 

257 logb = header['LOGBCI'] 

258 except KeyError: 

259 raise Ignore('FITS file does not have a LOGBCI field') 

260 

261 objid = header['OBJECT'] 

262 logb_string = np.format_float_positional(logb, 1, trim='0', sign=True) 

263 fig, _ = plot_bayes_factor( 

264 logb, title=rf'Coherence of {objid} $[\ln\,B = {logb_string}]$') 

265 outfile = io.BytesIO() 

266 fig.savefig(outfile, format='png') 

267 return outfile.getvalue() 

268 

269 

270@lvalert.handler('superevent', 

271 'mdc_superevent', 

272 shared=False) 

273def handle_plot_coherence(alert): 

274 """LVAlert handler to plot and upload a visualization of the coherence 

275 Bayes factor. 

276 

277 Notes 

278 ----- 

279 Under the hood, this just calls :meth:`plot_coherence`. 

280 

281 """ 

282 if alert['alert_type'] != 'log': 

283 return # not for us 

284 if not alert['data']['filename'].endswith('.fits'): 

285 return # not for us 

286 

287 graceid = alert['uid'] 

288 f = alert['data']['filename'] 

289 v = alert['data']['file_version'] 

290 fv = '{},{}'.format(f, v) 

291 

292 ( 

293 gracedb.download.s(fv, graceid) 

294 | 

295 plot_coherence.s() 

296 | 

297 gracedb.upload.s( 

298 f.replace('.fits', '.coherence.png'), graceid, 

299 message=( 

300 f'Bayes factor for coherence vs. incoherence of ' 

301 f'<a href="/api/superevents/{graceid}/files/{fv}">' 

302 f'{fv}</a>'), 

303 tags=['sky_loc'] 

304 ) 

305 ).delay()