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
I'm not sure if it would work if the axes were swapped, e.g. (dec, ra). In that case, one would end up calling ``swapaxes(-2, -1).swapaxes(-1,-2)`` which would leave the order unchanged. However, I do agree it can probably be simplified.