Rework AltAzShadowMaskBasisFunction

Keeps previous capability of allowing disjoint altitude or azimuth regions, while also respecting kinematic model constraints
Added more extensive unit test

rubin_scheduler/scheduler/basis_functions/mask_basis_funcs.py

@@ -172,14 +172,13 @@ def _calc_value(self, conditions, indx=None):
 
 
 class AltAzShadowMaskBasisFunction(BaseBasisFunction):
-    """Mask out range altitudes and azimuths, then extend the
-    mask so if observations are taken in pairs, the second in the pair will
-    not have moved into a masked region.
+    """Mask out a range of altitudes and azimuths, including
+     regions which will enter the mask within `shadow_minutes`.
 
-    Masks any alt/az regions as specified by the conditions object, or
-    values as provided by the kwargs. If both are provided,
-    the most restrictive of each will be used.
-    This mask is then extended using `shadow minutes`.
+    Masks any alt/az regions as specified by the conditions object in
+    `conditions.tel_az_limits` and `conditions.tel_alt_limits`,
+    as well as values as provided by `min_alt`, `max_alt`,
+    `min_az` and `max_az`.
 
     Parameters
     ----------
@@ -192,7 +191,6 @@ class AltAzShadowMaskBasisFunction(BaseBasisFunction):
     min_az : `float`
         Minimum azimuth value to apply to the mask. Default 0 (degrees).
         These azimuth values are absolute azimuth, not cumulative.
-        The limits in the telescope model are based on cumulative azimuth.
         The absolute and cumulative azimuth only diverge if the azimuth
         range is greater than 360 degrees.
     max_az : `float`
@@ -203,68 +201,113 @@ class AltAzShadowMaskBasisFunction(BaseBasisFunction):
         be chosen to be scheduled and when it might be observed.
         For pairs, the minimum pair time + some buffer is good.
         For sequences, try the length of the sequence + some buffer.
-    pad : `float`
-        Pad the conditions alt/az limits by this amount (degrees).
-        Default corresponds to approximately the radius of the fov.
+        Note that this just sets up the shadow *at* the shadow_minutes time
+        (and not all times up to shadow_minutes).
     """
 
     def __init__(
-        self, nside=None, min_alt=20.0, max_alt=82.0, min_az=0, max_az=360, shadow_minutes=40.0, pad=2.0
+        self,
+        nside=None,
+        min_alt=20.0,
+        max_alt=82.0,
+        min_az=0,
+        max_az=360,
+        shadow_minutes=40.0,
     ):
         super().__init__(nside=nside)
         self.min_alt = np.radians(min_alt)
         self.max_alt = np.radians(max_alt)
         self.min_az = np.radians(min_az)
         self.max_az = np.radians(max_az)
         self.shadow_time = shadow_minutes / 60.0 / 24.0  # To days
-        self.pad = np.radians(pad)
 
-    def _calc_value(self, conditions, indx=None):
-        # Mask everything to start
-        result = np.zeros(hp.nside2npix(self.nside), dtype=float) + np.nan
+        self.r_min_alt = IntRounded(self.min_alt)
+        self.r_max_alt = IntRounded(self.max_alt)
 
-        in_range_alt = np.zeros(hp.nside2npix(self.nside), dtype=int)
-        in_range_az = np.zeros(hp.nside2npix(self.nside), dtype=int)
+    def _calc_value(self, conditions, indx=None):
+        # Basis function value will be 0 except where masked (then np.nan)
+        result = np.zeros(hp.nside2npix(self.nside), dtype=float)
 
         # Compute the alt,az values in the future. Use the conditions object
         # so the results are cached and can be used by other surveys is
         # needed. Technically this could fail if the masked region is
         # very narrow or shadow time is very large.
         future_alt, future_az = conditions.future_alt_az(np.max(conditions.mjd + self.shadow_time))
+        r_future_alt = IntRounded(future_alt)
+        r_current_alt = IntRounded(conditions.alt)
+
+        # Check the basis function altitude limits, now and future
+        result[np.where(r_current_alt < self.r_min_alt)] = np.nan
+        result[np.where(r_current_alt > self.r_max_alt)] = np.nan
+        result[np.where(r_future_alt < self.r_min_alt)] = np.nan
+        result[np.where(r_future_alt > self.r_max_alt)] = np.nan
+        # Check the conditions objects altitude limits, now and future
+        if conditions.tel_alt_limits is not None:
+            combined = np.zeros(hp.nside2npix(self.nside), dtype=float)
+            for limits in conditions.tel_alt_limits:
+                # For conditions-based limits, must add pad
+                # And remember that discontiguous areas can be allowed
+                in_bounds = np.ones(hp.nside2npix(self.nside), dtype=float)
+                min_alt = IntRounded(limits[0] + conditions.altaz_limit_pad)
+                max_alt = IntRounded(limits[1] - conditions.altaz_limit_pad)
+                in_bounds[np.where(r_current_alt < min_alt)] = 0
+                in_bounds[np.where(r_current_alt > max_alt)] = 0
+                in_bounds[np.where(r_future_alt < min_alt)] = 0
+                in_bounds[np.where(r_future_alt > max_alt)] = 0
+                combined += in_bounds
+            result[np.where(combined == 0)] = np.nan
+        # And check against the kinematic hard limits.
+        # These are separate so they don't override user tel_alt_limits.
+        min_alt = IntRounded(conditions.kinematic_alt_limits[0] + conditions.altaz_limit_pad)
+        max_alt = IntRounded(conditions.kinematic_alt_limits[1] - conditions.altaz_limit_pad)
+        result[np.where(r_current_alt < min_alt)] = np.nan
+        result[np.where(r_current_alt > max_alt)] = np.nan
+        result[np.where(r_future_alt < min_alt)] = np.nan
+        result[np.where(r_future_alt > max_alt)] = np.nan
 
-        # Find minimum alt limits
-        min_alt = IntRounded(np.max([conditions.tel_alt_min + self.pad, self.min_alt]))
-        max_alt = IntRounded(np.min([conditions.tel_alt_max - self.pad, self.max_alt]))
-
-        # Check allowable altitude range against current and future alt values
-        good = np.where((IntRounded(conditions.alt) >= min_alt) & (IntRounded(conditions.alt) <= max_alt))[0]
-
-        in_range_alt[good] += 1
-        good = np.where((IntRounded(future_alt) >= min_alt) & (IntRounded(future_alt) <= max_alt))[0]
-        in_range_alt[good] += 1
-
-        # Check allowable azimuth range against azimuth values
         # note that % (mod) is not supported for IntRounded
         two_pi = 2 * np.pi
-        itwo_pi = IntRounded(two_pi)
-        if IntRounded(conditions.tel_az_max) - IntRounded(conditions.tel_az_min) >= itwo_pi:
-            in_range1 = np.zeros(len(conditions.az)) + 2
-        else:
-            azr = (conditions.tel_az_max - conditions.tel_az_min) % (two_pi) - 2 * self.pad
-            # Check current and future
-            in_range1 = np.where((conditions.az - conditions.tel_az_min - self.pad) % (two_pi) <= azr, 1, 0)
-            in_range1 += np.where((future_az - conditions.tel_az_min - self.pad) % (two_pi) <= azr, 1, 0)
-        if IntRounded(self.max_az) - IntRounded(self.min_az) >= itwo_pi:
-            in_range2 = np.zeros(len(conditions.az)) + 2
-        else:
-            azr = (self.max_az - self.min_az) % (two_pi)
-            in_range2 = np.where((conditions.az - self.min_az) % (two_pi) <= azr, 1, 0)
-            in_range2 += np.where((future_az - self.min_az) % (two_pi) <= azr, 1, 0)
-        in_range_az[np.where((in_range1 > 1) & (in_range2 > 1))] = 2
-
-        passed_all = np.where((in_range_alt > 1) & (in_range_az > 1))[0]
-        # Unmask the parts of the sky which are and will remain in range
-        result[passed_all] = 0
+        # Check the basis function azimuth limits, now and future
+        if np.abs(self.max_az - self.min_az) < two_pi:
+            az_range = (self.max_az - self.min_az) % (two_pi)
+            out_of_bounds = np.where((conditions.az - self.min_az) % (two_pi) > az_range)[0]
+            result[out_of_bounds] = np.nan
+            out_of_bounds = np.where((future_az - self.min_az) % (two_pi) > az_range)[0]
+            result[out_of_bounds] = np.nan
+        # Check the conditions objects azimuth limits, now and future
+        if conditions.tel_az_limits is not None:
+            combined = np.zeros(hp.nside2npix(self.nside), dtype=float)
+            for limits in conditions.tel_az_limits:
+                in_bounds = np.ones(hp.nside2npix(self.nside), dtype=float)
+                min_az = limits[0]
+                max_az = limits[1]
+                if np.abs(max_az - min_az) < two_pi:
+                    az_range = (max_az - min_az) % (two_pi) - 2 * conditions.altaz_limit_pad
+                    out_of_bounds = np.where(
+                        (conditions.az - min_az - conditions.altaz_limit_pad) % (two_pi) > az_range
+                    )[0]
+                    in_bounds[out_of_bounds] = 0
+                    out_of_bounds = np.where(
+                        (future_az - min_az - conditions.altaz_limit_pad) % (two_pi) > az_range
+                    )[0]
+                    in_bounds[out_of_bounds] = 0
+                combined += in_bounds
+            result[np.where(combined == 0)] = np.nan
+        # Check against the kinematic hard limits.
+        if np.abs(conditions.kinematic_az_limits[1] - conditions.kinematic_az_limits[0]) < two_pi:
+            az_range = (conditions.kinematic_az_limits[1] - conditions.kinematic_az_limits[0]) % (
+                two_pi
+            ) - conditions.altaz_limit_pad * 2
+            out_of_bounds = np.where(
+                (conditions.az - conditions.kinematic_az_limits[0] - conditions.altaz_limit_pad) % (two_pi)
+                > az_range
+            )[0]
+            result[out_of_bounds] = np.nan
+            out_of_bounds = np.where(
+                (future_az - conditions.kinematic_az_limits[0] - conditions.altaz_limit_pad) % (two_pi)
+                > az_range
+            )[0]
+            result[out_of_bounds] = np.nan
 
         return result
 

rubin_scheduler/scheduler/features/conditions.py

@@ -188,13 +188,15 @@ def __init__(
         scheduled_observations : `np.ndarray`, (M,)
             A list of MJD times when there are scheduled observations.
             Defaults to empty array.
-        tel_az_min : `float`
-        tel_az_max : `float
-        tel_alt_min : `float`
-        tel_alt_max : `float`
-            Valid altitude and azimuth ranges (radians).
-            Order matters for azimuth - 270 degrees to 90 degrees
-            is different than 90 degrees to 270 degrees.
+        tel_az_limits : `list` [[`float`, `float`]]
+            A list of lists giving valid azimuth ranges. e.g.,
+            [0, 2*np.pi] would mean all azimuth values are valid, while
+            [[0, np.pi/2], [3*np.pi/2, 2*np.pi]] would mean anywhere in
+            the south is invalid.  Radians.
+        tel_alt_limits : `list` [[`float`, `float`]]
+            A list of lists giving valid altitude ranges. Radians.
+        altaz_limit_pad : `float`
+            Pad to surround the tel_az_limits and tel_alt_limits with.
 
         Attributes (calculated on demand and cached)
         ------------------------------------------
@@ -338,10 +340,11 @@ def _init_attributes(self):
         self.cumulative_azimuth_rad = None
 
         # Telescope limits
-        self.tel_az_min = None
-        self.tel_az_max = None
-        self.tel_alt_min = None
-        self.tel_alt_max = None
+        self.tel_az_limits = None
+        self.tel_alt_limits = None
+        self.kinematic_alt_limits = None
+        self.kinematic_az_limits = None
+        self.altaz_limit_pad = None
 
         # Full sky cloud map
         self._cloud_map = None

rubin_scheduler/scheduler/model_observatory/model_observatory.py

@@ -126,6 +126,18 @@ class ModelObservatory:
         to the next night. Default "sun_n12_rising", e.g., sun at
         -12 degrees and rising. Other options are "sun_n18_rising"
         and "sunrise".
+    tel_az_limits : `list` [[`float`, `float`]]
+        A list of lists giving valid azimuth ranges. e.g.,
+        [0, 180] would mean all azimuth values are valid, while
+        [[0, 90], [270, 360]] or [270, 90] would mean anywhere in
+        the south is invalid.  Degrees.
+    tel_alt_limits : `list` [[`float`, `float`]]
+        A list of lists giving valid altitude ranges. Degrees.
+        For both the alt and az limits, if a particular alt (or az)
+        value is included in any limit, it is valid for all of them.
+        Altitude limits of [[20, 40], [40, 60]] would allow altitudes
+        betwen 20-40 and 40-60, but [[20, 40], [40, 60], [20, 86]]
+        will allow altitudes anywhere between 20-86 degrees.
     telescope : `str`
         Telescope name for rotation computations. Default "rubin".
     """
@@ -152,6 +164,8 @@ def __init__(
         wind_data=None,
         starting_time_key="sun_n12_setting",
         ending_time_key="sun_n12_rising",
+        tel_alt_limits=None,
+        tel_az_limits=None,
         telescope="rubin",
     ):
         if nside is None:
@@ -284,6 +298,28 @@ def __init__(
         else:
             self.observatory = kinem_model
 
+        # Pick up tel alt and az limits from kwargs
+        if tel_alt_limits is not None:
+            self.tel_alt_limits = np.radians(tel_alt_limits)
+        else:
+            self.tel_alt_limits = []
+        if tel_az_limits is not None:
+            self.tel_az_limits = np.radians(tel_az_limits)
+        else:
+            self.tel_az_limits = []
+        # Add the observatory alt/az limits to the user-defined limits
+        # But we do have to be careful that we're not overriding more
+        # restrictive limits that were already set.
+        # So we'll just keep these separate.
+        self.kinematic_tel_alt_limits = [
+            self.observatory.telalt_minpos_rad,
+            self.observatory.telalt_maxpos_rad,
+        ]
+        self.kinematic_tel_az_limits = [self.observatory.telaz_minpos_rad, self.observatory.telaz_maxpos_rad]
+        # Each of these limits will be treated as hard limits that we don't
+        # want pointings to stray into, so add a pad around the values
+        self.altaz_limit_pad = np.radians(2.0)
+
         # Let's make sure we're at an openable MJD
         good_mjd = False
         to_set_mjd = mjd
@@ -436,10 +472,11 @@ def return_conditions(self):
         self.conditions.mjd_start = self.mjd_start
 
         # Telescope limits
-        self.conditions.tel_az_min = self.observatory.telaz_minpos_rad
-        self.conditions.tel_az_max = self.observatory.telaz_maxpos_rad
-        self.conditions.tel_alt_min = self.observatory.telalt_minpos_rad
-        self.conditions.tel_alt_max = self.observatory.telalt_maxpos_rad
+        self.conditions.tel_az_limits = self.tel_az_limits
+        self.conditions.tel_alt_limits = self.tel_alt_limits
+        self.conditions.kinematic_alt_limits = self.kinematic_tel_alt_limits
+        self.conditions.kinematic_az_limits = self.kinematic_tel_az_limits
+        self.conditions.altaz_limit_pad = self.altaz_limit_pad
 
         # Planet positions from almanac
         self.conditions.planet_positions = self.almanac.get_planet_positions(self.mjd)

rubin_scheduler/scheduler/surveys/scripted_surveys.py

@@ -212,20 +212,36 @@ def _check_alts_ha(self, observation, conditions):
 
         # Also check the alt,az limits given by the conditions object
         count = in_range * 0
-        ir = np.where((alt[in_range] >= conditions.tel_az_min) & (alt[in_range] <= conditions.tel_az_max))[0]
-        count[ir] += 1
-        good = np.where(count > 0)[0]
-        in_range = in_range[good]
+        if conditions.tel_alt_limits is not None:
+            for limits in conditions.tel_alt_limits:
+                ir = np.where((alt[in_range] >= np.min(limits)) & (alt[in_range] <= np.max(limits)))[0]
+                count[ir] += 1
+            good = np.where(count > 0)[0]
+            in_range = in_range[good]
+        # Check against kinematic limits too
+        ir = np.where(
+            (alt[in_range] >= np.min(conditions.kinematic_alt_limits))
+            & (alt[in_range] <= np.max(conditions.kinematic_alt_limits))
+        )[0]
+        in_range = in_range[ir]
 
         count = in_range * 0
-        if (conditions.tel_az_max - conditions.tel_az_min) >= (2 * np.pi):
-            count += 1
-        else:
-            az_range = (conditions.tel_az_max - conditions.tel_az_min) % (2 * np.pi)
-            ir = np.where((az[in_range] - conditions.tel_az_min) % (2 * np.pi) <= az_range)[0]
-            count[ir] += 1
-        good = np.where(count > 0)[0]
-        in_range = in_range[good]
+        if conditions.tel_az_limits is not None:
+            for limits in conditions.tel_az_limits:
+                az_min = limits[0]
+                az_max = limits[1]
+                if np.abs(az_max - az_min) >= (2 * np.pi):
+                    count += 1
+                else:
+                    az_range = (az_max - az_min) % (2 * np.pi)
+                    ir = np.where((az[in_range] - az_min) % (2 * np.pi) <= az_range)[0]
+                    count[ir] += 1
+            good = np.where(count > 0)[0]
+            in_range = in_range[good]
+        if np.abs(conditions.kinematic_az_limits[1] - conditions.kinematic_az_limits[0]) < (2 * np.pi):
+            az_range = (conditions.kinematic_az_limits[1] - conditions.kinematic_az_limits[0]) % (2 * np.pi)
+            ir = np.where((az[in_range] - az_min) % (2 * np.pi) <= az_range)[0]
+            in_range = in_range[ir]
 
         # Check that filter needed is mounted
         good = np.isin(observation["filter"][in_range], conditions.mounted_filters)