Wrapping interpolation for all-sky CAR projection

[woodmd]

I'm trying to convert an all-sky CAR projection to HEALPix. Because reproject_to_healpix calls map_coordinates with mode='constant' I get a strip of null pixels at the edges of the map. It would be nice if this method could allow wrapping the interpolation in longitude (e.g. by calling with mode='wrap'). I realize wrapping the interpolation doesn't make sense for other projections so maybe you could provide a new method or option that would provide a special treatment for all-sky CAR projections. Note that another issue is how to stop it from wrapping the interpolation in latitude. There I guess the solution is to pad the array passed to map_coordinates with one row at low and high latitudes.

[woodmd]

A quick addendum... the same issue also arises for all-sky CAR to any WCS projection. Thinking about it a bit more I realize that my proposal to use mode='wrap' will not work since values are defined at pixel centers which leaves the lower/upper half of the pixels on the boundary undefined. I guess the solution would be to pad the array passed to map_coordinates with columns that are the average of the columns at the min and max longitude. Again since this is a rather special case maybe it makes sense to have a dedicated method for handling all-sky CAR projections.

Also it seems like there is a bug in the logic used in _reproject_celestial for defining the bounding box, e.g. https://github.com/astrofrog/reproject/blob/3a0287a8c3002f7da1ee891b06f7ad1df0657ee8/reproject/interpolation/core_celestial.py#L87-L94 Shouldn't the upper bound of the slice be (nx,ny) instead of (nx-1,ny-1)?

[astrofrog]

@woodmd - thanks for reporting this issue! Could you provide a minimal example that demonstrates the problem? (just to make it easier to work on a fix)

[woodmd]

Sure. Below is a script that reproduces the issue for both WCS and HPX reprojection. In both examples I'm starting from an all-sky CAR map where I've set the pixel values to the longitude coordinate. If the projection wraps correctly the projected maps should also have pixel value equal to longitude. In both examples you can see that the projection breaks down at lon=180 where the all-sky CAR projection wraps.

Note that for my own testing I put together a quick-and-dirty solution by hacking the existing methods in reproject: https://github.com/woodmd/gammapy/blob/3a5a10a7bec73ca2ba8255802f896a2ad1233294/gammapy/maps/reproject.py#L38-L110 Rather than changing the interpolation method I just pad the input array with mode='wrap'. It's probably not the most efficient implementation since it's passing the full data array to map_coordinates but it wasn't unreasonably slow in my testing with medium-size images (2k x 1k). The more elegant solution would probably be to stitch together together slices of the image from the lower/upper edge of the map. Another option might be to roll the longitude dimension by an offset that moves the projection inside the image boundaries.

import numpy as np
from astropy.io import fits
from astropy.wcs import WCS
from reproject.interpolation import reproject_interp
from reproject.healpix import reproject_to_healpix

nxpix = 36
nypix = 18
binsz = 10.0

pars_in = {
    'NAXIS': 2, 'NAXIS1': nxpix, 'NAXIS2': nypix,
    'CTYPE1': 'RA---CAR', 'CTYPE2': 'DEC--CAR',
    'CRPIX1': (nxpix+1.)/2., 'CRVAL1': 0.0, 'CUNIT1': 'deg', 'CDELT1': -binsz,
    'CRPIX2': (nypix+1.)/2., 'CRVAL2': 0.0, 'CUNIT2': 'deg', 'CDELT2': binsz,
}

pars_out = {
    'NAXIS': 2, 'NAXIS1': 20, 'NAXIS2': 20,
    'CTYPE1': 'RA---CAR', 'CTYPE2': 'DEC--CAR',
    'CRPIX1': (20+1.)/2., 'CRVAL1': 180.0, 'CUNIT1': 'deg', 'CDELT1': -1.0,
    'CRPIX2': (20+1.)/2., 'CRVAL2': 0.0, 'CUNIT2': 'deg', 'CDELT2': 1.0,
}

wcs_in = WCS(fits.Header(pars_in))
wcs_out = WCS(fits.Header(pars_out))

xpix, ypix = np.indices((nxpix,nypix))
x, y = wcs_in.wcs_pix2world(xpix,ypix,0)

# All-sky CAR to partial-sky CAR
data_wcs, _ = reproject_interp((x.T, wcs_in), wcs_out,
                           shape_out=(pars_out['NAXIS2'],pars_out['NAXIS1']))

# All-sky CAR to all-sky HPX
data_hpx, _ = reproject_to_healpix((x.T, wcs_in), 'icrs', nside=16, nested=False)

[astrofrog]

Thanks for the example! My time is limited in the next few weeks so if anyone else wants to take a look at this, please feel free to! Otherwise I'll get to it in a few weeks.