Some elegant way to produce a WCS from an image, given a list of pixel positions and a matched list of sky positions would be very nice.

A simple version, ignoring the distortion between focal plane and detector plane, could be hand rolled by the user following the steps outlined here once the Scale and Shift models in astropy become fittable; there is a pull request for that here.

This has been merged into the main branch; coming soon to a version near you.