Use multi-ordered sky maps in RAVEN pipeline; fixes #408

gwcelery/tasks/external_skymaps.py

@@ -10,6 +10,7 @@ from ligo.skymap.tool import ligo_skymap_combine
 import gcn
 import healpy as hp
 import lxml.etree
+from mhealpy import HealpixMap
 import re
 import ssl
 import urllib
@@ -32,21 +33,22 @@ def create_combined_skymap(se_id, ext_id):
     ext_skymap_filename = get_skymap_filename(ext_id)
     new_skymap_filename = re.findall(r'(.*).fits.gz', se_skymap_filename)[0]
 
-    #  FIXME: put download functions in canvas
-    se_skymap = gracedb.download(se_skymap_filename, se_id)
-    ext_skymap = gracedb.download(ext_skymap_filename, ext_id)
     message = 'Combined LVC-external sky map using {0} and {1}'.format(
         se_skymap_filename, ext_skymap_filename)
     message_png = (
         'Mollweide projection of <a href="/api/events/{graceid}/files/'
         '{filename}">{filename}</a>').format(
-            graceid=se_id, filename=new_skymap_filename + '-ext.fits.gz')
+            graceid=se_id, filename=new_skymap_filename + '-ext.multiorder.fits')
 
-    (
-        combine_skymaps.si(se_skymap, ext_skymap)
+    (   
+        _download_skymaps.si(
+            se_skymap_filename, ext_skymap_filename, se_id, ext_id
+        )
+        |
+        combine_skymaps.s()
         |
         group(
-            gracedb.upload.s(new_skymap_filename + '-ext.fits.gz', se_id,
+            gracedb.upload.s(new_skymap_filename + '-ext.multiorder.fits', se_id,
                              message, ['sky_loc', 'public']),
 
             skymaps.plot_allsky.s()
@@ -69,7 +71,7 @@ def get_skymap_filename(graceid):
     if 'S' in graceid:
         for message in reversed(gracedb_log):
             filename = message['filename']
-            if filename.endswith('.fits.gz'):
+            if filename.endswith('.multiorder.fits'):
                 return filename
     else:
         for message in reversed(gracedb_log):
@@ -80,6 +82,14 @@ def get_skymap_filename(graceid):
     raise ValueError('No skymap available for {0} yet.'.format(graceid))
 
 
+@app.task(shared=False)
+def _download_skymaps(se_filename, ext_filename, se_id, ext_id):
+    """Download both superevent and external sky map to be combined."""
+    se_skymap = gracedb.download(se_filename, se_id)
+    ext_skymap = gracedb.download(ext_filename, ext_id)
+    return se_skymap, ext_skymap
+
+
 @app.task(shared=False)
 def combine_skymaps(skymap1filebytes, skymap2filebytes):
     """This task combines the two input skymaps, in this case the external
@@ -90,8 +100,10 @@ def combine_skymaps(skymap1filebytes, skymap2filebytes):
             NamedTemporaryFile(content=skymap1filebytes) as skymap1file, \
             NamedTemporaryFile(content=skymap2filebytes) as skymap2file, \
             handling_system_exit():
-        ligo_skymap_combine.main([skymap1file.name,
-                                  skymap2file.name, combinedskymap.name])
+        sky1 = HealpixMap.read_map(skymap1file.name, density=True)
+        sky2 = HealpixMap.read_map(skymap2file.name, density=False)
+        comb_sky = sky1 * sky2
+        comb_sky.write_map(combinedskymap.name)
         return combinedskymap.read()
 
 

gwcelery/tasks/raven.py

@@ -47,6 +47,7 @@ def calculate_coincidence_far(superevent, exttrig, tl, th):
                    se_dict=superevent, ext_dict=exttrig,
                    grb_search=exttrig['search'],
                    se_fitsfile=se_skymap, ext_fitsfile=ext_skymap,
+                   se_moc=True, ext_moc=False,
                    incl_sky=True, gracedb=gracedb.client,
                    far_grb=exttrig['far'])
     else: