arbitrary dither pattern detailer

rubin_scheduler/scheduler/detailers/dither_detailer.py

@@ -4,12 +4,13 @@
     "CameraRotDetailer",
     "CameraSmallRotPerObservationListDetailer",
     "ComCamGridDitherDetailer",
+    "DeltaCoordDitherDetailer",
 )
 
 import numpy as np
 
 from rubin_scheduler.scheduler.detailers import BaseDetailer
-from rubin_scheduler.scheduler.utils import wrap_ra_dec
+from rubin_scheduler.scheduler.utils import ObservationArray, rotx, thetaphi2xyz, wrap_ra_dec, xyz2thetaphi
 from rubin_scheduler.utils import (
     _approx_altaz2pa,
     _approx_ra_dec2_alt_az,
@@ -84,6 +85,75 @@ def __call__(self, obs_array, conditions):
         return obs_array
 
 
+class DeltaCoordDitherDetailer(BaseDetailer):
+    """Dither pattern set by user
+
+    Parameters
+    ----------
+    delta_ra : `np.array`
+        Angular distances to move in the RA-direction (degrees).
+    delta_dec : `np.array`
+        Angular distances to move in the dec-direction (degree).
+    delta_rotskypos : `np.array`
+        Angular shifts to make in the rotskypos values (degrees). Default None
+        applies no rotational shift
+    """
+
+    def __init__(self, delta_ra, delta_dec, delta_rotskypos=None):
+
+        if delta_rotskypos is not None:
+            if np.size(delta_ra) != np.size(delta_rotskypos):
+                raise ValueError(
+                    "size of delta_ra (%i) is not equal to size of delta_rotskypos (%i)"
+                    % (np.size(delta_ra), np.size(delta_rotskypos))
+                )
+        if np.size(delta_ra) != np.size(delta_dec):
+            raise ValueError(
+                "Sizes of delta_ra (%i) and delta_dec (%i) do not match."
+                % (np.size(delta_ra), np.size(delta_dec))
+            )
+
+        self.survey_features = {}
+        self.delta_ra = np.radians(delta_ra)
+        self.delta_dec = np.radians(delta_dec)
+        if delta_rotskypos is None:
+            self.delta_rotskypos = delta_rotskypos
+        else:
+            self.delta_rotskypos = np.radians(delta_rotskypos)
+
+    def __call__(self, obs_array, conditions):
+
+        output_array_list = []
+        for obs in obs_array:
+            dithered_observations = ObservationArray(self.delta_ra.size)
+            for key in obs.dtype.names:
+                dithered_observations[key] = obs[key]
+            new_decs = self.delta_dec + 0
+            new_ras = self.delta_ra + np.pi / 2  # Adding 90 deg here for later rotation
+            # Rotate new positions to be centered on correct RA,dec.
+            # This way positions on poles should work fine.
+
+            #  Rotate ra and dec about the x-axis
+            x, y, z = thetaphi2xyz(new_ras, new_decs + np.pi / 2.0)
+            xp, yp, zp = rotx(obs["dec"], x, y, z)
+            theta, phi = xyz2thetaphi(xp, yp, zp)
+
+            new_decs = phi - np.pi / 2
+            new_ras = theta - np.pi / 2 + obs["RA"]  # Need to remove rotation from earlier
+
+            # Make sure coords are in proper range
+            new_ras, new_decs = wrap_ra_dec(new_ras, new_decs)
+
+            dithered_observations["RA"] = new_ras
+            dithered_observations["dec"] = new_decs
+            if self.delta_rotskypos is not None:
+                dithered_observations["rotSkyPos_desired"] += self.delta_rotskypos
+            output_array_list.append(dithered_observations)
+
+        result = np.concatenate(output_array_list)
+        return result
+
+
 class EuclidDitherDetailer(BaseDetailer):
     """Directional dithering for Euclid DDFs
 

rubin_scheduler/scheduler/surveys/base_survey.py

@@ -12,6 +12,7 @@
     HpInLsstFov,
     ObservationArray,
     comcam_tessellate,
+    rotx,
     thetaphi2xyz,
     xyz2thetaphi,
 )
@@ -430,16 +431,6 @@ def bf_short_labels(self):
         return label_map
 
 
-def rotx(theta, x, y, z):
-    """rotate the x,y,z points theta radians about x axis"""
-    sin_t = np.sin(theta)
-    cos_t = np.cos(theta)
-    xp = x
-    yp = y * cos_t + z * sin_t
-    zp = -y * sin_t + z * cos_t
-    return xp, yp, zp
-
-
 class BaseMarkovSurvey(BaseSurvey):
     """A Markov Decision Function survey object. Uses Basis functions
     to compute a final reward function and decide what to observe based

rubin_scheduler/scheduler/utils/utils.py

@@ -22,6 +22,7 @@
     "xyz2thetaphi",
     "mean_azimuth",
     "wrap_ra_dec",
+    "rotx",
 )
 
 import datetime
@@ -63,18 +64,38 @@ def smallest_signed_angle(a1, a2):
 
 
 def thetaphi2xyz(theta, phi):
+    """Convert theta,phi to x,y,z position on the unit sphere
+
+    Parameters
+    ----------
+    theta : `float`
+        Theta coordinate in radians (should be 0-2pi).
+    phi : `float`
+        Phi coordinate in radians (should run from 0-pi).
+    """
     x = np.sin(phi) * np.cos(theta)
     y = np.sin(phi) * np.sin(theta)
     z = np.cos(phi)
     return x, y, z
 
 
 def xyz2thetaphi(x, y, z):
+    """x,y,z position on unit sphere to theta,phi coords."""
     phi = np.arccos(z)
     theta = np.arctan2(y, x)
     return theta, phi
 
 
+def rotx(theta, x, y, z):
+    """rotate the x,y,z points theta radians about x axis"""
+    sin_t = np.sin(theta)
+    cos_t = np.cos(theta)
+    xp = x
+    yp = y * cos_t + z * sin_t
+    zp = -y * sin_t + z * cos_t
+    return xp, yp, zp
+
+
 def wrap_ra_dec(ra, dec):
     # XXX--from MAF, should put in general utils
     """

tests/scheduler/test_detailers.py

@@ -79,6 +79,62 @@ def test_start_field(self):
 
         assert len(obs_out) == 4
 
+    def test_delta_dither(self):
+        obs = ObservationArray(2)
+        obs["RA"] = np.radians(0)
+        obs["dec"] = np.radians(20)
+
+        det = detailers.DeltaCoordDitherDetailer(np.array([0.0, 0.0]), np.array([-10.0, 10.0]))
+
+        output = det(obs, [])
+
+        assert np.size(output) == 4
+        assert np.unique(output["RA"]) == np.unique(obs["RA"])
+
+        obs["RA"] = np.radians(90)
+        det = detailers.DeltaCoordDitherDetailer(np.array([0.0, 0.0]), np.array([-10.0, 10.0]))
+        output = det(obs, [])
+        assert np.size(output) == 4
+        assert np.unique(output["RA"]) == np.unique(obs["RA"])
+
+        ra_step = 5.0
+        det = detailers.DeltaCoordDitherDetailer(np.array([-ra_step, ra_step]), np.array([0.0, 0.0]))
+        output = det(obs, [])
+
+        # Make sure ra diff is larger
+        assert np.all(np.abs(output["RA"] - obs["RA"][0]) > np.radians(ra_step))
+
+        # Try having rotation as well
+        det = detailers.DeltaCoordDitherDetailer(
+            np.array([-ra_step, ra_step]), np.array([0.0, 0.0]), delta_rotskypos=np.array([5.0, 10.0])
+        )
+        output = det(obs, [])
+
+        assert np.size(output) == 4
+        assert np.size(np.unique(output["rotSkyPos_desired"])) == 2
+
+        # Make sure one-element works
+        obs = ObservationArray(1)
+        obs["RA"] = np.radians(0)
+        obs["dec"] = np.radians(20)
+        det = detailers.DeltaCoordDitherDetailer(np.array([0.0, 0.0]), np.array([-10.0, 10.0]))
+        output = det(obs, [])
+
+        assert np.size(output) == 2
+        assert np.unique(output["RA"]) == np.unique(obs["RA"])
+
+        # Test going to the pole
+        obs = ObservationArray(2)
+        obs["RA"] = np.radians(0)
+        obs["dec"] = np.radians(-90)
+
+        det = detailers.DeltaCoordDitherDetailer(np.array([-1.0, 1.0, -1, 1]), np.array([1.0, 1.0, -1, -1]))
+        output = det(obs, [])
+
+        # This should make a grid all at the same dec
+        assert np.size(np.unique(output["dec"])) == 1
+        assert output["dec"][0] > obs["dec"][0]
+
     def test_random_band(self):
         obs = ObservationArray(1)
         obs["band"] = "r"