Refactor and clean up GROWTH tiling schemes

- Move functions to the dorado.scheduling.tesselation module.
- Clean up existing GROWTH methods.
- Update all methods to take a single resolution parameter, area.
- Convert docstrings to Numpydoc format.
- Add example code.
- Add to Sphinx project documentation.
- Use Astropy quantities for angle arguments.
- Write output in ECSV format.
- Remove logger.

doc/conf.py

@@ -37,6 +37,7 @@
     'matplotlib.sphinxext.plot_directive',
     'sphinx.ext.autodoc',
     'sphinx.ext.autosectionlabel',
+    'sphinx.ext.extlinks',
     'sphinx.ext.intersphinx',
     'sphinx.ext.napoleon'
 ]
@@ -62,6 +63,14 @@
 autosummary_generate = True
 
 
+# -- Options for extlinks extension ------------------------------------------
+
+extlinks = {
+    'arxiv': ('https://arxiv.org/abs/%s', 'arXiv:'),
+    'doi': ('https://doi.org/%s', 'doi:')
+}
+
+
 # -- Options for intersphinx -------------------------------------------------
 
 intersphinx_mapping = {

dorado/scheduling/scripts/tile.py

@@ -6,106 +6,37 @@
 # SPDX-License-Identifier: NASA-1.3
 #
 """Create a tesselation."""
-import logging
-from astropy.coordinates import SkyCoord
 import astropy.units as u
+from astropy.table import QTable
 import numpy as np
 
 from ligo.skymap.tool import ArgumentParser, FileType
 
-log = logging.getLogger(__name__)
+from .. import tesselation
+
+METHODS = {'geodesic': tesselation.geodesic,
+           'golden-angle-spiral': tesselation.golden_angle_spiral,
+           'sinusoidal': tesselation.sinusoidal}
 
 
 def parser():
     p = ArgumentParser()
-    p.add_argument('fovsize', default=3.3,
-                   type=float,
-                   help='in degrees.')
-    p.add_argument('scale', default=0.8,
-                   type=float,
-                   help='scale for overlap.')
-    p.add_argument('output', metavar='tess.dat', type=FileType('wb'),
-                   help='Output filename')
+    p.add_argument('--area', default='50 deg2', type=u.Quantity,
+                   help='Average area per tile')
+    p.add_argument('--method', default='geodesic', help='Tiling algorithm',
+                   choices=tuple(METHODS.keys()))
+    p.add_argument('-o', '--output', metavar='OUTPUT.ecsv',
+                   type=FileType('wb'), help='Output filename')
     return p
 
 
-def tesselation_spiral_packing(FOV_size, scale=0.80):
-    """
-    Spiral-based spherical packing scheme for tiling
-    Used by GRANDMA follow-up during O3 (2004.04277)
-
-    :param FOV_size: FOV in degrees for packing
-    :param scale: scaling term for tile overlaps
-
-    :return: SkyCoord
-    """
-
-    FOV = FOV_size*FOV_size*scale
-
-    area_of_sphere = 4*np.pi*(180/np.pi)**2
-    n = int(np.ceil(area_of_sphere/FOV))
-
-    golden_angle = np.pi * (3 - np.sqrt(5))
-    theta = golden_angle * np.arange(n)
-    z = np.linspace(1 - 1.0 / n, 1.0 / n - 1, n)
-    radius = np.sqrt(1 - z * z)
-
-    coords = SkyCoord(radius, theta * u.rad, z,
-                      representation_type='cylindrical')
-    coords.representation_type = 'unitspherical'
-    return coords
-
-
-def tesselation_phi_theta_packing(FOV_size, scale=0.97):
-    """
-    Phi/Theta spherical packing scheme for tiling
-    Used by GROWTH DECam follow-up during O3 (1906.00806)
-
-    :param FOV_size: FOV in degrees for packing
-    :param scale: scaling term for tile overlaps
-
-    :return: SkyCoord
-    """
-
-    sphere_radius = 1.0
-    circle_radius = np.deg2rad(FOV_size/2.0) * scale
-    vertical_count = int((np.pi*sphere_radius)/(2*circle_radius))
-
-    phis = []
-    thetas = []
-
-    phi = -0.5*np.pi
-    phi_step = np.pi/vertical_count
-    while phi < 0.5*np.pi:
-        horizontal_count = int((2*np.pi*np.cos(phi)*sphere_radius) /
-                               (2*circle_radius))
-        if horizontal_count == 0:
-            horizontal_count = 1
-        theta = 0
-        theta_step = 2*np.pi/horizontal_count
-        while theta < 2*np.pi-1e-8:
-            phis.append(phi)
-            thetas.append(theta)
-            theta += theta_step
-        phi += phi_step
-    dec = np.array(np.rad2deg(phis))
-    ra = np.array(np.rad2deg(thetas))
-
-    coords = SkyCoord(ra * u.deg, dec*u.deg)
-    coords.representation_type = 'unitspherical'
-    return coords
-
-
 def main(args=None):
     args = parser().parse_args(args)
 
-    log.info('creating tesselation')
-    coords = tesselation_phi_theta_packing(args.fovsize, scale=args.scale)
-
-    fid = open(args.output.name, 'w')
-    for ii, coord in enumerate(coords):
-        fid.write('%d %.5f %.5f\n' % (ii, coord.ra.deg, coord.dec.deg))
-    fid.close()
+    method = METHODS[args.method]
+    coords = method(args.area)
+    table = QTable({'field_id': np.arange(len(coords)), 'center': coords})
+    table.write(args.output, format='ascii.ecsv')
 
 
 if __name__ == '__main__':

dorado/scheduling/tesselation/__init__.py

@@ -7,5 +7,7 @@
 #
 """Methods for tesselating the sky into survey tiles."""
 from ._geodesic import geodesic
+from ._spiral import golden_angle_spiral
+from ._sinusoidal import sinusoidal
 
-__all__ = ('geodesic',)
+__all__ = ('geodesic', 'golden_angle_spiral', 'sinusoidal')

dorado/scheduling/tesselation/_geodesic.py

@@ -9,6 +9,7 @@
 
 from anti_lib_progs.geodesic import get_poly, grid_to_points, make_grid
 from astropy.coordinates import SkyCoord
+from astropy import units as u
 import numpy as np
 
 
@@ -45,19 +46,20 @@ def solve_number_of_vertices(n, base, class_):
     return base_count * t + 2, t, b, c
 
 
-def geodesic(n, base='icosahedron', class_='I'):
+def geodesic(area, base='icosahedron', class_='I'):
     """Generate a geodesic polyhedron with the fewest vertices >= `n`.
 
     Parameters
     ----------
-    n : int
-        The target number of vertices.
+    area : :class:`astropy.units.Quantity`
+        The average area per tile in any Astropy solid angle units:
+        for example, :samp:`10 * astropy.units.deg**2` or
+        :samp:`0.1 * astropy.units.steradian`.
     base : {``'icosahedron'``,  ``'octahedron'``, ``'tetrahedron'``}
         The base polyhedron of the tesselation.
     class_ : {``'I'``, ``'II'``, ``'III'``}
-        The class of the geodesic polyhedron, which constrains how the faces of
-        the base polyhedron may be broken down. Class III permits the most
-        freedom.
+        The class of the geodesic polyhedron, which constrains the allowed
+        values of the number of points. Class III permits the most freedom.
 
     Returns
     -------
@@ -74,14 +76,15 @@ def geodesic(n, base='icosahedron', class_='I'):
     .. plot::
         :context: reset
 
+        from astropy import units as u
         from matplotlib import pyplot as plt
         import ligo.skymap.plot
         import numpy as np
 
         from dorado.scheduling import tesselation
 
         n_vertices_target = 1024
-        vertices = tesselation.geodesic(n_vertices_target)
+        vertices = tesselation.geodesic(4 * np.pi * u.sr / n_vertices_target)
         n_vertices = len(vertices)
 
         ax = plt.axes(projection='astro globe', center='0d 25d')
@@ -93,7 +96,8 @@ def geodesic(n, base='icosahedron', class_='I'):
     .. plot::
         :context: close-figs
 
-        vertices = tesselation.geodesic(n_vertices_target, class_='II')
+        vertices = tesselation.geodesic(4 * np.pi * u.sr / n_vertices_target,
+                                        class_='II')
         n_vertices = len(vertices)
 
         ax = plt.axes(projection='astro globe', center='0d 25d')
@@ -105,7 +109,8 @@ def geodesic(n, base='icosahedron', class_='I'):
     .. plot::
         :context: close-figs
 
-        vertices = tesselation.geodesic(n_vertices_target, class_='III')
+        vertices = tesselation.geodesic(4 * np.pi * u.sr / n_vertices_target,
+                                        class_='III')
         n_vertices = len(vertices)
 
         ax = plt.axes(projection='astro globe', center='0d 25d')
@@ -115,6 +120,8 @@ def geodesic(n, base='icosahedron', class_='I'):
         ax.grid()
 
     """
+    n = int(np.ceil(1 / area.to_value(u.spat)))
+
     # Adapted from
     # https://github.com/antiprism/antiprism_python/blob/master/anti_lib_progs/geodesic.py
     verts = []

dorado/scheduling/tesselation/_sinusoidal.py

@@ -0,0 +1,67 @@
+#
+# Copyright © 2020 United States Government as represented by the Administrator
+# of the National Aeronautics and Space Administration. No copyright is claimed
+# in the United States under Title 17, U.S. Code. All Other Rights Reserved.
+#
+# SPDX-License-Identifier: NASA-1.3
+#
+from astropy.coordinates import SkyCoord
+from astropy import units as u
+import numpy as np
+
+
+def sinusoidal(area):
+    """Generate a uniform grid on a sinusoidal equal area projection.
+
+    This is similar to what was used for GRANDMA follow-up in LIGO/Virgo
+    Observing Run 3 (O3). See :doi:`10.3847/2041-8213/ab3399`.
+
+    Parameters
+    ----------
+    area : :class:`astropy.units.Quantity`
+        The average area per tile in any Astropy solid angle units:
+        for example, :samp:`10 * astropy.units.deg**2` or
+        :samp:`0.1 * astropy.units.steradian`.
+
+    Returns
+    -------
+    coords : :class:`astropy.coordinates.SkyCoord`
+        The coordinates of the tiles.
+
+    See also
+    --------
+    <https://en.wikipedia.org/wiki/Sinusoidal_projection>
+
+    Example
+    -------
+
+    .. plot::
+
+        from astropy import units as u
+        from matplotlib import pyplot as plt
+        import ligo.skymap.plot
+
+        from dorado.scheduling import tesselation
+
+        vertices = tesselation.sinusoidal(100 * u.deg**2)
+
+        ax = plt.axes(projection='astro globe', center='0d 25d')
+        ax.set_title('Sinusoidal projection grid')
+        ax.plot_coord(vertices, '.')
+        ax.grid()
+    """
+    # Diameter of the field of view
+    diameter = 2 * np.sqrt(area.to_value(u.sr) / np.pi)
+
+    # Declinations of equal-declination strips
+    n_decs = int(np.ceil(np.pi / diameter)) + 1
+    decs = np.linspace(-0.5 * np.pi, 0.5 * np.pi, n_decs)
+
+    # Number of RA steps in each equal-declination strip
+    n_ras = np.ceil(2 * np.pi * np.cos(decs) / diameter).astype(int)
+    n_ras = np.maximum(1, n_ras)
+
+    ras = np.concatenate([np.linspace(0, 2 * np.pi, n, endpoint=False)
+                          for n in n_ras])
+    decs = np.concatenate([np.repeat(dec, n) for n, dec in zip(n_ras, decs)])
+    return SkyCoord(ras * u.rad, decs * u.rad)

dorado/scheduling/tesselation/_spiral.py

@@ -0,0 +1,58 @@
+#
+# Copyright © 2020 United States Government as represented by the Administrator
+# of the National Aeronautics and Space Administration. No copyright is claimed
+# in the United States under Title 17, U.S. Code. All Other Rights Reserved.
+#
+# SPDX-License-Identifier: NASA-1.3
+#
+from astropy.coordinates import SkyCoord
+from astropy import units as u
+import numpy as np
+
+GOLDEN_ANGLE = np.pi * (3 - np.sqrt(5)) * u.rad
+
+
+def golden_angle_spiral(area):
+    """Generate a tile grid from a spiral employing the golden angle.
+
+    This is a spiral-based spherical packing scheme that was used by GRANDMA
+    during LIGO/Virgo O3 (see :doi:`10.1093/mnras/staa1846`).
+
+    Parameters
+    ----------
+    area : :class:`astropy.units.Quantity`
+        The average area per tile in any Astropy solid angle units:
+        for example, :samp:`10 * astropy.units.deg**2` or
+        :samp:`0.1 * astropy.units.steradian`.
+
+    Returns
+    -------
+    coords : :class:`astropy.coordinates.SkyCoord`
+        The coordinates of the tiles.
+
+    See also
+    --------
+    <https://en.wikipedia.org/wiki/Golden_angle>
+
+    Example
+    -------
+
+    .. plot::
+
+        from astropy import units as u
+        from matplotlib import pyplot as plt
+        import ligo.skymap.plot
+
+        from dorado.scheduling import tesselation
+
+        vertices = tesselation.golden_angle_spiral(100 * u.deg**2)
+
+        ax = plt.axes(projection='astro globe', center='0d 25d')
+        ax.set_title('Golden angle spiral')
+        ax.plot_coord(vertices, '.')
+        ax.grid()
+    """
+    n = int(np.ceil(1 / area.to_value(u.spat)))
+    ra = GOLDEN_ANGLE * np.arange(n)
+    dec = np.arcsin(np.linspace(1 - 1.0 / n, 1.0 / n - 1, n)) * u.rad
+    return SkyCoord(ra, dec)