# See LICENSE.rst for BSD 3-clause license info
# -*- coding: utf-8 -*-
"""
desimodel.focalplane.gfa
========================
Guide-Focus-Alignment location routines
"""
import numpy as np
from matplotlib.path import Path
import astropy.table
from ..io import load_deviceloc, load_gfa, load_platescale, load_tiles
from ..footprint import find_points_in_tiles, find_points_radec, get_tile_radec
from .geometry import get_radius_deg, radec2xy, qs2xy
class GFALocations(object):
def __init__(self, gfatable=None, scale=1.0):
'''
Wrapper class for GFA locations
Options:
gfatable: table with GFA corners in X, Y and either GFA_LOC or PETAL
scale: scale factor for GFA size
Notes:
* GFA coordinates are identified by a positive integer GFA_LOC or PETAL
and four rows of corners specified in the X,Y tangent plane to the focal surface.
* Default loads the GFA locations from desimodel.io.load_gfa().
scale>1 enables selecting targets that are just off of the GFA to assist
with initial acquisition.
'''
from matplotlib.path import Path
if gfatable is None:
gfatable = load_gfa()
if 'GFA_LOC' in gfatable.dtype.names:
gfa_locations = gfatable['GFA_LOC']
elif 'PETAL' in gfatable.dtype.names:
gfa_locations = gfatable['PETAL']
else:
raise ValueError('gfatable must have GFA_LOC or PETAL column')
self.gfatable = gfatable
self.gfa_polygons = list()
self.gfa_locations = np.array(sorted(set(gfa_locations)))
self.max_radius_mm = 0.0
for loc in self.gfa_locations:
ii = (gfa_locations == loc)
xcorners = gfatable['X'][ii]
ycorners = gfatable['Y'][ii]
if scale != 1.0:
x0, y0 = np.mean(xcorners), np.mean(ycorners)
xcorners = x0 + scale*(xcorners-x0)
ycorners = y0 + scale*(ycorners-y0)
r = np.sqrt(xcorners**2 + ycorners**2)
self.max_radius_mm = max(self.max_radius_mm, np.max(r))
self.gfa_polygons.append(Path(list(zip(xcorners, ycorners))))
self.max_radius_deg = float(get_radius_deg(self.max_radius_mm, 0.0))
def xy_on_gfa(self, gfa_loc, x, y):
'''
Return boolean array of whether points (x,y) are on GFA at GFA_LOC
Args:
gfa_loc: integer GFA location [0-9] = PETAL_LOC for DESI
x, y: array of focal tangent plane locations in mm
Returns boolean array of whether (x,y) is on GFA at location GFA_LOC
'''
if gfa_loc not in self.gfa_locations:
raise ValueError('GFA_LOC {} not in {}'.format(gfa_loc, self.gfa_locations))
scalar_inputs = np.isscalar(x)
x = np.atleast_1d(x)
y = np.atleast_1d(y)
xy = np.vstack([x, y]).T
i = np.where(self.gfa_locations == gfa_loc)[0][0]
inside = self.gfa_polygons[i].contains_points(xy)
if scalar_inputs:
return inside[0]
else:
return inside
def qs_on_gfa(self, gfa_loc, q, s):
'''
Return boolean array of whether points (q,s) are on GFA at GFA_LOC
Args:
gfa_loc: integer GFA location ID
q: angular coordinate in degrees
s: radial coordinate along focal surface in mm
Returns boolean array of whether (q,s) is on GFA at location GFA_LOC
'''
x, y = qs2xy(q, s)
return self.xy_on_gfa(gfa_loc, x, y)
def targets_on_gfa(self, telra, teldec, targets):
'''
Returns subset of targets table with new GFA_LOC column
Args:
telra, teldec: Telescope pointing (RA,dec) in degrees
targets: table with columns RA, DEC
Returns gfa_targets table: subset of targets with new GFA_LOC column
'''
#- Trim to targets that might be on this pointing
ps = load_platescale()
rmax = ps['theta'][-1]
columns = targets.dtype.names
if ('RA' in columns) and ('DEC' in columns):
ra = targets['RA']
dec = targets['DEC']
elif ('TARGET_RA' in columns) and ('TARGET_DEC' in columns):
ra = targets['TARGET_RA']
dec = targets['TARGET_DEC']
else:
raise ValueError('input targets must have columns RA,DEC or TARGET_RA,TARGET_DEC; has columns {}'.format(targets.dtype.names))
ii = find_points_radec(telra, teldec, ra, dec, radius=rmax)
targets = targets[ii]
ra = ra[ii]
dec = dec[ii]
#- Identify targets on GFAs
x, y = radec2xy(telra, teldec, ra, dec)
gfa_loc = np.zeros(len(targets), dtype=int) - 1
keep = np.zeros(len(targets), dtype=bool)
for loc in self.gfa_locations:
ii = self.xy_on_gfa(loc, x, y)
keep |= ii
gfa_loc[ii] = loc
results = astropy.table.Table(targets[keep], copy=False)
results['GFA_LOC'] = gfa_loc[keep].astype(np.int16)
assert np.all(results['GFA_LOC'] >= 0)
assert np.all(np.in1d(results['GFA_LOC'], self.gfa_locations))
return results
#-------------------------------------------------------------------------
[docs]def on_gfa(telra, teldec, ra, dec, scale=1.0):
"""Checks if a target is on any of the 10 GFAs given `telra`, `teldec`
and an array of `ra` and `dec` pointings, as well as a parameter for
degrees of tolerance one would like to allow.
Parameters
----------
telra : :class:`float`
The telescope's arbitrary RA pointing.
teldec : :class:`float`
The telescope's arbitrary Dec pointing.
ra : array-like
An array of RA values for locations in the sky.
dec : array-like
An array of Dec values for locations in the sky.
scale : :class:`float`, optional
Scale factor for GFA size to allow slightly off the edge targets
Returns
-------
array
An array of the same length as input `ra`, giving the GFA_LOC that
each (ra, dec) is on. -1 means not on a GFA.
"""
gfa = GFALocations(scale=scale)
x, y = radec2xy(telra, teldec, ra, dec)
gfaloc = np.full(len(x), -1, dtype=int)
for loc in gfa.gfa_locations:
ii = gfa.xy_on_gfa(loc, x, y)
gfaloc[ii] = loc
return gfaloc
[docs]def on_tile_gfa(tileid, targets, scale=1.0):
"""This function takes a tileid, a table of targets, and an optional
buffer_arcsec parameter to return the indices of targets lying on the GFA
as well as the GFA locations from 0-9.
Parameters
----------
tileid : :class:`int`
DESI tile ID, used to lookup telescope (RA, dec).
targets : Table
Table with columns RA, DEC.
scale : :class:`float`, optional
Scale factor for GFA size to allow slightly off the edge targets
Returns
-------
Table
Subset of input targets with new GFA_LOC column
"""
gfa = GFALocations(scale=scale)
telra, teldec = get_tile_radec(tileid)
return gfa.targets_on_gfa(telra, teldec, targets)
[docs]def get_gfa_targets(targets, rfluxlim=1000, tiles=None, scale=1.0):
"""This function takes a table of targets, as well as optional parameters
including a minimum flux in the r-band, a list of tiles, and a buffer in arcseconds
and returns a table of targets on the GFA satisfying a minimum flux_r
Parameters
----------
targets : Table
Table columns RA, DEC, FLUX_R.
rfluxlim : :class:`float`, optional
r-band flux limit; default 1000 = rmag 15.
tiles : Table, optional
Table of tiles, default to :func:`desimodel.io.load_tiles`.
scale : :class:`float`, optional
Scale factor for GFA size to allow slightly off the edge targets
Returns
-------
Table
A subset of input `targets` with additional columns:
TILEID: (integer) DESI tile ID; GFA_LOC: (integer) GFA location [0-9].
Notes
-----
* The same target could be repeated with different TILEID, GFA_LOC.
* The function returns an empty Table if no targets are on any GFAs or
of sufficient brightness.
"""
# Checks if the flux_r meets a minimum threshold
targets = astropy.table.Table(targets[targets['FLUX_R'] > rfluxlim])
if len(targets) == 0:
targets['TILEID'] = np.zeros(0, dtype=np.int32)
targets['GFA_LOC'] = np.zeros(0, dtype=np.int8)
return targets
if(tiles is None):
tiles = load_tiles()
gfa = GFALocations(scale=scale)
points = find_points_in_tiles(tiles, targets['RA'], targets['DEC'],
radius=gfa.max_radius_deg)
target_tables = list()
for i, pointlist in enumerate(points):
if pointlist:
tileid = tiles[i]['TILEID']
tiletargets = on_tile_gfa(tileid, targets, scale=scale)
if len(tiletargets) > 0:
tiletargets['TILEID'] = tileid
target_tables.append(tiletargets)
if len(target_tables)>1:
result = astropy.table.vstack(target_tables)
else:
#- Create blank table with correct dtype
result = astropy.table.Table(dtype=targets.dtype)
result.add_column(astropy.table.Column(name='TILEID', dtype='i4'))
return result