From 2a8a19bbb48c5b1ab19fe4d9f9dbedda320e798a Mon Sep 17 00:00:00 2001 From: Leo Singer Date: Tue, 23 Mar 2021 10:33:52 -0400 Subject: [PATCH] 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 | 9 ++ dorado/scheduling/scripts/tile.py | 101 +++---------------- dorado/scheduling/tesselation/__init__.py | 4 +- dorado/scheduling/tesselation/_geodesic.py | 25 +++-- dorado/scheduling/tesselation/_sinusoidal.py | 67 ++++++++++++ dorado/scheduling/tesselation/_spiral.py | 58 +++++++++++ 6 files changed, 169 insertions(+), 95 deletions(-) create mode 100644 dorado/scheduling/tesselation/_sinusoidal.py create mode 100644 dorado/scheduling/tesselation/_spiral.py diff --git a/doc/conf.py b/doc/conf.py index c1c2d1e..8d9530b 100644 --- a/doc/conf.py +++ b/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 = { diff --git a/dorado/scheduling/scripts/tile.py b/dorado/scheduling/scripts/tile.py index 241b681..3e51012 100644 --- a/dorado/scheduling/scripts/tile.py +++ b/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__': diff --git a/dorado/scheduling/tesselation/__init__.py b/dorado/scheduling/tesselation/__init__.py index 565391c..69daded 100644 --- a/dorado/scheduling/tesselation/__init__.py +++ b/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') diff --git a/dorado/scheduling/tesselation/_geodesic.py b/dorado/scheduling/tesselation/_geodesic.py index fc869be..bee247e 100644 --- a/dorado/scheduling/tesselation/_geodesic.py +++ b/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,6 +76,7 @@ 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 @@ -81,7 +84,7 @@ def geodesic(n, base='icosahedron', class_='I'): 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 = [] diff --git a/dorado/scheduling/tesselation/_sinusoidal.py b/dorado/scheduling/tesselation/_sinusoidal.py new file mode 100644 index 0000000..7fdc9a9 --- /dev/null +++ b/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 + -------- + + + 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) diff --git a/dorado/scheduling/tesselation/_spiral.py b/dorado/scheduling/tesselation/_spiral.py new file mode 100644 index 0000000..40987c1 --- /dev/null +++ b/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 + -------- + + + 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)