3

I've got a stationary all-sky camera and would like to fit the WCS solution (map pixels x,y to alt,az). If I identify a few stars, I can get an initial solution thusly

import numpy as np
from astropy.modeling import models, fitting
from astropy.modeling.projections import Sky2Pix_ZEA, AffineTransformation2D
from scipy.optimize import minimize    

# Stars of known x,y on chip and alt,az in sky.
star_obs = np.array([('Achernar', 3441.0, 2918.0, 49.947050461141515, 215.1280625253878),
                    ('Achernar', 3576.0, 3018.0, 43.715460044327585, 217.98734214492922),
                    ('Betelgeuse', 2123.0, 971.0, 43.319872170968644, 40.984431336638984),
                    ('Betelgeuse', 2330.0, 956.0, 47.122656796538564, 32.091831385823845),
                    ('Betelgeuse', 2677.0, 949.0, 51.30177229534061, 14.886238412655885),
                    ('Canopus', 2221.0, 2671.0, 55.59320568854928, 141.12216403321605)],
                    dtype=[('star_name', 'S10'), ('x', '<f8'), ('y', '<f8'),
                    ('alt', '<f8'), ('az', '<f8')])    

class sky2pix(object):
    def __init__(self, x, y, alt, az):
        projection = Sky2Pix_ZEA()
        self.affine = AffineTransformation2D()
        self.x = x
        self.y = y
        self.projx, self.projy = projection(az, alt)

    def __call__(self, x0):
        self.affine.translation.value = x0[0:2]
        self.affine.matrix.value = x0[2:]
        newx, newy = self.affine(self.projx, self.projy)
        residuals = np.sum((newx - self.x)**2 + (newy-self.y)**2)
        return residuals

fun = sky2pix(star_obs['x'], star_obs['y'], star_obs['alt'], star_obs['az'])
x0 = np.array([np.median(star_obs['x']), np.median(star_obs['y']), 1., 0., 0., 1.])
fit_result = minimize(fun, x0)

I'd like to make use of astropy's ability to make compound models and just say proj_aff = Sky2Pix_ZEA + AffineTransformation2D and then feed that into a astropy.modeling.fitting routine, but I can't tell how to deal with AffineTransformation2D returning 2 outputs.

I.P. Freeley
  • 885
  • 1
  • 9
  • 19
  • 1
    Have you tried using [astrometry.net](http://astrometry.net)? This doesn't directly answer your question, but it might nevertheless solve your problem. – apdnu Apr 21 '17 at 20:32
  • Yes. Turns out astrometry.net is not well suited to fisheye lenses, so it can only converge if I give it a small cut-out near the zenith. – I.P. Freeley Apr 21 '17 at 20:56

1 Answers1

0

Currently, none of the fitters in astropy.modeling can deal with more than one output (a limitation that I also ran into), but this was not clear from the documentation. Moreover, the models Sky2Pix_ZEA and AffineTransformation2D currently have the attribute fittable=False, for some reason.

I have now opened an issue for supporting multiple outputs and documented the limitation in the meantime, but I'm not sure how soon someone will be able to work on this (almost certainly not in time for the upcoming v2.0 release).

Note that it is possible to write your own fitters by sub-classing Fitter, as mentioned here, if that's of interest. A colleague has written a couple of specialized ones in only tens of lines, though there are a few layers of interfacing with a statistic function and optimizer to understand (from the existing code).

James
  • 1
  • 2