Native API

[nden]

This issue concerns the GWCS Native API and is related to #561. We are considering streamlining the API and

Two use cases are listed below for consideration. The proposed behavior is in bold and is very close to the current behavior but makes some assumptions. The strict approach is tempting in terms of software engineering and keeping a unambiguous interface but may not meet expectations. Any comments from power users?

@Cadair @astrofrog @chris-simpson @jehturner @karllark @drlaw1558

Use case: transforms do not support units

type of inputs	strict approach	accommodating approach	astropy.wcs
numerical	Return numbers	Return numbers	Return numbers
quantities	Raise an error	Convert inputs to numbers using units of the input system; return numbers (current behavior); OR return quantities using the units of the output system; print a warning	Return numbers
HLO.	Raise an error	Convert inputs to numbers using units of the input system; return numbers (current behavior); OR return quantities using the units of the output system; print a warning	Raise an error

Use case: transforms support units

type of inputs	strict approach	accommodating approach
numerical	Raise an error	Return quantities assuming units in the input frame; print a warning (current behavior)
quantities	Return quantities	Return quantities (current behavior)
HLO.	Raise an error	Return quantities (current behavior)

[drlaw1558]

I'm not sure that I fully follow the issue, but my inclination would be to go with the accommodating approach.

As an example, I'm thinking of the JWST GWCS used in datamodels. At the moment I can transform integer values in the detector frame to float values in RA/DEC and vice versa (which is super useful) without having to deal with astropy.units. It would be very annoying to have to set up an astropy quantity with both values and units just to evaluate this transform.

[chris-simpson]

Perhaps I'm not clear about what you mean by "native API", but I thought this was clarified in #561, that the output_level keyword would control what gets returned (by default being the same level as the input), at least for numbers vs HLOs (there wasn't really a differentiation between Quantity and HLO since I think it's only SkyCoord where they differ?) To gel with the philosophy of #561, I think the code has to be accommodating, but I think a transform that supports units should (by default) return a numerical answer (in the units of the output frame) if the inputs are numerical.

[nden]

@drlaw1558 Thanks for your thoughts! The motivation for this discussion is a proposed interface change in #561. My belief is it is not necessary and we have an agreement within STScI developers to proceed as specified above.

@chris-simpson I've updated the issue to clarify what I mean by "native API".

All use cases in the proposed interface are unambiguous except for one - transforms support units and inputs do not. Most of the time this will be the case of supplying pixel inputs to transforms with units. Requiring quantities, i.e. adding *u.pix to inputs will be annoying to end users. Because we propose to return quantities in this case, the invert function is unambiguous if using the return values.

[chris-simpson]

Within the confines of what's being discussed here, I think the output should always be at the same "level" as the input. Numerical inputs return numbers, Quantity inputs return Quantities, and HLO inputs return HLOs.

If the transform is unit-less, then input Quantities/HLOs should be converted to the units of the start frame before evaluating the transform. Warnings should be suppressable. Personally, I wouldn't even bother with a warning in this case because, to my way of thinking, this seems like the "natural" way for the API to operate given the pain of adding units to a Model -- defining the units in the CoordinateFrame would seem to make this unambiguous and unnecessary -- but I understand the desire to produce one.

[nden]

Within the confines of what's being discussed here, I think the output should always be at the same "level" as the input.

I can see this working, essentially honoring the type of the input, not the type of the transform. And it is closer to the current behavior. The only exception is HLO input.

An HLO input to a celestial->detector transform cannot return an HLO. At most it can return Quantities with units of pixels but this will be annoying. If we agree on this exception, i.e. return Quantities if input is HLO, I think this would work.
Alternatively, HLO input can be excluded from the native API and allowed in the APE 14 one only.

[chris-simpson]

Just to check we're on the same page, I view "HLO" as meaning "HLO where possible, otherwise Quantity" -- for a standard spectral datacube, I would consider "HLO" to mean(SkyCoord, wavelength) for example (which is one HLO and one Quantity), whereas "Quantity" means (RA, Dec, wavelength) (or whatever order is defined). So if your CF doesn't have axes that can be represented as HLOs (e.g., a pixel frame) then "HLO" means "Quantity", and that means u.pix should be returned (I know that's a bit painful, but I don't think we should make an exceptions specifically for those units).

[nden]

We are almost on the same page. HLO in this case is (SkyCood, SpectralCoord). Quantities will be number * unit, not necessarily RA, DEC objects. Is this what you are thinking?

[chris-simpson]

Ah yes, I forgot about SpectralCoord, thanks (DRAGONS doesn't use it because it's newer than our initial spectroscopy code). Is there a TimeCoord if the data being represented is a time series? I don't think so but, even if there is, there are coordinate systems where some axes can't be represented as HLOs, e.g., a datacube which is (pixel, pixel, SpectralCoord). So there can be hybrid HLO/Quantity inputs and outputs and in a fully generalized sense one could perform a transform between one frame where everything can be represented by HLOs and one where only some (or none) can.

[nden]

Time object is the HLO for time. There's also StokesCoordinate.

[chris-simpson]

Sorry to belabour the point but I still don't think you can separate HLOs and Quantities, because some axes (such as pixel axes) don't have an HLO representation.

Take a standard imaging WCS, from (x,y) pixels to (RA,Dec). The pixel coordinates can only be represented as Quantity objects but the sky coordinates can be an HLO, which is presumably what the user wants as an output. (Although you could just return two Quantity objects and let the use turn them into a SkyCoord themselves.) And the backward transform can't return an HLO if the input is a SkyCoord because there's no pixel HLO.

Or consider an "image" of reduced data from a MOS instrument where multiple spectra are dispersed along one axis. That axis can be represented as a SpectralCoord HLO but the other axis may be simply "object number" (a Quantity with u.dimensionless) which can't.

[nden]

Sorry to belabour the point but I still don't think you can separate HLOs and Quantities, because some axes (such as pixel axes) don't have an HLO representation.

Yes, this is well understood and I made a similar point above.
However, the practical fact is that astropy's APE 14 does not behave consistently with Quantities. For spatial and temporal WCSs these are not considered HLOs, while for spectral they are. I asked for clarification in astropy/astropy#18828.

That said, I think if we accept this distinction, HLOs should not be dealt with by the Native API, and should work with the APE 14 API only.

Internally, at STScI, there are developers who think the strict API should be implemented. The most unrestrictive one will work with HLOs. I am looking for a middle ground which will be least surprising to end users and will be intuitive in most use cases. The rest will be covered by documentation. I'm arguing that's the accommodating approach column in the table above.
For reference I added a column with astropy.wcs behaviour. I would be fine matching that approach as well.

None of this affects the STScI pipelines because they don't work with units internally. The data analysis tools are not affected either, because they use the APE 14 API.

What are the Gemini's use cases?

[chris-simpson]

All our transforms are unit-free. This seemed to make sense to me because putting units in there as well as having units in the CoordinateFrames looked like it could lead to ambiguity. The downside is that when get_transform() is used to extract the transform between two frames then you only get the transform so you need to be aware of the units. In practice, however, we're generally only doing these sorts of manipulations for pixel-to-pixel distortion transforms. But this is one of the things I liked about #561 that it would allow for the correct propagation of units in subsections of the full transform. Of course, it's fairly straightforward to make a new gWCS object from

[(gwcs.frame1, gwcs.get_transform(frame1, frame2)],
 (gwcs.frame2, None)]

but then one is obliged to use pixel_to_world and world_to_pixel to handle units even though the input and output frames might not be pixel and world.

DRAGONS has some "standards" for the CoordinateFrames -- for example, all spectral axes are in nanometres and the spectral axis is always the first one -- to make life easier when coding and faster when executing. (Early on, I did consider writing a bunch of helper functions that would, for example, allow you to extract only the wavelength axis from a pixel-to-world transform by inspecting the output frame and identifying which was the spectral axis, but it was just simpler to define our own standards.) All this means there are very few cases where DRAGONS uses gWCS with units (I think your PR had all those cases) because either we're using pixels (where adding u.pix is a pain) or we know what we're getting. The cases where we do are because we want a heavyweight object like a SkyCoord (which generally we don't require and which will inevitably slow down the execution), although of course it's possible to make these objects from the unit-free output.

So I think we could work with anything as long as it's clear. We're increasingly finding ourselves working with data structures where one axis is discrete (and dimensionless) rather than continuous (e.g., an echelle spectrum where one axis is wavelength and the other is the spectral order, or a cube where the third axis is "readout number") but I don't think this has any significant ramifications.

[nden]

It sounds like DRAGONS' use of gwcs is similar to STScI pipelines.

There's an incomplete PR implementing a discrete coordinate frame. We should make sure it covers all known cases.

Can you elaborate on the statement below. I didn't quite understand what #561 will enable
But this is one of the things I liked about #561 that it would allow for the correct propagation of units in subsections of the full transform.