tickets/SP-2518: Enhance skyproj to support additional coverage plots

schedview/collect/visits.py

@@ -143,7 +143,10 @@ def read_ddf_visits(*args, **kwargs) -> pd.DataFrame:
     ddf_field_names = tuple(ddf_locations().keys())
     ddf_visits = []
     for ddf_name in ddf_field_names:
-        this_field_visits = all_visits.query(f"{target_column_name}.str.contains('{ddf_name}')")
+        matches_target_column = all_visits[target_column_name].str.contains(ddf_name, case=False)
+        matches_note_column = all_visits["scheduler_note"].str.contains(f"DD:{ddf_name}", case=False)
+        this_field_visits = all_visits.loc[matches_target_column | matches_note_column, :]
+
         # Add a column with just the DDF field name and nothing else,
         # so that if the same DDF appears in different places
         # in the target_name value they still get identified as the

schedview/plot/colors.py

@@ -4,12 +4,12 @@
 # Follow RTN-045 for mapping filters to plot colors
 
 PLOT_BAND_COLORS = {
-    "u": "#0c71ff",
-    "g": "#49be61",
-    "r": "#c61c00",
-    "i": "#ffc200",
-    "z": "#f341a2",
-    "y": "#5d0000",
+    "u": "#1600ea",
+    "g": "#31de1f",
+    "r": "#b52626",
+    "i": "#370201",
+    "z": "#ba52ff",
+    "y": "#61a2b3",
 }
 
 # Use
@@ -19,51 +19,50 @@
 #    colorblind_safe=True,
 #    grid_size=256
 # )
-# to generate extra colors, attempting to get distinguishable ones, even
-# for the color blind.
+# to generate extra colors, attempting to get distinguishable ones.
 EXTRA_COLORS = [
-    "#0000b6",
-    "#04e9ff",
-    "#006200",
-    "#750065",
-    "#ff3ce6",
-    "#d75d47",
-    "#bc3c7f",
-    "#04f0bd",
-    "#6302ac",
-    "#b408b9",
-    "#81b000",
-    "#ff79b9",
-    "#640500",
-    "#65a9ff",
-    "#b20000",
-    "#00f468",
-    "#009236",
-    "#0135ff",
-    "#2ac879",
-    "#d43cff",
-    "#ff864e",
-    "#95085e",
-    "#ff95ff",
-    "#f721b0",
-    "#9401d0",
-    "#ad0e8a",
-    "#370091",
-    "#7c04ff",
-    "#7a0890",
-    "#bf54b1",
-    "#f94f00",
-    "#b45900",
-    "#d783ff",
-    "#0893d7",
-    "#b4c900",
-    "#862113",
-    "#8f5bd0",
-    "#fe75d6",
-    "#1aaf49",
-    "#006bbb",
-    "#0c03fb",
-    "#7b9000",
+    "#007000",
+    "#ff8c72",
+    "#006169",
+    "#000073",
+    "#9461ff",
+    "#002e00",
+    "#d8a2ff",
+    "#26af00",
+    "#dc1a25",
+    "#36164c",
+    "#10d9aa",
+    "#48936f",
+    "#0113ff",
+    "#84010e",
+    "#9d7cb5",
+    "#01593c",
+    "#7c3ca8",
+    "#e2b200",
+    "#015cff",
+    "#924e47",
+    "#00b2fd",
+    "#81c9d3",
+    "#ff5701",
+    "#2b8683",
+    "#00457c",
+    "#c87168",
+    "#0ed767",
+    "#0080aa",
+    "#7c7f31",
+    "#3702a6",
+    "#364801",
+    "#983e00",
+    "#5dae9d",
+    "#540005",
+    "#79a762",
+    "#c073f8",
+    "#0d7849",
+    "#280162",
+    "#73a9c9",
+    "#013ed6",
+    "#5a3c71",
+    "#9804f1",
 ]
 
 

schedview/plot/survey_skyproj.py

@@ -1,9 +1,11 @@
 import datetime
 from functools import partial
+from typing import Callable, Literal, Tuple, TypeAlias
 
 import astropy.units as u
 import colorcet
 import healpy as hp
+import matplotlib as mpl
 import matplotlib.pyplot as plt
 import numpy as np
 import pandas as pd
@@ -15,11 +17,12 @@
 
 # Import ModelObservatory to make sphinx happy
 from rubin_scheduler.scheduler.model_observatory import ModelObservatory  # noqa F401
+from rubin_sim import maf
 
 import schedview.compute.astro
 from schedview import band_column
 from schedview.compute.camera import LsstCameraFootprintPerimeter
-from schedview.compute.maf import compute_hpix_metric_in_bands
+from schedview.compute.maf import compute_hpix_metric_in_bands, compute_metric
 
 DEFAULT_MAP_KWARGS = {"cmap": colorcet.cm.blues}
 
@@ -30,6 +33,12 @@
     extent=[0.0, 360.0, -90, 89],
 )
 
+HorizonOptions: TypeAlias = Tuple[
+    Literal["morning west"], Literal["morning east"], Literal["evening west"], Literal["evening east"]
+]
+
+DEFAULT_HORIZONS: HorizonOptions = ("morning west", "morning east", "evening west", "evening east")
+
 
 def compute_circle_points(
     center_ra,
@@ -93,7 +102,17 @@ def compute_circle_points(
     return circle
 
 
-def map_healpix(map_hpix, model_observatory, night_events, axes=None, sky_map_factory=None, **kwargs):
+def map_healpix(
+    map_hpix,
+    model_observatory,
+    night_events,
+    axes=None,
+    sky_map_factory=None,
+    horizons: HorizonOptions | None = DEFAULT_HORIZONS,
+    horizon_zd: float = 60.0,
+    colorbar: bool = True,
+    **kwargs,
+):
     """Plots visits over a healpix map, with astronomical annotations.
 
     Parameters
@@ -109,6 +128,8 @@ def map_healpix(map_hpix, model_observatory, night_events, axes=None, sky_map_fa
         A matplotlib set of axes.
     sky_map_factory : `skyproj.skyproj.LaeaSkyproj` or `skyproj.McBrydeSkyproj`
         Factory to make sky_map instances
+    colorbar : `bool`
+        Show a colorbar. Defaults to True.
     **kwargs
         Additional keyword arguments are passed to
         `skyproj.skyproj.LaeaSkyproj.draw_hpxmap`.
@@ -152,24 +173,42 @@ def map_healpix(map_hpix, model_observatory, night_events, axes=None, sky_map_fa
     sun_moon_positions = model_observatory.almanac.get_sun_moon_positions(mjd)
     sun_ra = np.degrees(sun_moon_positions["sun_RA"].item())
     sun_decl = np.degrees(sun_moon_positions["sun_dec"].item())
-    sky_map.ax.scatter(sun_ra, sun_decl, color="yellow")
+    sky_map.ax.scatter(sun_ra, sun_decl, color="brown")
 
     moon_ra = np.degrees(sun_moon_positions["moon_RA"])
     moon_decl = np.degrees(sun_moon_positions["moon_dec"])
     sky_map.ax.scatter(moon_ra, moon_decl, color="orange")
 
     # Night limit horizons
-    latitude = model_observatory.site.latitude
-    zd = 70
-    evening = compute_circle_points(night_events.loc["sun_n12_setting", "LST"], latitude, zd, 180, 360)
-    sky_map.ax.plot(evening.ra, evening.decl, color="red", linestyle="dashed")
-    evening = compute_circle_points(night_events.loc["sun_n12_setting", "LST"], latitude, zd, 0, 180)
-    sky_map.ax.plot(evening.ra, evening.decl, color="red", linestyle="dotted")
-
-    morning = compute_circle_points(night_events.loc["sun_n12_rising", "LST"], latitude, zd, 0, 180)
-    sky_map.ax.plot(morning.ra, morning.decl, color="red", linestyle="dashed")
-    morning = compute_circle_points(night_events.loc["sun_n12_rising", "LST"], latitude, zd, 180, 360)
-    sky_map.ax.plot(morning.ra, morning.decl, color="red", linestyle="dotted")
+    if horizons is not None and len(horizons) > 0:
+        assert isinstance(horizons, tuple)
+
+        latitude = model_observatory.site.latitude
+        if "morning west" in horizons:
+            evening = compute_circle_points(
+                night_events.loc["sun_n12_setting", "LST"], latitude, horizon_zd, 180, 360
+            )
+            sky_map.ax.plot(evening.ra, evening.decl, color="red", linestyle="dashed")
+        if "morning east" in horizons:
+            evening = compute_circle_points(
+                night_events.loc["sun_n12_setting", "LST"], latitude, horizon_zd, 0, 180
+            )
+            sky_map.ax.plot(evening.ra, evening.decl, color="red", linestyle="dotted")
+
+        if "evening east" in horizons:
+            morning = compute_circle_points(
+                night_events.loc["sun_n12_rising", "LST"], latitude, horizon_zd, 0, 180
+            )
+            sky_map.ax.plot(morning.ra, morning.decl, color="red", linestyle="dashed")
+        if "evening west" in horizons:
+            morning = compute_circle_points(
+                night_events.loc["sun_n12_rising", "LST"], latitude, horizon_zd, 180, 360
+            )
+            sky_map.ax.plot(morning.ra, morning.decl, color="red", linestyle="dotted")
+
+    if colorbar:
+        sky_map.draw_colorbar(location="bottom", pad=0.14)
+
     return sky_map
 
 
@@ -364,3 +403,148 @@ def create_metric_map_grid(metric, visits, observatory, nside=32, **kwargs) -> F
     night_events = schedview.compute.astro.night_events(day_obs_date)
     fig = create_hpix_map_grid(metric_hpix, observatory, night_events, **kwargs)
     return fig
+
+
+def map_count_healpix(
+    *args,
+    num_colors: int = 6,
+    cmap: str | mpl.colors.Colormap = "cividis_r",
+    clip: bool = True,
+    colorbar: bool = True,
+    **kwargs,
+) -> skyproj.skyproj.LaeaSkyproj:
+    """Plot a healpix map appropriate when values are low integers greater
+    than zero.
+
+    Parameters
+    ----------
+    *args
+        Positional arguments forwarded to :func:`map_healpix`.  Usually
+        ``map_hpix``, ``model_observatory``, ``night_events`` and an optional
+        ``axes`` instance.
+    num_colors : `int`, optional
+        Maximum number of integers.
+        Defaults to 6.
+    cmap : 'str' or `matplotlib.colors.Colormap`, optional
+        Colormap name or instance to use for the discrete bins. If a string,
+        the colormap is obtained via ``plt.get_cmap(cmap, num_colors)``.
+    clip : `bool`, optional
+        If ``True`` (default), values above the highest bin are clipped
+        and the colorbar shows a “+” suffix for the final tick label.
+    **kwargs
+        Additional keyword arguments forwarded to :func:`map_healpix`.
+
+    Returns
+    -------
+    sky_map : `skyproj.skyproj.LaeaSkyproj`
+        The sky projection object with the healpix map drawn and a colorbar
+        added.  The map is displayed with discrete colors corresponding to
+        the specified number of bins.
+    """
+    if isinstance(cmap, str):
+        cmap = plt.get_cmap(cmap, num_colors)
+
+    norm = mpl.colors.BoundaryNorm(np.arange(num_colors + 1) + 0.5, ncolors=num_colors, clip=clip)
+
+    sky_map = map_healpix(*args, cmap=cmap, norm=norm, colorbar=False, **kwargs)
+
+    if colorbar:
+        ticks = 1 + np.arange(num_colors)
+        cbar = sky_map.draw_colorbar(ticks=ticks, location="bottom", pad=0.1)
+
+        if clip:
+            cbar_ticks_labels = [t.get_text() for t in cbar.ax.get_xticklabels()]
+            cbar_ticks_labels[-1] = cbar_ticks_labels[-1] + "+"
+            cbar.set_ticks(ticks, labels=cbar_ticks_labels)
+
+    return sky_map
+
+
+def map_hpix_in_laea_and_mcbryde(
+    map_function: Callable, *args, lon_0: float = 0.0, figsize: tuple = (15, 10), **kwargs
+) -> Figure:
+    """Map a healpix array in both LAEA and McBryde projections, side by side.
+
+    Parameters
+    ----------
+    map_function: `Callable`
+        A function that takes sky_map_function and axes arguments, and plots
+        the healpix map on these axes. The ``axes`` agument must take an
+        instance of `matplotlib.axes.Axes`, and ``sky_map_factory`` must be
+        a dictionary with keys ``laea`` and ``mcbryde`` with values that
+        are functions that return instances of `skyproj.LaeaSkyproj` and
+        `skyproj.McBrydeSkyproj`, respectively.
+    lon_0: `float`
+        The central longitude for the McBryde projection.
+    figsize: `tuple`
+        Figure size, passed to ``plt.sublots``
+    *args:
+        Additional positional arguments are passed as positional arguments
+        to ``map_function``.
+    **kwargs:
+        Additional keyword arguments are passed as keyword arguments to
+        ``map_funtion``.
+
+    Returns
+    -------
+    fig : `Figure`
+        A matplotlib figure with the maps on two axes.
+    """
+    map_factories = {
+        "laea": DEFAULT_SKY_MAP_FACTORY,
+        "mcbryde": partial(skyproj.McBrydeSkyproj, lon_0=lon_0),
+    }
+
+    fig, axes = plt.subplots(1, 2, figsize=figsize, gridspec_kw={"width_ratios": [1, 2]})
+    for plot_index, proj in enumerate(map_factories):
+        these_axes = axes[plot_index]
+        map_function(*args, axes=these_axes, sky_map_factory=map_factories[proj], **kwargs)
+    return fig
+
+
+def map_metric_in_laea_and_mcbryde(
+    visits: pd.DataFrame,
+    metric: maf.metrics.base_metric.BaseMetric,
+    map_function: Callable,
+    *args,
+    nside: int = 32,
+    **kwargs,
+) -> Figure:
+    """Map a healpix array in both LAEA and McBryde projections, side by side.
+
+    Parameters
+    ----------
+    visits: `pandas.DataFrame`
+        A table of visits with any columns necessary for computation of the
+        supplied ``metric``
+    metric: `rubin_sim.maf.metrics.base_metric.BaseMetric`
+        The MAF metric to compute.
+    map_function: `Callable`
+        A function that takes sky_map_function and axes arguments, and plots
+        the healpix map on these axes. The ``axes`` agument must take an
+        instance of `matplotlib.axes.Axes`, and ``sky_map_factory`` must be
+        a dictionary with keys ``laea`` and ``mcbryde`` with values that
+        are functions that return instances of `skyproj.LaeaSkyproj` and
+        `skyproj.McBrydeSkyproj`, respectively.
+    *args:
+        Additional positional arguments are passed as positional arguments
+        to ``map_function``.
+    nside: `int`
+        The nside of the healpix array on which to compute the metric.
+    **kwargs:
+        Additional keyword arguments are passed as keyword arguments to
+        ``map_funtion``.
+
+    Returns
+    -------
+    fig : `Figure`
+        A matplotlib figure with the maps on two axes.
+    """
+    metric_bundle = maf.MetricBundle(
+        metric,
+        maf.HealpixSlicer(nside=nside, verbose=False),
+    )
+    compute_metric(visits, metric_bundle)
+
+    fig = map_hpix_in_laea_and_mcbryde(map_function, metric_bundle.metric_values, *args, **kwargs)
+    return fig