[AstroPy] Reprojecting each slice in a data cube based on another cube's WCS?

[Eric L. N. Jensen]

Hi all,

I'm working on a comparison of disk models to some ALMA synthesis imaging data, and in order to be able to compare models to data on a pixel-by-pixel basis, I'd like to be able to reproject each of the planes in a model data cube, based on the world coordinate system (WCS) information in the header of the observed data. I thought I might be able to do this with the python wrapper to Montage (using mProject), but it doesn't appear that mProject will handle 3D data (though I have posted a query to the Montage users' list as well to see if I'm missing something there). 

My files are FITS files with NAXIS = 3, where the first two axes are RA and Dec and the third axis is frequency or velocity, i.e. each slice of the third dimension is a 2-D image at a different frequency. 

The CRVAL/CDELT/CRPIX values of the third axis are already identical in both models and data, so no reprojection/resampling along the third axis is necessary - essentially I'm looking for an efficient way to spatially re-project each 2D image plane -  ideally without disassembling all the slices, reprojecting them individually, and reassembling them into a single file. 

I'm sure that python has various lower-level interpolation/resampling algorithms that I can use if necessary, but I'm hoping not to have to do the FITS header keyword parsing / recalculating / reprojecting work myself if this is already possible with existing high-level tools.  (The field of view is small, so it's really just a linear resampling problem, no complicated wide-field projection issues.)

Thanks in advance for your help with this,

Eric