Merge pull request #114 from lsst/tickets/SP-1667

tickets/SP-1667: Merge the Times Square OR4 nightsum with the nightly nightsum

schedview/__init__.py

@@ -1,2 +1,3 @@
+from .dayobs import DayObs
 from .sphere import *
 from .version import __version__

schedview/collect/__init__.py

@@ -6,6 +6,7 @@
     "load_bright_stars",
     "read_ddf_visits",
     "read_multiple_opsims",
+    "read_night_visits",
 ]
 
 from .footprint import get_footprint

schedview/collect/consdb.py

@@ -1,6 +1,7 @@
 import numpy as np
 import pandas as pd
 import rubin_sim.maf
+from rubin_scheduler.scheduler.utils import ObservationArray
 from rubin_scheduler.utils.consdb import load_consdb_visits
 
 
@@ -11,7 +12,12 @@ def read_consdb(
     **kwargs,
 ) -> pd.DataFrame:
     consdb_visits = load_consdb_visits(instrument, *args, **kwargs)
-    visit_records: np.recarray = consdb_visits.opsim.to_records()
+    # If the return is empty, to_records won't be able to determine the
+    # right dtypes, so use ObservationArray()[0:0] instead to make the
+    # empty recarray with the correct dtypes.
+    visit_records: np.recarray = (
+        consdb_visits.opsim.to_records() if len(consdb_visits.opsim) > 0 else ObservationArray()[0:0]
+    )
     for stacker in stackers:
         visit_records = stacker.run(visit_records)
     visits = pd.DataFrame(visit_records)

schedview/collect/visits.py

@@ -0,0 +1,78 @@
+import pandas as pd
+from lsst.resources import ResourcePath
+from rubin_scheduler.utils import ddf_locations
+from rubin_scheduler.utils.consdb import KNOWN_INSTRUMENTS
+from rubin_sim import maf
+
+from schedview import DayObs
+
+from .consdb import read_consdb
+from .opsim import read_opsim
+
+# Use old-style format, because f-strings are not reusable
+OPSIMDB_TEMPLATE = (
+    "/sdf/group/rubin/web_data/sim-data/sims_featureScheduler_runs{sim_version}/baseline/"
+    + "baseline_v{sim_version}_10yrs.db"
+)
+
+NIGHT_STACKERS = [
+    maf.HourAngleStacker(),
+    maf.stackers.ObservationStartDatetime64Stacker(),
+    maf.stackers.TeffStacker(),
+    maf.stackers.OverheadStacker(),
+    maf.stackers.DayObsISOStacker(),
+]
+
+DDF_STACKERS = [
+    maf.stackers.ObservationStartDatetime64Stacker(),
+    maf.stackers.TeffStacker(),
+    maf.stackers.DayObsISOStacker(),
+]
+
+
+def read_visits(
+    day_obs: str | int | DayObs,
+    visit_source: str,
+    stackers: list[maf.stackers.base_stacker.BaseStacker] = [],
+    num_nights: int = 1,
+) -> pd.DataFrame:
+
+    if visit_source in KNOWN_INSTRUMENTS:
+        visits = read_consdb(
+            visit_source,
+            stackers=stackers,
+            day_obs=DayObs.from_date(day_obs).date.isoformat(),
+            num_nights=num_nights,
+        )
+    else:
+        baseline_opsim_rp = ResourcePath(OPSIMDB_TEMPLATE.format(sim_version=visit_source))
+        mjd: int = DayObs.from_date(day_obs).mjd
+        visits = read_opsim(
+            baseline_opsim_rp,
+            constraint=f"FLOOR(observationStartMJD-0.5)<={mjd}"
+            + f" AND FLOOR(observationStartMJD-0.5)>({mjd-num_nights})",
+            stackers=stackers,
+        )
+    return visits
+
+
+def read_ddf_visits(*args, **kwargs) -> pd.DataFrame:
+
+    if "stackers" not in kwargs:
+        kwargs["stackers"] = DDF_STACKERS
+
+    all_visits = read_visits(*args, **kwargs)
+
+    ddf_field_names = tuple(ddf_locations().keys())
+    # Different versions of the schedule include a DD: prefix, or not.
+    # Catch them all.
+    ddf_field_names = ddf_field_names + tuple([f"DD:{n}" for n in ddf_field_names])
+
+    # Figure out which column has the target names.
+    target_column_name = "target_name" if "target_name" in all_visits.columns else "target"
+    if target_column_name not in all_visits.columns:
+        raise ValueError("Cannot find a column in visits with which to identify DDF fields.")
+
+    ddf_visits = all_visits.loc[all_visits[target_column_name].isin(ddf_field_names)]
+
+    return ddf_visits

schedview/dayobs.py

@@ -0,0 +1,195 @@
+import copy
+import datetime
+from dataclasses import dataclass
+from functools import cached_property
+from typing import Literal, Self
+
+import dateutil.parser
+from astropy.time import Time
+
+DAYOBS_TZ = datetime.timezone(datetime.timedelta(hours=-12))
+MJD_EPOCH = datetime.date(1858, 11, 17)
+ONE_DAY = datetime.timedelta(days=1)
+
+IntDateFormat = Literal["mjd", "yyyymmdd", "auto"]
+
+
+# Use a frozen dataclass so that we can used cached properties without
+# having to worry about updating them.
+@dataclass(frozen=True)
+class DayObs:
+    """Represent a day of observation, dayobs, as defined in SITCOMTN-032.
+    That is, the date in the -12hr timezone, such that the "date" of an
+    entire night of observing is the calendar date in Chile of the evening in
+    which the night begins.
+
+    Parameters
+    ----------
+    date : `datetime.date`
+        The calendar date.
+    int_format : `str`
+        ``yyyymmdd`` if an integer representation is a mapping of decimal
+        digits to year, month, and day; ``mjd`` if the integer representation
+        is the Modified Julian Date. The default is ``mjd``.
+    """
+
+    # Use the standard library date as a lowest common denominator
+    # that cannot represent more or less than we want it to.
+    date: datetime.date
+    int_format: IntDateFormat = "auto"
+
+    @classmethod
+    def from_date(cls, arg: datetime.date | int | str | Self, int_format: IntDateFormat = "auto") -> Self:
+        """Create a representation of day_obs for a given date.
+        This factory takes a representation of the date (already in the
+        -12hr timezone as defined in SITCOMTN-032), not a date and time.
+
+        Parameters
+        ----------
+        arg : `datetime.date` | `int` | `str`
+            The representation of dayobs, as defined in SITCOMTN-032.
+            If an integer, it will be interpreted according to the
+            ``int_format`` argument.
+            If it is "yesterday", "today", or "tomorrow", it uses the day_obs
+            reletive to the current time.
+            A reasonable attempt is made to interpret other strings as
+            dates.
+        int_format : `str` (optional)
+            One of ``mjd``, in which case integers are interpreted as Modified
+            Julian Dates,
+            ``yyyymmdd``, in which case integers encode year, month, and day
+            into decimal digits,
+            or ``auto``, in which case 8 digit decimals are interpretd as
+            yyyymmdd and others as mjd.
+
+        Returns
+        -------
+        day_obs_converter: `DayObsConverter`
+            A new instance of the converter.
+        """
+
+        # If the argument is convertable to an int, do.
+        if isinstance(arg, str) and int_format in ("yyyymmdd", "auto"):
+            try:
+                arg = datetime.datetime.strptime(arg, "%Y%m%d").date()
+            except ValueError:
+                pass
+
+        match arg:
+            case DayObs():
+                return copy.deepcopy(arg)
+            case datetime.date():
+                date = arg
+            case int():
+                if int_format == "yyyymmdd" or (
+                    len(str(arg)) == len(str("YYYYMMDD")) and int_format == "auto"
+                ):
+                    # Digits encode YYYYMMDD
+                    date = datetime.datetime.strptime(str(arg), "%Y%m%d").date()
+                else:
+                    if int_format not in ("auto", "mjd"):
+                        raise ValueError("Invalid integer format.")
+
+                    date = MJD_EPOCH + datetime.timedelta(days=arg)
+            case "yesterday":
+                date = (datetime.datetime.now(tz=DAYOBS_TZ) - ONE_DAY).date()
+            case "today" | "tonight":
+                date = datetime.datetime.now(tz=DAYOBS_TZ).date()
+            case "tomorrow":
+                date = (datetime.datetime.now(tz=DAYOBS_TZ) + ONE_DAY).date()
+            case str():
+                dayobs_datetime = dateutil.parser.parse(arg)
+                if dayobs_datetime.hour != 0 or dayobs_datetime.minute != 0 or dayobs_datetime.second != 0:
+                    raise ValueError("The argument to from_dayobs must be a date, not a date and time.")
+                if dayobs_datetime.tzinfo is not None and dayobs_datetime.tzinfo != DAYOBS_TZ:
+                    raise ValueError("The argument to from_dayobs must naive or in the dayobs timezone")
+                date = datetime.date(dayobs_datetime.year, dayobs_datetime.month, dayobs_datetime.day)
+            case _:
+                raise NotImplementedError()
+
+        return cls(date, int_format)
+
+    @classmethod
+    def from_time(
+        cls, arg: datetime.datetime | str | float | Time, int_format: IntDateFormat = "mjd"
+    ) -> Self:
+        """Create a representation of the dayobs that includes a given time.
+
+        Parameters
+        ----------
+        arg : datetime.datetime | str | float | Time
+            A time in the date to return the day_obs of.
+            Floats are interpreted as Modified Julian Dates (in UTC).
+            "now" is the time now.
+            Representations without timezones are assumed to be in UTC.
+        int_format : `str`
+            If `mjd`, represent the date as an MJD when cast to an integer.
+            If `yyyymmdd`, encode year month and day into decimal digits
+            instead.
+            `mjd` by default.
+
+        Returns
+        -------
+        Returns
+        -------
+        day_obs_converter: `DayObsConverter`
+            A new instance of the converter.
+        """
+
+        match arg:
+            case datetime.datetime():
+                dayobs_datetime: datetime.datetime = arg
+            case "now":
+                dayobs_datetime: datetime.datetime = datetime.datetime.now(tz=datetime.timezone.utc)
+            case str():
+                dayobs_datetime: datetime.datetime = dateutil.parser.parse(arg)
+            case float():
+                # Interpret floats as UTC MJDs
+                maybe_datetime = Time(arg, format="mjd").datetime
+                if isinstance(maybe_datetime, datetime.datetime):
+                    dayobs_datetime: datetime.datetime = maybe_datetime
+                else:
+                    raise ValueError(f"{cls.__name__} currently only accepts scalars.")
+            case Time():
+                maybe_datetime = arg.datetime
+                if isinstance(maybe_datetime, datetime.datetime):
+                    dayobs_datetime: datetime.datetime = maybe_datetime
+                else:
+                    raise ValueError(f"{cls.__name__} currently only accepts scalars.")
+            case _:
+                raise NotImplementedError()
+
+        assert isinstance(dayobs_datetime, datetime.datetime)
+
+        # If dayobs_datetime is not timezone aware, assume UTC
+        if dayobs_datetime.tzinfo is None:
+            dayobs_datetime = dayobs_datetime.replace(tzinfo=datetime.timezone.utc)
+
+        dayobs_date = dayobs_datetime.astimezone(DAYOBS_TZ).date()
+        return cls(dayobs_date, int_format)
+
+    @cached_property
+    def yyyymmdd(self) -> int:
+        return self.date.day + 100 * (self.date.month + 100 * self.date.year)
+
+    @cached_property
+    def mjd(self) -> int:
+        return (self.date - MJD_EPOCH).days
+
+    @cached_property
+    def start(self) -> Time:
+        start_datetime = datetime.datetime(self.date.year, self.date.month, self.date.day, tzinfo=DAYOBS_TZ)
+        return Time(start_datetime)
+
+    @cached_property
+    def end(self) -> Time:
+        end_datetime = (
+            datetime.datetime(self.date.year, self.date.month, self.date.day, tzinfo=DAYOBS_TZ) + ONE_DAY
+        )
+        return Time(end_datetime)
+
+    def __int__(self) -> int:
+        return self.mjd if self.int_format in ("auto", "mjd") else self.yyyymmdd
+
+    def __str__(self):
+        return self.date.isoformat()

schedview/plot/survey.py

@@ -449,7 +449,7 @@ def create_metric_visit_map_grid(
     if len(metric_visits):
         metric_hpix = compute_hpix_metric_in_bands(metric_visits, metric, nside=nside)
     else:
-        metric_hpix = np.zeros(hp.nside2npix(nside))
+        metric_hpix = {b: np.zeros(hp.nside2npix(nside)) for b in visits["filter"].unique()}
 
     if len(visits):
         if use_matplotlib:

tests/test_dayobs.py

@@ -0,0 +1,56 @@
+import datetime
+import unittest
+from collections import namedtuple
+
+from astropy.time import Time
+
+from schedview import DayObs
+
+DayObsTestData = namedtuple("DayObsTestData", ("date", "yyyymmdd", "iso_date", "mjd", "iso_times"))
+
+
+class TestDayObs(unittest.TestCase):
+    test_values = (
+        DayObsTestData(
+            datetime.date(2024, 9, 3),
+            20240903,
+            "2024-09-03",
+            60556,
+            ("2024-09-03T12:00:00Z", "2024-09-04T05:00:00Z"),
+        ),
+        DayObsTestData(
+            datetime.date(2024, 10, 31),
+            20241031,
+            "2024-10-31",
+            60614,
+            ("2024-10-31T12:00:00Z", "2024-11-01T05:00:00Z"),
+        ),
+    )
+
+    @staticmethod
+    def assert_equal(dayobs, test_data):
+        assert dayobs.date == test_data.date
+        assert dayobs.yyyymmdd == test_data.yyyymmdd
+        assert dayobs.date.isoformat() == test_data.iso_date
+        assert dayobs.mjd == test_data.mjd
+        assert str(dayobs) == test_data.iso_date
+        if dayobs.int_format in ("mjd", "auto"):
+            assert int(dayobs) == test_data.mjd
+        else:
+            assert int(dayobs) == test_data.yyyymmdd
+
+    def test_dayobs(self):
+        for d in self.test_values:
+            self.assert_equal(DayObs.from_date(d.date), d)
+            self.assert_equal(DayObs.from_date(d.yyyymmdd), d)
+            self.assert_equal(DayObs.from_date(d.iso_date), d)
+            self.assert_equal(DayObs.from_date(d.mjd), d)
+
+            for iso_time in d.iso_times:
+                self.assert_equal(DayObs.from_time(iso_time), d)
+
+                t = Time(iso_time)
+                self.assert_equal(DayObs.from_time(t), d)
+                self.assert_equal(DayObs.from_time(t.datetime), d)
+                self.assert_equal(DayObs.from_time(t.iso), d)
+                self.assert_equal(DayObs.from_time(t.mjd), d)