From a11b49711e71d1427eef632d4b644f90abae8019 Mon Sep 17 00:00:00 2001 From: Nadia Dencheva Date: Thu, 17 Nov 2016 17:39:16 -0500 Subject: [PATCH 1/5] save work --- docs/gwcs/using_wcs.rst | 32 +++++++++++++++++++++ docs/index.rst | 62 +++++++++++++++++++---------------------- gwcs/tests/test_wcs.py | 37 ++++++++++++++++++------ gwcs/utils.py | 50 +++++++++++++++++++++++++++++++++ gwcs/wcs.py | 22 ++++++++++++--- 5 files changed, 157 insertions(+), 46 deletions(-) create mode 100644 docs/gwcs/using_wcs.rst diff --git a/docs/gwcs/using_wcs.rst b/docs/gwcs/using_wcs.rst new file mode 100644 index 00000000..d8c19d2f --- /dev/null +++ b/docs/gwcs/using_wcs.rst @@ -0,0 +1,32 @@ +Using the WCS object +==================== + +Let's use the WCS object created in `Getting Started`_ in order to show the interface. + +To see what frames are defined: + + >>> print(wcsobj.available_frames) + ['detector', 'focal', 'icrs'] + +Some methods allow managing the transforms in a more detailed manner. + +Transforms between frames can be retrieved and evaluated separately. + + >>> distortion = wcsobj.get_transform('detector', 'focal') + >>> distortion(1, 2) + (13.4, 0.) + +Transforms in the pipeline can be replaced by new transforms. + + >>> new_transform = Shift(1) & Shift(1.5) | distortion + >>> wcsobj.set_transform('detector', 'focal', new_transform) + >>> wcsobj(1, 2) + (5.257230028926096, -72.53171157138964) + +A transform can be inserted before or after a frame in the pipeline. + + >>> scale = Scale(2) & Scale(1) + >>> w.insert_transform('icrs', scale, after=False) + >>> w(1, 2) + (10.514460057852192, -72.53171157138964) + diff --git a/docs/index.rst b/docs/index.rst index 830585e4..6e0fa7c2 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -1,5 +1,5 @@ GWCS Documentation -==================== +================== `GWCS `__ is a package for constructing and managing the World Coordinate System (WCS) of astronomical data. @@ -26,17 +26,17 @@ Installation `gwcs `__ requires: -- `numpy `__ 1.6 or later +- `numpy `__ 1.7 or later - `astropy `__ 1.1 or later -- `asdf `__ +- `asdf `__ 1.2.1 Getting Started --------------- -The simplest way to initialize a WCS object is to pass a `forward_transform` and an `output_frame` +The simplest way to initialize a WCS object is to pass a ``forward_transform`` and an ``output_frame`` to `~gwcs.wcs.WCS`. As an example, consider a typical basic FITS WCS of an image. The following imports are generally useful: @@ -45,49 +45,42 @@ The following imports are generally useful: >>> from astropy import coordinates as coord >>> from astropy import units as u >>> from gwcs import wcs - >>> from gwcs import coordinate_frames + >>> from gwcs import coordinate_frames as cf -The `forward_transform` is constructed as a combined model using `astropy.modeling`. The frames -can be strings or subclasses of `~gwcs.coordinate_frames.CoordinateFrame`. +The ``forward_transform`` is constructed as a combined model using `~astropy.modeling`. The frames +can be strings or subclasses of `~gwcs.coordinate_frames.CoordinateFrame`. There are additional benefits +having them as objects. >>> transform = Shift(-10.5) & Shift(-13.2) | Rotation2D(0.0023) | \ Scale(.01) & Scale(.04) | Pix2Sky_TAN() | RotateNative2Celestial(5.6, -72.05, 180) - >>> sky_frame = coordinate_frames.CelestialFrame(reference_frame=coord.ICRS(), name='icrs') - >>> w = wcs.WCS(forward_transform=transform, output_frame=sky_frame) + >>> sky_frame = cf.CelestialFrame(reference_frame=coord.ICRS(), name='icrs') + >>> wcsobj = wcs.WCS(forward_transform=transform, output_frame=sky_frame) To convert a pixel (x, y) = (1, 2) to sky coordinates, call the WCS object as a function: - >>> sky = w(1, 2) - >>> print(sky) + >>> ra, dec = wcsobj(x, y) + >>> print(ra, dec) (5.284139265842845, -72.49775640633503) -The :meth:`~gwcs.wcs.WCS.invert` method evaluates the :meth:`~gwcs.wcs.WCS.backward_transform` -if available, otherwise applies an iterative method to calculate the reverse coordinates. - - >>> w.invert(*sky) - (1.000000000009388, 2.0000000000112728) - -Some methods allow managing the transforms in a more detailed manner. +It is possible to get the result as a `~astropy.coordinates.SkyCoord` object`: -Transforms between frames can be retrieved and evaluated separately. + >>> sky_coord = wcsobj(x, y, output="numericals_plus") + >>> print(sky_coord) + - >>> dist = w.get_transform('detector', 'focal') - >>> dist(1, 2) - (0.16, 0.8) +This result can now be transformed to any other Celestial coordinate system supported by +`astropy.coordinates`. -Transforms in the pipeline can be replaced by new transforms. + >>> sky_coord.transform_to("galactic") + - >>> new_transform = Shift(1) & Shift(1.5) | distortion - >>> w.set_transform('detector', 'focal', new_transform) - >>> w(1, 2) - (5.257230028926096, -72.53171157138964) - -A transform can be inserted before or after a frame in the pipeline. +The :meth:`~gwcs.wcs.WCS.invert` method evaluates the :meth:`~gwcs.wcs.WCS.backward_transform` +if available, otherwise applies an iterative method to calculate the reverse coordinates. - >>> scale = Scale(2) & Scale(1) - >>> w.insert_transform('icrs', scale, after=False) - >>> w(1, 2) - (10.514460057852192, -72.53171157138964) + >>> wcsobj.invert(ra, dec) + (1.000000000009388, 2.0000000000112728) Using `gwcs` @@ -97,7 +90,8 @@ Using `gwcs` .. toctree:: :maxdepth: 2 - gwcs/create_wcs + gwcs/using_wcs.rst + gwcs/create_wcs.rst gwcs/selector_model.rst gwcs/wcstools.rst diff --git a/gwcs/tests/test_wcs.py b/gwcs/tests/test_wcs.py index 0fd0ca42..65dccd63 100644 --- a/gwcs/tests/test_wcs.py +++ b/gwcs/tests/test_wcs.py @@ -14,6 +14,7 @@ from .. import wcs from ..wcstools import * from .. import coordinate_frames as cf +from .. import utils from ..utils import CoordinateFrameError @@ -24,13 +25,13 @@ icrs = cf.CelestialFrame(reference_frame=coord.ICRS(), name='icrs') detector = cf.Frame2D(name='detector', axes_order=(0, 1)) focal = cf.Frame2D(name='focal', axes_order=(0, 1), unit=(u.m, u.m)) +spec = cf.SpectralFrame(name='wave', unit=[u.m, ], axes_order=(2, ), axes_names=('lambda', )) pipe = [(detector, m1), (focal, m2), (icrs, None) ] - # Test initializing a WCS def test_create_wcs(): @@ -124,12 +125,26 @@ def test_return_coordinates(): x = 1 y = 2.3 numerical_result = (26.8, -0.6) - + # Celestial frame num_plus_output = w(x, y, output='numericals_plus') assert_allclose(w(x, y), numerical_result) - assert_allclose(num_plus_output.ra.value, numerical_result[0]) - assert_allclose(num_plus_output.dec.value, numerical_result[1]) + assert_allclose(utils._get_values(w.unit, num_plus_output), numerical_result) + assert_allclose(w.invert(num_plus_output), (x, y)) assert isinstance(num_plus_output, coord.SkyCoord) + # Spectral frame + poly = models.Polynomial1D(1, c0=1, c1=2) + w = wcs.WCS(forward_transform=poly, output_frame=spec) + numerical_result = poly(y) + num_plus_output = w(y, output='numericals_plus') + assert_allclose(utils._get_values(w.unit, num_plus_output), numerical_result) + assert isinstance(num_plus_output, u.Quantity) + # CompositeFrame - [celestial, spectral] + output_frame = cf.CompositeFrame(frames=[icrs, spec]) + transform = m1 & poly + w = wcs.WCS(forward_transform=transform, output_frame=output_frame) + numerical_result = transform(x, y, y) + num_plus_output = w(x, y, y, output='numericals_plus') + assert_allclose(utils._get_values(w.unit, *num_plus_output), numerical_result) def test_from_fiducial_sky(): @@ -228,7 +243,7 @@ def setup_class(self): n2c = models.RotateNative2Celestial(phi, lon, theta, name='sky_rotation') tan = models.Pix2Sky_TAN(name='tangent_projection') - sky_cs = cf.CelestialFrame(reference_frame=coord.ICRS()) + sky_cs = cf.CelestialFrame(reference_frame=coord.ICRS(), name='sky') wcs_forward = wcslin | tan | n2c pipeline = [('detector', distortion), ('focal', wcs_forward), @@ -253,7 +268,7 @@ def test_distortion(self): def test_wcslinear(self): ra, dec = self.fitsw.wcs_pix2world(self.xv, self.yv, 1) - sky = self.wcs.get_transform('focal', 'CelestialFrame')(self.xv, self.yv) + sky = self.wcs.get_transform('focal', 'sky')(self.xv, self.yv) assert_allclose(ra, sky[0]) assert_allclose(dec, sky[1]) @@ -277,9 +292,15 @@ def test_footprint(self): assert_allclose(footprint, fits_footprint) def test_inverse(self): - sky_coord = self.wcs(1, 2) + sky_coord = self.wcs(1, 2, output="numericals_plus") with pytest.raises(NotImplementedError): - self.wcs.invert(sky_coord[0], sky_coord[1]) + self.wcs.invert(sky_coord) + + def test_back_coordinates(self): + sky_coord = self.wcs(1, 2, output="numericals_plus") + sky2foc = self.wcs.get_transform('sky', 'focal') + res = self.wcs.transform('sky', 'focal', sky_coord) + assert_allclose(res, self.wcs.get_transform('detector', 'focal')(1, 2)) def test_units(self): assert(self.wcs.unit == (u.degree, u.degree)) diff --git a/gwcs/utils.py b/gwcs/utils.py index 5618ea20..b5a93454 100644 --- a/gwcs/utils.py +++ b/gwcs/utils.py @@ -13,6 +13,7 @@ from astropy.modeling import core, projections from astropy.io import fits from astropy import coordinates as coords +from astropy import units as u # these ctype values do not include yzLN and yzLT pairs @@ -102,6 +103,39 @@ def _get_slice(axis_domain): return slice(x, y, step) +def _get_values(units, *args): + """ + Return the values of SkyCoord or Quantity objects. + + Parameters + ---------- + units : str or `~astropy.units.Unit` + Units of the wcs object. + The input values are converted to ``units`` before the values are returned. + """ + val = [] + values = [] + print('args', args) + for arg in args: + print('arg', arg) + if isinstance(arg, coords.SkyCoord): + try: + print('arg1', arg) + lon = arg.data.lon + lat = arg.data.lat + except AttributeError: + lon = arg.spherical.lon + lat = arg.spherical.lat + val.extend([lon, lat]) + elif isinstance(arg, u.Quantity): + val.append(arg) + else: + raise TypeError("Unsupported coordinate type {}".format(arg)) + for va, un in zip(val, units): + values.append(va.to(un).value) + return values + + def _compute_lon_pole(skycoord, projection): """ Compute the longitude of the celestial pole of a standard frame in the @@ -422,6 +456,22 @@ def create_projection_transform(projcode): return projklass(**projparams) +def isnumerical(val): + """ + Determine if a value is numerical (number or np.array of numbers). + """ + dtypes = ['uint64', 'float64', 'int8', 'int64', 'int16', 'uint16', 'uint8', + 'float32', 'int32', 'uint32'] + isnum = True + if isinstance(val, coords.SkyCoord): + isnum = False + elif isinstance(val, u.Quantity): + isnum = False + elif isinstance(val, np.ndarray) and val.dtype not in dtypes: + isnum = False + return isnum + + # ######### axis separability ######### # Functions to determine axis separability # The interface will change most likely diff --git a/gwcs/wcs.py b/gwcs/wcs.py index 86670c28..3279c315 100644 --- a/gwcs/wcs.py +++ b/gwcs/wcs.py @@ -225,7 +225,10 @@ def __call__(self, *args, **kwargs): result = self.forward_transform(*args) output = kwargs.pop('output', None) if output == 'numericals_plus': - result = self.output_frame.coordinates(*result) + if self.output_frame.naxes == 1: + result = self.output_frame.coordinates(result) + else: + result = self.output_frame.coordinates(*result) elif output is not None and output != "numericals": raise ValueError("Type of output unrecognized {0}".format(output)) return result @@ -239,18 +242,23 @@ def invert(self, *args, **kwargs): Parameters ---------- - args : float or array like + args : float, array like, `~astropy.coordinates.SkyCoord` or `~astropy.units.Unit` coordinates to be inverted kwargs : dict keyword arguments to be passed to the iterative invert method. """ + if not utils.isnumerical(args[0]): + args = utils._get_values(self.unit, *args) try: result = self.backward_transform(*args) except (NotImplementedError, KeyError): result = self._invert(*args, **kwargs) output = kwargs.pop('output', None) if output == 'numericals_plus': - return self.input_frame.coordinates(*result) + if self.input_frame.naxes == 1: + return self.input_frame.coordinates(result) + else: + return self.input_frame.coordinates(*result) else: return result @@ -280,12 +288,18 @@ def transform(self, from_frame, to_frame, *args, **kwargs): `~astropy.units.Quantity` object. """ transform = self.get_transform(from_frame, to_frame) + if not utils.isnumerical(args[0]): + args = utils._get_values(self.unit, *args) + result = transform(*args) output = kwargs.pop("output", None) if output == "numericals_plus": to_frame_name, to_frame_obj = self._get_frame_name(to_frame) if to_frame_obj is not None: - result = to_frame_obj.coordinates(*result) + if to_frame_obj.naxes == 1: + result = to_frame_obj.coordinates(result) + else: + result = to_frame_obj.coordinates(*result) else: raise TypeError("Coordinate objects could not be created because" "frame {0} is not defined.".format(to_frame_name)) From cb510e39f5c3d1b8ddd1fd7109e92b138f47d013 Mon Sep 17 00:00:00 2001 From: Nadia Dencheva Date: Mon, 28 Nov 2016 15:15:34 -0500 Subject: [PATCH 2/5] remove extarnwous print --- gwcs/coordinate_frames.py | 1 - 1 file changed, 1 deletion(-) diff --git a/gwcs/coordinate_frames.py b/gwcs/coordinate_frames.py index 0ae8c65c..66e9f396 100644 --- a/gwcs/coordinate_frames.py +++ b/gwcs/coordinate_frames.py @@ -321,7 +321,6 @@ def coordinates(self, *args): coo = [] for frame in self.frames: fargs = [args[i] for i in frame.axes_order] - print(frame, fargs, frame.axes_order) coo.append(frame.coordinates(*fargs)) return coo From 29c89dc23a9198fdc98d75be71b220e98479c372 Mon Sep 17 00:00:00 2001 From: Nadia Dencheva Date: Tue, 20 Dec 2016 09:36:51 -0500 Subject: [PATCH 3/5] update examples --- docs/index.rst | 119 ++++++++++++++++++++++++++++++++++++------------- 1 file changed, 87 insertions(+), 32 deletions(-) diff --git a/docs/index.rst b/docs/index.rst index 6e0fa7c2..185f019e 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -1,7 +1,7 @@ GWCS Documentation -================== +==================== -`GWCS `__ is a package for constructing and managing +`GWCS `__ is a package for managing the World Coordinate System (WCS) of astronomical data. Introduction @@ -26,20 +26,40 @@ Installation `gwcs `__ requires: -- `numpy `__ 1.7 or later +- `numpy `__ 1.67 or later -- `astropy `__ 1.1 or later +- `astropy `__ 1.2 or later -- `asdf `__ 1.2.1 +- `asdf `__ +To install from source:: + + git clone https://github.com/spacetelescope/gwcs.git + cd gwcs + python setup.py install + +To install the latest release:: + + pip install gwcs + +GWCS is also available as part of astroconda. Getting Started --------------- -The simplest way to initialize a WCS object is to pass a ``forward_transform`` and an ``output_frame`` -to `~gwcs.wcs.WCS`. As an example, consider a typical basic FITS WCS of an image. +The WCS data model represents a pipeline of transformations from some +initial coordinate frame to a standard coordinate frame. +It is implemented as a list of steps executed in order. Each step defines a +starting coordinate frame and the transform to the next frame in the pipeline. +The last step has no transform, only a frame which is the output frame of +the total transform. As a minimum a WCS object has an input_frame (defaults to "detector"), +an output_frame and the transform between them. + +As an example, consider a typical WCS of an image with some distortion. + The following imports are generally useful: + >>> import numpy as np >>> from astropy.modeling.models import (Shift, Scale, Rotation2D, Pix2Sky_TAN, RotateNative2Celestial, Mapping, Polynomial2D) >>> from astropy import coordinates as coord @@ -47,40 +67,76 @@ The following imports are generally useful: >>> from gwcs import wcs >>> from gwcs import coordinate_frames as cf -The ``forward_transform`` is constructed as a combined model using `~astropy.modeling`. The frames -can be strings or subclasses of `~gwcs.coordinate_frames.CoordinateFrame`. There are additional benefits -having them as objects. +The `forward_transform` is constructed as a combined model using `astropy.modeling`. +The frames are subclasses of `~gwcs.coordinate_frames.CoordinateFrame` (although strings are +acceptable too). - >>> transform = Shift(-10.5) & Shift(-13.2) | Rotation2D(0.0023) | \ - Scale(.01) & Scale(.04) | Pix2Sky_TAN() | RotateNative2Celestial(5.6, -72.05, 180) +First create polynomial transform to represent distortion: + + >>> polyx = Polynomial2D(4) + >>> polyx.parameters = np.random.randn(15) + >>> polyy = Polynomial2D(4) + >>> polyy.parameters = np.random.randn(15) + >>> distortion = (Mapping((0, 1, 0, 1)) | polyx & polyy).rename("distortion") + +Next create a transform from distortion fdree coordinates to ICRS. + + >>> dist2sky = (Shift(-10.5) & Shift(-13.2) | Rotation2D(0.0023) | \ + Scale(.01) & Scale(.04) | Pix2Sky_TAN() | RotateNative2Celestial(5.6, -72.05, 180)).rename("focal2sky") + +Create three coordinate frames, associated with the detector, a celestial frame, and an intermediate one. + + >>> detector_frame = cf.Frame2D(name="detector", axes_names=("x", "y"), unit=(u.pix, u.pix)) >>> sky_frame = cf.CelestialFrame(reference_frame=coord.ICRS(), name='icrs') - >>> wcsobj = wcs.WCS(forward_transform=transform, output_frame=sky_frame) + >>> focal_frame = cf.Frame2D(name="focal_frame", unit=(u.arcsec, u.arcsec)) + +Create the WCS pipeline and initialize the WCS: + + >>> pipeline = [(detector_frame, distortion), + (focal_frame, dist2sky), + (sky_frame, None) + ] + >>> wcsobj = wcs.WCS(pipeline) + >>> print(wcsobj) + From Transform + ----------- ---------- + detector distortion + focal_frame focal2sky + icrs None To convert a pixel (x, y) = (1, 2) to sky coordinates, call the WCS object as a function: - >>> ra, dec = wcsobj(x, y) - >>> print(ra, dec) - (5.284139265842845, -72.49775640633503) + >>> sky = wcsobj(1, 2) + >>> print(sky) + (5.759024831907874, -72.22626601919619) + +#The :meth:`~gwcs.wcs.WCS.invert` method evaluates the :meth:`~gwcs.wcs.WCS.backward_transform` +#if available, otherwise applies an iterative method to calculate the reverse coordinates. +# +# >>> w.invert(*sky) +# (1.000000000009388, 2.0000000000112728) +# +#Some methods allow managing the transforms in a more detailed manner. -It is possible to get the result as a `~astropy.coordinates.SkyCoord` object`: +Transforms between frames can be retrieved and evaluated separately. - >>> sky_coord = wcsobj(x, y, output="numericals_plus") - >>> print(sky_coord) - + >>> dist = w.get_transform('detector_frame', 'focal_frame') + >>> dist(1, 2) + (15.354214305518118, 8.791536957201615) -This result can now be transformed to any other Celestial coordinate system supported by -`astropy.coordinates`. +Transforms in the pipeline can be replaced by new transforms. - >>> sky_coord.transform_to("galactic") - + >>> new_transform = distortion | Shift(1) & Shift(1.5) + >>> wcsobj.set_transform('detector', 'focal', new_transform) + >>> wcsobj(1, 2) + (5.791157736884894, -72.16623599444335) -The :meth:`~gwcs.wcs.WCS.invert` method evaluates the :meth:`~gwcs.wcs.WCS.backward_transform` -if available, otherwise applies an iterative method to calculate the reverse coordinates. +A transform can be inserted before or after a frame in the pipeline. - >>> wcsobj.invert(ra, dec) - (1.000000000009388, 2.0000000000112728) + >>> scale = Scale(2) & Scale(1) + >>> wcsobj.insert_transform('icrs', scale, after=False) + >>> wcsobj(1, 2) + (11.582315473769787, -72.16623599444335) Using `gwcs` @@ -90,8 +146,7 @@ Using `gwcs` .. toctree:: :maxdepth: 2 - gwcs/using_wcs.rst - gwcs/create_wcs.rst + gwcs/create_wcs gwcs/selector_model.rst gwcs/wcstools.rst From 03771d57c13e1b7df81640068848558324e97eeb Mon Sep 17 00:00:00 2001 From: Nadia Dencheva Date: Wed, 21 Dec 2016 16:05:50 -0500 Subject: [PATCH 4/5] update documentation and tests --- docs/gwcs/create_wcs.rst | 113 ----------------------------------- docs/gwcs/selector_model.rst | 16 ++--- docs/gwcs/using_wcs.rst | 66 +++++++++++++++++--- docs/index.rst | 85 ++++++++------------------ gwcs/wcs.py | 6 ++ 5 files changed, 96 insertions(+), 190 deletions(-) delete mode 100644 docs/gwcs/create_wcs.rst diff --git a/docs/gwcs/create_wcs.rst b/docs/gwcs/create_wcs.rst deleted file mode 100644 index 6efaab53..00000000 --- a/docs/gwcs/create_wcs.rst +++ /dev/null @@ -1,113 +0,0 @@ -An imaging example -================== - - -This is a step by step example of creating a WCS for an imaging observation. -The forward transformation is from detector pixels to sky coordinates in ICRS. -There's an intermediate coordinate frame, called 'focal', after applying the distortion. - -The following packages must be imported. - - - >>> import numpy as np - >>> from astropy.modeling.models import (Shift, Scale, Rotation2D, - Pix2Sky_TAN, RotateNative2Celestial, Mapping, Polynomial2D, AffineTransformation2D) - >>> from astropy import coordinates as coord - >>> from astropy import units as u - >>> from gwcs import wcs - >>> from gwcs import coordinate_frames as cf - - -Create the transform -~~~~~~~~~~~~~~~~~~~~ - -First create a model for the distortion transform. Let's assume the distortion -in each direction is represented by a 2nd degree Polynomial model of `x` and `y`. - - >>> dist_x = Polynomial2D(2, c0_0=0.0013, c1_0=0.5, c2_0=0, c0_1=0.823, c0_2=1.4, c1_1=1.7, name='x_distortion') - >>> dist_y = Polynomial2D(2, c0_0=0.03, c1_0=0.25, c2_0=1.2, c0_1=0.3, c0_2=0.4, c1_1=0.7, name='y_distortion') - >>> distortion = dist_x & dist_y - -The last line above joins the two distortion models in one model which now takes -4 inputs (x, y, x, y), two for each model. In order for this to work a -:class:`~astropy.modeling.mappings.Mapping` model must be prepended to the ``distortion`` transform. - - >>> distortion_mapping = Mapping((0, 1, 0, 1), name='distortion_mapping') - >>> distortion = distortion_mapping | dist_x & dist_y - -Next we create the transform from focal plane to sky. For this example, suppose the WCS is in a FITS -header represented through the usual FITS WCS keywords - a point on the detector (CRPIX1/2) corresponds -to a point on a projection plane tangent to the celestial sphere (CRVAL1/2). The FITS keywords are: .. - -:: - - CRPIX1 = 2048.0 - CRPIX2 = 1024.0 - CRVAL1 = 5.63056810618 - CRVAL2 = -72.0545718428 - LONPOLE = 180 - PC1_1 = 1.29058668e-05 - PC1_2 = 5.95320246e-06 - PC2_1 = 5.02215196e-06 - PC2_2 = -1.26450104e-05 - CTYPE1 = 'RA---TAN' - CTYPE2 = 'DEC---TAN' - -.. note:: FITS WCS keywords are given simply as an example and because it's the most often - used way to represent this information. However, this code is not limited to FITS WCS. - -The WCS information above represents a serial compound model consisting of a shift in the focal plane -by CRPIX, a rotation by the PC matrix, a tangent deprojection and a sky rotation. Using -`astropy.modeling `__ this can be written as - - >>> shift = Shift(CRPIX1, name="x_shift") & Shift(CRPIX2, name="y_shift") - >>> plane_rotation = AffineTransformation2D(matrix=np.array([[PC1_1, PC1_2], [PC2_1, PC2_2]])) - >>> tangent = Pix2Sky_TAN() - >>> sky_rotation = RotateNative2Celestial(CRVAL1, CRVAL2, LONPOLE) - -Chaining these models into one compound models creates the total transformation from focal plane to sky. - - >>> focal2sky = shift | plane_rotation | tangent | sky_rotation - - -Create the coordinate frames -~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - - >>> outframe = cf.CelestialFrame(reference_frame=coord.ICRS(), name='icrs') - >>> outframe.unit - [Unit("deg"), Unit("deg")] - >>> inframe = cf.Frame2D(name='detector') - >>> print(inframe.unit) - [Unit("pix"), Unit("pix")] - >>> focal = cf.Frame2D(name='focal') - - -Create the WCS object -~~~~~~~~~~~~~~~~~~~~~ - -In this case it is convenient to initialize the WCS object with a list of tuples, -where each tuple (step_frame, step_transform) represents a transform "step_transform" -from frame "step_frame" to the next frame in the WCS pipeline. -The transform in the last step is always None to indicate end of the pipeline. - - >>> pipeline = [(inframe, distortion), (focal, focal2sky), (outframe, None)] - >>> w = wcs.WCS(pipeline) - >>> w(1, 2) - (5.736718396223817, -72.057214400243) - -Frame objects allow to extend the functionality by using `astropy.coordinates` and `astropy.units`. -Frames are available as attributes of the WCS object. - - >>> w.available_frames - ['detector', 'focal', 'icrs'] - >>> w.icrs - , unit=[Unit("deg"), Unit("deg")], name=icrs)> - >>> w.icrs.coordinates(1, 2) - - >>> w.icrs.transform_to('galactic', 1, 2) - - - - diff --git a/docs/gwcs/selector_model.rst b/docs/gwcs/selector_model.rst index abac721b..772a5e2f 100644 --- a/docs/gwcs/selector_model.rst +++ b/docs/gwcs/selector_model.rst @@ -4,7 +4,7 @@ Discontinuous WCS - An IFU Example There are a couple of models in GWCS which help with managing discontinuous WCSs. Given (x, y) pixel indices, `~gwcs.selector.LabelMapperArray` returns a label (int or str) which uniquely identifies a location in a frame. `~gwcs.selector.RegionsSelector` -maps different transforms to different locations in the frame of `~gwcs.selector.RegionsSelector.label_mapper`. +maps different transforms to different locations in the frame of `~gwcs.selector._LabelMapper`. An example use case is an IFU observation. The locations on the detector where slices are projected are labeled with integer numbers. Given an (x, y) pixel, the `~gwcs.selector.LabelMapperArray` @@ -71,7 +71,7 @@ Create the pixel to world transform for the entire IFU: The WCS object now can evaluate simultaneously the transforms of all slices >>> wifu = wcs.WCS(forward_transform=regions_transform, output_frame=cframe, input_frame=det) - >>> x, y = mask.data.shape + >>> x, y = mask.mapper.shape >>> x, y = np.mgrid[:x, :y] >>> r, d, l = wifu(x, y) @@ -91,19 +91,15 @@ is functionally equivalent to the above commands bu returns coordinate objects: >>> wifu.available_frames ['detector', 'world'] - >>> wifu.world.coordinates(10, 200) + >>> wifu.output_frame.coordinates(10, 200) (, ) Frames provide additional information: - >>> print output_frame.axes_type + >>> print(wifu.output_frame.axes_type) [u'SPATIAL', u'SPECTRAL', u'SPATIAL'] - >>> print output_frame.axes_names + >>> print(wifu.output_frame.axes_names) [u'ra', 'lambda', u'dec'] - >>> output_frame.unit + >>> print(wifu.output_frame.unit) [Unit("deg"), Unit("micron"), Unit("deg")] - - - - diff --git a/docs/gwcs/using_wcs.rst b/docs/gwcs/using_wcs.rst index d8c19d2f..5c7f40d6 100644 --- a/docs/gwcs/using_wcs.rst +++ b/docs/gwcs/using_wcs.rst @@ -1,12 +1,57 @@ Using the WCS object ==================== -Let's use the WCS object created in `Getting Started`_ in order to show the interface. +Let's expand the WCS created in :ref:`getting-started` by adding a distortion. + +First create polynomial transform to represent distortion: + + >>> polyx = Polynomial2D(4) + >>> polyx.parameters = np.random.randn(15) + >>> polyy = Polynomial2D(4) + >>> polyy.parameters = np.random.randn(15) + >>> distortion = (Mapping((0, 1, 0, 1)) | polyx & polyy).rename("distortion") + +Create an intermediate frame. The distortion transforms positions on the +detector into this frame. + + >>> focal_frame = cf.Frame2D(name="focal_frame", unit=(u.arcsec, u.arcsec)) + +Create the WCS pipeline and initialize the WCS: + + >>> pipeline = [(detector_frame, distortion), + (focal_frame, det2sky), + (sky_frame, None) + ] + >>> wcsobj = wcs.WCS(pipeline) + >>> print(wcsobj) + From Transform + ----------- ---------- + detector distortion + focal_frame focal2sky + icrs None To see what frames are defined: >>> print(wcsobj.available_frames) - ['detector', 'focal', 'icrs'] + ['detector', 'focal_frame', 'icrs'] + >>> wcsobj.input_frame + + >>> wcsobj.output_frame + )> + +Because the ``output_frame`` is a `~gwcs.coordinate_frames.CoordinateFrame` object we can get +the result of the WCS transform as an `astropy.coordinates.SkyCoord` object and transform +them to other standard coordinate frames supported by `astropy.coordinates`. + + >>> skycoord = wcsobj(1, 2, output="numericals_plus") + >>> print(skycoord) + + >>> print(skycoord.transform_to("galactic")) + Some methods allow managing the transforms in a more detailed manner. @@ -14,19 +59,24 @@ Transforms between frames can be retrieved and evaluated separately. >>> distortion = wcsobj.get_transform('detector', 'focal') >>> distortion(1, 2) - (13.4, 0.) + (4.807433286098964, 4.924746607074259) Transforms in the pipeline can be replaced by new transforms. >>> new_transform = Shift(1) & Shift(1.5) | distortion - >>> wcsobj.set_transform('detector', 'focal', new_transform) + >>> wcsobj.set_transform('detector', 'focal_frame', new_transform) >>> wcsobj(1, 2) - (5.257230028926096, -72.53171157138964) + (7.641677379945592, -71.18890415491595) A transform can be inserted before or after a frame in the pipeline. >>> scale = Scale(2) & Scale(1) - >>> w.insert_transform('icrs', scale, after=False) - >>> w(1, 2) - (10.514460057852192, -72.53171157138964) + >>> wcsobj.insert_transform('icrs', scale, after=False) + >>> wcsobj(1, 2) + (15.283354759891184, -71.18890415491595) + +The WCS object has an attribute ``domain`` which describes the range of +acceptable values for each input axis. + >>> wcsobj.domain = [{'lower': 0, 'upper': 2048, 'includes_lower': True, 'includes_upper': False}, + {'lower': 0, 'upper': 1000, 'includes_lower': True, 'includes_upper': False}] diff --git a/docs/index.rst b/docs/index.rst index 185f019e..4f6ce5d2 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -1,5 +1,5 @@ GWCS Documentation -==================== +================== `GWCS `__ is a package for managing the World Coordinate System (WCS) of astronomical data. @@ -16,9 +16,8 @@ using the flexible framework of compound models in `~astropy.modeling`. In the case of a celestial output frame `~astropy.coordinates` provides further transformations between standard coordinate frames. Spectral output coordinates are instances of `~astropy.units.Quantity` and are -transformed to other units with the tools in that package. -The goal is to provide a flexible toolkit which is easily extendable by adding new -transforms and frames. +transformed to other units with the tools in that package. The goal is to provide +a flexible toolkit which is easily extendable by adding new transforms and frames. Installation @@ -26,7 +25,7 @@ Installation `gwcs `__ requires: -- `numpy `__ 1.67 or later +- `numpy `__ 1.7 or later - `astropy `__ 1.2 or later @@ -44,18 +43,20 @@ To install the latest release:: GWCS is also available as part of astroconda. +.. _getting-started: + Getting Started ---------------- +=============== The WCS data model represents a pipeline of transformations from some initial coordinate frame to a standard coordinate frame. It is implemented as a list of steps executed in order. Each step defines a starting coordinate frame and the transform to the next frame in the pipeline. The last step has no transform, only a frame which is the output frame of -the total transform. As a minimum a WCS object has an input_frame (defaults to "detector"), -an output_frame and the transform between them. +the total transform. As a minimum a WCS object has an ``input_frame`` (defaults to "detector"), +an ``output_frame`` and the transform between them. -As an example, consider a typical WCS of an image with some distortion. +As an example, consider a typical WCS of an image without distortion. The following imports are generally useful: @@ -67,76 +68,41 @@ The following imports are generally useful: >>> from gwcs import wcs >>> from gwcs import coordinate_frames as cf -The `forward_transform` is constructed as a combined model using `astropy.modeling`. +The ``forward_transform`` is constructed as a combined model using `astropy.modeling`. The frames are subclasses of `~gwcs.coordinate_frames.CoordinateFrame` (although strings are acceptable too). -First create polynomial transform to represent distortion: - - >>> polyx = Polynomial2D(4) - >>> polyx.parameters = np.random.randn(15) - >>> polyy = Polynomial2D(4) - >>> polyy.parameters = np.random.randn(15) - >>> distortion = (Mapping((0, 1, 0, 1)) | polyx & polyy).rename("distortion") - -Next create a transform from distortion fdree coordinates to ICRS. +Create a transform to convert detector coordinates to ICRS. - >>> dist2sky = (Shift(-10.5) & Shift(-13.2) | Rotation2D(0.0023) | \ - Scale(.01) & Scale(.04) | Pix2Sky_TAN() | RotateNative2Celestial(5.6, -72.05, 180)).rename("focal2sky") + >>> det2sky = (Shift(-10.5) & Shift(-13.2) | Rotation2D(0.0023) | \ + Scale(.01) & Scale(.04) | Pix2Sky_TAN() | \ + RotateNative2Celestial(5.6, -72.05, 180)).rename("detector2sky") -Create three coordinate frames, associated with the detector, a celestial frame, and an intermediate one. +Create a coordinate frame associated with the detector and a celestial frame. >>> detector_frame = cf.Frame2D(name="detector", axes_names=("x", "y"), unit=(u.pix, u.pix)) >>> sky_frame = cf.CelestialFrame(reference_frame=coord.ICRS(), name='icrs') - >>> focal_frame = cf.Frame2D(name="focal_frame", unit=(u.arcsec, u.arcsec)) -Create the WCS pipeline and initialize the WCS: +Initialize the WCS: - >>> pipeline = [(detector_frame, distortion), - (focal_frame, dist2sky), - (sky_frame, None) - ] - >>> wcsobj = wcs.WCS(pipeline) + >>> wcsobj = wcs.WCS(forward_transform=det2sky, input_frame=detector_frame, output_frame=sky_frame) >>> print(wcsobj) From Transform ----------- ---------- - detector distortion - focal_frame focal2sky + detector detector2sky icrs None To convert a pixel (x, y) = (1, 2) to sky coordinates, call the WCS object as a function: >>> sky = wcsobj(1, 2) >>> print(sky) - (5.759024831907874, -72.22626601919619) - -#The :meth:`~gwcs.wcs.WCS.invert` method evaluates the :meth:`~gwcs.wcs.WCS.backward_transform` -#if available, otherwise applies an iterative method to calculate the reverse coordinates. -# -# >>> w.invert(*sky) -# (1.000000000009388, 2.0000000000112728) -# -#Some methods allow managing the transforms in a more detailed manner. - -Transforms between frames can be retrieved and evaluated separately. + (5.284139265842838, -72.49775640633504) - >>> dist = w.get_transform('detector_frame', 'focal_frame') - >>> dist(1, 2) - (15.354214305518118, 8.791536957201615) +The :meth:`~gwcs.wcs.WCS.invert` method evaluates the :meth:`~gwcs.wcs.WCS.backward_transform` +if available, otherwise applies an iterative method to calculate the reverse coordinates. -Transforms in the pipeline can be replaced by new transforms. - - >>> new_transform = distortion | Shift(1) & Shift(1.5) - >>> wcsobj.set_transform('detector', 'focal', new_transform) - >>> wcsobj(1, 2) - (5.791157736884894, -72.16623599444335) - -A transform can be inserted before or after a frame in the pipeline. - - >>> scale = Scale(2) & Scale(1) - >>> wcsobj.insert_transform('icrs', scale, after=False) - >>> wcsobj(1, 2) - (11.582315473769787, -72.16623599444335) + >>> wcsobj.invert(*sky) + (1.000, 2.000) Using `gwcs` @@ -146,11 +112,12 @@ Using `gwcs` .. toctree:: :maxdepth: 2 - gwcs/create_wcs + gwcs/using_wcs.rst gwcs/selector_model.rst gwcs/wcstools.rst + See also -------- diff --git a/gwcs/wcs.py b/gwcs/wcs.py index 3279c315..4d560987 100644 --- a/gwcs/wcs.py +++ b/gwcs/wcs.py @@ -391,6 +391,9 @@ def pipeline(self): @property def domain(self): + """ + Return the range of acceptable values for each input axis. + """ frames = self.available_frames transform_meta = self.get_transform(frames[0], frames[1]).meta if 'domain' in transform_meta: @@ -400,6 +403,9 @@ def domain(self): @domain.setter def domain(self, value): + """ + Set the range of acceptable values for each input axis. + """ self._validate_domain(value) frames = self.available_frames transform = self.get_transform(frames[0], frames[1]) From e6f93b6df65053655cdb35b02b6ddee954655ff3 Mon Sep 17 00:00:00 2001 From: Nadia Dencheva Date: Wed, 21 Dec 2016 16:31:07 -0500 Subject: [PATCH 5/5] pychecker corrections --- gwcs/coordinate_frames.py | 1 - gwcs/region.py | 18 ++++++++---------- gwcs/selector.py | 18 ++++++------------ gwcs/utils.py | 2 +- gwcs/wcs.py | 1 - 5 files changed, 15 insertions(+), 25 deletions(-) diff --git a/gwcs/coordinate_frames.py b/gwcs/coordinate_frames.py index 66e9f396..f156f629 100644 --- a/gwcs/coordinate_frames.py +++ b/gwcs/coordinate_frames.py @@ -9,7 +9,6 @@ from astropy import utils as astutil from astropy import coordinates as coord from astropy.extern import six -from . import utils as gwutils __all__ = ['Frame2D', 'CelestialFrame', 'SpectralFrame', 'CompositeFrame', diff --git a/gwcs/region.py b/gwcs/region.py index afd789ec..1ccad4df 100644 --- a/gwcs/region.py +++ b/gwcs/region.py @@ -58,6 +58,7 @@ def scan(self, mask): Pixels which are not included in any region are marked with 0 or "". """ + class Polygon(Region): """ @@ -84,7 +85,7 @@ def __init__(self, rid, vertices, coord_frame="detector"): self._vertices = np.asarray(vertices) self._bbox = self._get_bounding_box() - self._scan_line_range = list(range(self._bbox[1], self._bbox[3]+self._bbox[1]+1)) + self._scan_line_range = list(range(self._bbox[1], self._bbox[3] + self._bbox[1] + 1)) # constructs a Global Edge Table (GET) in bbox coordinates self._GET = self._construct_ordered_GET() @@ -164,12 +165,12 @@ def scan(self, data): while y <= scline: AET = self.update_AET(y, AET) scan_line = Edge('scan_line', start=[self._bbox[0], y], - stop=[self._bbox[0]+self._bbox[2], y]) + stop=[self._bbox[0] + self._bbox[2], y]) x = [np.ceil(e.compute_AET_entry(scan_line)[1]) for e in AET if e is not None] xnew = np.asarray(np.sort(x), dtype=np.int) for i, j in zip(xnew[::2], xnew[1::2]): data[y][i:j+1] = self._rid - y = y+1 + y = y + 1 return data def update_AET(self, y, AET): @@ -195,12 +196,9 @@ def update_AET(self, y, AET): def __contains__(self, px): """even-odd algorithm or smth else better sould be used""" - minx = self._vertices[:, 0].min() - maxx = self._vertices[:, 0].max() - miny = self._vertices[:, 1].min() - maxy = self._vertices[:, 1].max() - return px[0] >= self._bbox[0] and px[0] <= self._bbox[0]+self._bbox[2] and \ - px[1] >= self._bbox[1] and px[1] <= self._bbox[1]+self._bbox[3] + return px[0] >= self._bbox[0] and px[0] <= self._bbox[0] + self._bbox[2] and \ + px[1] >= self._bbox[1] and px[1] <= self._bbox[1] + self._bbox[3] + class Edge(object): @@ -316,7 +314,7 @@ def intersection(self, edge): v = edge._stop - edge._start w = self._start - edge._start D = np.cross(u, v) - return np.cross(v, w)/D * u + self._start + return np.cross(v, w) / D * u + self._start def is_parallel(self, edge): u = self._stop - self._start diff --git a/gwcs/selector.py b/gwcs/selector.py index f19ed290..655a3da6 100644 --- a/gwcs/selector.py +++ b/gwcs/selector.py @@ -67,14 +67,8 @@ """ from __future__ import absolute_import, division, unicode_literals, print_function -import copy -import numbers import numpy as np -import warnings - -from astropy.extern import six from astropy.modeling.core import Model -from astropy.modeling.parameters import Parameter from . import region from .utils import RegionError, _toindex @@ -372,15 +366,15 @@ def _has_overlapping(ranges): Test a list of tuple representing ranges of values has no overlapping ranges. """ d = dict(ranges) - start = ranges[:,0] - end = ranges[:,1] + start = ranges[:, 0] + end = ranges[:, 1] start.sort() l = [] for v in start: l.append([v, d[v]]) l = np.array(l) - start = np.roll(l[: ,0], -1) - end = l[: ,1] + start = np.roll(l[:, 0], -1) + end = l[:, 1] if any((end - start)[:-1] > 0) or any(start[-1] > end): return True else: @@ -404,7 +398,7 @@ def _find_range(self, value_range, value): Index of the tuple which defines a range holding the input value. None, if the input value is not within any available range. """ - a, b = value_range[:,0], value_range[:,1] + a, b = value_range[:, 0], value_range[:, 1] ind = np.logical_and(value >= a, value <= b).nonzero()[0] if ind.size > 1: raise ValueError("There are overlapping ranges.") @@ -484,7 +478,7 @@ def __init__(self, inputs, outputs, selector, label_mapper, undefined_transform_ self._outputs = outputs self.label_mapper = label_mapper self._undefined_transform_value = undefined_transform_value - self._selector = selector #copy.deepcopy(selector) + self._selector = selector # copy.deepcopy(selector) if " " in selector.keys() or 0 in selector.keys(): raise ValueError('"0" and " " are not allowed as keys.') diff --git a/gwcs/utils.py b/gwcs/utils.py index b5a93454..316801d1 100644 --- a/gwcs/utils.py +++ b/gwcs/utils.py @@ -739,7 +739,7 @@ def separable_axes(wcsobj, start_frame=None, end_frame=None): raise ValueError("A starting frame is needed to determine separability of axes.") sep = is_separable(transform) - return [sep[ax] for ax in self.axes_order] + return [sep[ax] for ax in end_frame.axes_order] _operators = {'&': _cstack, '|': _cdot, '+': _arith_oper, '-': _arith_oper, diff --git a/gwcs/wcs.py b/gwcs/wcs.py index 4d560987..bda5e65c 100644 --- a/gwcs/wcs.py +++ b/gwcs/wcs.py @@ -474,4 +474,3 @@ def footprint(self, domain=None, center=True): vertices += .5 result = self.__call__(*vertices) return np.asarray(result) -