MOC calculation (#154)

This PR uses MOCs.

gwemopt/args.py

@@ -37,7 +37,6 @@ def parse_args(args):
     parser.add_argument("--doPerturbativeTiling", action="store_true", default=False)
     parser.add_argument("--doOrderByObservability", action="store_true", default=False)
 
-    parser.add_argument("--doUseCatalog", action="store_true", default=False)
     parser.add_argument("--catalogDir", help="catalog directory", default=CATALOG_DIR)
     parser.add_argument(
         "--catalog", help="Galaxy catalog name (e.g GLADE)", default=None
@@ -54,7 +53,6 @@ def parse_args(args):
         type=int,
         default=2000,
     )
-    parser.add_argument("--doObservability", action="store_true", default=False)
 
     parser.add_argument("--doEfficiency", action="store_true", default=False)
     parser.add_argument(
@@ -83,7 +81,6 @@ def parse_args(args):
 
     parser.add_argument("--Ndet", default=1, type=int)
     parser.add_argument("--nside", default=256, type=int)
-    parser.add_argument("--DScale", default=1.0, type=float)
     parser.add_argument("--Tobs", default="0.0,1.0")
 
     parser.add_argument("--mindiff", default=0.0, type=float)
@@ -111,8 +108,6 @@ def parse_args(args):
     parser.add_argument("--doSingleExposure", action="store_true", default=False)
     parser.add_argument("--filters", default="r,g,r")
     parser.add_argument("--doAlternatingFilters", action="store_true", default=False)
-    parser.add_argument("--doRASlices", action="store_true", default=False)
-    parser.add_argument("--nside_down", default=2, type=int)
     parser.add_argument("--max_filter_sets", default=4, type=int)
     parser.add_argument("--iterativeOverlap", default=0.0, type=float)
     parser.add_argument("--maximumOverlap", default=1.0, type=float)
@@ -137,10 +132,8 @@ def parse_args(args):
     parser.add_argument("--doBlocks", action="store_true", default=False)
     parser.add_argument("--Nblocks", default=4, type=int)
 
-    parser.add_argument("--doRASlice", action="store_true", default=False)
-    parser.add_argument("--raslice", default="0.0,24.0")
-
     parser.add_argument("--absmag", default=-15.0, type=float)
+    parser.add_argument("--raslice", default=None)
 
     parser.add_argument(
         "--galactic_limit", help="Galactic limit.", default=0.0, type=float
@@ -155,10 +148,9 @@ def parse_args(args):
     parser.add_argument("--true_inclination", default=0.0, type=float)
 
     parser.add_argument("--doTrueLocation", action="store_true", default=False)
-    parser.add_argument("--observability_thresh", default=0.05, type=float)
-    parser.add_argument("--doObservabilityExit", action="store_true", default=False)
 
     parser.add_argument("--inclination", action="store_true", default=False)
+    parser.add_argument("--projection", default="astro mollweide")
 
     parser.add_argument("--solverType", default="heuristic")
     parser.add_argument("--milpSolver", default="PULP_CBC_CMD")

gwemopt/catalogs/get.py

@@ -2,8 +2,10 @@
 
 import healpy as hp
 import numpy as np
+from astropy.table import Table
 from astroquery.vizier import Vizier
 from ligo.skymap import distance
+from ligo.skymap.bayestar import derasterize
 from scipy.stats import norm
 
 from gwemopt.catalogs.clu import CluCatalog
@@ -63,7 +65,7 @@ def get_catalog(params, map_struct, export_catalog: bool = True):
 
     n, cl = params["powerlaw_n"], params["powerlaw_cl"]
 
-    prob_scaled = copy.deepcopy(map_struct["prob"])
+    prob_scaled = copy.deepcopy(map_struct["skymap_raster"]["PROB"])
     prob_sorted = np.sort(prob_scaled)[::-1]
     prob_indexes = np.argsort(prob_scaled)[::-1]
     prob_cumsum = np.cumsum(prob_sorted)
@@ -72,34 +74,42 @@ def get_catalog(params, map_struct, export_catalog: bool = True):
     prob_scaled = prob_scaled**n
 
     ipix = hp.ang2pix(
-        map_struct["nside"],
+        params["nside"],
         np.array(cat_df["ra"]),
         np.array(cat_df["dec"]),
         lonlat=True,
     )
 
-    if "distnorm" in map_struct:
-        if map_struct["distnorm"] is not None:
+    if "DISTNORM" in map_struct:
+        if map_struct["skymap_raster"]["DISTNORM"] is not None:
             # creat an mask to cut at 3 sigma in distance
             mask = np.zeros(len(cat_df["distmpc"]))
 
-            # calculate the moments from distmu, distsigma and distnorm
-            mom_mean, mom_std, mom_norm = distance.parameters_to_moments(
-                map_struct["distmu"], map_struct["distsigma"]
-            )
-
             condition_indexer = np.where(
-                (cat_df["distmpc"] < (mom_mean[ipix] + (3 * mom_std[ipix])))
-                & (cat_df["distmpc"] > (mom_mean[ipix] - (3 * mom_std[ipix])))
+                (
+                    cat_df["distmpc"]
+                    < (
+                        map_struct["skymap_raster"]["DISTMEAN"][ipix]
+                        + (3 * map_struct["skymap_raster"]["DISTSTD"][ipix])
+                    )
+                )
+                & (
+                    cat_df["distmpc"]
+                    > (
+                        map_struct["skymap_raster"]["DISTMEAN"][ipix]
+                        - (3 * map_struct["skymap_raster"]["DISTSTD"][ipix])
+                    )
+                )
             )
             mask[condition_indexer] = 1
 
             s_loc = (
                 prob_scaled[ipix]
                 * (
-                    map_struct["distnorm"][ipix]
+                    map_struct["skymap_raster"]["DISTNORM"][ipix]
                     * norm(
-                        map_struct["distmu"][ipix], map_struct["distsigma"][ipix]
+                        map_struct["skymap_raster"]["DISTMU"][ipix],
+                        map_struct["skymap_raster"]["DISTSIGMA"][ipix],
                     ).pdf(cat_df["distmpc"])
                 )
                 ** params["powerlaw_dist_exp"]
@@ -202,7 +212,7 @@ def get_catalog(params, map_struct, export_catalog: bool = True):
         s_mass = s_loc * (1 + (alpha_mass * beta_mass * s_mass))
 
     s = np.array(s_loc * Slum * sdet)
-    prob = np.zeros(map_struct["prob"].shape)
+    prob = np.zeros(map_struct["skymap_raster"]["PROB"].shape)
     if params["galaxy_grade"] == "Sloc":
         for j in range(len(ipix)):
             prob[ipix[j]] += s_loc[j]
@@ -220,11 +230,12 @@ def get_catalog(params, map_struct, export_catalog: bool = True):
             "You are trying to use a galaxy grade that is not implemented yet."
         )
 
+    prob[np.isnan(prob)] = 0.0
     prob = prob / np.sum(prob)
 
-    map_struct["prob_catalog"] = prob
-    if params["doUseCatalog"]:
-        map_struct["prob"] = prob
+    map_struct["skymap_raster"]["PROB"] = prob
+    skymap = map_struct["skymap_raster"].copy()
+    map_struct["skymap"] = derasterize(skymap)
 
     cat_df["grade"] = grade
     cat_df["S"] = s

gwemopt/chipgaps/__init__.py

@@ -2,5 +2,5 @@
 Module for chip gaps.
 """
 
-from gwemopt.chipgaps.decam import get_decam_quadrant_ipix
-from gwemopt.chipgaps.ztf import get_ztf_quadrant_ipix
+from gwemopt.chipgaps.decam import get_decam_quadrant_moc
+from gwemopt.chipgaps.ztf import get_ztf_quadrant_moc

gwemopt/chipgaps/decam.py

@@ -24,6 +24,7 @@
 from astropy import units as u
 from astropy.coordinates import ICRS, SkyCoord
 from astropy.wcs import WCS
+from mocpy import MOC
 
 
 class DECamtile:
@@ -210,7 +211,7 @@ def get_decam_ccds(ra, dec, save_footprint=False):
     return np.transpose(offsets, (2, 0, 1))
 
 
-def get_decam_quadrant_ipix(nside, ra, dec):
+def get_decam_quadrant_moc(ra, dec):
     ccd_coords = get_decam_ccds(0, 0)
 
     skyoffset_frames = SkyCoord(ra, dec, unit=u.deg).skyoffset_frame()
@@ -219,13 +220,15 @@ def get_decam_quadrant_ipix(nside, ra, dec):
         unit=u.deg,
         frame=skyoffset_frames[:, np.newaxis, np.newaxis],
     ).transform_to(ICRS)
-    ccd_xyz = np.moveaxis(ccd_coords_icrs.cartesian.xyz.value, 0, -1)[0]
 
-    ipixs = []
-    for xyz in ccd_xyz:
-        ipix = hp.query_polygon(nside, xyz)
-        ipixs.append(ipix.tolist())
-    return ipixs
+    moc = None
+    for subfield_id, coords in enumerate(ccd_coords_icrs[0]):
+        if moc is None:
+            moc = MOC.from_polygon_skycoord(coords)
+        else:
+            moc = moc + MOC.from_polygon_skycoord(coords)
+
+    return moc
 
 
 def ccd_xy_to_radec(alpha_p, delta_p, ccd_centers_xy):

gwemopt/chipgaps/ztf.py

@@ -25,6 +25,7 @@
 from astropy import units as u
 from astropy.coordinates import ICRS, SkyCoord
 from astropy.wcs import WCS
+from mocpy import MOC
 
 
 class ZTFtile:
@@ -273,7 +274,7 @@ def get_ztf_quadrants():
     return np.transpose(offsets, (2, 0, 1))
 
 
-def get_ztf_quadrant_ipix(nside, ra, dec, subfield_ids=None):
+def get_ztf_quadrant_moc(ra, dec, subfield_ids=None):
     quadrant_coords = get_ztf_quadrants()
 
     skyoffset_frames = SkyCoord(ra, dec, unit=u.deg).skyoffset_frame()
@@ -282,13 +283,16 @@ def get_ztf_quadrant_ipix(nside, ra, dec, subfield_ids=None):
         unit=u.deg,
         frame=skyoffset_frames[:, np.newaxis, np.newaxis],
     ).transform_to(ICRS)
-    quadrant_xyz = np.moveaxis(quadrant_coords_icrs.cartesian.xyz.value, 0, -1)[0]
 
-    ipixs = []
-    for subfield_id, xyz in enumerate(quadrant_xyz):
+    moc = None
+    for subfield_id, coords in enumerate(quadrant_coords_icrs[0]):
         if subfield_ids is not None:
             if subfield_id not in subfield_ids:
                 continue
-        ipix = hp.query_polygon(nside, xyz)
-        ipixs.append(ipix.tolist())
-    return ipixs
+
+        if moc is None:
+            moc = MOC.from_polygon_skycoord(coords)
+        else:
+            moc = moc + MOC.from_polygon_skycoord(coords)
+
+    return moc