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