Generalized WCS

[nden]

This PR is a proposal for a WCS data object.

Please take a look and comment on it. It would be very useful to get feedback about the proposed functionality in view of different instrumentation/use cases.

I am adding a few people who might be interested in case they are not watching this repository. I'm sure I'll miss many.

[nden]

@astrofrog @eteq @perrygreenfield @wkerzendorf @keflavich @timj @mdboom

[timj]

Thanks. I'll also ask @dsberry to take a look.

[nden]

@timj Yes, please. I couldn't find his github tag.

[keflavich]

I would choose a different name for this method or describe it more completely; I'm not entirely clear what scan is supposed to do. Maybe create_mask?

Should the returned mask be a boolean array, or is it, e.g., an array with 0's and, say, 5's?

[nden]

The name is influenced by the algorithm I used which scans the region line by line to find the pixels included. But you are right, it may not be a good name in the context of WCS.

The returned mask is a byte array with region labels and 0's for pixels which are not in any region. The region labels may be int or str.

[keflavich]

OK, that makes some sense. What is the input mask? Is it like the example below (line 107)? That one seems like it should be the output of scan.

[nden]

Yes, line 107 is the output mask. The scan method is used by the create_regions_mask function which takes the shape of the science data array and the JSON file with regions definitions. It creates the byte array with 0's initially and for each region uses its scan method to find which science array pixels it contains and sets those values in the mask to the region's label.

[keflavich]

OK, I think I get it:

R = Region(3, coordsys)
assert R.__contains__(1,1) # just to show that 1,1 is in the region, nothing else is...
mask = [[0,0,0],[0,0,0],[0,0,0]]
newmask = R.scan(mask)
newmask
[[0,0,0],[0,3,0],[0,0,0]]

So if, say, I wanted a mask that selects all pixels in a given region, I'd do:

region3_mask = (R.scan(blank_mask)==3)

[cdeil]

Is the input mask modified and returned?
Or is the returned mask a new mask or copy of the input mask and only the shape of the input mask is relevant?

Shouldn't the mask have a world coordinate system attached?
Or is the region defined in pixel coordinates?

I think this is the equivalent of the region.get_mask method in pyregion?

mask = region.get_mask(shape=(300,300))

Since it's part of the API, +1 to a more explicit name for this method.

[nden]

Yep, that's it. And the create_mask function loops over all regions.
Once the mask is created (which is a one time operation) to get the indexes of all pixels from region 3:

mask==3

scan is a relatively expensive operation but it's done once.

A prototype implementation is here:

https://github.com/nden/code-experiments/blob/master/generalized_wcs_api/prototype/region.py

[nden]

@cdeil The input mask is created in create_regions_mask which then loops over all regions. The individual region's scan method modifies the mask (no copy).The region has coordinate system attached to it. All regions in a mask should have the same coordinate system. The coordinate system is not limited to cartesian.

[keflavich]

On a quick look, I have no strong objections except that I don't really understand what scan does.

In terms of how specutils is using these, I think we want the "x axis" (called "dispersion","wavelength","frequency", etc) to have a CoordinateSystem object with a SpectralFrame, but it would be an nddata object. Does this sound right?

In terms of initializing SpectralFrame, I think units needs to be a required parameter, but date_obs should not be (since it is possible to include data averaged over multiple dates for a single spectrum).

(minor: s/standart/standard/)

[keflavich]

Also, I'm pretty excited to see the region API & definitions. I think a Region needs to have write methods to output, e.g., ds9, CASA, GAIA, and other region types. Converting between regions and tables is probably the most unnecessarily complicated thing I do on a daily basis. However, as you suggest, there is overlap with photutils, and I do want to know what a region will do with partial pixels (i.e., what does "scan" or "make_mask" do?)

[timj]

Are rest_freq and rest_wave independent values? Or alternative views of the same thing?

[keflavich]

Good point; I think there should just be one defined and the other should be a property. Also, consider renaming to reference_frequency or reference_wavelength - a reference is required for transformation to z or velocity, but a rest frequency isn't required for the frame definition. At least, that's how I've interpreted the use of these terms.

[timj]

Rest Frequency is common usage. It's the FITS terminology as well. I think it would be confusing to use a different term for this than the commonly accepted name.

[keflavich]

OK, fair enough. In GBTIDL, at least, I've seen both being used - rest frequency referring to the vacuum frequency of the targeted line, reference frequency referring to where the spectrograph had been tuned and therefore where the default zero-velocity point will be. But since that's the FITS convention (RESTFRQ, right?), we should stick with it.

[nden]

Ok , I'll make rest_frequency an attribute and rest_wavelength a property. But I think it should be an attribute and not a parameter to a transform_to_some_other_coord_system() function because the class is also meant to be an info container

[cdeil]

Can you state somewhere at the top how this relates to the existing astropy.wcs module?
Do you plan to make this "generalized WCS" object a replacement for astropy.wcs.WCS or add it independently?

[nden]

The "generalized WCS" has been discussed a lot but I don't think there's a decision about that. It may be too early to decide this. As for the relation with astropy.wcs the FITSWCS class uses it.

[wkerzendorf]

There's a couple of (relatively big) ideas muddled together in this PR and I suggest we discuss them separately in separate PRs:

1. Coordinate Systems/ Reference Frames
   I believe that's an important point that is raised here: After extending numbers with units and making a quantity - the next step is that we have a quantity or multiple quantities that have a reference frame. Currently Time and Coordinates are such constructs and I think we should use them here instead of making new classes. However, we should also think about a common API for reference frames (like Time and Coordinates that can be used by things like WCS.
2. Regions
   As @nden mentioned this is really something that should be discussed with the photutils -people as they are working on something similar (again separating this out into a different PR will help). One thing that I see being discussed here is the way to serialize it: I think the API/implementation and data formats are separate questions. We should make a class that is very useful for storing and working with regions ... and then we should think how to load and store them in different file formats. We can then even make up our own format (using XML/JSON/YAML, ....).
3. WCS
   World Coordinate Systems are nothing but transformations going from one coordinate system to another. The Model class is the perfect container for storing such transformations. In addition, it can be fit (e.g if we know pairs of x,y and ra, dec we can fit a transformation). I disagree with a FITSWCS class. Again, FITS is a storage format - some of our WCS models will be able to write down their coefficients in FITS Keywords (standardized Keywords) and some will be able to be constructed from keywords.
   Like with regions we need to define an API that is useful for doing transformations between n-dimensional arrays to physical units and then we can think about storing them.

[keflavich]

@wkerzendorf - I disagree that FITSWCS should not be a class. wcs right now has internals that are similar to / built upon FITS keywords (wcs.wcs.cd, wcs.wcs.cdelt, etc.) that I often refer to because I'm very familiar with the FITS standard. As long as FITSWCS is a subclass of WCS as in this PR, I think it's helpful to have a class that's directly tied to the FITS format. I think it's particularly useful if there are WCS formats that cannot be written as FITS keywords to have that information within python - I want to know if I'm working with an object that I will not be able to save to a FITS file.

Since FITSWCS should be essentially an interface to the underlying WCS, I think it is appropriate to have as part of the API; I prefer to think about "transformations between n-dimensional arrays" as the "behind the scenes work" and FITSWCS as the user API. So in that sense, I may be quibbling about interpretation of what is API and what is "just code", but nonetheless I prefer to have a FITS-based interface.

[wkerzendorf]

@keflavich The API we're building is in its essence a transform between pixel and world coordinates (for the lack of a better name) and vice versa. FITS is not a language for transformation, but a storage format for data. In the case of WCS there are certain conventions how to encode certain transformations in FITS header keywords. But, you can store any transformation you like in FITS header keywords - but would need to invent your own convention.

So my suggestion is to separate storage on disk and transformation. IMHO we should be implementing some of the transformation and then programming functions/methods to store these into different formats including fits. For example (as you suggested on specutils):

class LinearWCS(Poly1DModel):

    @classmethod
    def from_fits_header(cls, fits_header):
        return cls(fits_header['crval1'], fits_header['cdelt1'])

    @classmethod
    def from_mycustom_format(cls, fname):
        offset, slope = map(float, file(fname).read().strip().split(','))
        return cls(offset, slope)



    def __init__(offset, slope):
        self.degree = 1
        self.c_0 = offset
        self.c_1 = slope

   def to_fits(self, fits_header):
        fits_header['crval1'] = self.offset
        fits_header['cdelt1'] = self.slope

In this case you can see if you can write to fits by seeing if there's a to_fits. But we're conceptually separating I/O and transformation.

[astrofrog]

@wkerzendorf - I think the point about FITSWCS is that various people have been referring to the current WCS standard as FITS-WCS (see e.g. Astropy paper). I think @mdboom and @perrygreenfield have been using the term for a while. However, it doesn't mean we can't find a better name.

[dsberry]

@wkerzendorf This is exactly how AST functions - it has its own rich model for representing WCS (a "FrameSet"), and also has a separate set of "Channel" classes for exporting and importing that model in various external formats. For instance, the FitsChan class can import a set of FITSWCS headers to create an AST FrameSet, and can - sometimes - export a FrameSet to create a set of FITSWCS headers. But because AST has a much richer collection of transformations than FITSWCS allows, this export operation will fail if there is no equivalent FITSWCS construct. In this case a flag is returned by the FitsChan export method indicating that the conversion was not possible. As an alternative, the FitsChan class can be told to use a set of non-standard FITS keywords that can represent any AST FrameSet. Plus there are other channel classes that convert to and from things like STC-S, etc.

It all depends on whether you want to lock yourself into the limitations imposed by the FITS-WCS standard.

[dsberry]

Don't you also need something to indicate how the time is represented - is it an MJD, a Julian Epoch, A Besselian epoch, or what?

[eteq]

@nden - a quick question before digging too deeply: how do you envision this relating to astropy.coordinates? I think we really don't want two have two (three, if you count FITS-WCS) different ways of representing celestial coordinates and doing the relevant transforms.

[astrofrog]

This might sound a bit crazy, but we could consider renaming astropy.coordinates to astropy.celestial and then we'd have astropy.time, astropy.celestial, and in future astropy.spectrum/spectral which could then all be used to define coordinate systems.

[dsberry]

On 2 July 2013 22:37, Thomas Robitaille notifications@github.com wrote:

This might sound a bit crazy, but we could consider renaming
astropy.coordinates to astropy.celestial and then we'd have astropy.time,
astropy.celestial, and in future astropy.spectrum/spectral which could
then all be used to define coordinate systems.

—
Reply to this email directly or view it on GitHubhttps://github.com//pull/10#issuecomment-20379398
.

Would they all inherit from a common base class? There are many things that
all coordinate systems of any type have in common. And then you could
create another sub-class to hold compound coordinate systems that
encapsulate instances of the other types - e.g. a compound coordinate class
for a cube could hold a celestial coordinate system and a spectral
coordinate system.

One thing I wished we done in AST at an early stage, but which it is now
too late to do anything about, is to separate the mathematical concept of a
coordinate system from the physical domain that the coordinate system
describes. If I was doing AST again from scratch, I'd have an abstract
base class called "Frame", with concrete sub-classes like CartesianFrame
PolarFrame, SphericalFrame, CylindricalFrame, CompoundFrame, etc. These
would simply define the mathematical properties of a specific type of
coordinate system (e.g. the number of axes, metric, geometry, etc) without
associating them with a specific physical domain like the sky. I would then
have a separate set of classes that represent physical domains, such as
SkyDomain, SpectralDomain, TimeDomain, FocalPlaneDomain, PixelDomain,
CompoundFrame and (most importantly) a GenericDomain to cover anything not
covered by the others. Each of these Domains would "know" about all the
common coordinate systems used to describe its physical domain
(sky,time,etc) and how to convert between them. It would represent just one
of these coordinate systems at any one time, and would encapsulate an
instance of an appropriate Frame class to do so.

This would have made it easier to provide (for example) a class to describe
3D positions within the Solar System ("SolarSystemDomain"), for which there
are many commonly used coordinate systems, some of which are Cartesian and
some spherical.

David

[astrofrog]

It looks like one is missing, as there are only three below? Should there be one for arbitrary Frames which are neither time, spectral, or spatial?

[astrofrog]

@nden - sorry for the late comments! I think it would be worth trying to implement any revisions before the coordination meeting next week since I think this is one of the most important things to focus on post-0.3. Thanks for working on this! :)

[astrofrog]

@nden - I meant to say, this API is looking very promising to me! I think we should actually start implementing some of the code to have some examples to play around with (maybe you have already, behind the scenes?).

[astrofrog]

I think it would be nice to see more examples of using the WCS - for instance, given a WCS that maps from Equatorial to pixel coordinates, can I easily get a WCS that maps from Galactic to pixel coordinates?

[astrofrog]

One request I have - it would be nice to be able to define 'Transform' objects similarly to Matplotlib in the sense that they can be easily combined. For instance, if I have a transformation from FK4 to FK5, I could do soemthing like:

WCS_new = WCS_orig + FK4toFK5

where WCS_orig`` maps from FK4 to pixel, andWCS_new`therefore maps from``FK5 to pixel``` or something like this?

EDIT - maybe this isn't necessary if one can simply call the WCS object with any coordinate object, and the WCS figures out how to get it into pixel coordinates?

[eteq]

A comment here that also applies anywhere this sort of thing is used: these sort of all-caps abbreviations I think are a bad idea. I know historically they come from situations where bytes were more precious, but here we can make a break with the past and instead use something more readable. So for example, spectral_ref_frame might instead be: ['wavelength', 'frequency', 'energy', 'wave number', 'air wavelength', 'radio velocity', 'optical velocity', 'optical z', 'relativistic v', 'relativistic beta']. (If these are actually needed to be attributes in some cases, replace the spaces with underscores)

[nden]

These came from WCS Paper III. While I'm not opposed to changing them we have to consider that they are in common use now and will have to support them in some way.

[eteq]

Oh, I agree we have to accept them as FITS/WCS header keywords, but we do not have to actually encourage their use in code. In fact, I want to make every effort to discouraging it - this isn't fortran, so we're not limited to 72 characters any more 😉

This reflects me general attitude that we are not aiming to reproduce WCS, but rather make a better and more pythonic take on WCS.