Source code for superphot_pipeline.astrometry.transformation
"""Define a class to apply sky-to-frame transformations."""
import numpy
[docs]class Transformation:
"""A clas that applies transformations generated by anmatch."""
[docs] def __init__(self, trans_fname=None):
"""See set_transformation()."""
self.coefficients = dict(x=[], y=[])
if trans_fname is not None:
self.read_anmatch_transformation(trans_fname)
[docs] def read_anmatch_transformation(self, trans_fname):
"""
Parse the transformation to apply from the given file.
Args:
trans_fname: The name of the transformation to parse for the
transformation to apply.
Returns:
None
"""
with open(trans_fname) as trans_lines:
for line in trans_lines:
if line[0] == '#':
if line[1:].strip().startswith('2MASS:'):
#ra is perfectly reasonable for sky position
#pylint: disable=invalid-name
ra, dec = map(float, line[1:].strip().split()[1:3])
#pylint: enable=invalid-name
self.center = dict(RA=ra * numpy.pi / 180.0,
Dec=dec * numpy.pi / 180.0)
else:
key, value = map(str.strip, line.split('='))
if key == 'type':
assert value == 'polynomial'
elif key == 'order':
self.order = int(value)
elif key == 'offset':
assert value == '0, 0'
elif key == 'scale':
assert value == '1'
elif key == 'basisshift':
assert value == '0, 0'
elif (
key[0] == 'd'
and
key[1] in ['x', 'y']
and
key[2:] == 'fit'
):
self.coefficients[key[1]] = [
float(word.strip()) for word in value.split(',')
]
assert(len(self.coefficients['x'])
==
(self.order + 1) * (self.order + 2) // 2)
assert(len(self.coefficients['y'])
==
(self.order + 1) * (self.order + 2) // 2)
[docs] def gnomonic_projection(self, sources, projected):
"""
Project the given sky position to a tangent plane (gnomonic projection).
Args:
sources: The the sky position to project (should have `'RA'` and
`'Dec'` keys coordinates in degrees.
projected: A numpy array with `'xi'` and `'eta'` fields to fill
with the projected coordinates (in degrees).
Returns:
None
"""
degree_to_rad = numpy.pi / 180.0
ra_diff = (sources['RA'] * degree_to_rad - self.center['RA'])
cos_ra_diff = numpy.cos(ra_diff)
cos_source_dec = numpy.cos(sources['Dec'] * degree_to_rad)
cos_center_dec = numpy.cos(self.center['Dec'])
sin_source_dec = numpy.sin(sources['Dec'] * degree_to_rad)
sin_center_dec = numpy.sin(self.center['Dec'])
denominator = (
sin_center_dec * sin_source_dec
+
cos_center_dec * cos_source_dec * cos_ra_diff
) * degree_to_rad
projected['xi'] = (cos_source_dec * numpy.sin(ra_diff)) / denominator
projected['eta'] = (
cos_center_dec * sin_source_dec
-
sin_center_dec * cos_source_dec * cos_ra_diff
) / denominator
[docs] def evaluate_polynomial(self, projected, transformed):
"""
Apply the currently defined polynomial transfomation to (xi, eta).
Args:
projected: The projected (c.f. gnomonic_projection()) coordinates
to transform. Should have `'xi'` and `'eta'` keys containing the
projected coordinates to transform in degrees.
transformed: A numpy array with `'x'` and `'y'` fields to fill
with the transformed coordinates (in pixels).
Returns:
None
"""
transformed['x'] = transformed['y'] = 0.0
xi_to_order = 1.0
order_factorial = 1.0
coef_index = 0
for order in range(self.order + 1):
xi_term = numpy.copy(xi_to_order)
eta_term = 1.0
xi_factor = order_factorial
eta_factor = 1.0
for eta_power in range(order + 1):
for out_var in ['x', 'y']:
transformed[out_var] += (
self.coefficients[out_var][coef_index]
*
xi_term
*
eta_term
/
xi_factor
/
eta_factor
)
coef_index += 1
if eta_power < order:
xi_term /= projected['xi']
eta_term *= projected['eta']
xi_factor /= order - eta_power
eta_factor *= eta_power + 1
order_factorial *= order + 1
xi_to_order *= projected['xi']
[docs] def __call__(self, sources, save_intermediate=False):
"""
Return the x and y coordinates to which the given (ra, dec) project to.
Args:
sources: The sky position to project to frame coordinates. Should
have keys named `'RA'` and `'Dec'` in degrees.
save_intermediate: Should the intermediate coordinates be saved.
Those are the coordinates after apllying the sky projection and
before applying the polynomial transformation.
Returns:
projected: The position in pixels relative to the lower left
frame of the corner to which the given sky position projects
under the currently defined transformation, as a numpy array
with fields named `'x'` and `'y'`. If `save_intermediate` is
True, the result also contains fields called `'xi'` and `'eta'`
containing the intermediate coordinates.
"""
try:
projected_shape = sources['RA'].shape
except TypeError:
projected_shape = (1,)
#pylint false positive
#pylint: disable=no-member
projected_dtype = [('x', numpy.float64), ('y', numpy.float64)]
intermediate_dtype = [('xi', numpy.float64), ('eta', numpy.float64)]
#pylint: enable=no-member
if save_intermediate:
projected_dtype.extend(intermediate_dtype)
projected = numpy.empty(projected_shape, dtype=projected_dtype)
intermediate = (projected if save_intermediate
else numpy.empty(projected_shape, intermediate_dtype))
self.gnomonic_projection(sources, intermediate)
self.evaluate_polynomial(intermediate, projected)
return projected