Add WCS.celestial examples to the docs

[keflavich]

Sometimes users will encounter this error:

NotImplementedError: Currently only data with a 2-d celestial WCS can be reprojected using flux-conserving algorithm

when trying to reproject 2D data with extra dimensions in the WCS. The solution is very easy, just use wcs.celestial, i.e., the celestial subset of the WCS, instead of the full 3D or 4D WCS. There's presently no example of this in the reproject docs (at least as far as I could see by searching), so we should add one.

[keflavich]

A typical example use case is like this:

fh1 = fits.open(file1)
fh2 = fits.open(file2)
w1 = wcs.WCS(fh1[0].header)
w2 = wcs.WCS(fh2[0].header)
repr,_ = reproject_interp(input_data=(fh1[0].data.squeeze(), w1.celestial), 
                          output_projection=w2.celestial,
                          shape_out=fh2[0].data.squeeze().shape)
hdu = fits.PrimaryHDU(data=repr, header=w2.celestial.to_header()) 
# or
hdu = fits.PrimaryHDU(data=repr, header=fh2[0].header)
hdu.writeto(file3)

[keflavich]

the example I put up needs some explanation of why we reproject to w2.celestial and specify shape instead of reprojecting to fh2[0].header - i.e., the .celestial is dropping the extra axes.

Also, why use header=w2.celestial vs header=fh2[0].header? The former will preserve the data shape but will abandon all header parameters other than WCS parameters. The latter will keep the header exactly untouched, and will therefore add dummy dimensions to the data if there were any in the original data (if the original data had naxis* = [500,500,1,1], then the second version of hdu.write will have that shape)