Implementation of wcs_to_celestial_frame

This is an alternative to https://github.com/astropy/astropy/pull/2681 and uses a function-based implementation. We still need to figure out how to make this extensible, but it is a lot cleaner for dealing with the defaults, etc. so far for FITS-WCS.

* [x] extensibility
* [x] unit tests
From experimenting, it seemed that doing ``wcs = wcs.sub([WCSSUB_CELESTIAL])`` does seem to always re-arrange the axes, but should I not rely on that?