GWCS slicing

[keflavich]

Is slicing of GWCS implemented, or can it be?

A specific example use case is slicing a WCS for visualization purposes with wcsaxes:

hdu = fits.open('filename.fits')[('SCI', 1)]
ww = WCS(hdu.header)
slcs = slice(50, 100), slice(100, 150)
data = hdu.data[slcs]
ax = pl.subplot(projection=ww[slcs])
ax.imshow(data)

see also:

astropy/astropy#13691

[keflavich]

This can be done:

SlicedLowLevelWCS(gwcs1, slices=[slice(35, 45), slice(55, 65)])

but could you add a __getitem__ method as per the instructions here? https://docs.astropy.org/en/latest/wcs/wcsapi.html#slicing-of-wcs-objects

[mcara]

What is the purpose of slicing? Is it changing the origin of the image (pixel) coordinate system (that is, shift pixel coordinates)?

[Cadair]

For me (and ndcube) the main purpose of slicing is when you subset an array[1] you still want to have a valid WCS for that smaller array. This is conceptually achieved most of the time by offsetting the reference pixel, however, this is not what SlicedLowLevelWCS does it wraps all the methods and properties, and adjusts the inputs based on the stored slices, for example here's what it does for axis_correlation_matrix: https://github.com/astropy/astropy/blob/main/astropy/wcs/wcsapi/wrappers/sliced_wcs.py#L308-L310

The advantage of this approach is that it works for all APE 14 compliant WCSes, the major disadvantages are that the type of the sliced WCS changes, i.e. you don't have a gwcs.WCS object any more you have a SlicedLowLevelWCS object or a HighLevelWCSWrapper object (if you want high level API), and that you can no longer serialise it in the same way (I haven't written an ASDF schema / converter for SlicedLowLevelWCS [yet]).

What I would like to see is that "we" all come to an agreement about how we want to handle slicing of WCSes to minimise the friction these issues cause. For astropy.wcs.WCS there's an even smaller set of slices where another astropy.wcs.WCS object can be returned, but this currently doesn't happen if you use SlicedLowLevelWCS directly (which I would like to address).

Obviously doing "native" gWCS slicing is harder in the completely generalised sense of any model with any frames, but I have been thinking for a while, it might be nice to provide helpers for the simpler cases such as a distortionless celestial transform etc. Which could then offer a more predictable slicing interface.

In either case, as pointed out by @keflavich if gwcs.WCS implemented a __slice__ method then it could use SlicedLowLevelWCS and HighLevelWCSWrapper maybe with some other API trickery to make it 🦆 more like a gwcs.WCS object.

---

I actually think that the best short term way to address @keflavich's original use case is in astropy core, by making SlicedLowLevelWCS better (so the canonical way to slice any APE 14 compliant WCS is with SlicedLowLevelWCS) and to add some serialisation options (both FITS and gWCS/ASDF).

However, it might be worth starting to think about if we want to start work on a WCS APE aimed at facilitating modifications of WCSes in a generic way exposing optional API to the implementations to allow a more native experience when it's supported.

---

[1] You have to limit the set of valid slices as you can't use steps with WCS as the interpretation is ambiguous.

[keflavich]

@Cadair nailed it. Short version:

1. Subsets (cutouts)
2. Resampling (this is harder / can be ambiguous)

[havok2063]

Has there been any traction on this? I'm interested in slicing gwcs for creating cutouts.

[izkgao]

@havok2063

This is a temporary solution.

I wrapped up the astropy.nddata.Cutout2D and replaced the wcs attribute with the sliced gwcs. The way I sliced gwcs is to add offsets.

from astropy.modeling import models as astmodels
from astropy.nddata import Cutout2D
from astropy.wcs import WCS
from copy import deepcopy

def slice_gwcs(wcs, slices):
    # slices order (y, x)
    wcs = deepcopy(wcs)
    xmin, xmax = slices[1].start, slices[1].stop
    ymin, ymax = slices[0].start, slices[0].stop
    shape = (ymax - ymin, xmax - xmin)

    offsets = astmodels.Shift(xmin, name='offset1') & astmodels.Shift(ymin, name='offset2')

    wcs.insert_transform('detector', offsets, after=True)
    wcs.bounding_box = ((0, shape[1] - 1), (0, shape[0] - 1))
    wcs.pixel_shape = shape[::-1]
    wcs.array_shape = shape
    return wcs

def Cutout2D_GWCS(data,
                  position,
                  size,
                  wcs=None,
                  mode='trim',
                  fill_value=np.nan,
                  copy=False):
    # Get astropy fits wcs
    fitswcs = WCS(wcs.to_fits_sip())
    # Get cutout
    cutout = Cutout2D(data=data,
                      position=position,
                      size=size,
                      wcs=fitswcs,
                      mode=mode,
                      fill_value=fill_value,
                      copy=copy)

    # Slice gwcs
    cutout.wcs = slice_gwcs(wcs, cutout.slices_original)
    return cutout

Here is an example code including some simple testing results.

from astropy import units as u
from astropy.coordinates import SkyCoord
from astropy.io import fits
from astropy.wcs import WCS
from jwst import datamodels
from jwst.assign_wcs.util import update_fits_wcsinfo

# Read image
# Set `memmap` to True for large image
# model = datamodels.open('your_data.fits', memmap=True)
wcs = model.meta.wcs

# Set up a random coordinate and size
pos = SkyCoord(ra=150.1, dec=2.26, unit='deg')
size = 15 * u.arcsec

# Get cutout
cutout = Cutout2D_GWCS(data=model.data, position=pos, size=size,
                       wcs=wcs, mode='partial', copy=True)

# Check shape
print(cutout.shape)
# (500, 500)

# Check orginal origin of the cutout
print(cutout.origin_original)
# (57252, 63672)

# Get the (ra, dec) using origin wcs
print(wcs.forward_transform(57252, 63672), wcs.forward_transform(57252+499, 63672+499))
# (150.10207691461721, 2.2579230085389055) (150.09791530128845, 2.2620812905624015)

# Get the (ra, dec) using sliced wcs
print(cutout.wcs.forward_transform(0, 0), cutout.wcs.forward_transform(499, 499))
# (150.10207691461721, 2.2579230085389055) (150.09791530128845, 2.2620812905624015)

# Create a datamodel and save
dat = datamodels.ImageModel(cutout.data)
dat.meta.wcs = cutout.wcs
update_fits_wcsinfo(dat)
dat.save('cutout.fits')

# Read fits
hdu = fits.open('cutout.fits')
fwcs = WCS(hdu[1].header)
print(fwcs.pixel_to_world_values(0, 0))
# (array(150.10207691), array(2.25792301))
hdu.close()

[havok2063]

@izkgao yes this is very helpful, thanks!