Skip to content
This repository was archived by the owner on May 26, 2022. It is now read-only.
Open
Show file tree
Hide file tree
Changes from 4 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
140 changes: 140 additions & 0 deletions generalized_wcs/coordinate_systems_api.md
Original file line number Diff line number Diff line change
@@ -0,0 +1,140 @@
Coordinate Systems
------------------

The proposed API for coordinate systems is based on the
[STC document](http://www.ivoa.net/documents/PR/STC/STC-20050315.html)
and [WCS Paper III](http://fits.gsfc.nasa.gov/fits_wcs.html)
and is a very high level description of coordinate system classes.

A coordinate system consists of one or more coordinate frames which
describe a reference position and a reference frame.

The base class of all frames is `CoordinateFrame`. It is kept intentionally
very simple so that it can be easily extended.

class CoordinateFrame(object):
"""
Base class for CoordinateFrames

Parameters
----------
system : string
type of the frame
num_axes : int
axes_names : list of strings
units : list of units
"""
def __init__(self, system, num_axes, axes_names=None, units=None):
""" Initialize a frame"""

def transform_to(self, other):
"""
Transform from the current reference system to other if
the system attribute of the two matches
"""

Four basic frames are defined as subclasses of `CoordinateFrame`:

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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?


class TimeFrame(CoordinateFrame):
"""
Time Frame

Parameters
----------
time_scale : time scale

reference_position :an instance of ReferencePosition

reference_direction : (optional)
"""

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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?

standard_time_scale = ["TT", "TDT", "ET", "TAI", "IAT",
"UTC", "GPS", "TDB", "TEB", "TCG",
"TCB", "LST"]

def __init__(self, time_scale, reference_position, reference_direction=None, units="s"):
assert time_scale in standard_time_scale, "Unrecognized time scale"

super("Time", numaxes=1, axes_names=['time'], units=[units])
self._time_scale =time_scale
self._reference_position = reference_position

class SkyFrame(CoordinateFrame):
"""
Space Frame Representation
"""
standart_ref_frame = ["FK4", "FK5", "ICRS", '...']

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You need the observer location for things like horizon coordinates. Plus an equinox for FK4/FK5/Ecliptic, etc. Plus an epoch of observation to allow conversion between systems that move with respect to each other (like FK4 and FK5)

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Also, there is a subtlety to do with the "transform_to" method - since in general sky coordinates systems may move with respect to each other over time, you need something to indicate which reference frame is to be considered "fixed". For instance, say you have two skyframes both representing Helioecliptic longitude and latitude (i.e. the Sun defines zero longitude) but they have different epoch of observations. If you use "transform_to" to transform a position from one SkyFrame to the other, does the position change? Since both SkyFrames represent HelioEcliptic, you could argue that the position should not change. But on the other hand, that position will change when viewed against the background of distant stars (i.e. if you convert the position from the first skyframe, to ICRS, and then convert back to the second SkyFrame). The AST library addresses this issue by means of the "AlignSystem" attribute attached to the base Frame class (see http://starlink.jach.hawaii.edu/docs/sun211.htx/node454.html).

Copy link
Copy Markdown
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@dsberry Thanks for the link. There's a lot to think about and revise. Thanks for all comments.

def __init__(self, reference_frame, reference_position, obstime, equinox=None,
projection="", axes_names=["", ""], units=["",""], offset_center=None):

super(SkyFrame, self).__init__('Space', num_axes=2, axes_names=axes_names, units=units)
self._equinox = equinox
self._projection = projection
self._offset_center = offset_center
self._reference_frame = reference_frame
self._reference_position = reference_position

class SpectralFrame(CoordinateFrame):
"""
Represents Spectral Frame

Parameters
----------
reference_frame : string
one of spectral_ref_frame
reference_position : an instance of ReferencePosition

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Need to have an observer location otherwise you can't handle TOPOCENTRIC.


"""
spectral_ref_frame = ["WAVE", "FREQ", "ENER", "WAVEN", "AWAV", "VRAD",
"VOPT", "ZOPT", "VELO", "BETA"]
def __init__(self, reference_frame, reference_pos, date_obs, rest_freq=None, rest_wave=None,
units=""):

super(SpectralFrame, self).__init__('Spec', num_axes=1, axes_names=[reference_frame], units=[""])
self._reference_frame = reference_frame
self._reference_position = reference_position
self.rest_freq = rest_freq
self.rest_wave = rest_wave

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Copy Markdown
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

self.dateobs = dateobs

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

There needs to be a standard of rest in here. LSRK, LSRD, Topo, Helio, Barycentric, etc. They are all listed in FITS Paper 3 and are critical for mm/submm observing. Technically planet frames also have to be there...

Copy link
Copy Markdown
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

See ReferencePosition below, which is basically Table 1 in the STC document. Is standart_of_rest better name?

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ah. Right so is that what STC calls the standard of rest? I had assumed Reference Position was the point on the sky that the spectrum referred to. AST calls that RefPos. You need to know which direction you are looking in so that you can work out where that is relative to the standard of rest. So if ReferencePosition is the standard of rest where do you specify the position on the sky of the spectrum itself?

A `CoordinateSystem` has one or more `CoordinateFrame` objects:

class CoordinateSystem(object):
"""
A coordinate system has one or more frames.

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You're getting this use of "coordinate system" from the STC, right? This is yet another example of "coordinate system" being an overloaded term: the STC (I believe) means it to be something like "multiple coordinate frames mixed together to create a higher dimensional space" - that is a "system of coordinates" in the same sense of the word as a "system of equations". But that's not the common usage of "coordinate system" meaning something like equatorial J2000 ir whatever. So I think we would want to try to use different terminology here. (some possible, admittedly not great, ideas: MultiFrame, SystemOfFrames, or MultiDimCoordinates)

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In common parlance, a "coordinate system" is surely any collection of axes
that can be used to identify points in some specified domain. I would have
thought this definition applies to all domains of any type or
dimensionality. So sure - you could use equatorial J2000 RA and Dec axes to
form a "coordinate system" to describe point on the 2D celestial sphere.
But it would be just as valid to describe a set of (ra,dec,frequency) axes
as a "coordinate system" that describes positions within the domain of a a
spectral cube.

David

On 7 February 2014 07:11, Erik Tollerud notifications@github.com wrote:

In generalized_wcs/coordinate_systems_api.md:

  •    def transform_to(self, other, *args):
    
  •        #use affine transform
    
  • class DetectorFrame(CartesianFrame):
  •    def **init**(self, reference_pixel, units=[u.pixel, u.pixel], name=None, reference_position="Local",
    
  •                 axes_names=None):
    
  •        super(DetectorFrame, self).**init**(2, reference_position=reference_position, units=units, name=name, axes_names=axes_names)
    
  •        self._reference_pixel =
    
    +A CoordinateSystem has one or more CoordinateFrame objects:
    +
  • class CoordinateSystem(object):
  •    """
    
  •    A coordinate system has one or more frames.
    

You're getting this use of "coordinate system" from the STC, right? This
is yet another example of "coordinate system" being an overloaded term: the
STC (I believe) means it to be something like "multiple coordinate frames
mixed together to create a higher dimensional space" - that is a "system of
coordinates" in the same sense of the word as a "system of equations". But
that's not the common usage of "coordinate system" meaning something
like equatorial J2000 ir whatever. So I think we would want to try to use
different terminology here. (some possible, admittedly not great, ideas:
MultiFrame, SystemOfFrames, or MultiDimCoordinates)

Reply to this email directly or view it on GitHubhttps://github.com//pull/10/files#r9531363
.

Copy link
Copy Markdown
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I agree with David here. I'm not sure what the common usage of "coordinate system" is among astronomers, it must be field dependent. But I'm pretty sure "coordinate system" is used with spectral cubes as well as 1D spectra. The "coordinate system" in "WCS" is not meant to be limited to space only.

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@dsberry @nden : I'm saying that in common parlance, "coordinate system" does not always mean the same thing. If I go up and down the hallway and ask different astronomers what "coordinate systems" are (I actually just did that experiment), I get answers ranging "You mean like Cartesian or polar?" to "The way you describe objects in a catalog" to "It depends on the context".

Put more in the context of this proposal's language, some people of "coordinate system" as meaning the frame, while others think of it as a collection of frames (what you mean here), while yet others are thinking of specific points in a frame ("this catalog are in the ICRS coordinate system").

All I'm saying is that not worth worrying about which is "right", but rather just pick a phrasing that doesn't have this language ambiguity. That said, it may just be there's no other better word for it...


Parameters
----------
name : string
a user defined name
frames : list
list of frames [Time, Sky, Spectral]
"""
def __init__(self, frames, name="CompositeSystem"):
self.frames = frames

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

So is this a way to define multi-dimensional coordinate Frames such as you may use for a spectral cube? Why does the list of frames not include CoordinateFrame and CoordinateSystem? That is, you may want to include basic CoordinateFrames in your multi-dimensional system, and you may want to combine composite systems together to form a composite system with higher dimensionality.

Copy link
Copy Markdown
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You are right, this is too restrictive.Will change it.


class ReferencePosition(object):
"""
Encapsulates a reference position for a CoordFrame

Table 1 in the STC document.

"""
standart_reference_positions = ["GEOCENTER", "BARYCENTER", "HELIOCENTER",
"TOPOCENTER", "LSR", "LSRK", "LSRD",
"GALACTIC_CENTER", "MOON", "LOCAL_GROUP_CENTER",
"EMBARYCENTER", "RELOCATABLE", "UNKNOWNRefPos"]

def __init__(self, frame, name):
if self._validate(frame, name):
self._name = name
else:
raise ReferencePositionError

def _validate_ref_pos(self, frame):
"""
validates that the reference position is allowed for the frame
"""
return True
163 changes: 163 additions & 0 deletions generalized_wcs/region_api.md
Original file line number Diff line number Diff line change
@@ -0,0 +1,163 @@
Regions
-------

This file describes the region api.

Regions are necessary to describe multiple WCSs in an observation,
for example IFUs. This proposal follows the Region definitions in the
[STC document](http://www.ivoa.net/documents/PR/STC/STC-20050315.html).

* Note : This has some overlap with the aperture definitions in photutils.
Should we consider merging this work?

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Eventually, yes, but don't let photutils affect the API design here, we can still change the photutils API.


A Region is defined in the context of a coordinate system. Subclasses should
define `__contains__` method and a `scan` method.

class Region(object):
"""
Base class for regions.

Parameters
-------------
rid : int or string
region ID
coordinate_system : astropy.wcs.CoordinateSystem instance or a string
in the context of WCS this would be an instance of wcs.CoordinateSysem
"""
def __init__(self, rid, coordinate_system):
self._coordinate_system = coordinate_system
self._rid = rid

def __contains__(self, x, y):
"""
Determines if a pixel is within a region.

Parameters
----------
x,y : float
x , y values of a pixel

Returns
-------
True or False

Subclasses must define this method.
"""

def scan(self, mask):
"""
Sets mask values to region id for all pixels within the region.
Subclasses must define this method.

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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?

Copy link
Copy Markdown
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Copy Markdown
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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)

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Copy Markdown
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

Copy link
Copy Markdown
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@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.


Parameters
----------
mask : ndarray
a byte array with the shape of the observation to be used as a mask

Returns
-------
mask : array where the value of the elements is the region ID or 0 (for
pixels which are not included in any region).
"""

An example of a Polygon region class

class Polygon(Region):
"""
Represents a 2D polygon region with multiple vertices

Parameters
----------
rid : string
polygon id
vertices : list of (x,y) tuples or lists
The list is ordered in such a way that when traversed in a
counterclockwise direction, the enclosed area is the polygon.
The last vertex must coincide with the first vertex, minimum
4 vertices are needed to define a triangle
coord_system : string
coordinate system

"""
def __init__(self, rid, vertices, coord_system="Cartesian"):
self._rid = rid
self._vertices = vertices
self.coord_system = coord_system

def __contains__(self, x, y):

def scan(mask):




`create_regions_mask` is a function which creates a byte array with regions labels

def create_regions_mask(mask_shape, regions_def, regions_schema=None):
"""
Given a JSON file with regions definitions and a schema, this function

- creates a byte array of shape mask_shape
- reads in and validates the regions definitions
- scans each region and marks each pixel in the output array with the region ID
- returns the mask

Given 2 regions on a detector of size 10x10 px, the mask may look like this:

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is a detail, but we have to think about what happens when two regions overlap or at least overlap a common pixel.

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@astrofrog: presumably that's why it's a byte array (I assume that means 1 bit per region)? Or maybe you're thinking of some later complications due to overlap? The only thing is, a byte array could become very (redundantly) large with, say, 1000 MOS slits! If you assume the science array is a 32-bit float, you could end up with the mask being 30x the size of the science array, which isn't reasonable just to allow for a few overlaps. Hmmm... Did I understand that correctly, Nadia?

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Actually, I think Nadia must have meant an array of bytes. I saw "byte array" and thought "bit array". So @astrofrog's point stands...


0 0 0 0 0 0 0 0 0 0
1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 2 2
2 2 2 2 2 2 2 2 2 2
2 2 2 2 2 2 2 2 2 2
2 2 0 0 0 0 2 2 2 2
2 0 0 0 0 0 0 0 2 2

"""

JSON is used for regions definitions. The schema for a Set of Regions and a Polygon region is

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The JSON schema is used for region serialization to string?
In that case can you also show an example or two for a region JSON string and not only a schema example?

Are we inventing our own region format or is there some reference for this?
At the moment I think ds9 regions are very much in use ... do we want to support them somehow, maybe by a from_ds9 parser or something?

Copy link
Copy Markdown
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It is used for a serialization to a file. An example of a regions file defining two polygons is here:

https://github.com/nden/code-experiments/blob/master/generalized_wcs_api/prototype/polygon_region_example.json

The schema follows the definition of regions in the STC document which really defines their format. It is simply the implementation of this format in JSON. I could have used XML or ds9 definitions but for various reasons I chose JSON. It doesn't really matter. As I said in the beginning of wcs_api this data model has a second part to it which I am not addressing here, namely its serialization to a file. JSON seems as a plausible choice as anything else but I don't want to divert the discussion to serialization because it's really a separate issue.

Regions are used in many contexts and I think eventually we probably would have to support multiple serializations including ds9 regions.

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

+1 to a way to convert to/from ds9 regions, but this is lower priority IMHO.

defined below. It will be expanded with other kinds of regions.

{"$schema": "http://json-schema.org/draft-03/schema#",
"title": "A set of Regions",
"type": "array",
"items":{
"type": "object",
"title": "Polygon Region",
"description": "A representation of a rectangular region on a detector",
"properties": {
"id": {
"type": ["integer", "string"],
"required": true,
"unique": true
},
"coordinate_system": {
"type": "object",
"properties":{
"name": {
"type": "string",
"enum": ["Cartesian", "Sky"]
}
},
"required": true
},
"vertices": {
"type": "array",
"description": "Array of vertices describing the Polygon",
"items": {
"type": "array",
"description": "Pairs of (x, y) coordinates, representing a vertex",
"items": {
"type": "integer"
},
"minItems": 2,
"maxItems": 2
},
"minItems": 4,
"required": true
}
}
}
}
43 changes: 43 additions & 0 deletions generalized_wcs/selector.md
Original file line number Diff line number Diff line change
@@ -0,0 +1,43 @@
SelectorModel
-------------

The `SelectorModel` provides mappings between a transform and some other quantity.
The obvious case is an IFU detector where each region maps to a different transform.
However, the goal is for this to be general and allow for other use cases, for example
multiple orders, slitless spectroscopy.

For this reason the base class defines only the mapping and additional functionality
is left for specific classes.

The `__call__` method is defined by subclasses and implements the real work of
figuring out how to match input coordinates with a transform.

class SelectorModel(object):
"""
Parameters
----------
labels : a list of strings or objects
transforms : a list of transforms
transforms match labels
"""
def __init__(self, labels, transforms):
self._selector = dict(labels, transforms)

def __call__(self, label):
raise NotImplementedError

An example of a selector of regions of an IFU is below. It has an additional attribute
`regions_map` which is a byte array with region labels. The actual definitions of the
regions are read in from a JSON file, see the [Region API](https://github.com/nden/astropy-api/blob/generalized_wcs/generalized_wcs/region_api.md)

class RegionsSelector(SelectorModel):
def __init__(self, labels, transforms, regions):
super(RegionsSelector, self).__init__(labels, transforms)
#create_regions_mask_from_json is defined in the regions API.
self.regions_map = regions.create_regions_map_from_json(mask_shape, regions, schema)

def __call__(self, x, y):
"""
Transforms x and y using self._selector and self.regions_map

"""
Loading