Convenience tool for up/downsampled spatial WCS

[keflavich]

One use case for reproject is to resample data onto a new grid to enable direct comparison between images. For objects on different parts of the sky, this comparison can be done by keeping the WCS centers the same (or specifying a WCS center) and changing the pixel size - which is different from, but somewhat related to, what find_optimal_celestial_wcs does.

Suggested feature is a tool to make a WCS with a spec like:

def resampled_wcs(original_wcs=None, pixel_scale=None, center=None, rotation=None, frame=None, shape=None):
    """
    Return a WCS with properties inherited from the original WCS unless they are overridden.
    """
    new_wcs = original_wcs.copy() # if original_wcs is not None - or load the WCS from the header, preserve size
    .... etc ...

where for resampling you would do something like:

data = fits.open(filename)
target_hdr = resampled_wcs(data[0].header, pixel_scale=0.5*u.arcsec)
reproj, _ = reproject_interp(data, target_hdr)
result = fits.PrimaryHDU(data=reproj, header=target_hdr)

cc @privong

[Cadair]

We have some functionality in ndcube which can do this kind of thing for any ape 14 WCS.

[keflavich]

could you link it here?

[Cadair]

https://github.com/sunpy/ndcube/blob/main/ndcube/wcs/wrappers/resampled_wcs.py#L7

I couldn't dig it up on my phone 🙈

These are things which were planning on being upstreamed to astropy but were instead taken into ndcube to incubate them for a bit.

[keflavich]

Thanks. That looks useful, but reproject doesn't support APE14 now, so we can't use it directly.

[keflavich]

Uh, rather, reproject does support APE14... #166... but... my example doesn't.

[Cadair]

Why doesn't your example?

[keflavich]

from ndcube.wcs.wrappers import ResampledLowLevelWCS
from spectral_cube import SpectralCube
import reproject


# read whatever random cube I have laying around
cube = SpectralCube.read('G000.00+00.00_H2CO_2pol.fits')[:25,:20,:21]
ww = ResampledLowLevelWCS(cube.wcs, [2,2,1])

reproject.reproject_interp(cube[0].hdu, ww)

gives

Traceback (most recent call last):
  File "<ipython-input-24-b22d4859483b>", line 1, in <module>
    reproject.reproject_interp(cube[0].hdu, ww)
  File "/home/adam/repos/astropy/astropy/utils/decorators.py", line 536, in wrapper
    return function(*args, **kwargs)
  File "/home/adam/anaconda3/envs/python3.9/lib/python3.9/site-packages/reproject/interpolation/high_level.py", line 79, in reproject_interp
    wcs_out, shape_out = parse_output_projection(output_projection, shape_out=shape_out,
  File "/home/adam/anaconda3/envs/python3.9/lib/python3.9/site-packages/reproject/utils.py", line 129, in parse_output_projection
    raise TypeError('output_projection should either be a Header, a WCS '
TypeError: output_projection should either be a Header, a WCS object, or a filename

[keflavich]

It should work:

reproject/reproject/utils.py

Line 117 in 43b0d8a

elif isinstance(output_projection, BaseHighLevelWCS):

but doesn't.

[keflavich]

Oh. reproject wants a high-level WCS, but we only have a low-level WCS. What do we do to wrap this to get a high-level WCS?

[keflavich]

aha, got it:

from ndcube.wcs.wrappers import ResampledLowLevelWCS
from spectral_cube import SpectralCube
import reproject

reproject.reproject_interp(cube[0].hdu, ww)
from astropy.wcs.wcsapi import HighLevelWCSWrapper

# read whatever random cube I have laying around
cube = SpectralCube.read('G000.00+00.00_H2CO_2pol.fits')[:25,:20,:21]
ww = ResampledLowLevelWCS(cube.wcs.celesital, 2)

reproject.reproject_interp(cube[0].hdu, HighLevelWCSWrapper(ww), shape_out=(10,10))

[keflavich]

spectral-cube's reproject only supports FITS WCS right now though

[keflavich]

and ResampledLowLevelWCS doesn't support serialization right now (no to_header methhod)