Add interpolation routines

@cdeil @keflavich - here are some initial reprojection routines based on my PR to astropy. This now includes a cube reprojection routine if the third axis is independent.

@keflavich - I was going to implement your bad pixel preservation logic from FITS_tools, but map_coordinates seems to do the right thing by default (preserve nans and infs). Do you think this is because scipy implemented it in the mean time? Otherwise, any feature from FITS_tools that I'm missing here?

This needs:
* [ ] Tests
* [ ] Docs
* [x] A full n-d reprojection function
For example, if `wcs_in.wcs.lng == 0` and `wcs_in.wcs.lat == 1`, you'd end up with `array.swapaxes(-2, -2).swapaxes(-1,-1)`, which is equivalent to "do nothing".