tickets/SP-3294: add a maf stacker to add rise and set times for each field

rubin_sim/maf/stackers/__init__.py

@@ -8,6 +8,7 @@
 from .mo_stackers import *
 from .n_follow_stacker import *
 from .neo_dist_stacker import *
+from .riseset_stacker import *
 from .sdss_stackers import *
 from .sn_stacker import *
 from .teff_stacker import *

rubin_sim/maf/stackers/riseset_stacker.py

@@ -0,0 +1,161 @@
+__all__ = ("compute_gen_oblique_ascension", "RiseSetStacker")
+
+import numpy as np
+from rubin_scheduler.utils import Site, calc_lmst
+
+from .base_stacker import BaseStacker
+
+SIDEREAL_DAY = 0.9972696  # solar days per sidereal day
+
+
+def compute_gen_oblique_ascension(ra, dec, lat, alt):
+    """Compute the generalized oblique ascension.
+
+    The generalized oblique ascension is the local mean sidereal time at which
+    a point with coordinates (ra, dec) crosses altitude alt at a site with
+    geographic latitude lat.  It equals ra minus the ascensional difference at
+    that altitude.
+
+    Parameters
+    ----------
+    ra : `float` or `numpy.ndarray`
+        Right ascension, in radians.
+    dec : `float` or `numpy.ndarray`
+        Declination, in radians.
+    lat : `float`
+        Geographic latitude of the observer, in radians.
+    alt : `float`
+        Altitude of the crossing, in radians.
+
+    Returns
+    -------
+    oblique_ascension : `float` or `numpy.ndarray`
+        Generalized oblique ascension in radians.  NaN where the crossing
+        altitude is never reached (object is circumpolar with respect to alt,
+        or never rises to alt).
+
+    Notes
+    -----
+    This computes a "generalized" oblique ascension: the traditional oblique
+    ascension corresponds to the special case alt = 0 (rising or setting on
+    the horizon).  Right ascension is a special case of oblique ascension
+    where latitude is zero: at the equator every point rises and sets at
+    HA = +/-90 deg, so the LMST at rising is always RA - 90 deg, i.e. the
+    zero-latitude oblique ascension is a fixed offset from RA.
+
+    The key intermediate quantity is the ascensional difference, the hour
+    angle at which the point crosses altitude alt:
+
+        cos(ascensional_difference) =
+            (sin(alt) - sin(dec) * sin(lat)) / (cos(dec) * cos(lat))
+
+    When alt = 0 this reduces to the classical formula
+    cos(D) = -tan(dec) * tan(lat).
+    """
+    cos_asc_diff = (np.sin(alt) - np.sin(dec) * np.sin(lat)) / (np.cos(dec) * np.cos(lat))
+
+    # Values outside [-1, 1] mean the altitude is never reached.
+    normal = (cos_asc_diff >= -1.0) & (cos_asc_diff <= 1.0)
+    ascensional_difference = np.where(normal, np.arccos(np.clip(cos_asc_diff, -1.0, 1.0)), np.nan)
+
+    return ra - ascensional_difference
+
+
+class RiseSetStacker(BaseStacker):
+    """Add rise_mjd and set_mjd columns for each observation.
+
+    rise_mjd is the most recent MJD before observationStartMJD at which the
+    field rose above alt_limit (from below).
+
+    set_mjd is the next MJD after observationStartMJD at which the field will
+    fall below alt_limit.
+
+    Both are NaN for fields that are circumpolar relative to alt_limit
+    (always above) or that never rise above alt_limit.
+
+    The calculation is purely geometric: it ignores refraction and uses
+    sidereal (not solar) time for rise/set crossing times.
+
+    Parameters
+    ----------
+    mjd_col : `str`, optional
+        Column name for observation start MJD. Default 'observationStartMJD'.
+    ra_col : `str`, optional
+        Column name for RA. Default 'fieldRA'.
+    dec_col : `str`, optional
+        Column name for Dec. Default 'fieldDec'.
+    degrees : `bool`, optional
+        If True, ra/dec are in degrees. Default True.
+    site : `str` or `rubin_scheduler.utils.Site`, optional
+        Observatory name or Site object. Default 'LSST'.
+    alt_limit : `float`, optional
+        Limiting altitude in degrees. Default 20.0.
+    """
+
+    cols_added = ["rise_mjd", "set_mjd"]
+
+    def __init__(
+        self,
+        mjd_col="observationStartMJD",
+        ra_col="fieldRA",
+        dec_col="fieldDec",
+        degrees=True,
+        site="LSST",
+        alt_limit=20.0,
+    ):
+        self.mjd_col = mjd_col
+        self.ra_col = ra_col
+        self.dec_col = dec_col
+        self.degrees = degrees
+        self.alt_limit = alt_limit
+        self.cols_req = [mjd_col, ra_col, dec_col]
+        self.units = ["MJD", "MJD"]
+        self.cols_added_dtypes = [float, float]
+
+        if isinstance(site, str):
+            self.site = Site(name=site)
+        else:
+            self.site = site
+
+    def _run(self, sim_data, cols_present=False):
+        if cols_present:
+            return sim_data
+
+        mjd = sim_data[self.mjd_col]
+        ra = sim_data[self.ra_col]
+        dec = sim_data[self.dec_col]
+
+        if self.degrees:
+            ra = np.radians(ra)
+            dec = np.radians(dec)
+
+        lat = self.site.latitude_rad
+        lon = self.site.longitude_rad
+        alt = np.radians(self.alt_limit)
+
+        # Compute LMST (returned in hours by calc_lmst) then convert to
+        # radians.
+        lmst_rad = calc_lmst(mjd, lon) / 12.0 * np.pi
+
+        # Oblique ascension at rise (LMST when the field crosses alt from
+        # below) and at set (LMST when it crosses from above).  The
+        # ascensional difference embedded in compute_gen_oblique_ascension
+        # is the HA offset from the meridian to the crossing; negating it
+        # gives the set LMST.
+        oa_rise = compute_gen_oblique_ascension(ra, dec, lat, alt)
+        oa_set = 2.0 * ra - oa_rise  # ra + ascensional_difference
+
+        normal = np.isfinite(oa_rise)
+
+        # HA elapsed since last rise / remaining until next set,
+        # both in [0, 2*pi).
+        ha_since_rise = (lmst_rad - oa_rise) % (2.0 * np.pi)
+        ha_until_set = (oa_set - lmst_rad) % (2.0 * np.pi)
+
+        rise_mjd = mjd - ha_since_rise / (2.0 * np.pi) * SIDEREAL_DAY
+        set_mjd = mjd + ha_until_set / (2.0 * np.pi) * SIDEREAL_DAY
+
+        sim_data["rise_mjd"] = np.where(normal, rise_mjd, np.nan)
+        sim_data["set_mjd"] = np.where(normal, set_mjd, np.nan)
+
+        return sim_data

tests/maf/test_stackers.py

@@ -6,7 +6,7 @@
 import numpy as np
 import pandas as pd
 from astropy import units as u
-from astropy.coordinates import SkyCoord
+from astropy.coordinates import AltAz, EarthLocation, SkyCoord
 from astropy.time import Time
 from rubin_scheduler.utils import Site, _alt_az_pa_from_ra_dec, calc_lmst
 
@@ -419,5 +419,128 @@ def testHelp(self):
             stackers.BaseStacker.help(doc=True)
 
 
+class TestRiseSetStacker(unittest.TestCase):
+    """Tests for RiseSetStacker."""
+
+    def setUp(self):
+        self.alt_limit = 20.0
+        self.stacker = stackers.RiseSetStacker(alt_limit=self.alt_limit)
+        site = self.stacker.site
+        self.location = EarthLocation(
+            lat=site.latitude * u.deg,
+            lon=site.longitude * u.deg,
+            height=site.height * u.m,
+        )
+
+    def _get_alt(self, ra_deg, dec_deg, mjd):
+        """Return altitude in degrees from astropy.
+
+        pressure=0 disables the atmospheric refraction correction, matching
+        the stacker's purely geometric calculation.
+        """
+        coord = SkyCoord(ra=ra_deg * u.deg, dec=dec_deg * u.deg)
+        frame = AltAz(
+            obstime=Time(mjd, format="mjd", scale="utc"),
+            location=self.location,
+            pressure=0 * u.hPa,
+        )
+        return coord.transform_to(frame).alt.deg
+
+    def test_rise_set_above_horizon(self):
+        """rise_mjd/set_mjd bracket the observation and match alt_limit.
+
+        For several fields that are above alt_limit at the time of the
+        (simulated) observation, verify that:
+          - rise_mjd <= observationStartMJD < set_mjd
+          - The altitude at rise_mjd and set_mjd is alt_limit within 0.5 deg.
+
+        0.5 deg in altitude corresponds to roughly 2-3 minutes of time near the
+        horizon for these fields, comfortably within the 1-minute requirement.
+
+        All fields are tested together in a single stacker.run() call to
+        exercise the vectorised code path.
+        """
+        # (RA deg, Dec deg) pairs observable from Rubin (lat ~-30.24 deg)
+        # with a normal rise and set above 20 deg.  Fields must have
+        # -79.8 deg < Dec < +39.8 deg to avoid circumpolar/never-rises.
+        test_cases = [
+            (60.0, -40.0),  # moderate southern declination
+            (180.0, -25.0),  # near-equatorial
+            (300.0, -55.0),  # deep southern sky
+        ]
+        mjd_base = 51545.0  # near J2000; precession negligible at this epoch
+
+        # For each field, scan forward in steps of ~7 min to find an MJD where
+        # the field is at least 1 deg above alt_limit (so we are not right at a
+        # transition boundary).
+        dtype = [("observationStartMJD", float), ("fieldRA", float), ("fieldDec", float)]
+        rows = []
+        for ra_deg, dec_deg in test_cases:
+            for delta in np.arange(0.0, 1.0, 0.005):
+                if self._get_alt(ra_deg, dec_deg, mjd_base + delta) > self.alt_limit + 1.0:
+                    rows.append((mjd_base + delta, ra_deg, dec_deg))
+                    break
+        self.assertEqual(
+            len(rows),
+            len(test_cases),
+            "Could not find an above-horizon time for every test field",
+        )
+
+        sim_data = np.array(rows, dtype=dtype)
+        result = self.stacker.run(sim_data)
+
+        alt_tol = 0.5  # degrees; corresponds to ~2-3 min near the horizon
+
+        for i, (test_mjd, ra_deg, dec_deg) in enumerate(rows):
+            with self.subTest(ra=ra_deg, dec=dec_deg):
+                rise_mjd = result["rise_mjd"][i]
+                set_mjd = result["set_mjd"][i]
+
+                # Results must not be NaN for a field above the horizon.
+                self.assertFalse(np.isnan(rise_mjd), "rise_mjd should not be NaN")
+                self.assertFalse(np.isnan(set_mjd), "set_mjd should not be NaN")
+
+                # The observation must fall within the [rise, set) window.
+                self.assertLessEqual(rise_mjd, test_mjd)
+                self.assertGreater(set_mjd, test_mjd)
+
+                # The altitude at each crossing must equal alt_limit.
+                self.assertAlmostEqual(
+                    self._get_alt(ra_deg, dec_deg, rise_mjd),
+                    self.alt_limit,
+                    delta=alt_tol,
+                    msg=f"Altitude at rise_mjd wrong for RA={ra_deg}, Dec={dec_deg}",
+                )
+                self.assertAlmostEqual(
+                    self._get_alt(ra_deg, dec_deg, set_mjd),
+                    self.alt_limit,
+                    delta=alt_tol,
+                    msg=f"Altitude at set_mjd wrong for RA={ra_deg}, Dec={dec_deg}",
+                )
+
+    def test_rise_set_nan_cases(self):
+        """NaN is returned for circumpolar and never-above-limit fields.
+
+        At Rubin (lat ~-30.24 deg) with alt_limit = 20 deg:
+
+        - Dec = -80 deg: the lower-transit altitude is ~20.2 deg, so the field
+          is always above the limit (circumpolar w.r.t. alt_limit).
+        - Dec = +60 deg: the upper-transit altitude is ~-0.2 deg, so the field
+          never rises above the limit (or even above the horizon).
+        """
+        dtype = [("observationStartMJD", float), ("fieldRA", float), ("fieldDec", float)]
+        mjd = 51545.0
+
+        # Circumpolar field: should return NaN for both columns.
+        result = self.stacker.run(np.array([(mjd, 0.0, -80.0)], dtype=dtype))
+        self.assertTrue(np.isnan(result["rise_mjd"][0]), "Circumpolar field: rise_mjd should be NaN")
+        self.assertTrue(np.isnan(result["set_mjd"][0]), "Circumpolar field: set_mjd should be NaN")
+
+        # Field that never rises above alt_limit: should also return NaN.
+        result = self.stacker.run(np.array([(mjd, 0.0, 60.0)], dtype=dtype))
+        self.assertTrue(np.isnan(result["rise_mjd"][0]), "Never-rises field: rise_mjd should be NaN")
+        self.assertTrue(np.isnan(result["set_mjd"][0]), "Never-rises field: set_mjd should be NaN")
+
+
 if __name__ == "__main__":
     unittest.main()