make certain plots optional

gwemopt/args.py

@@ -25,8 +25,12 @@ def parse_args(args):
         "--timeallocationType", help="time allocation type", default="powerlaw"
     )
 
-    parser.add_argument("--doPlots", action="store_true", default=False)
-    parser.add_argument("--doMovie", action="store_true", default=False)
+    parser.add_argument(
+        "--plots",
+        help="comma delimited list of plot types (skymap,tiles,coverage,schedule,efficiency)",
+        default="",
+    )
+    parser.add_argument("--movie", action="store_true", default=False)
     parser.add_argument("--doTiles", action="store_true", default=False)
     parser.add_argument("--tilesType", help="tiling type", default="moc")
     parser.add_argument("--doMindifFilt", action="store_true", default=False)

gwemopt/moc.py

@@ -81,8 +81,6 @@ def create_moc(params, map_struct=None):
                         continue
                     moc_struct[index] = moclists[ii]
             else:
-                moc = get_decam_quadrant_moc(tesselation[:, 1], tesselation[:, 2])
-
                 for ii, tess in tqdm(enumerate(tesselation), total=len(tesselation)):
                     index, ra, dec = tess[0], tess[1], tess[2]
                     if (

gwemopt/params.py

@@ -11,6 +11,7 @@
 import gwemopt
 from gwemopt.paths import CONFIG_DIR, REFS_DIR, TESSELATION_DIR
 from gwemopt.tiles import TILE_TYPES
+from gwemopt.utils import tesselation
 
 
 def params_struct(opts):
@@ -51,9 +52,9 @@ def params_struct(opts):
             )
             if not os.path.isfile(tessfile):
                 if params["config"][telescope]["FOV_type"] == "circle":
-                    gwemopt.tiles.tesselation_spiral(params["config"][telescope])
+                    tesselation.tesselation_spiral(params["config"][telescope])
                 elif params["config"][telescope]["FOV_type"] == "square":
-                    gwemopt.tiles.tesselation_packing(params["config"][telescope])
+                    tesselation.tesselation_packing(params["config"][telescope])
             if opts.tilesType == "galaxy":
                 params["config"][telescope]["tesselation"] = np.empty((3,))
             else:
@@ -166,4 +167,8 @@ def params_struct(opts):
         opts.parallelBackend if hasattr(opts, "parallelBackend") else "threading"
     )
 
+    params["movie"] = opts.movie if hasattr(opts, "movie") else False
+
+    params["plots"] = opts.plots.split(",") if hasattr(opts, "plots") else []
+
     return params

gwemopt/plotting/__init__.py

@@ -1,5 +1,5 @@
 from gwemopt.plotting.movie import make_movie
-from gwemopt.plotting.plot_coverage import make_coverage_plots
+from gwemopt.plotting.plot_coverage import make_coverage_movie, make_coverage_plots
 from gwemopt.plotting.plot_efficiency import make_efficiency_plots
 from gwemopt.plotting.plot_inclination import plot_inclination
 from gwemopt.plotting.plot_schedule import make_schedule_plots

gwemopt/plotting/plot_coverage.py

@@ -264,62 +264,67 @@ def plot_coverage_scaled(params, map_struct, coverage_struct, plot_sun_moon, max
     plt.savefig(plot_name, dpi=200)
     plt.close()
 
-    if params["doMovie"]:
-        print("Creating movie from schedule...")
 
-        idx = np.isfinite(coverage_struct["data"][:, 2])
-        mjd_min = np.min(coverage_struct["data"][idx, 2])
-        mjd_max = np.max(coverage_struct["data"][idx, 2])
-        mjd_N = 100
+def plot_coverage_movie(params, map_struct, coverage_struct, plot_sun_moon, max_time):
 
-        mjds = np.linspace(mjd_min, mjd_max, num=mjd_N)
-        moviedir = params["outputDir"].joinpath("movie")
-        moviedir.mkdir(exist_ok=True, parents=True)
+    hdu = map_struct["hdu"]
+    columns = [col.name for col in hdu.columns]
+
+    print("Creating movie from schedule...")
+
+    idx = np.isfinite(coverage_struct["data"][:, 2])
+    mjd_min = np.min(coverage_struct["data"][idx, 2])
+    mjd_max = np.max(coverage_struct["data"][idx, 2])
+    mjd_N = 100
+
+    mjds = np.linspace(mjd_min, mjd_max, num=mjd_N)
+    moviedir = params["outputDir"].joinpath("movie")
+    moviedir.mkdir(exist_ok=True, parents=True)
 
-        def make_plot(jj, params, map_struct, coverage_struct):
-            hdu = map_struct["hdu"]
-            columns = [col.name for col in hdu.columns]
+    def make_plot(jj, params, map_struct, coverage_struct):
+        hdu = map_struct["hdu"]
+        columns = [col.name for col in hdu.columns]
 
-            mjd = mjds[jj]
-            plot_name = moviedir.joinpath(f"coverage-{jj:04d}.png")
-            title = f"Coverage Map: {mjd:.2f}"
+        mjd = mjds[jj]
+        plot_name = moviedir.joinpath(f"coverage-{jj:04d}.png")
+        title = f"Coverage Map: {mjd:.2f}"
 
-            fig = plt.figure(figsize=(8, 6), dpi=100)
-            args = {"projection": params["projection"]}
-            if args["projection"] == "astro globe":
-                args["center"] = map_struct["center"]
-            ax = plt.axes([0.05, 0.05, 0.9, 0.9], **args)
-            ax.imshow_hpx(hdu, field=columns.index("PROB"), cmap="cylon")
-            ax.grid()
+        fig = plt.figure(figsize=(8, 6), dpi=100)
+        args = {"projection": params["projection"]}
+        if args["projection"] == "astro globe":
+            args["center"] = map_struct["center"]
+        ax = plt.axes([0.05, 0.05, 0.9, 0.9], **args)
+        ax.imshow_hpx(hdu, field=columns.index("PROB"), cmap="cylon")
+        ax.grid()
 
-            idx = np.where(coverage_struct["data"][:, 2] <= mjd)[0]
-            for ii in idx:
-                moc = coverage_struct["moc"][ii]
-                data = coverage_struct["data"][ii]
+        idx = np.where(coverage_struct["data"][:, 2] <= mjd)[0]
+        for ii in idx:
+            moc = coverage_struct["moc"][ii]
+            data = coverage_struct["data"][ii]
 
-                alpha = data[4] / max_time
-                if alpha > 1:
-                    alpha = 1.0
+            alpha = data[4] / max_time
+            if alpha > 1:
+                alpha = 1.0
 
-                moc.fill(
-                    ax=ax,
-                    wcs=ax.wcs,
-                    alpha=alpha,
-                    fill=True,
-                    color="black",
-                    linewidth=1,
-                )
+            moc.fill(
+                ax=ax,
+                wcs=ax.wcs,
+                alpha=alpha,
+                fill=True,
+                color="black",
+                linewidth=1,
+            )
 
-            plt.savefig(plot_name, dpi=200)
-            plt.close()
+        plt.savefig(plot_name, dpi=200)
+        plt.close()
 
-        for jj in tqdm(range(len(mjds))):
-            make_plot(jj, params, map_struct, coverage_struct)
+    for jj in tqdm(range(len(mjds))):
+        make_plot(jj, params, map_struct, coverage_struct)
 
-        moviefiles = moviedir.joinpath("coverage-%04d.png")
-        filename = params["outputDir"].joinpath("coverage.mpg")
+    moviefiles = moviedir.joinpath("coverage-%04d.png")
+    filename = params["outputDir"].joinpath("coverage.mpg")
 
-        make_movie(moviefiles, filename)
+    make_movie(moviefiles, filename)
 
 
 def make_coverage_plots(
@@ -337,3 +342,16 @@ def make_coverage_plots(
     )
 
     plot_coverage_scaled(params, map_struct, coverage_struct, plot_sun_moon, max_time)
+
+
+def make_coverage_movie(
+    params, map_struct, coverage_struct, plot_sun_moon: bool = True
+):
+
+    idx = np.isfinite(coverage_struct["data"][:, 4])
+    if not idx.size:
+        return
+
+    max_time = np.max(coverage_struct["data"][idx, 4])
+
+    plot_coverage_movie(params, map_struct, coverage_struct, plot_sun_moon, max_time)

gwemopt/run.py

@@ -15,6 +15,7 @@
 from gwemopt.params import params_struct
 from gwemopt.paths import DEFAULT_BASE_OUTPUT_DIR
 from gwemopt.plotting import (
+    make_coverage_movie,
     make_coverage_plots,
     make_efficiency_plots,
     make_tile_plots,
@@ -64,7 +65,7 @@ def run(args=None):
     else:
         catalog_struct = None
 
-    if args.doPlots:
+    if "skymap" in params["plots"]:
         print("Plotting skymap...")
         plot_skymap(params, map_struct)
         if args.inclination:
@@ -97,7 +98,7 @@ def run(args=None):
         else:
             raise ValueError(f"Unknown tilesType: {params['tilesType']}")
 
-        if args.doPlots:
+        if "tiles" in params["plots"]:
             print("Plotting tiles struct...")
             make_tile_plots(params, map_struct, tile_structs)
 
@@ -120,10 +121,21 @@ def run(args=None):
         print("Summary of coverage...")
         summary(params, map_struct, coverage_struct, catalog_struct=catalog_struct)
 
-        print("Plotting coverage...")
-        make_coverage_plots(
-            params, map_struct, coverage_struct, catalog_struct=catalog_struct
-        )
+        if "coverage" in params["plots"]:
+            print("Plotting coverage...")
+            make_coverage_plots(
+                params,
+                map_struct,
+                coverage_struct,
+                catalog_struct=catalog_struct,
+            )
+
+        if params["movie"]:
+            make_coverage_movie(
+                params,
+                map_struct,
+                coverage_struct,
+            )
 
     if args.doEfficiency:
         if args.doSchedule or args.doCoverage:
@@ -155,7 +167,7 @@ def run(args=None):
                         f'{efficiency_structs[key]["3D"]*100:.2f}% '
                     )
 
-            if args.doPlots:
+            if "efficiency" in params["plots"]:
                 print("Plotting efficiency...")
                 make_efficiency_plots(params, map_struct, efficiency_structs)
         else:

gwemopt/scheduler.py

@@ -639,7 +639,7 @@ def scheduler(params, config_struct, tile_struct):
     else:
         raise ValueError(f'Unknown solverType {params["solverType"]}')
 
-    if params["doPlots"]:
+    if "schedule" in params["plots"]:
         from gwemopt.plotting import make_schedule_plots
 
         make_schedule_plots(params, exposurelist, keys)

gwemopt/tests/test_coverage.py

@@ -43,7 +43,6 @@ def test_coverage():
             "-e",
             str(test_skymap),
             "--doTiles",
-            "--doPlots",
             "--doCoverage",
             "--coverageFiles",
             ",".join(schedule_list),

gwemopt/tests/test_schedule.py

@@ -70,7 +70,6 @@ def test_scheduler():
                 "-e",
                 str(test_skymap),
                 "--doTiles",
-                "--doPlots",
                 "--doSchedule",
                 "--timeallocationType",
                 "powerlaw",

gwemopt/tiles.py

@@ -1098,65 +1098,3 @@ def compute_tiles_map(
         )
 
     return vals
-
-
-def tesselation_spiral(config_struct, scale=0.80):
-    if config_struct["FOV_type"] == "square":
-        FOV = config_struct["FOV"] * config_struct["FOV"] * scale
-    elif config_struct["FOV_type"] == "circle":
-        FOV = np.pi * config_struct["FOV"] * config_struct["FOV"] * scale
-
-    area_of_sphere = 4 * np.pi * (180 / np.pi) ** 2
-    n = int(np.ceil(area_of_sphere / FOV))
-    print("Using %d points to tile the sphere..." % n)
-
-    golden_angle = np.pi * (3 - np.sqrt(5))
-    theta = golden_angle * np.arange(n)
-    z = np.linspace(1 - 1.0 / n, 1.0 / n - 1, n)
-    radius = np.sqrt(1 - z * z)
-
-    points = np.zeros((n, 3))
-    points[:, 0] = radius * np.cos(theta)
-    points[:, 1] = radius * np.sin(theta)
-    points[:, 2] = z
-
-    ra, dec = hp.pixelfunc.vec2ang(points, lonlat=True)
-    fid = open(config_struct["tesselationFile"], "w")
-    for ii in range(len(ra)):
-        fid.write("%d %.5f %.5f\n" % (ii, ra[ii], dec[ii]))
-    fid.close()
-
-
-def tesselation_packing(config_struct, scale=0.97):
-    sphere_radius = 1.0
-    if config_struct["FOV_type"] == "square":
-        circle_radius = np.deg2rad(config_struct["FOV"] / 2.0) * scale
-    elif config_struct["FOV_type"] == "circle":
-        circle_radius = np.deg2rad(config_struct["FOV"]) * scale
-    vertical_count = int((np.pi * sphere_radius) / (2 * circle_radius))
-
-    phis = []
-    thetas = []
-
-    phi = -0.5 * np.pi
-    phi_step = np.pi / vertical_count
-    while phi < 0.5 * np.pi:
-        horizontal_count = int(
-            (2 * np.pi * np.cos(phi) * sphere_radius) / (2 * circle_radius)
-        )
-        if horizontal_count == 0:
-            horizontal_count = 1
-        theta = 0
-        theta_step = 2 * np.pi / horizontal_count
-        while theta < 2 * np.pi - 1e-8:
-            phis.append(phi)
-            thetas.append(theta)
-            theta += theta_step
-        phi += phi_step
-    dec = np.array(np.rad2deg(phis))
-    ra = np.array(np.rad2deg(thetas))
-
-    fid = open(config_struct["tesselationFile"], "w")
-    for ii in range(len(ra)):
-        fid.write("%d %.5f %.5f\n" % (ii, ra[ii], dec[ii]))
-    fid.close()