From 127a4f94ae26c7983af927d8586184bcc356ba6d Mon Sep 17 00:00:00 2001 From: Eric Neilsen Date: Thu, 7 Mar 2024 09:34:42 -0800 Subject: [PATCH 1/8] initial create_hpix_visit_map_grid and supporting code --- schedview/plot/mplutils.py | 395 ++++++++++++++++++++++++++++++++++ schedview/plot/survey_mpl.py | 135 ++++++++++++ tests/test_plot_mplutils.py | 95 ++++++++ tests/test_plot_survey_mpl.py | 50 +++++ 4 files changed, 675 insertions(+) create mode 100644 schedview/plot/mplutils.py create mode 100644 schedview/plot/survey_mpl.py create mode 100644 tests/test_plot_mplutils.py create mode 100644 tests/test_plot_survey_mpl.py diff --git a/schedview/plot/mplutils.py b/schedview/plot/mplutils.py new file mode 100644 index 00000000..8826d08c --- /dev/null +++ b/schedview/plot/mplutils.py @@ -0,0 +1,395 @@ +import astropy +import astropy.units as u +import cartopy +import cartopy.crs as ccrs +import colorcet +import healpy as hp +import matplotlib as mpl +import numpy as np +import pandas as pd +from astropy.coordinates import SkyCoord + +DEFAULT_MERIDIANS = np.arange(-180, 210, 30.0) +DEFAULT_PARALLELS = np.arange(-60, 90.0, 30.0) +EARTH_DIAMETER = 2 * astropy.constants.R_earth.to(u.m).value +PC = ccrs.PlateCarree() +LAEA_SOUTH = ccrs.LambertAzimuthalEqualArea(central_longitude=0, central_latitude=-89.99) + + +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 sky_map( + axes, + parallels=DEFAULT_PARALLELS, + meridians=DEFAULT_MERIDIANS, + full_sky=True, + west_right=True, + label_parallels=True, + label_meridians=True, + south=True, + grid_color="gray", +): + """Turn a set of matplotlib axes into a sky map. + + Parameters + ---------- + axes : `cartopy.mpl.geoaxes.GeoAxes` + A matplotlib set of axes with a cartopy transform specifying the + projection. + parallels : `numpy.ndarray` + An array specifying the locations of graticule parallels, in + degrees. Defaults to every 30 degrees. + meridians : `numpy.ndarray` + An array specifying the locations of graticule meridians, in + degrees. Defaults to every 30 degrees. + full_sky : `bool` + Show the full sky, not just where there is data. Defaults to True. + west_right : `bool` + Show west to the right (the surface of a sphere as seen from the + inside, as in standing on the earth looking at the sky) rather than + to the left (the surface of a sphere as seen from the outside, as + in looking down on a the Earth). Defaults to True. + label_parallels : `bool` + Show labels on the parallel graticules. Defaults to True. + label_meridians : `bool` + Show labels on the meridian graticules. Defaults to True. + south : `bool` + If True, place the south equatorial pole at the center (if True). + If False place the north equatorial pole at the center. + Defaults to True. + grid_color : `str` + The color for the graticules. Defaults to `gray`. + """ + if full_sky: + # Cartopy was designed for geography, so it scales to units of + # meters on the Earth. + axes.set_xlim(-1 * EARTH_DIAMETER, EARTH_DIAMETER) + axes.set_ylim(-1 * EARTH_DIAMETER, EARTH_DIAMETER) + + if west_right: + axes.set_xlim(reversed(axes.get_xlim())) + + gl = axes.gridlines(crs=PC, draw_labels=False, color=grid_color) + gl.xlocator = mpl.ticker.FixedLocator(meridians) + gl.ylocator = mpl.ticker.FixedLocator(parallels) + gl.xformatter = cartopy.mpl.gridliner.LONGITUDE_FORMATTER + gl.yformatter = cartopy.mpl.gridliner.LATITUDE_FORMATTER + + meridian_label = np.array([0, 90, 180, 270]) + if label_meridians: + if south: + meridian_label_dec = 89 + m = meridian_label[0] + axes.text( + m, + meridian_label_dec, + r"$%d^\circ$" % m, + horizontalalignment="center", + verticalalignment="bottom", + color=grid_color, + transform=PC, + ) + m = meridian_label[1] + axes.text( + m, + meridian_label_dec, + r"$%d^\circ$" % m, + horizontalalignment="right", + verticalalignment="center", + color=grid_color, + transform=PC, + ) + m = meridian_label[2] + axes.text( + m, + meridian_label_dec, + r"$%d^\circ$" % m, + horizontalalignment="center", + verticalalignment="top", + color=grid_color, + transform=PC, + ) + m = meridian_label[3] + axes.text( + m, + meridian_label_dec, + r"$%d^\circ$" % m, + horizontalalignment="left", + verticalalignment="center", + color=grid_color, + transform=PC, + ) + else: + meridian_label_dec = -89 + m = meridian_label[0] + axes.text( + m, + meridian_label_dec, + r"$%d^\circ$" % m, + horizontalalignment="center", + verticalalignment="top", + color=grid_color, + transform=PC, + ) + m = meridian_label[1] + axes.text( + m, + meridian_label_dec, + r"$%d^\circ$" % m, + horizontalalignment="right", + verticalalignment="center", + color=grid_color, + transform=PC, + ) + m = meridian_label[2] + axes.text( + m, + meridian_label_dec, + r"$%d^\circ$" % m, + horizontalalignment="center", + verticalalignment="bottom", + color=grid_color, + transform=PC, + ) + m = meridian_label[3] + axes.text( + m, + meridian_label_dec, + r"$%d^\circ$" % m, + horizontalalignment="left", + verticalalignment="center", + color=grid_color, + transform=PC, + ) + + if label_parallels: + for p in parallels: + if p < 60: + axes.text( + 270 - 90, + p, + r"$%d^\circ$" % p, + horizontalalignment="left", + verticalalignment="top", + color=grid_color, + weight="bold", + transform=PC, + ) + + +def map_healpix(axes, hp_map, resolution=None, scale_limits=None, **in_kwargs): + """Show healpixels on a set of GeoAxes. + + Parameters + ---------- + axes : `cartopy.mpl.geoaxes.GeoAxes` + A matplotlib set of axes with a cartopy transform specifying the + projection. + hp_map : `numpy.ndarray` or `numpy.ma.MaskedArray` + An array of healpixes values. + resolution : `float` or `None` + The resolution onto which to sample healpixel points (deg/pixel). + None sets it according to the nside of the healpix map. + Defaults to None. + scale_limits : `list` [`float`] or `None` + The minimum and maximum values to map to color values. + None sets it dynamically using `astropy.visualization.ZScaleInterval` + Defaults to None. + **kwargs + Additional keyword arguments are passed to + `cartopy.mpl.geoaxes.GeoAxes.imshow`. + """ + kwargs = {"cmap": colorcet.cm.blues} + kwargs.update(in_kwargs) + + nside = hp.npix2nside(hp_map.size) + if resolution is None: + resolution = hp.max_pixrad(nside, degrees=True) / 4 + + ra, decl = np.meshgrid(np.arange(-180, 180, resolution), np.arange(-90, 90, resolution)) + + if scale_limits is None: + try: + unmasked_hpix = hp_map[~hp_map.mask] + except AttributeError: + # this array doesn't have a mask + unmasked_hpix = hp_map + + scale_limits = astropy.visualization.ZScaleInterval().get_limits(unmasked_hpix) + + im = hp_map[hp.ang2pix(nside, ra, decl, lonlat=True)] + + axes.imshow( + im, extent=(-180, 180, 90, -90), transform=PC, vmin=scale_limits[0], vmax=scale_limits[1], **kwargs + ) + + +def plot_ecliptic(axes, **kwargs): + """Plot the ecliptic plane on a set of GeoAxes. + + Parameters + ---------- + axes : `cartopy.mpl.geoaxes.GeoAxes` + A matplotlib set of axes with a cartopy transform specifying the + projection. + **kwargs + Additional keyword arguments are passed to + `cartopy.mpl.geoaxes.GeoAxes.plot`. + """ + 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) + axes.plot(points_on_ecliptic.ra, points_on_ecliptic.decl, **kwargs) + + +def plot_galactic_plane(axes, **kwargs): + """Plot the galactic plane on a set of GeoAxes. + + Parameters + ---------- + axes : `cartopy.mpl.geoaxes.GeoAxes` + A matplotlib set of axes with a cartopy transform specifying the + projection. + **kwargs + Additional keyword arguments are passed to + `cartopy.mpl.geoaxes.GeoAxes.plot`. + """ + 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) + axes.plot(points_on_galactic_plane.ra, points_on_galactic_plane.decl, **kwargs) + + +def plot_sun(axes, model_observatory, night_events, **kwargs): + """Plot the sun on a set of GeoAxes. + + Parameters + ---------- + axes : `cartopy.mpl.geoaxes.GeoAxes` + A matplotlib set of axes with a cartopy transform specifying the + projection. + model_observatory : `ModelObservatory` + The model observatory (which supplies the location and sun position). + night_events : `pd.DataFrame` + A table of almanac events for the night, as generated by + `schedview.compute.astro.night_events`. + **kwargs + Additional keyword arguments are passed to + `cartopy.mpl.geoaxes.GeoAxes.scatter`. + """ + 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()) + axes.scatter(sun_ra, sun_decl, **kwargs) + + +def plot_moon(axes, model_observatory, night_events, **kwargs): + """Plot the moon at the start, middle, and end of night + on a set of GeoAxes. + + Parameters + ---------- + axes : `cartopy.mpl.geoaxes.GeoAxes` + A matplotlib set of axes with a cartopy transform specifying the + projection. + model_observatory : `ModelObservatory` + The model observatory (which supplies the location and sun position). + night_events : `pd.DataFrame` + A table of almanac events for the night, as generated by + `schedview.compute.astro.night_events`. + **kwargs + Additional keyword arguments are passed to + `cartopy.mpl.geoaxes.GeoAxes.scatter`. + """ + mjd = night_events.loc[["sunset", "night_middle", "sunrise"], "MJD"] + sun_moon_positions = model_observatory.almanac.get_sun_moon_positions(mjd) + sun_ra = np.degrees(sun_moon_positions["moon_RA"]) + sun_decl = np.degrees(sun_moon_positions["moon_dec"]) + axes.scatter(sun_ra, sun_decl, **kwargs) + + +def plot_horizons(axes, night_events, latitude, zd=90, **kwargs): + """Plot the limits of what can be observed at a given zenith distance + on a night. + + Parameters + ---------- + axes : `cartopy.mpl.geoaxes.GeoAxes` + A matplotlib set of axes with a cartopy transform specifying the + projection. + night_events : `pd.DataFrame` + A table of almanac events for the night, as generated by + `schedview.compute.astro.night_events`. + latitute : `float` + The latitude of the observatory, in degrees. + zd : `float` + The highest observable zenith distance, in degrees. + Defaults to 90. + **kwargs + Additional keyword arguments are passed to + `cartopy.mpl.geoaxes.GeoAxes.scatter`. + """ + evening = compute_circle_points(night_events.loc["sun_n12_setting", "LST"], latitude, zd, 180, 380) + axes.plot(evening.ra, evening.decl, **kwargs) + + morning = compute_circle_points(night_events.loc["sun_n12_rising", "LST"], latitude, zd, 0, 180) + axes.plot(morning.ra, morning.decl, **kwargs) diff --git a/schedview/plot/survey_mpl.py b/schedview/plot/survey_mpl.py new file mode 100644 index 00000000..3649c6e7 --- /dev/null +++ b/schedview/plot/survey_mpl.py @@ -0,0 +1,135 @@ +import matplotlib as mpl +import matplotlib.pyplot as plt +import numpy as np + +# Import ModelObservatory to make sphinx happy +from rubin_scheduler.scheduler.model_observatory import ModelObservatory # noqa F401 + +from schedview.compute.camera import LsstCameraFootprintPerimeter + +from .mplutils import ( + LAEA_SOUTH, + PC, + map_healpix, + plot_ecliptic, + plot_galactic_plane, + plot_moon, + plot_sun, + sky_map, +) + + +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 : `cartopy.mpl.geoaxes.GeoAxes` + A matplotlib set of axes with a cartopy transform specifying the + projection. + **kwargs + Additional keyword arguments are passed to + `map_healpix`. + + Returns + ------- + axes : `cartopy.mpl.geoaxes.GeoAxes` + A matplotlib set of axes with a cartopy transform specifying the + projection, map, visits, and astronomical annotations. + """ + + if axes is None: + fig = plt.figure() + axes = fig.add_subplot(1, 1, 1, projection=LAEA_SOUTH) + + map_healpix(axes, map_hpix, **kwargs) + plot_ecliptic(axes, transform=PC, color="green") + plot_galactic_plane(axes, transform=PC, color="blue") + plot_sun(axes, model_observatory, night_events, transform=PC, color="yellow") + plot_moon(axes, model_observatory, night_events, transform=PC, color="orange") + + camera_perimeter = LsstCameraFootprintPerimeter() + ras, decls = camera_perimeter(visits.fieldRA, visits.fieldDec, visits.rotSkyPos) + + # Keep matplotlib & cartopy from going the wrong way around. + for these_ras in ras: + if these_ras.max() - these_ras.min() > 180: + low_ras = these_ras < 180 + these_ras[low_ras] = these_ras[low_ras] + 360 + + vertices = [np.stack(c).T for c in zip(ras, decls)] + pointing_collection = mpl.collections.PolyCollection(vertices) + pointing_collection.set(transform=PC, edgecolor="black", linewidth=1, facecolor="None") + + axes.add_collection(pointing_collection) + sky_map(axes) + return axes + + +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 + `map_healpix`. + + 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, projection=LAEA_SOUTH) + axes.set_title(band, loc="left") + map_visits_over_healpix( + visits_in_band, hpix_maps[band], model_observatory, night_events, axes=axes, **kwargs + ) + + return fig diff --git a/tests/test_plot_mplutils.py b/tests/test_plot_mplutils.py new file mode 100644 index 00000000..ad30a6f3 --- /dev/null +++ b/tests/test_plot_mplutils.py @@ -0,0 +1,95 @@ +import datetime + +import astropy +import astropy.units as u +import healpy as hp +import matplotlib.pyplot as plt +import numpy as np +from astropy.coordinates import SkyCoord +from rubin_scheduler.scheduler.model_observatory import ModelObservatory + +import schedview.compute.astro +from schedview.plot.mplutils import ( + LAEA_SOUTH, + compute_circle_points, + map_healpix, + plot_ecliptic, + plot_galactic_plane, + plot_horizons, + plot_moon, + plot_sun, + sky_map, +) + +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_sky_map(): + figure = plt.figure() + axes = figure.add_subplot(projection=LAEA_SOUTH) + sky_map(axes) + + +def test_map_healpix(): + figure = plt.figure() + axes = figure.add_subplot(projection=LAEA_SOUTH) + hp_map = np.random.default_rng(RANDOM_SEED).uniform(0, 1, hp.nside2npix(4)) + map_healpix(axes, hp_map) + + +def test_plot_ecliptic(): + figure = plt.figure() + axes = figure.add_subplot(projection=LAEA_SOUTH) + plot_ecliptic(axes) + + +def test_plot_galactic_plane(): + figure = plt.figure() + axes = figure.add_subplot(projection=LAEA_SOUTH) + plot_galactic_plane(axes) + + +def test_plot_sun(): + figure = plt.figure() + axes = figure.add_subplot(projection=LAEA_SOUTH) + model_observatory = ModelObservatory(init_load_length=1) + astropy.utils.iers.conf.iers_degraded_accuracy = "ignore" + night_events = schedview.compute.astro.night_events(datetime.date(2025, 1, 1)) + plot_sun(axes, model_observatory, night_events) + + +def test_plot_moon(): + figure = plt.figure() + axes = figure.add_subplot(projection=LAEA_SOUTH) + model_observatory = ModelObservatory(init_load_length=1) + astropy.utils.iers.conf.iers_degraded_accuracy = "ignore" + night_events = schedview.compute.astro.night_events(datetime.date(2025, 1, 1)) + plot_moon(axes, model_observatory, night_events) + + +def test_plot_horizons(): + figure = plt.figure() + axes = figure.add_subplot(projection=LAEA_SOUTH) + astropy.utils.iers.conf.iers_degraded_accuracy = "ignore" + night_events = schedview.compute.astro.night_events(datetime.date(2025, 1, 1)) + plot_horizons(axes, night_events, -32, 80) diff --git a/tests/test_plot_survey_mpl.py b/tests/test_plot_survey_mpl.py new file mode 100644 index 00000000..3022ea43 --- /dev/null +++ b/tests/test_plot_survey_mpl.py @@ -0,0 +1,50 @@ +import datetime + +import astropy +import healpy as hp +import numpy as np +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_mpl import create_hpix_visit_map_grid, map_visits_over_healpix + +RANDOM_SEED = 6563 + + +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) From c600897b6ff751516fed6e8b4ec8861a37176ca4 Mon Sep 17 00:00:00 2001 From: Eric Neilsen Date: Fri, 8 Mar 2024 14:31:46 -0800 Subject: [PATCH 2/8] add cartopy and matplotlib dependencies --- pyproject.toml | 2 ++ requirements.txt | 2 ++ 2 files changed, 4 insertions(+) diff --git a/pyproject.toml b/pyproject.toml index 9e878044..d61319cf 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -23,11 +23,13 @@ urls = {documentation = "https://schedview.lsst.io", repository = "https://githu dynamic = [ "version" ] dependencies = [ "astropy >= 5.3", + "cartopy", "colorcet", "bokeh >= 3.1.1", "healpy", "holoviews", "hvplot", + "matplotlib", "numpy", "pandas", "panel >= 1.1.0", diff --git a/requirements.txt b/requirements.txt index 57eef49a..f6581031 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,9 +1,11 @@ astropy bokeh +cartopy colorcet healpy holoviews hvplot +matplotlib numpy pandas panel From 1e0ca3a39351a17364e2042631da3ad6e52058d5 Mon Sep 17 00:00:00 2001 From: Eric Neilsen Date: Mon, 11 Mar 2024 09:03:13 -0700 Subject: [PATCH 3/8] skyproj implementation of visits-over-hpix grid --- pyproject.toml | 1 + requirements.txt | 1 + schedview/plot/survey_skyproj.py | 217 ++++++++++++++++++++++++++++++ tests/test_plot_survey_skyproj.py | 76 +++++++++++ 4 files changed, 295 insertions(+) create mode 100644 schedview/plot/survey_skyproj.py create mode 100644 tests/test_plot_survey_skyproj.py diff --git a/pyproject.toml b/pyproject.toml index d61319cf..cc1edea7 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -35,6 +35,7 @@ dependencies = [ "panel >= 1.1.0", "param", "pytz", + "skyproj", "rubin-scheduler", "lsst.resources", "uranography >= 1.1.0 ", diff --git a/requirements.txt b/requirements.txt index f6581031..7ec85d93 100644 --- a/requirements.txt +++ b/requirements.txt @@ -11,6 +11,7 @@ pandas panel param pytz +skyproj uranography rubin-scheduler rubin-sim diff --git a/schedview/plot/survey_skyproj.py b/schedview/plot/survey_skyproj.py new file mode 100644 index 00000000..a320b1c2 --- /dev/null +++ b/schedview/plot/survey_skyproj.py @@ -0,0 +1,217 @@ +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 with a cartopy transform specifying the + projection. + **kwargs + Additional keyword arguments are passed to + `skyproj.skyproj.LaeaSkyproj.draw_hpxmap`. + + Returns + ------- + axes : `cartopy.mpl.geoaxes.GeoAxes` + A matplotlib set of axes with a cartopy transform specifying the + projection, map, visits, and astronomical annotations. + """ + + 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, 380) + 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 axes + + +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") + map_visits_over_healpix( + visits_in_band, hpix_maps[band], model_observatory, night_events, axes=axes, **kwargs + ) + + return fig diff --git a/tests/test_plot_survey_skyproj.py b/tests/test_plot_survey_skyproj.py new file mode 100644 index 00000000..7b93fe04 --- /dev/null +++ b/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) From 00d139dc3a0a46acf93edabdbf2736e9ce90c27e Mon Sep 17 00:00:00 2001 From: Eric Neilsen Date: Mon, 11 Mar 2024 09:04:49 -0700 Subject: [PATCH 4/8] docstring fix --- schedview/plot/survey_skyproj.py | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/schedview/plot/survey_skyproj.py b/schedview/plot/survey_skyproj.py index a320b1c2..8da319cd 100644 --- a/schedview/plot/survey_skyproj.py +++ b/schedview/plot/survey_skyproj.py @@ -94,17 +94,15 @@ def map_visits_over_healpix(visits, map_hpix, model_observatory, night_events, a 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 with a cartopy transform specifying the - projection. + A matplotlib set of axes. **kwargs Additional keyword arguments are passed to `skyproj.skyproj.LaeaSkyproj.draw_hpxmap`. Returns ------- - axes : `cartopy.mpl.geoaxes.GeoAxes` - A matplotlib set of axes with a cartopy transform specifying the - projection, map, visits, and astronomical annotations. + axes : `matplotlib.axes._axes.Axes` + A matplotlib set of axes. """ if axes is None: From 0a7ce183c114de0adc7f0d30e6f38e9ccee0374e Mon Sep 17 00:00:00 2001 From: Eric Neilsen Date: Mon, 11 Mar 2024 09:24:38 -0700 Subject: [PATCH 5/8] remove cartopy-based hpix/visit map grid --- schedview/plot/mplutils.py | 395 ---------------------------------- schedview/plot/survey_mpl.py | 135 ------------ tests/test_plot_mplutils.py | 95 -------- tests/test_plot_survey_mpl.py | 50 ----- 4 files changed, 675 deletions(-) delete mode 100644 schedview/plot/mplutils.py delete mode 100644 schedview/plot/survey_mpl.py delete mode 100644 tests/test_plot_mplutils.py delete mode 100644 tests/test_plot_survey_mpl.py diff --git a/schedview/plot/mplutils.py b/schedview/plot/mplutils.py deleted file mode 100644 index 8826d08c..00000000 --- a/schedview/plot/mplutils.py +++ /dev/null @@ -1,395 +0,0 @@ -import astropy -import astropy.units as u -import cartopy -import cartopy.crs as ccrs -import colorcet -import healpy as hp -import matplotlib as mpl -import numpy as np -import pandas as pd -from astropy.coordinates import SkyCoord - -DEFAULT_MERIDIANS = np.arange(-180, 210, 30.0) -DEFAULT_PARALLELS = np.arange(-60, 90.0, 30.0) -EARTH_DIAMETER = 2 * astropy.constants.R_earth.to(u.m).value -PC = ccrs.PlateCarree() -LAEA_SOUTH = ccrs.LambertAzimuthalEqualArea(central_longitude=0, central_latitude=-89.99) - - -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 sky_map( - axes, - parallels=DEFAULT_PARALLELS, - meridians=DEFAULT_MERIDIANS, - full_sky=True, - west_right=True, - label_parallels=True, - label_meridians=True, - south=True, - grid_color="gray", -): - """Turn a set of matplotlib axes into a sky map. - - Parameters - ---------- - axes : `cartopy.mpl.geoaxes.GeoAxes` - A matplotlib set of axes with a cartopy transform specifying the - projection. - parallels : `numpy.ndarray` - An array specifying the locations of graticule parallels, in - degrees. Defaults to every 30 degrees. - meridians : `numpy.ndarray` - An array specifying the locations of graticule meridians, in - degrees. Defaults to every 30 degrees. - full_sky : `bool` - Show the full sky, not just where there is data. Defaults to True. - west_right : `bool` - Show west to the right (the surface of a sphere as seen from the - inside, as in standing on the earth looking at the sky) rather than - to the left (the surface of a sphere as seen from the outside, as - in looking down on a the Earth). Defaults to True. - label_parallels : `bool` - Show labels on the parallel graticules. Defaults to True. - label_meridians : `bool` - Show labels on the meridian graticules. Defaults to True. - south : `bool` - If True, place the south equatorial pole at the center (if True). - If False place the north equatorial pole at the center. - Defaults to True. - grid_color : `str` - The color for the graticules. Defaults to `gray`. - """ - if full_sky: - # Cartopy was designed for geography, so it scales to units of - # meters on the Earth. - axes.set_xlim(-1 * EARTH_DIAMETER, EARTH_DIAMETER) - axes.set_ylim(-1 * EARTH_DIAMETER, EARTH_DIAMETER) - - if west_right: - axes.set_xlim(reversed(axes.get_xlim())) - - gl = axes.gridlines(crs=PC, draw_labels=False, color=grid_color) - gl.xlocator = mpl.ticker.FixedLocator(meridians) - gl.ylocator = mpl.ticker.FixedLocator(parallels) - gl.xformatter = cartopy.mpl.gridliner.LONGITUDE_FORMATTER - gl.yformatter = cartopy.mpl.gridliner.LATITUDE_FORMATTER - - meridian_label = np.array([0, 90, 180, 270]) - if label_meridians: - if south: - meridian_label_dec = 89 - m = meridian_label[0] - axes.text( - m, - meridian_label_dec, - r"$%d^\circ$" % m, - horizontalalignment="center", - verticalalignment="bottom", - color=grid_color, - transform=PC, - ) - m = meridian_label[1] - axes.text( - m, - meridian_label_dec, - r"$%d^\circ$" % m, - horizontalalignment="right", - verticalalignment="center", - color=grid_color, - transform=PC, - ) - m = meridian_label[2] - axes.text( - m, - meridian_label_dec, - r"$%d^\circ$" % m, - horizontalalignment="center", - verticalalignment="top", - color=grid_color, - transform=PC, - ) - m = meridian_label[3] - axes.text( - m, - meridian_label_dec, - r"$%d^\circ$" % m, - horizontalalignment="left", - verticalalignment="center", - color=grid_color, - transform=PC, - ) - else: - meridian_label_dec = -89 - m = meridian_label[0] - axes.text( - m, - meridian_label_dec, - r"$%d^\circ$" % m, - horizontalalignment="center", - verticalalignment="top", - color=grid_color, - transform=PC, - ) - m = meridian_label[1] - axes.text( - m, - meridian_label_dec, - r"$%d^\circ$" % m, - horizontalalignment="right", - verticalalignment="center", - color=grid_color, - transform=PC, - ) - m = meridian_label[2] - axes.text( - m, - meridian_label_dec, - r"$%d^\circ$" % m, - horizontalalignment="center", - verticalalignment="bottom", - color=grid_color, - transform=PC, - ) - m = meridian_label[3] - axes.text( - m, - meridian_label_dec, - r"$%d^\circ$" % m, - horizontalalignment="left", - verticalalignment="center", - color=grid_color, - transform=PC, - ) - - if label_parallels: - for p in parallels: - if p < 60: - axes.text( - 270 - 90, - p, - r"$%d^\circ$" % p, - horizontalalignment="left", - verticalalignment="top", - color=grid_color, - weight="bold", - transform=PC, - ) - - -def map_healpix(axes, hp_map, resolution=None, scale_limits=None, **in_kwargs): - """Show healpixels on a set of GeoAxes. - - Parameters - ---------- - axes : `cartopy.mpl.geoaxes.GeoAxes` - A matplotlib set of axes with a cartopy transform specifying the - projection. - hp_map : `numpy.ndarray` or `numpy.ma.MaskedArray` - An array of healpixes values. - resolution : `float` or `None` - The resolution onto which to sample healpixel points (deg/pixel). - None sets it according to the nside of the healpix map. - Defaults to None. - scale_limits : `list` [`float`] or `None` - The minimum and maximum values to map to color values. - None sets it dynamically using `astropy.visualization.ZScaleInterval` - Defaults to None. - **kwargs - Additional keyword arguments are passed to - `cartopy.mpl.geoaxes.GeoAxes.imshow`. - """ - kwargs = {"cmap": colorcet.cm.blues} - kwargs.update(in_kwargs) - - nside = hp.npix2nside(hp_map.size) - if resolution is None: - resolution = hp.max_pixrad(nside, degrees=True) / 4 - - ra, decl = np.meshgrid(np.arange(-180, 180, resolution), np.arange(-90, 90, resolution)) - - if scale_limits is None: - try: - unmasked_hpix = hp_map[~hp_map.mask] - except AttributeError: - # this array doesn't have a mask - unmasked_hpix = hp_map - - scale_limits = astropy.visualization.ZScaleInterval().get_limits(unmasked_hpix) - - im = hp_map[hp.ang2pix(nside, ra, decl, lonlat=True)] - - axes.imshow( - im, extent=(-180, 180, 90, -90), transform=PC, vmin=scale_limits[0], vmax=scale_limits[1], **kwargs - ) - - -def plot_ecliptic(axes, **kwargs): - """Plot the ecliptic plane on a set of GeoAxes. - - Parameters - ---------- - axes : `cartopy.mpl.geoaxes.GeoAxes` - A matplotlib set of axes with a cartopy transform specifying the - projection. - **kwargs - Additional keyword arguments are passed to - `cartopy.mpl.geoaxes.GeoAxes.plot`. - """ - 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) - axes.plot(points_on_ecliptic.ra, points_on_ecliptic.decl, **kwargs) - - -def plot_galactic_plane(axes, **kwargs): - """Plot the galactic plane on a set of GeoAxes. - - Parameters - ---------- - axes : `cartopy.mpl.geoaxes.GeoAxes` - A matplotlib set of axes with a cartopy transform specifying the - projection. - **kwargs - Additional keyword arguments are passed to - `cartopy.mpl.geoaxes.GeoAxes.plot`. - """ - 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) - axes.plot(points_on_galactic_plane.ra, points_on_galactic_plane.decl, **kwargs) - - -def plot_sun(axes, model_observatory, night_events, **kwargs): - """Plot the sun on a set of GeoAxes. - - Parameters - ---------- - axes : `cartopy.mpl.geoaxes.GeoAxes` - A matplotlib set of axes with a cartopy transform specifying the - projection. - model_observatory : `ModelObservatory` - The model observatory (which supplies the location and sun position). - night_events : `pd.DataFrame` - A table of almanac events for the night, as generated by - `schedview.compute.astro.night_events`. - **kwargs - Additional keyword arguments are passed to - `cartopy.mpl.geoaxes.GeoAxes.scatter`. - """ - 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()) - axes.scatter(sun_ra, sun_decl, **kwargs) - - -def plot_moon(axes, model_observatory, night_events, **kwargs): - """Plot the moon at the start, middle, and end of night - on a set of GeoAxes. - - Parameters - ---------- - axes : `cartopy.mpl.geoaxes.GeoAxes` - A matplotlib set of axes with a cartopy transform specifying the - projection. - model_observatory : `ModelObservatory` - The model observatory (which supplies the location and sun position). - night_events : `pd.DataFrame` - A table of almanac events for the night, as generated by - `schedview.compute.astro.night_events`. - **kwargs - Additional keyword arguments are passed to - `cartopy.mpl.geoaxes.GeoAxes.scatter`. - """ - mjd = night_events.loc[["sunset", "night_middle", "sunrise"], "MJD"] - sun_moon_positions = model_observatory.almanac.get_sun_moon_positions(mjd) - sun_ra = np.degrees(sun_moon_positions["moon_RA"]) - sun_decl = np.degrees(sun_moon_positions["moon_dec"]) - axes.scatter(sun_ra, sun_decl, **kwargs) - - -def plot_horizons(axes, night_events, latitude, zd=90, **kwargs): - """Plot the limits of what can be observed at a given zenith distance - on a night. - - Parameters - ---------- - axes : `cartopy.mpl.geoaxes.GeoAxes` - A matplotlib set of axes with a cartopy transform specifying the - projection. - night_events : `pd.DataFrame` - A table of almanac events for the night, as generated by - `schedview.compute.astro.night_events`. - latitute : `float` - The latitude of the observatory, in degrees. - zd : `float` - The highest observable zenith distance, in degrees. - Defaults to 90. - **kwargs - Additional keyword arguments are passed to - `cartopy.mpl.geoaxes.GeoAxes.scatter`. - """ - evening = compute_circle_points(night_events.loc["sun_n12_setting", "LST"], latitude, zd, 180, 380) - axes.plot(evening.ra, evening.decl, **kwargs) - - morning = compute_circle_points(night_events.loc["sun_n12_rising", "LST"], latitude, zd, 0, 180) - axes.plot(morning.ra, morning.decl, **kwargs) diff --git a/schedview/plot/survey_mpl.py b/schedview/plot/survey_mpl.py deleted file mode 100644 index 3649c6e7..00000000 --- a/schedview/plot/survey_mpl.py +++ /dev/null @@ -1,135 +0,0 @@ -import matplotlib as mpl -import matplotlib.pyplot as plt -import numpy as np - -# Import ModelObservatory to make sphinx happy -from rubin_scheduler.scheduler.model_observatory import ModelObservatory # noqa F401 - -from schedview.compute.camera import LsstCameraFootprintPerimeter - -from .mplutils import ( - LAEA_SOUTH, - PC, - map_healpix, - plot_ecliptic, - plot_galactic_plane, - plot_moon, - plot_sun, - sky_map, -) - - -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 : `cartopy.mpl.geoaxes.GeoAxes` - A matplotlib set of axes with a cartopy transform specifying the - projection. - **kwargs - Additional keyword arguments are passed to - `map_healpix`. - - Returns - ------- - axes : `cartopy.mpl.geoaxes.GeoAxes` - A matplotlib set of axes with a cartopy transform specifying the - projection, map, visits, and astronomical annotations. - """ - - if axes is None: - fig = plt.figure() - axes = fig.add_subplot(1, 1, 1, projection=LAEA_SOUTH) - - map_healpix(axes, map_hpix, **kwargs) - plot_ecliptic(axes, transform=PC, color="green") - plot_galactic_plane(axes, transform=PC, color="blue") - plot_sun(axes, model_observatory, night_events, transform=PC, color="yellow") - plot_moon(axes, model_observatory, night_events, transform=PC, color="orange") - - camera_perimeter = LsstCameraFootprintPerimeter() - ras, decls = camera_perimeter(visits.fieldRA, visits.fieldDec, visits.rotSkyPos) - - # Keep matplotlib & cartopy from going the wrong way around. - for these_ras in ras: - if these_ras.max() - these_ras.min() > 180: - low_ras = these_ras < 180 - these_ras[low_ras] = these_ras[low_ras] + 360 - - vertices = [np.stack(c).T for c in zip(ras, decls)] - pointing_collection = mpl.collections.PolyCollection(vertices) - pointing_collection.set(transform=PC, edgecolor="black", linewidth=1, facecolor="None") - - axes.add_collection(pointing_collection) - sky_map(axes) - return axes - - -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 - `map_healpix`. - - 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, projection=LAEA_SOUTH) - axes.set_title(band, loc="left") - map_visits_over_healpix( - visits_in_band, hpix_maps[band], model_observatory, night_events, axes=axes, **kwargs - ) - - return fig diff --git a/tests/test_plot_mplutils.py b/tests/test_plot_mplutils.py deleted file mode 100644 index ad30a6f3..00000000 --- a/tests/test_plot_mplutils.py +++ /dev/null @@ -1,95 +0,0 @@ -import datetime - -import astropy -import astropy.units as u -import healpy as hp -import matplotlib.pyplot as plt -import numpy as np -from astropy.coordinates import SkyCoord -from rubin_scheduler.scheduler.model_observatory import ModelObservatory - -import schedview.compute.astro -from schedview.plot.mplutils import ( - LAEA_SOUTH, - compute_circle_points, - map_healpix, - plot_ecliptic, - plot_galactic_plane, - plot_horizons, - plot_moon, - plot_sun, - sky_map, -) - -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_sky_map(): - figure = plt.figure() - axes = figure.add_subplot(projection=LAEA_SOUTH) - sky_map(axes) - - -def test_map_healpix(): - figure = plt.figure() - axes = figure.add_subplot(projection=LAEA_SOUTH) - hp_map = np.random.default_rng(RANDOM_SEED).uniform(0, 1, hp.nside2npix(4)) - map_healpix(axes, hp_map) - - -def test_plot_ecliptic(): - figure = plt.figure() - axes = figure.add_subplot(projection=LAEA_SOUTH) - plot_ecliptic(axes) - - -def test_plot_galactic_plane(): - figure = plt.figure() - axes = figure.add_subplot(projection=LAEA_SOUTH) - plot_galactic_plane(axes) - - -def test_plot_sun(): - figure = plt.figure() - axes = figure.add_subplot(projection=LAEA_SOUTH) - model_observatory = ModelObservatory(init_load_length=1) - astropy.utils.iers.conf.iers_degraded_accuracy = "ignore" - night_events = schedview.compute.astro.night_events(datetime.date(2025, 1, 1)) - plot_sun(axes, model_observatory, night_events) - - -def test_plot_moon(): - figure = plt.figure() - axes = figure.add_subplot(projection=LAEA_SOUTH) - model_observatory = ModelObservatory(init_load_length=1) - astropy.utils.iers.conf.iers_degraded_accuracy = "ignore" - night_events = schedview.compute.astro.night_events(datetime.date(2025, 1, 1)) - plot_moon(axes, model_observatory, night_events) - - -def test_plot_horizons(): - figure = plt.figure() - axes = figure.add_subplot(projection=LAEA_SOUTH) - astropy.utils.iers.conf.iers_degraded_accuracy = "ignore" - night_events = schedview.compute.astro.night_events(datetime.date(2025, 1, 1)) - plot_horizons(axes, night_events, -32, 80) diff --git a/tests/test_plot_survey_mpl.py b/tests/test_plot_survey_mpl.py deleted file mode 100644 index 3022ea43..00000000 --- a/tests/test_plot_survey_mpl.py +++ /dev/null @@ -1,50 +0,0 @@ -import datetime - -import astropy -import healpy as hp -import numpy as np -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_mpl import create_hpix_visit_map_grid, map_visits_over_healpix - -RANDOM_SEED = 6563 - - -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) From 3f5fe1265c616101e587f242080eb27c610575e5 Mon Sep 17 00:00:00 2001 From: Eric Neilsen Date: Mon, 11 Mar 2024 09:24:54 -0700 Subject: [PATCH 6/8] remove cartopy dependency --- pyproject.toml | 1 - requirements.txt | 1 - 2 files changed, 2 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index cc1edea7..ea0f55eb 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -23,7 +23,6 @@ urls = {documentation = "https://schedview.lsst.io", repository = "https://githu dynamic = [ "version" ] dependencies = [ "astropy >= 5.3", - "cartopy", "colorcet", "bokeh >= 3.1.1", "healpy", diff --git a/requirements.txt b/requirements.txt index 7ec85d93..01d7b0de 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,6 +1,5 @@ astropy bokeh -cartopy colorcet healpy holoviews From b0ff7db5c79a80c8d7f2cfe343e245f55278a3b6 Mon Sep 17 00:00:00 2001 From: Eric Neilsen Date: Mon, 11 Mar 2024 12:15:33 -0700 Subject: [PATCH 7/8] give the hpix-visit maps their titles back --- schedview/plot/survey_skyproj.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/schedview/plot/survey_skyproj.py b/schedview/plot/survey_skyproj.py index 8da319cd..09792258 100644 --- a/schedview/plot/survey_skyproj.py +++ b/schedview/plot/survey_skyproj.py @@ -155,7 +155,7 @@ def map_visits_over_healpix(visits, map_hpix, model_observatory, night_events, a for visit_ras, visit_decls in zip(ras, decls): sky_map.draw_polygon(visit_ras, visit_decls, edgecolor="black", linewidth=0.2) - return axes + return sky_map.ax def create_hpix_visit_map_grid( @@ -208,8 +208,9 @@ def create_hpix_visit_map_grid( 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") - map_visits_over_healpix( + 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 From db2d7f6c4053ff6688a4d113c7119f754bed12af Mon Sep 17 00:00:00 2001 From: Eric Neilsen Date: Mon, 11 Mar 2024 12:17:27 -0700 Subject: [PATCH 8/8] fix dotted-dashed bondary in visit-hpix map --- schedview/plot/survey_skyproj.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/schedview/plot/survey_skyproj.py b/schedview/plot/survey_skyproj.py index 09792258..7a58e085 100644 --- a/schedview/plot/survey_skyproj.py +++ b/schedview/plot/survey_skyproj.py @@ -139,7 +139,7 @@ def map_visits_over_healpix(visits, map_hpix, model_observatory, night_events, a # Night limit horizons latitude = model_observatory.site.latitude zd = 70 - evening = compute_circle_points(night_events.loc["sun_n12_setting", "LST"], latitude, zd, 180, 380) + 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")