diff --git a/generalized_wcs/README.md b/generalized_wcs/README.md new file mode 100644 index 0000000..0245ab4 --- /dev/null +++ b/generalized_wcs/README.md @@ -0,0 +1,11 @@ +Proposal for a WCS package + +*wcs_api.md* contains the main class. + +*coordinate_systems.md* defines the coordinate systems. + +*selector.md* contains a definition of a SelectorModel (a subclass of `astropy.modeling.Model`, to be used when a WCS requires more than one transform based on a selection criteria, e.g. for IFUs. + +*region_api.md* explains the role of regions in WCS. + +*example.md* shows how an instrument scientist would use the API. diff --git a/generalized_wcs/coordinate_systems_api.md b/generalized_wcs/coordinate_systems_api.md new file mode 100644 index 0000000..ace9cc0 --- /dev/null +++ b/generalized_wcs/coordinate_systems_api.md @@ -0,0 +1,207 @@ +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 high level description of coordinate system classes. + +[Prototype implementation (incomplete)] ( https://github.com/nden/code-experiments/blob/master/generalized_wcs_api/prototype/coordinate_systems.py) + +Note: Some of this may be replaced (or merged) with the high level coordinate class +being worked on although it may be that it is only a space coordinate class. + +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 + ---------- + num_axes : int + number of axes + ref_system : str + reference system (see subclasses) + ref_pos : str or ``ReferencePosition`` + reference position or standard of rest + units : list of units + unit for each axis + axes_names : list of str + names of the axes in this frame + name : str + name (alias) for this frame + """ + def __init__(self, num_axes, ref_system=None, ref_pos=None, units=None, axes_names=None, name=None): + """ Initialize a frame""" + + def transform_to(self, other): + """ + Transform from the current reference system to other + """ + +The following basic frames are defined as subclasses of `CoordinateFrame`: + + class TimeFrame(CoordinateFrame): + """ + Time Frame + + Parameters + ---------- + scale : str + time scale, one of ``standard_time_scale`` + format : str + time representation + obslat : float or ``astropy.coordinates.Latitude`` + observer's latitude + obslon : float or ``astropy.coordinates.Longitude`` + observer's longitude + units : str or ``astropy.units.Unit`` + unit for each axis + axes_names : list of str + names of the axes in this frame + """ + + standard_time_scale = time.Time.SCALES + + time_format = time.Time.FORMATS + + def __init__(self, scale, format, obslat=0., obslon=0., units="s", axes_names="Time", name="Time"): + + + def transform_to(self, val, tosys, val2=None, precision=None, in_subfmt=None, + out_subfmt=None, copy=False): + """Transforms to a different scale/representation""" + + class SkyFrame(CoordinateFrame): + """ + Sky Frame Representation + + Parameters + ---------- + reference_system : str + one of standard_ref_frame + reference_position : str, ReferencePosition instance + obstime : time.Time instance + observation time + equinox : time.Time + equinox + ##TODO? List of supported projections per frame + projection : str + projection type + axes_names : iterable or str + axes labels + units : str or units.Unit instance or iterable of those + units on axes + distance : Quantity + distance from origin + see ``~astropy.coordinates.distance` + obslat : float or ~astropy.coordinates.Latitute` + observer's latitude + obslon : float or ~astropy.coordinates.Longitude + observer's longitude + """ + standart_ref_frame = ["FK4", "FK5", "ICRS", '...'] + + def __init__(self, reference_system, reference_position, obstime, equinox=None, + projection="", axes_names=["", ""], units=["",""], distance=None, + obslat=None, obslon=None, name=None): + + def transform_to(self, lat, lon, other, distance=None): + """ + Transform from the current reference system to other. + """ + + + class SpectralFrame(CoordinateFrame): + """ + Represents Spectral Frame + + Parameters + ---------- + reference_frame : str + one of spectral_ref_frame + reference_position : ``ReferencePosition`` + obstime : time.Time instance + units : str or units.Unit instance + units for the spectral frame + axes_names : str + spectral axis name + If None, reference_system is used + rest_freq : float or Quantity (default None) + rest frequency + obslat : float or instance of `~astropy.coordinates.Latitude` + observer's latitude + obslon : float or instanceof `~astropy.coordinates.Longitude` + observer's longitude + """ + + spectral_ref_frame = ["WAVE", "FREQ", "ENER", "WAVEN", "AWAV", "VRAD", + "VOPT", "ZOPT", "VELO", "BETA"] + + standard_reference_position = ["GEOCENTER", "BARYCENTER", "HELIOCENTER", + "TOPOCENTER", "LSR", "LSRK", "LSRD", + "GALACTIC_CENTER", "MOON", "LOCAL_GROUP_CENTER"] + + def __init__(self, reference_system, reference_position, units, axes_names=None, + obstime=None, rest_freq=None, obslat=None, obslon=None, name=None): + + def transform_to(self, coord, other): + """ Transforms to other spectral coordinate system""" + + class CartesianFrame(CoordinateFrame): + def __init__(self, num_axes, reference_position, units, projection=None, name=None, axes_names=None): + ##TODO: reference position + super(CartesianFrame, self).__init__(num_axes, ref_pos=reference_position, units=units, + axes_names=axes_names, name=name) + + 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. + + Parameters + ---------- + name : string + a user defined name + frames : list + list of frames [Time, Sky, Spectral] + """ + def __init__(self, frames, name="CompositeSystem"): + self.frames = frames + + + 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 + diff --git a/generalized_wcs/design.md b/generalized_wcs/design.md new file mode 100644 index 0000000..af2ad6b --- /dev/null +++ b/generalized_wcs/design.md @@ -0,0 +1,7 @@ +The goal is to provide a framework for WCS transformations which is general and flexible enough that it is easy to extend by adding new transformations or coordinate classes. + +The WCS class is represented as a directed graph whose nodes are coordinate systems and edges are transformations +between them. This allows one to specify multiple coordinate systems and the transformations between them. +(The [prototype](https://github.com/nden/code-experiments/blob/master/generalized_wcs_api/prototype/wcs.py) uses [networkx](http://networkx.github.io).) + +Further advantage of this approach is that this can be extended by the transformations in the celestial coordinates package `astropy.coordinates` which is also represented by a graph and thus providing access to its claasses and transformations. diff --git a/generalized_wcs/examples.md b/generalized_wcs/examples.md new file mode 100644 index 0000000..eecffb5 --- /dev/null +++ b/generalized_wcs/examples.md @@ -0,0 +1,99 @@ +These are examples of how an instrument scientist can use the package to construct a WCS for an instrument. All examples assume the WCS information was already read in from somewhere (either a FITS header or a JSON file). This does not deal with WCS serialization issues. + + +JSON files used in the examples: + +[dist_info](https://github.com/nden/code-experiments/blob/master/generalized_wcs_api/prototype/jwst_example/reference_files/spec_regions.json) - coonstains models for distortion corrections + + +[spec_info](https://github.com/nden/code-experiments/blob/master/generalized_wcs_api/prototype/jwst_example/reference_files/spec_wcs.json) - contsains spectral models + + +[regions_def](https://github.com/nden/code-experiments/blob/master/generalized_wcs_api/prototype/jwst_example/reference_files/regions_miri.json) - constains definitions of IFU regions as polygon vertices in pixel space + +[wcs_regions](https://github.com/nden/code-experiments/blob/master/generalized_wcs_api/prototype/jwst_example/reference_files/wcs_regions.json) - contains basic WCS for each IFU region as offsets from a primary region. The WCS for the primary region is in the FITS header + +[spec_regions](https://github.com/nden/code-experiments/blob/master/generalized_wcs_api/prototype/jwst_example/reference_files/spec_regions.json) - contains spectral models for all regions + + + + class ImagingWCS(WCS): + """ + wcs_info is a dict of basic FITS WCS keywords from the FITS header + + """ + def __init__(wcs_info, distortion_info): + dist = self.create_distortion_transformstion(distortion_info) + linear_trans = self.create_linear_wcs_transformation(wcs_info) + coord_sys = self.create_coordinate_system(wcs_info) + foc2sky = self.create_foc2sky(wcs_info) + trans = modeling.SerialCompositeModel([dist, linear_trans, foc2sky_trans]) + super(ImagingWCS, self).__init__(trans, coord_sys) + + def create_linear_wcs_transformation(self, wcs_info): + """Creates a composite transformation of a shift, rotation and scaling""" + return linear_trans + + def create_distortion_transformation(self, dist_info): + """ + Creates a distortion transformation. + + An example of a JSON representation. + """ + return dist_trans + + + def create_foc2sky(self, wcs_info): + """ Create a transformation from focal plane to sky - deprojection + celestial rotation""" + return foc2sky + + def create_coordinate_system(self, wcs_info): + """Construts an instance of wsc.CoordinateSystem""" + return coord_sys + + +Spectral WCS: + + class SpectralWCS(WCS): + + def __init__(self, spec_info): + trans = self.create_spectral_transform(spec_info) + coord_sys = self.create_coordinate_system(wcs_info) + super(SpectralWCS, self).__init__(trans, coord_sys) + + def create_coordinate_system(self, wcs_info): + """Construts an instance of wsc.CoordinateSystem""" + return coord_sys + + def create_spectral_transform(spec_info): + """ Constructs a transform from teh model(s) in sepc_info""" + return spec_model + +An examplel of a composite spectral transform + + class CompositeSpectralTransform(object): + + def __init__(self, wcs_info, spec_info, dist_info=None): + self.spatial_transform = ImagingTransform(wcs_info, dist_info) + self.spectral_transform = SpectralTransform(spec_info) + + def __call__(self, x, y): + alpha, delta = self.spatial_transform(x, y) + wave = self.spectral_transform(x) + return alpha, delta, wave + + +IFU Example: + + class IFUWCS(WCS): + + def __init__(self, regions_def, wcs_info, wcs_regions, spec_regions, dist_info): + trans = selector.RegionsSelector(mask_shape, regions_def, wcs_info, wcs_regions, spec_regions, dist_info ) + coord_sys = self.create_coordinate_system(wcs_info) + super(IFUWCS, self).__init__(trans, coord_sys) + + def create_coordinate_system(self, wcs_info): + """ + Construts an instance of wcs.CoordinateSystem which in this case is a composite systems of a SkyFrame and SpectralFrame. + """ + return coord_sys diff --git a/generalized_wcs/region_api.md b/generalized_wcs/region_api.md new file mode 100644 index 0000000..e8223c3 --- /dev/null +++ b/generalized_wcs/region_api.md @@ -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? + +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. + + 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: + + 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 +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 + } + } + } + } diff --git a/generalized_wcs/selector.md b/generalized_wcs/selector.md new file mode 100644 index 0000000..1897292 --- /dev/null +++ b/generalized_wcs/selector.md @@ -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 + + """ \ No newline at end of file diff --git a/generalized_wcs/wcs_api.md b/generalized_wcs/wcs_api.md new file mode 100644 index 0000000..e41d4bc --- /dev/null +++ b/generalized_wcs/wcs_api.md @@ -0,0 +1,193 @@ +WCS data model - (WIP) +-------------- + + +A [prototype](https://github.com/nden/code-experiments/tree/master/generalized_wcs_api/prototype) implementation and an +[example](https://github.com/nden/code-experiments/tree/master/generalized_wcs_api/prototype/jwst_example) of using this +proposal to construct JWST WCS. + +User interface +-------------- + +Create a WCS object which transforms from `input_coordinate_system` to `output_coordinate_system` using `transform`: + + wcsobj = WCS(transform, output_coordinate_system, input_coordinate_system) + +``input_cooridnate_system`` defaults to an instance of [DetectorFrame](https://github.com/nden/code-experiments/blob/master/generalized_wcs_api/prototype/coordinate_systems.py#L276) with units of `astropy.units.pixel`. To create a WCS which transforms from detector to sky coordinates: + + sky = coordinate_systems.SkyFrame('FK5', 'BARYCENTER', obstime, equinox, projection='TAN', + axes_names=['RA', 'DEC'], units=['deg', 'deg']) + wcsobj = WCS(transform, output_coordinate_system=sky, reference_pixel=CRPIX) + +Transform pixel coordinates to world coordinates: + + ra, dec = wcsobj(x, y) + +Transform from world coordinates to pixels: + + x, y = wcsobj.invert(ra, dec) + +`transform` may be an instance of [SelectorModel](https://github.com/nden/astropy-api/blob/generalized_wcs/generalized_wcs/selector.md) which maps transforms to other quantities, +for example regions on detector but not limited to this. In this case: + + wcsobj(x,y) + +transforms the coordinates using the correct transform. +This is handled by the `SelectorModel` class. + +One can choose a specific transform, for example the transform for region 3 of an IFU: + + wcs_region_3 = wcsobj.select(label=3) + +If `transform` is an instance of `SelectorModel`, the `invert` method maps to the `invert()` +methods of all individual transforms. + +It is possible to transform the output coordinates to a different coordinate system using other packages in astropy + + wcsobj.output_coordinate_system.transform_to(ra, dec, other_system) + +One can add additional coordinate systems and transformations between them + + focal = coordinate_systems.DetectorFrame(name="focal") + wcsobj.add_transform(fromsys=self.input_coordinate_system, tosys=focal, transform=pix2foc) + +and later use their names to do a transformation + + undistorted_coordinates = wcsobj.transform(fromsys=`detector`, tosys='focal', x, y) + +To see all coordinate systems defined by a WCS object: + + wcsobj.coordinate_systems + +The special case of FITS WCS is handled through the `FITSWCS` class which inherits from `WCS` +and uses the current `astropy.wcs` to create the coordinate system object and the transforms. + + fwcs = FITSWCS(fits_file, ext, key=' ') + +In more detail +-------------- + +* [WCS API](https://github.com/nden/astropy-api/blob/generalized_wcs/generalized_wcs/wcs_api.md#wcs-api) +* [Selector Object](https://github.com/nden/astropy-api/blob/generalized_wcs/generalized_wcs/selector.md) +* [Regions API](https://github.com/nden/astropy-api/blob/generalized_wcs/generalized_wcs/region_api.md) +* [Coordinate Systems API](https://github.com/nden/astropy-api/blob/generalized_wcs/generalized_wcs/coordinate_systems_api.md) + +WCS API +------- + +The WCS class is represented as a directed graph whose nodes are coordinate systems and edges are transformations +between them. This allows one to specify multiple coordinate systems and the transformations between them. +(The [prototype](https://github.com/nden/code-experiments/blob/master/generalized_wcs_api/prototype/wcs.py) uses `networkx`.) + + class WCS(nx.DiGraph): + """ + WCS Data Model + + Parameters + ---------- + transform : astropy.modeling.Model or a callable + a callable which performs the transformation + + output_coordinate_system : astropy.wcs.CoordinateSystem + output system + input_coordinate_system : astropy.wcs.CoordinateSystem + input coordinate system + default is `wcs.DetectorFrame` + """ + + def __init__(self, transform, output_coordinate_system, input_coordinate_system=None): + self.transform = transform + self.output_coordinate_system = output_coordinate_system + self.units = self.output_coordinate_system.units + + def select(self, label): + """ + If the transform is an instance of `SelectorModel` pick one of the available transforms + + Parameters + ---------- + label : int or str + + """ + if isinstance(self.transform, SelectorModel): + return self.transform[label] + else: + raise Error + + def add_transform(self, fromsys, tosys, transform): + """ + Add new coordinate systems and transforms between them to this WCS + """ + ##TODO: check that name is unique + self.add_edge(fromsys, tosys, transform=transform) + + + def transform_to(self, fromsys, tosys, *args): + """ + Perform a transformation between any two systems registered with this WCS. + + Parameters + ---------- + fromsys : str + name of `fromsys` coordinate system + tosys : str + name of `tosys` coordinate system + """ + + @property + def coordinate_systems(self): + return self.nodes() + + + def __call__(self, args): + """ + Performs the forward transformation pix --> world. + """ + return self.transform(x, y) + + def invert(self, args): + try: + inverse_transform = self.transform.inverse() + return inverse_transform(*args) + except NotImplementedError: + return self.transform.invert(*args) + + +FITS WCS uses `~astropy.wcs`. This takes advantage of knowledge about FITS +specific keywords and parsing the FITS header. + + class FITSWCS(WCS): + """ + Implements FITS WCS + + Uses `astropy.wcs` + + Parameters + ---------- + fits_file : string or fits.HDUList object + ext : int + extension number + key : str + alternate WCS key + + """ + def __init__(self, fits_file, ext=0, key=""): + """ + Creates an `~astropy.wcs.WCS` object and uses its `all_pix2world` + method as a transform. + Also constructs a `CoordSystem` object from it. + If available, adds additional coordinate systems and transformations between them + based on `wcs.pix2foc` and similar methods. + + """ + + def __call__(self, args): + """ + The forward transform uses `all_pix2world` method of the object + """ + + def invert(self, args): + """ + The inverse transform uses the `all_sky2pix` method. + """ +