Support flat GW sky maps as well

gwcelery/tasks/external_skymaps.py

@@ -24,14 +24,16 @@ from ..import _version
 
 
 @app.task(shared=False)
-def create_combined_skymap(se_id, ext_id):
+def create_combined_skymap(se_id, ext_id, gw_moc=True):
     """Creates and uploads the combined LVC-Fermi skymap.
 
     This also uploads the external trigger skymap to the external trigger
     GraceDB page.
     """
-    se_skymap_filename = get_skymap_filename(se_id)
+    se_skymap_filename = get_skymap_filename(se_id, moc=gw_moc)
     ext_skymap_filename = get_skymap_filename(ext_id)
+    new_filename = \
+        'combined-ext.multiorder.fits' if gw_moc else 'combined-ext.fits.gz'
 
     message = 'Combined LVC-external sky map using {0} and {1}'.format(
         se_skymap_filename, ext_skymap_filename)
@@ -39,7 +41,7 @@ def create_combined_skymap(se_id, ext_id):
         'Mollweide projection of <a href="/api/events/{graceid}/files/'
         '{filename}">{filename}</a>').format(
             graceid=ext_id,
-            filename='combined-ext.multiorder.fits')
+            filename=new_filename)
 
     (
         _download_skymaps.si(
@@ -49,7 +51,7 @@ def create_combined_skymap(se_id, ext_id):
         combine_skymaps.s()
         |
         group(
-            gracedb.upload.s('combined-ext.multiorder.fits', ext_id,
+            gracedb.upload.s(new_filename, ext_id,
                              message, ['sky_loc', 'ext_coinc']),
 
             skymaps.plot_allsky.s()
@@ -68,18 +70,21 @@ def create_combined_skymap(se_id, ext_id):
 
 @app.task(autoretry_for=(ValueError,), retry_backoff=10,
           retry_backoff_max=600)
-def get_skymap_filename(graceid):
+def get_skymap_filename(graceid, moc=True):
     """Get the skymap fits filename.
 
     If not available, will try again 10 seconds later, then 20, then 40, etc.
     until up to 10 minutes after initial attempt.
     """
     gracedb_log = gracedb.get_log(graceid)
+    ending_string = 'multiorder.fits' if moc else 'fits.gz'
     if 'S' in graceid:
+        # Try first to get a multiordered sky map
         for message in reversed(gracedb_log):
             filename = message['filename']
-            if filename.endswith('.multiorder.fits'):
+            if filename.endswith(ending_string):
                 return filename
+        # Try next to get a flattened sky map
     else:
         for message in reversed(gracedb_log):
             filename = message['filename']
@@ -108,38 +113,46 @@ def combine_skymaps(skymapfilesbytes):
             NamedTemporaryFile(content=skymap1filebytes) as skymap1file, \
             NamedTemporaryFile(content=skymap2filebytes) as skymap2file, \
             handling_system_exit():
-        #  FIXME: Use method that regrids the combined sky map e.g. mhealpy
-        #  Load sky maps
-        gw_sky = read_sky_map(skymap1file.name, moc=True)
-        ext_sky, ext_header = read_sky_map(skymap2file.name, moc=False)
-
-        #  Find ra/dec of each GW pixel
-        level, ipix = ah.uniq_to_level_ipix(gw_sky["UNIQ"])
-        nsides = ah.level_to_nside(level)
-        areas = ah.nside_to_pixel_area(nsides)
-        ra_gw, dec_gw = ah.healpix_to_lonlat(ipix, nsides, order='nested')
-
-        #  Find corresponding external sky map indicies
-        ext_nside = ah.npix_to_nside(len(ext_sky))
-        ext_ind = ah.lonlat_to_healpix(
-            ra_gw, dec_gw, ext_nside,
-            order='nested' if ext_header['nest'] else 'ring')
-
-        #  Weight GW prob density by external sky map probabilities
-        gw_sky['PROBDENSITY'] *= ext_sky[ext_ind]
-        gw_sky['PROBDENSITY'] /= np.sum(gw_sky['PROBDENSITY'] * areas).value
-
-        #  Modify GW sky map with new data
-        distmean, diststd = parameters_to_marginal_moments(
-            gw_sky['PROBDENSITY'] * areas.value,
-            gw_sky['DISTMU'], gw_sky['DISTSIGMA'])
-        gw_sky.meta['distmean'], gw_sky.meta['diststd'] = distmean, diststd
-        gw_sky.meta['instruments'].update(ext_header['instruments'])
-        gw_sky.meta['HISTORY'].extend([
-            '', 'The values were reweighted by using data from {}'.format(
-                list(ext_header['instruments'])[0])])
-
-        write_sky_map(combinedskymap.name, gw_sky)
+
+        # Determine whether GW sky map is multiordered or flat
+        gw_moc = '.multiorder.fits' in skymap1file.name
+
+        if gw_moc:
+            #  FIXME: Use method that regrids the combined sky map e.g. mhealpy
+            #  Load sky maps
+            gw_sky = read_sky_map(skymap1file.name, moc=True)
+            ext_sky, ext_header = read_sky_map(skymap2file.name, moc=False)
+
+            #  Find ra/dec of each GW pixel
+            level, ipix = ah.uniq_to_level_ipix(gw_sky["UNIQ"])
+            nsides = ah.level_to_nside(level)
+            areas = ah.nside_to_pixel_area(nsides)
+            ra_gw, dec_gw = ah.healpix_to_lonlat(ipix, nsides, order='nested')
+
+            #  Find corresponding external sky map indicies
+            ext_nside = ah.npix_to_nside(len(ext_sky))
+            ext_ind = ah.lonlat_to_healpix(
+                ra_gw, dec_gw, ext_nside,
+                order='nested' if ext_header['nest'] else 'ring')
+
+            #  Weight GW prob density by external sky map probabilities
+            gw_sky['PROBDENSITY'] *= ext_sky[ext_ind]
+            gw_sky['PROBDENSITY'] /= \
+                np.sum(gw_sky['PROBDENSITY'] * areas).value
+
+            #  Modify GW sky map with new data
+            distmean, diststd = parameters_to_marginal_moments(
+                gw_sky['PROBDENSITY'] * areas.value,
+                gw_sky['DISTMU'], gw_sky['DISTSIGMA'])
+            gw_sky.meta['distmean'], gw_sky.meta['diststd'] = distmean, diststd
+            gw_sky.meta['instruments'].update(ext_header['instruments'])
+            gw_sky.meta['HISTORY'].extend([
+                '', 'The values were reweighted by using data from {}'.format(
+                    list(ext_header['instruments'])[0])])
+            write_sky_map(combinedskymap.name, gw_sky)
+        else:
+            ligo_skymap_combine.main([skymap1file.name,
+                                      skymap2file.name, combinedskymap.name])
         return combinedskymap.read()
 
 

gwcelery/tasks/external_triggers.py

@@ -313,7 +313,7 @@ def handle_grb_igwn_alert(alert):
                                         superevent, tl, th, gw_group),
 
                 external_skymaps.create_combined_skymap.si(
-                    superevent_id, ext_id
+                    superevent_id, ext_id, gw_moc='CBC'==gw_group
                 )
             ).delay()
         elif 'EM_COINC' in alert['object']['labels']: