Add method to reweight GW sky map using external; no longer regrid

gwcelery/tasks/external_skymaps.py

 """Create and upload external sky maps."""
 from astropy import units as u
 from astropy.coordinates import ICRS, SkyCoord
+import astropy_healpix as ah
 from astropy_healpix import HEALPix, pixel_resolution_to_nside
 from celery import group
 #  import astropy.utils.data
 import numpy as np
-from ligo.skymap.io import fits
+from ligo.skymap.io import fits, read_sky_map, write_sky_map
 import gcn
 import healpy as hp
 import lxml.etree
@@ -55,6 +56,8 @@ def create_combined_skymap(se_id, ext_id):
             gracedb.upload.s('combined-ext.png', se_id,
                              message_png, ['sky_loc', 'ext_coinc', 'public'])
         )
+        |
+        gracedb.create_label.si('COMBINEDSKYMAP_READY', se_id)
     ).delay()
 
 
@@ -100,10 +103,29 @@ def combine_skymaps(skymapfilesbytes):
             NamedTemporaryFile(content=skymap1filebytes) as skymap1file, \
             NamedTemporaryFile(content=skymap2filebytes) as skymap2file, \
             handling_system_exit():
-        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, column_names='PROBDENSITY')
+        #  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
+        write_sky_map(combinedskymap.name, gw_skymap)
         return combinedskymap.read()