Merge pull request #77 from lsst/tickets/PREOPS-4975

tickets/PREOPS-4975: provide matplotlib implementation of hpix and visit maps

pyproject.toml

@@ -28,11 +28,13 @@ dependencies = [
     "healpy",
     "holoviews",
     "hvplot",
+    "matplotlib",
     "numpy",
     "pandas",
     "panel >= 1.1.0",
     "param",
     "pytz",
+    "skyproj",
     "rubin-scheduler",
     "lsst.resources",
     "uranography >= 1.1.0 ",

requirements.txt

@@ -4,11 +4,13 @@ colorcet
 healpy
 holoviews
 hvplot
+matplotlib
 numpy
 pandas
 panel
 param
 pytz
+skyproj
 uranography
 rubin-scheduler
 rubin-sim

schedview/plot/survey_skyproj.py

@@ -0,0 +1,216 @@
+import astropy.units as u
+import colorcet
+import matplotlib.pyplot as plt
+import numpy as np
+import pandas as pd
+import skyproj
+from astropy.coordinates import SkyCoord
+
+# Import ModelObservatory to make sphinx happy
+from rubin_scheduler.scheduler.model_observatory import ModelObservatory  # noqa F401
+
+from schedview.compute.camera import LsstCameraFootprintPerimeter
+
+DEFAULT_MAP_KWARGS = {"cmap": colorcet.cm.blues}
+
+
+def compute_circle_points(
+    center_ra,
+    center_decl,
+    radius=90.0,
+    start_bear=0,
+    end_bear=360,
+    step=1,
+):
+    """Create points along a circle or arc on a sphere
+
+    Parameters
+    ----------
+    center_ra : `float`
+        R.A. of the center of the circle (deg.).
+    center_decl : `float`
+        Decl. of the center of the circle (deg.).
+    radius : float, optional
+        Radius of the circle (deg.), by default 90.0
+    start_bear : int, optional
+        Bearing (E. of N.) of the start of the circle (deg.), by default 0
+    end_bear : int, optional
+        Bearing (E. of N.) of the end of the circle (deg.), by default 360
+    step : int, optional
+        Spacing of the points along the circle (deg.), by default 1
+
+    Returns
+    -------
+    circle : `pandas.DataFrame`
+        DataFrame with points in the circle.
+    """
+    ras = []
+    decls = []
+
+    bearing_angles = np.arange(start_bear, end_bear + step, step) * u.deg
+    center_coords = SkyCoord(center_ra * u.deg, center_decl * u.deg)
+    circle_coords = center_coords.directional_offset_by(bearing_angles, radius * u.deg)
+    bearings = bearing_angles.value.tolist()
+    ras = circle_coords.ra.deg.tolist()
+    decls = circle_coords.dec.deg.tolist()
+
+    # de-normalize RA so there are no large jumps, which can
+    # confuse cartopy + matplotlib
+
+    previous_ra = ras[0]
+    for ra_index in range(1, len(ras)):
+        if ras[ra_index] - previous_ra > 180:
+            ras[ra_index] -= 360
+        if ras[ra_index] - previous_ra < -180:
+            ras[ra_index] += 360
+        previous_ra = ras[ra_index]
+
+    circle = pd.DataFrame(
+        {
+            "bearing": bearings,
+            "ra": ras,
+            "decl": decls,
+        }
+    ).set_index("bearing")
+
+    return circle
+
+
+def map_visits_over_healpix(visits, map_hpix, model_observatory, night_events, axes=None, **kwargs):
+    """Plots visits over a healpix map, with astronomical annotations.
+
+    Parameters
+    ---------
+    visits : `pd.DataFrame`
+        A DataFrame of visits, with columns ``fieldRA`` and ``fieldDec``,
+        with the coordinates (in degrees), ``observationStartMJD`` with the
+        MJD of the start of each visit, and ``filter`` with the band
+        (as a string).
+    map_hpix : `numpy.ndarray`
+        The healpixel array to show.
+    model_observatory : `ModelObservatory`
+        The model observatory.
+    night_events : `pandas.DataFrame`
+        A table of almanac events for the night, as generated by
+        `schedview.compute.astro.night_events`.
+    axes : `matplotlib.axes._axes.Axes`
+        A matplotlib set of axes.
+    **kwargs
+        Additional keyword arguments are passed to
+        `skyproj.skyproj.LaeaSkyproj.draw_hpxmap`.
+
+    Returns
+    -------
+    axes : `matplotlib.axes._axes.Axes`
+        A matplotlib set of axes.
+    """
+
+    if axes is None:
+        fig = plt.figure()
+        axes = fig.add_subplot(1, 1, 1)
+
+    sky_map = skyproj.LaeaSkyproj(ax=axes, lat_0=-np.degrees(np.arccos(np.finfo(float).resolution)), lon_0=0)
+
+    hpxmap_kwargs = {"cmap": colorcet.cm.blues}
+    hpxmap_kwargs.update(kwargs)
+    sky_map.draw_hpxmap(map_hpix, zoom=False, **hpxmap_kwargs)
+
+    # Ecliptic
+    ecliptic_pole = SkyCoord(lon=0 * u.degree, lat=90 * u.degree, frame="geocentricmeanecliptic").icrs
+    points_on_ecliptic = compute_circle_points(ecliptic_pole.ra.deg, ecliptic_pole.dec.deg)
+    sky_map.draw_polygon(points_on_ecliptic.ra, points_on_ecliptic.decl, edgecolor="green")
+
+    # Galactic plane
+    galactic_pole = SkyCoord(l=0 * u.degree, b=90 * u.degree, frame="galactic").icrs
+    points_on_galactic_plane = compute_circle_points(galactic_pole.ra.deg, galactic_pole.dec.deg)
+    sky_map.draw_polygon(points_on_galactic_plane.ra, points_on_galactic_plane.decl, edgecolor="blue")
+
+    # Sun and moon
+    mjd = night_events.loc["night_middle", "MJD"]
+    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.scatter(sun_ra, sun_decl, color="yellow")
+
+    moon_ra = np.degrees(sun_moon_positions["moon_RA"])
+    moon_decl = np.degrees(sun_moon_positions["moon_dec"])
+    sky_map.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.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.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.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.plot(morning.ra, morning.decl, color="red", linestyle="dotted")
+
+    camera_perimeter = LsstCameraFootprintPerimeter()
+    ras, decls = camera_perimeter(visits.fieldRA, visits.fieldDec, visits.rotSkyPos)
+
+    for visit_ras, visit_decls in zip(ras, decls):
+        sky_map.draw_polygon(visit_ras, visit_decls, edgecolor="black", linewidth=0.2)
+
+    return sky_map.ax
+
+
+def create_hpix_visit_map_grid(
+    visits, hpix_maps, model_observatory, night_events, fig=None, num_rows=2, **kwargs
+):
+    """Plot an array of visits over a healpix maps, with astronomical
+    annotations.
+
+    Parameters
+    ---------
+    visits : `pd.DataFrame`
+        A DataFrame of visits, with columns ``fieldRA`` and ``fieldDec``,
+        with the coordinates (in degrees), ``observationStartMJD`` with the
+        MJD of the start of each visit, and ``filter`` with the band
+        (as a string).
+    hpix_maps : `dict` [`str`, `numpy.ndarray`]
+        The healpixel array to show.
+    model_observatory : `ModelObservatory`
+        The model observatory.
+    night_events : `pandas.DataFrame`
+        A table of almanac events for the night, as generated by
+        `schedview.compute.astro.night_events`.
+    fig : `matplotlib.figure.Figure` or `None`
+        A matplotlib figure in which to plot the array of maps.
+        If None, a new figure is created.
+        Defaults to None.
+    nrows : `int`
+        The number of rows in the array of plots.
+        Defaults to 2.
+    **kwargs
+        Additional keyword arguments are passed to
+        `skyproj.skyproj.LaeaSkyproj.draw_hpxmap`.
+
+    Returns
+    -------
+    fig : `matplotlib.figure.Figure`
+        The matplotlib figure with the array of maps.
+    """
+
+    # Python's "//" operator rounds down, but we want to round up, so
+    # apply it to the negative, then negate the result.
+    num_columns = -(-len(hpix_maps) // num_rows)
+
+    if fig is None:
+        plot_height = 5 * num_rows
+        plot_width = 5.1 * num_columns
+        fig = plt.figure(figsize=(plot_width, plot_height))
+
+    for band_idx, band in enumerate(hpix_maps.keys()):
+        visits_in_band = visits.query(f"filter == '{band}'")
+        axes = fig.add_subplot(num_rows, num_columns, band_idx + 1)
+        axes.set_title(band, loc="left")
+        new_axes = map_visits_over_healpix(
+            visits_in_band, hpix_maps[band], model_observatory, night_events, axes=axes, **kwargs
+        )
+        new_axes.set_title(band, loc="left")
+
+    return fig

tests/test_plot_survey_skyproj.py

@@ -0,0 +1,76 @@
+import datetime
+
+import astropy
+import astropy.units as u
+import healpy as hp
+import numpy as np
+from astropy.coordinates import SkyCoord
+from astropy.time import Time
+from lsst.resources import ResourcePath
+from rubin_scheduler.scheduler.model_observatory import ModelObservatory
+
+import schedview
+import schedview.collect
+from schedview.plot.survey_skyproj import (
+    compute_circle_points,
+    create_hpix_visit_map_grid,
+    map_visits_over_healpix,
+)
+
+RANDOM_SEED = 6563
+
+
+def test_compute_circle_points():
+    center_ra, center_decl = 32.2, 10.3
+    radius = 12.2
+    circle_points = compute_circle_points(
+        center_ra,
+        center_decl,
+        radius=radius,
+        start_bear=0,
+        end_bear=360,
+        step=1,
+    )
+
+    center_coords = SkyCoord(center_ra, center_decl, unit=u.deg)
+    circle_coords = SkyCoord(circle_points.ra, circle_points.decl, unit=u.deg)
+    angles = center_coords.separation(circle_coords)
+    np.testing.assert_almost_equal(angles.deg, radius)
+    bearings = center_coords.position_angle(circle_coords)
+    np.testing.assert_almost_equal(circle_points.reset_index().bearing, bearings.deg)
+
+
+def test_map_visits_over_healpix():
+    hp_map = np.random.default_rng(RANDOM_SEED).uniform(0, 1, hp.nside2npix(4))
+
+    visits_path = ResourcePath("resource://schedview/data/sample_opsim.db")
+    visits = schedview.collect.read_opsim(visits_path)
+    visits_mjd = visits["observationStartMJD"].median()
+    time_datetime = Time(visits_mjd - 0.5, format="mjd").datetime
+
+    model_observatory = ModelObservatory(init_load_length=1)
+    astropy.utils.iers.conf.iers_degraded_accuracy = "ignore"
+    night_events = schedview.compute.astro.night_events(
+        datetime.date(time_datetime.year, time_datetime.month, time_datetime.day)
+    )
+
+    map_visits_over_healpix(visits, hp_map, model_observatory, night_events)
+
+
+def test_create_hpix_visit_map_grid():
+    hpix_maps = {}
+    for band in "ugrizy":
+        hpix_maps[band] = np.random.default_rng(RANDOM_SEED).uniform(0, 1, hp.nside2npix(4))
+
+    visits_path = ResourcePath("resource://schedview/data/sample_opsim.db")
+    visits = schedview.collect.read_opsim(visits_path)
+    visits_mjd = visits["observationStartMJD"].median()
+    time_datetime = Time(visits_mjd - 0.5, format="mjd").datetime
+
+    model_observatory = ModelObservatory(init_load_length=1)
+    astropy.utils.iers.conf.iers_degraded_accuracy = "ignore"
+    night_events = schedview.compute.astro.night_events(
+        datetime.date(time_datetime.year, time_datetime.month, time_datetime.day)
+    )
+
+    create_hpix_visit_map_grid(visits, hpix_maps, model_observatory, night_events)