Dithering utils - WIP

pocs/tests/utils/test_dither.py

@@ -0,0 +1,175 @@
+import pytest
+import astropy.units as u
+from astropy.coordinates import SkyCoord
+from astropy.coordinates import Angle
+
+from pocs.utils import dither
+
+
+def test_dice9_SkyCoord():
+    base = SkyCoord("16h52m42.2s -38d37m12s")
+
+    positions = dither.get_dither_positions(base_position=base,
+                                            n_positions=12,
+                                            pattern=dither.dice9,
+                                            pattern_offset=30 * u.arcminute)
+
+    assert isinstance(positions, SkyCoord)
+    assert len(positions) == 12
+    # postion 0 should be the base position
+    assert positions[0].separation(base) < Angle(1e12 * u.degree)
+    # With no random offset positions 9, 10, 11 should be the same as 0, 1, 2
+    assert positions[0:3].to_string() == positions[9:12].to_string()
+    # Position 1 should be 30 arcminute offset from base, in declination direction only
+    assert base.spherical_offsets_to(
+        positions[1])[0].radian == pytest.approx(Angle(0 * u.degree).radian)
+    assert base.spherical_offsets_to(
+        positions[1])[1].radian == pytest.approx(Angle(0.5 * u.degree).radian)
+    # Position 3 should be 30 arcminute offset from base in RA only.
+    assert base.spherical_offsets_to(
+        positions[3])[0].radian == pytest.approx(Angle(0.5 * u.degree).radian)
+    assert base.spherical_offsets_to(
+        positions[3])[1].radian == pytest.approx(Angle(0 * u.degree).radian)
+
+
+def test_dice9_string():
+    base = "16h52m42.2s -38d37m12s"
+
+    positions = dither.get_dither_positions(base_position=base,
+                                            n_positions=12,
+                                            pattern=dither.dice9,
+                                            pattern_offset=30 * u.arcminute)
+
+    base = SkyCoord(base)
+
+    assert isinstance(positions, SkyCoord)
+    assert len(positions) == 12
+    # postion 0 should be the base position
+    assert positions[0].separation(base) < Angle(1e12 * u.degree)
+    # With no random offset positions 9, 10, 11 should be the same as 0, 1, 2
+    assert positions[0:3].to_string() == positions[9:12].to_string()
+    # Position 1 should be 30 arcminute offset from base, in declination direction only
+    assert base.spherical_offsets_to(
+        positions[1])[0].radian == pytest.approx(Angle(0 * u.degree).radian)
+    assert base.spherical_offsets_to(
+        positions[1])[1].radian == pytest.approx(Angle(0.5 * u.degree).radian)
+    # Position 3 should be 30 arcminute offset from base in RA only.
+    assert base.spherical_offsets_to(
+        positions[3])[0].radian == pytest.approx(Angle(0.5 * u.degree).radian)
+    assert base.spherical_offsets_to(
+        positions[3])[1].radian == pytest.approx(Angle(0 * u.degree).radian)
+
+
+def test_dice9_bad_base_position():
+    with pytest.raises(ValueError):
+        dither.get_dither_positions(base_position=42,
+                                    n_positions=42,
+                                    pattern=dither.dice9,
+                                    pattern_offset=300 * u.arcsecond)
+
+
+def test_dice9_random():
+    base = SkyCoord("16h52m42.2s -38d37m12s")
+
+    positions = dither.get_dither_positions(base_position=base,
+                                            n_positions=12,
+                                            pattern=dither.dice9,
+                                            pattern_offset=30 * u.arcminute,
+                                            random_offset=30 * u.arcsecond)
+
+    assert isinstance(positions, SkyCoord)
+    assert len(positions) == 12
+    # postion 0 should be the base position
+    assert positions[0].separation(base) < Angle(30 * 2**0.5 * u.arcsecond)
+
+    angle_0 = Angle(0 * u.degree).radian
+    angle_05 = Angle(0.5 * u.degree).radian
+    angle_30 = Angle(30 * u.arcsecond).radian
+    position_1_offset = base.spherical_offsets_to(positions[1])
+    position_3_offset = base.spherical_offsets_to(positions[3])
+
+    # Position 1 should be 30 arcminute offset from base, in declination direction only
+    assert position_1_offset[0].radian == pytest.approx(angle_0, abs=angle_30)
+    assert position_1_offset[1].radian == pytest.approx(angle_05, abs=angle_30)
+
+    # Position 3 should be 30 arcminute offset from base in RA only.
+    assert position_3_offset[0].radian == pytest.approx(angle_05, abs=angle_30)
+    assert position_3_offset[1].radian == pytest.approx(angle_0, abs=angle_30)
+
+
+def test_random():
+    base = SkyCoord("16h52m42.2s -38d37m12s")
+
+    positions = dither.get_dither_positions(base_position=base,
+                                            n_positions=12,
+                                            random_offset=30 * u.arcsecond)
+    assert isinstance(positions, SkyCoord)
+    assert len(positions) == 12
+
+    angle_0 = Angle(0 * u.degree).radian
+    angle_30 = Angle(30 * u.arcsecond).radian
+    position_0_offset = base.spherical_offsets_to(positions[0])
+    position_1_offset = base.spherical_offsets_to(positions[1])
+
+    assert position_0_offset[0].radian == pytest.approx(angle_0, abs=angle_30)
+    assert position_0_offset[1].radian == pytest.approx(angle_0, abs=angle_30)
+
+    assert position_1_offset[0].radian == pytest.approx(angle_0, abs=angle_30)
+    assert position_1_offset[1].radian == pytest.approx(angle_0, abs=angle_30)
+
+
+def test_dice5():
+    base = SkyCoord("16h52m42.2s -38d37m12s")
+
+    positions = dither.get_dither_positions(base_position=base,
+                                            n_positions=12,
+                                            pattern=dither.dice5,
+                                            pattern_offset=30 * u.arcminute)
+
+    assert isinstance(positions, SkyCoord)
+    assert len(positions) == 12
+    # postion 0 should be the base position
+    assert positions[0].separation(base) < Angle(1e12 * u.degree)
+    # With no random offset positions 5, 6, 7 should be the same as 0, 1, 2
+    assert positions[0:3].to_string() == positions[5:8].to_string()
+    # Position 1 should be 30 arcminute offset from base, in RA and dec
+    assert base.spherical_offsets_to(
+        positions[1])[0].radian == pytest.approx(Angle(0.5 * u.degree).radian)
+    assert base.spherical_offsets_to(
+        positions[1])[1].radian == pytest.approx(Angle(0.5 * u.degree).radian)
+    # Position 3 should be 30 arcminute offset from base in RA and dec
+    assert base.spherical_offsets_to(positions[3])[0].radian == pytest.approx(
+        Angle(-0.5 * u.degree).radian)
+    assert base.spherical_offsets_to(positions[3])[1].radian == pytest.approx(
+        Angle(-0.5 * u.degree).radian)
+
+
+def test_custom_pattern():
+    base = SkyCoord("16h52m42.2s -38d37m12s")
+    cross = ((0, 0),
+             (0, 1),
+             (1, 0),
+             (0, -1),
+             (-1, 0))
+
+    positions = dither.get_dither_positions(base_position=base,
+                                            n_positions=12,
+                                            pattern=cross,
+                                            pattern_offset=1800 * u.arcsecond)
+
+    assert isinstance(positions, SkyCoord)
+    assert len(positions) == 12
+    # postion 0 should be the base position
+    assert positions[0].separation(base) < Angle(1e12 * u.degree)
+    # With no random offset positions 5, 6, 7 should be the same as 0, 1, 2
+    assert positions[0:3].to_string() == positions[5:8].to_string()
+    # Position 3 should be 30 arcminute offset from base, in declination direction only
+    assert base.spherical_offsets_to(
+        positions[3])[0].radian == pytest.approx(Angle(0 * u.degree).radian)
+    assert base.spherical_offsets_to(positions[3])[1].radian == pytest.approx(
+        Angle(-0.5 * u.degree).radian)
+    # Position 4 should be 30 arcminute offset from base in RA only.
+    assert base.spherical_offsets_to(positions[4])[0].radian == pytest.approx(
+        Angle(-0.5 * u.degree).radian)
+    assert base.spherical_offsets_to(
+        positions[4])[1].radian == pytest.approx(Angle(0 * u.degree).radian)

pocs/utils/dither.py

@@ -0,0 +1,95 @@
+import numpy as np
+import astropy.units as u
+from astropy.coordinates import SkyCoord
+from astropy.coordinates import SkyOffsetFrame
+from astropy.coordinates import ICRS
+
+# Pattern for dice 9 3x3 grid (sequence of (RA offset, dec offset) pairs)
+dice9 = ((0, 0),
+         (0, 1),
+         (1, 1),
+         (1, 0),
+         (1, -1),
+         (0, -1),
+         (-1, -1),
+         (-1, 0),
+         (-1, 1))
+
+
+# Pattern for dice 5 grid (sequence of (RA offset, dec offset) pairs)
+dice5 = ((0, 0),
+         (1, 1),
+         (1, -1),
+         (-1, -1),
+         (-1, 1))
+
+
+def get_dither_positions(base_position,
+                         n_positions,
+                         pattern=None,
+                         pattern_offset=30 * u.arcminute,
+                         random_offset=None):
+    """Create a a dithering patter for a given position.
+
+    Given a base position creates a SkyCoord list of dithered sky positions,
+    applying a dither pattern and/or random dither offsets.
+
+    Args:
+        base_position (SkyCoord or compatible): base position for the dither pattern,
+            either a SkyCoord or an object that can be converted to one by the SkyCoord
+            constructor (e.g. string).
+        n_positions (int): number of dithered sky positions to generate.
+        pattern (sequence of 2-tuples, optional): sequence of (RA offset, dec offset)
+            tuples, in units of the pattern_offset. If given pattern_offset must also
+            be specified.
+        pattern_offset (Quantity, optional): scale for the dither pattern. Should
+            be a Quantity with angular units, if a numeric type is passed instead
+            it will be assumed to be in arceconds. If pattern offset is given pattern
+            must be given too. Default 30 arcminutes.
+        random_offset (Quantity, optional): scale of the random offset to apply
+            to both RA and dec. Should be a Quantity with angular units, if numeric
+            type passed instead it will be assumed to be in arcseconds.
+
+    Returns:
+        SkyCoord: list of n_positions dithered sky positions
+
+    Raises:
+        ValueError: Raised if the `base_position` is not a valid `astropy.coordinates.SkyCoord`.
+    """
+    if not isinstance(base_position, SkyCoord):
+        try:
+            base_position = SkyCoord(base_position)
+        except ValueError:
+            raise ValueError(f"Base position '{base_position}' cannot be converted to a SkyCoord")
+
+    # Use provided pattern if given.
+    if pattern:
+        if not isinstance(pattern_offset, u.Quantity):
+            pattern_offset = pattern_offset * u.arcsec
+
+        # Get n_positions from the pattern
+        ra_offsets = [pattern[count % len(pattern)][0] for count in range(n_positions)]
+        dec_offsets = [pattern[count % len(pattern)][1] for count in range(n_positions)]
+
+        # Apply offsets to positions
+        ra_offsets *= pattern_offset
+        dec_offsets *= pattern_offset
+
+    else:
+        ra_offsets = np.zeros(n_positions) * u.arcsec
+        dec_offsets = np.zeros(n_positions) * u.arcsec
+
+    if random_offset:
+        if not isinstance(random_offset, u.Quantity):
+            random_offset = random_offset * u.arcsec
+
+        # Apply random offsets
+        ra_offsets += np.random.uniform(low=-1, high=+1, size=ra_offsets.shape) * random_offset
+        dec_offsets += np.random.uniform(low=-1, high=+1, size=dec_offsets.shape) * random_offset
+
+    offsets = SkyOffsetFrame(lon=ra_offsets, lat=dec_offsets, origin=base_position)
+    positions = offsets.transform_to(ICRS)
+
+    dither_coords = SkyCoord(positions)
+
+    return dither_coords