Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

SCSB-111: Improve documentation #483

Merged
merged 10 commits into from
Dec 18, 2023
7 changes: 4 additions & 3 deletions CHANGES.rst
Original file line number Diff line number Diff line change
@@ -1,10 +1,11 @@
0.20.0 (unreleased)
-------------------

other
^^^^^
- Improve documentation (part 1) [#483]
nden marked this conversation as resolved.
Show resolved Hide resolved

- Replace ``pkg_resources`` with ``importlib.metadata``. [#478]
- other
^^^^^
- Replace ``pkg_resources`` with ``importlib.metadata``. [#478]
nden marked this conversation as resolved.
Show resolved Hide resolved

0.19.0 (2023-09-15)
-------------------
Expand Down
101 changes: 101 additions & 0 deletions docs/gwcs/fits_analog.rst
Original file line number Diff line number Diff line change
@@ -0,0 +1,101 @@
.. _fits_equivalent_example
nden marked this conversation as resolved.
Show resolved Hide resolved

FITS Equivalent WCS Example
===========================

The following example shows how to construct a GWCS object equivalent to
a FITS imaging WCS without distortion, defined in this FITS imaging header::

WCSAXES = 2 / Number of coordinate axes
WCSNAME = '47 Tuc ' / Coordinate system title
CRPIX1 = 2048.0 / Pixel coordinate of reference point
CRPIX2 = 1024.0 / Pixel coordinate of reference point
PC1_1 = 1.290551569736E-05 / Coordinate transformation matrix element
PC1_2 = 5.9525007864732E-06 / Coordinate transformation matrix element
PC2_1 = 5.0226382102765E-06 / Coordinate transformation matrix element
PC2_2 = -1.2644844123757E-05 / Coordinate transformation matrix element
CDELT1 = 1.0 / [deg] Coordinate increment at reference point
CDELT2 = 1.0 / [deg] Coordinate increment at reference point
CUNIT1 = 'deg' / Units of coordinate increment and value
CUNIT2 = 'deg' / Units of coordinate increment and value
CTYPE1 = 'RA---TAN' / TAN (gnomonic) projection + SIP distortions
CTYPE2 = 'DEC--TAN' / TAN (gnomonic) projection + SIP distortions
CRVAL1 = 5.63056810618 / [deg] Coordinate value at reference point
CRVAL2 = -72.05457184279 / [deg] Coordinate value at reference point
LONPOLE = 180.0 / [deg] Native longitude of celestial pole
LATPOLE = -72.05457184279 / [deg] Native latitude of celestial pole
RADESYS = 'ICRS' / Equatorial coordinate system


The following imports are generally useful:
nden marked this conversation as resolved.
Show resolved Hide resolved

>>> import numpy as np
>>> from astropy.modeling import models
>>> from astropy import coordinates as coord
>>> from astropy import units as u
>>> 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`` are subclasses of `~gwcs.coordinate_frames.CoordinateFrame`. Although strings are
acceptable as ``coordinate_frames`` it is recommended this is used only in testing/debugging.
nden marked this conversation as resolved.
Show resolved Hide resolved

Using the `~astropy.modeling` package create a combined model to transform
detector coordinates to ICRS following the FITS WCS standard convention.

nden marked this conversation as resolved.
Show resolved Hide resolved
First, create a transform which shifts the input ``x`` and ``y`` coordinates by ``CRPIX``. We subtract 1 from the CRPIX values because the first pixel is considered pixel ``1`` in FITS WCS:

>>> shift_by_crpix = models.Shift(-(2048 - 1)*u.pix) & models.Shift(-(1024 - 1)*u.pix)

Create a transform which rotates the inputs using the ``PC matrix``.

>>> matrix = np.array([[1.290551569736E-05, 5.9525007864732E-06],
... [5.0226382102765E-06 , -1.2644844123757E-05]])
>>> rotation = models.AffineTransformation2D(matrix * u.deg,
... translation=[0, 0] * u.deg)
>>> rotation.input_units_equivalencies = {"x": u.pixel_scale(1*u.deg/u.pix),
... "y": u.pixel_scale(1*u.deg/u.pix)}
>>> rotation.inverse = models.AffineTransformation2D(np.linalg.inv(matrix) * u.pix,
... translation=[0, 0] * u.pix)
>>> rotation.inverse.input_units_equivalencies = {"x": u.pixel_scale(1*u.pix/u.deg),
... "y": u.pixel_scale(1*u.pix/u.deg)}

Create a tangent projection and a rotation on the sky using ``CRVAL``.

>>> tan = models.Pix2Sky_TAN()
>>> celestial_rotation = models.RotateNative2Celestial(5.63056810618*u.deg, -72.05457184279*u.deg, 180*u.deg)

>>> det2sky = shift_by_crpix | rotation | tan | celestial_rotation
>>> det2sky.name = "linear_transform"

Create a ``detector`` coordinate frame and a ``celestial`` ICRS frame.
Copy link
Member

Choose a reason for hiding this comment

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

Why do we need this detector frame? Can't we just bunch all transforms together from detector all the way to the sky? Also, what is a "detector frame"? is it just the value of the name argument in the Frame2D initializer?

Copy link
Collaborator

Choose a reason for hiding this comment

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

Generally the frames are information containers - they hold the units on the axes, the names of the axes, the axes type and physical type. You could create a WCS without an input frame but some functionality will not be available. For example, if the transforms are defined with units, you will need to make sure the inputs are with units, so any time you pass input in the case above, the inputs need to have u.pix attached to them. This can be annoying.
In addition the input frame does not always represent a detector and imo it's good to have the information on what inputs are expected.

Copy link
Member

Choose a reason for hiding this comment

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

@nden comment should be added to the docs. My comment was intended to prompt a discussion/explanation about organization of the (gwcs) pipelines. For example, why do some pipelines have multiple steps (JWST WCS comes to mind)?

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

I wonder if a paragraphs or two before this example outlining the general structure of a WCS object would be helpful. Any objection @nden?

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

What if the current "Getting Started" is replaced by an earlier section:

Basic Structure of a GWCS Object
--------------------------------

The key concept to be aware of is that a GWCS Object consists of a pipeline
of steps; each step contains a transform (i.e., an Astropy model) that 
converts the input coordinates of the step to the output coordinates of
the step. Furthermore, each step has an optional coordinate frame associated
with the step. The coordinate frame represents the input coordinate frame, not
the output coordinates. Most typically, the first step coordinate frame is
the detector pixel coordinates (the default). Since no step has a coordinate 
frame for the output coordinates, it is necessary to append a step with no
transform to the end of the pipeline to represent the output coordinate frame.
For imaging, this frame typically references one of the Astropy standard
Sky Coordinate Frames of Reference. The GWCS frames also serve to hold the
units on the axes, the names of the axes and the physical type of the axis
(e.g., wavelength).

Since it is often useful to obtain coordinates in an intermediate frame of
reference, GWCS allows the pipeline to consist of more than one transform.
For example, for spectrographs, it is useful to have access to coordinates
in the slit plane, and in such a case, the first step would transform from
the detector to the slit plane, and the second step from the slit plane to
sky coordinates and a wavelength. Constructed this way, it is possible to
extract from the GWCS the needed transforms between identified frames of
reference.  

The GWCS object can be saved to the ASDF format using the
`asdf <https://asdf.readthedocs.io/en/latest/>`__ package and validated 
using `ASDF Standard <https://asdf-standard.readthedocs.io/en/latest/>`__

There are two ways to save the GWCS object to a files:

- `Save a WCS object as a pure ASDF file`_ 

@mcara does this address your comment?

Copy link
Member

Choose a reason for hiding this comment

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

Yes! This is a great introduction.

Copy link
Collaborator

Choose a reason for hiding this comment

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

I like this too.

Furthermore, each step has an optional coordinate frame associated with the step.

Is the frame optional? I thought it's required.

a GWCS Object consists of a pipeline of steps

Suggest adding linear -> consists of a linear pipeline of steps
so that it's clear it's not an arbitrary graph.


>>> 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',
... unit=(u.deg, u.deg))

This WCS pipeline has only one step - from ``detector`` to ``sky``:

>>> pipeline = [(detector_frame, det2sky),
... (sky_frame, None)
... ]
>>> wcsobj = wcs.WCS(pipeline)
>>> print(wcsobj)
From Transform
-------- ----------------
detector linear_transform
icrs None

To convert a pixel (x, y) = (1, 2) to sky coordinates, call the WCS object as a function:
Copy link
Member

Choose a reason for hiding this comment

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

Now we have a complete WCS object. In the next example, we will use it to convert pixel coordinates (1, 2) to sky...

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

I've modified it to something similar.


>>> sky = wcsobj(1*u.pix, 2*u.pix, with_units=True)
>>> print(sky)
<SkyCoord (ICRS): (ra, dec) in deg
(5.52515954, -72.05190935)>

The :meth:`~gwcs.wcs.WCS.invert` method evaluates the :meth:`~gwcs.wcs.WCS.backward_transform`
nden marked this conversation as resolved.
Show resolved Hide resolved
if available, otherwise applies an iterative method to calculate the reverse coordinates.

>>> wcsobj.invert(sky)
(<Quantity 1. pix>, <Quantity 2. pix>)
Binary file modified docs/gwcs/ifu-regions.png
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
24 changes: 12 additions & 12 deletions docs/gwcs/ifu.rst
Original file line number Diff line number Diff line change
Expand Up @@ -36,6 +36,16 @@ First, import the usual packages.
>>> from gwcs import wcs, selector
>>> from gwcs import coordinate_frames as cf

Next, create the appropriate mapper object corresponding to the figure above:

>>> y, x = np.mgrid[:1000, :500]
>>> fmask = (((-x + 0.01 * y + 0.00002 * y**2)/ 500) * 13 - 0.5) + 14
>>> mask = fmask.astype(np.int8)
>>> mask[(mask % 2) == 1] = 0
>>> mask[mask > 13] = 0
>>> mask = mask // 2
>>> labelmapper = selector.LabelMapperArray(mask)

nden marked this conversation as resolved.
Show resolved Hide resolved
The output frame is common for all slits and is a composite frame with two subframes,
`~gwcs.coordinate_frames.CelestialFrame` and `~gwcs.coordinate_frames.SpectralFrame`.

Expand All @@ -59,28 +69,18 @@ In this example the mask is an array with the size of the detector where each it
corresponds to a pixel on the detector and its value is the slice number (label) this pixel
belongs to.

Assuming the array is stored in
`ASDF <https://asdf-standard.readthedocs.io/en/latest>`__ format, create the mask:

.. doctest-skip-all

>>> import asdf
>>> f = asdf.open('mask.asdf')
>>> data = f.tree['mask']
>>> mask = selector.LabelMapperArray(data)

Create the pixel to world transform for the entire IFU:

>>> regions_transform = selector.RegionsSelector(inputs=['x','y'],
... outputs=['ra', 'dec', 'lam'],
... selector=transforms,
... label_mapper=mask,
... label_mapper=labelmapper,
... undefined_transform_value=np.nan)

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)
>>> y, x = mask.mapper.shape
>>> y, x = labelmapper.mapper.shape
>>> y, x = np.mgrid[:y, :x]
>>> r, d, l = wifu(x, y)

Expand Down
4 changes: 4 additions & 0 deletions docs/gwcs/using_wcs.rst
Original file line number Diff line number Diff line change
Expand Up @@ -109,6 +109,9 @@ Inverse Transformations
Often, it is useful to be able to compute inverse transformation that converts
coordinates from the output frame back to the coordinates in the input frame.

Note. the ``backward_transform`` attribute is equivalent to
``forward_transform.inverse``.

In this section, for illustration purpose, we will be using the same 2D imaging
WCS from ``imaging_wcs_wdist.asdf`` created in :ref:`imaging_example` whose
forward transformation converts image coordinates to world coordinates and
Expand Down Expand Up @@ -175,3 +178,4 @@ point far away from the image for which numerical inverse fails.
[5.31656943e-06 2.72052603e-10]
[6.81557583e-06 1.06560533e-06]
[3.96365344e-04 6.41822468e-05]]

146 changes: 60 additions & 86 deletions docs/index.rst
Original file line number Diff line number Diff line change
Expand Up @@ -124,87 +124,87 @@ There are two ways to save the WCS to a file:
- `Save a WCS object as an ASDF extension in a FITS file`_
nden marked this conversation as resolved.
Show resolved Hide resolved


A step by step example of constructing an imaging GWCS object.
A step-by-step example of constructing an imaging GWCS object.
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

The following example shows how to construct a GWCS object equivalent to
a FITS imaging WCS without distortion, defined in this FITS imaging header::

WCSAXES = 2 / Number of coordinate axes
WCSNAME = '47 Tuc ' / Coordinate system title
CRPIX1 = 2048.0 / Pixel coordinate of reference point
CRPIX2 = 1024.0 / Pixel coordinate of reference point
PC1_1 = 1.290551569736E-05 / Coordinate transformation matrix element
PC1_2 = 5.9525007864732E-06 / Coordinate transformation matrix element
PC2_1 = 5.0226382102765E-06 / Coordinate transformation matrix element
PC2_2 = -1.2644844123757E-05 / Coordinate transformation matrix element
CDELT1 = 1.0 / [deg] Coordinate increment at reference point
CDELT2 = 1.0 / [deg] Coordinate increment at reference point
CUNIT1 = 'deg' / Units of coordinate increment and value
CUNIT2 = 'deg' / Units of coordinate increment and value
CTYPE1 = 'RA---TAN' / TAN (gnomonic) projection + SIP distortions
CTYPE2 = 'DEC--TAN' / TAN (gnomonic) projection + SIP distortions
CRVAL1 = 5.63056810618 / [deg] Coordinate value at reference point
CRVAL2 = -72.05457184279 / [deg] Coordinate value at reference point
LONPOLE = 180.0 / [deg] Native longitude of celestial pole
LATPOLE = -72.05457184279 / [deg] Native latitude of celestial pole
RADESYS = 'ICRS' / Equatorial coordinate system
The following example shows how to construct a GWCS object that maps
input pixel coordinates to sky coordinates. This example
involves 4 sequential transformations:

- Adjusting pixel coordinates such that the center of the array has
(0, 0) value (typical of most WCS definitions, but any pixel may
be the reference that is tied to the sky reference, even the (0, 0)
pixel, or even pixels outside of the detector).
- Scaling pixels such that the center pixel of the array has the expected
angular scale. (I.e., applying the plate scale)
- Projecting the resultant coordinates onto the sky using the tangent
projection. If the field of view is small, the inaccuracies resulting
leaving this out will be small; however, this is generally applied.
- Transforming the center pixel to the appropriate celestial coordinate
with the approprate orientation on the sky. For simplicity's sake,
we assume the detector array is already oriented with north up, and
that the array has the appropriate parity as the sky coordinates.


The detector has a 1000 pixel by 1000 pixel array.

For simplicity, no units will be used, but instead will be implicit.
Copy link
Member

Choose a reason for hiding this comment

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

shouldn't this go into the example in fits_analog.rst?

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

Again, do you mean duplicating (with needed changes) it there rather than moving?


The following imports are generally useful:

.. doctest-skip::

>>> import numpy as np
>>> from astropy.modeling import models
>>> from astropy import coordinates as coord
>>> from astropy import units as u
>>> 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`` are subclasses of `~gwcs.coordinate_frames.CoordinateFrame`. Although strings are
acceptable as ``coordinate_frames`` it is recommended this is used only in testing/debugging.
In the following transformation definitions, angular units are in degrees by
default.

.. doctest-skip::

Using the `~astropy.modeling` package create a combined model to transform
detector coordinates to ICRS following the FITS WCS standard convention.
>>> pixelshift = models.Shift(-500) & models.Shift(-500)
>>> pixelscale = models.Scale(0.1 / 3600.) & models.Scale(0.1 / 3600.) # 0.1 arcsec/pixel
>>> tangent_projection = models.Pix2Sky_TAN()
>>> celestial_rotation = models.RotateNative2Celestial(30., 45., 180.)

First, create a transform which shifts the input ``x`` and ``y`` coordinates by ``CRPIX``. We subtract 1 from the CRPIX values because the first pixel is considered pixel ``1`` in FITS WCS:
For the last transformation, the 3 arguments are, respectively:
nden marked this conversation as resolved.
Show resolved Hide resolved

>>> shift_by_crpix = models.Shift(-(2048 - 1)*u.pix) & models.Shift(-(1024 - 1)*u.pix)
- Celestial longitude (i.e., RA) of the fiducial point (e.g., (0, 0) in the input
spherical coordinates).
In this case we put the detector center at 30 degrees (RA = 2 hours)
- Celestial latitude (i.e., Dec) of the fiducial point. Here Dec = 45 degrees.
- Longitude of celestial pole in input coordinate system. With north up, this
always corresponds to a value of 180.

Create a transform which rotates the inputs using the ``PC matrix``.
The more general case where the detector is not aligned with north, would have
a rotation transform after the pixelship and pixelscale transformations to
Copy link
Member

Choose a reason for hiding this comment

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

Suggested change
a rotation transform after the pixelship and pixelscale transformations to
a rotation transform after the pixelshift and pixelscale transformations to

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

So modified, though I really like how pixelship sounds.

align the detector coordinates with north up.

>>> matrix = np.array([[1.290551569736E-05, 5.9525007864732E-06],
... [5.0226382102765E-06 , -1.2644844123757E-05]])
>>> rotation = models.AffineTransformation2D(matrix * u.deg,
... translation=[0, 0] * u.deg)
>>> rotation.input_units_equivalencies = {"x": u.pixel_scale(1*u.deg/u.pix),
... "y": u.pixel_scale(1*u.deg/u.pix)}
>>> rotation.inverse = models.AffineTransformation2D(np.linalg.inv(matrix) * u.pix,
... translation=[0, 0] * u.pix)
>>> rotation.inverse.input_units_equivalencies = {"x": u.pixel_scale(1*u.pix/u.deg),
... "y": u.pixel_scale(1*u.pix/u.deg)}
The net transformation from pixel coordinates to celestial coordinates then
becomes:

Create a tangent projection and a rotation on the sky using ``CRVAL``.
.. doctest-skip::

>>> tan = models.Pix2Sky_TAN()
>>> celestial_rotation = models.RotateNative2Celestial(5.63056810618*u.deg, -72.05457184279*u.deg, 180*u.deg)
>>> det2sky = pixelshift | pixelscale | tangent_projection | celestial_rotation

>>> det2sky = shift_by_crpix | rotation | tan | celestial_rotation
>>> det2sky.name = "linear_transform"
The remaining elements to defining the WCS are he input and output
frames of reference. While the GWCS scheme allows intermediate frames
of reference, this example doesn't have any. The output frame is
expressed with no associated transform

Create a ``detector`` coordinate frame and a ``celestial`` ICRS frame.
.. doctest-skip::

>>> 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',
... unit=(u.deg, u.deg))

This WCS pipeline has only one step - from ``detector`` to ``sky``:

>>> pipeline = [(detector_frame, det2sky),
... (sky_frame, None)
... ]
>>> wcsobj = wcs.WCS(pipeline)
>>> wcsobj = wcs.WCS([(detector_frame, det2sky),
... (sky_frame, None)
... ]
>>> print(wcsobj)
From Transform
-------- ----------------
Expand All @@ -213,14 +213,18 @@ This WCS pipeline has only one step - from ``detector`` to ``sky``:

To convert a pixel (x, y) = (1, 2) to sky coordinates, call the WCS object as a function:

>>> sky = wcsobj(1*u.pix, 2*u.pix, with_units=True)
.. doctest-skip::

>>> sky = wcsobj(1, 2)
>>> print(sky)
<SkyCoord (ICRS): (ra, dec) in deg
(5.52515954, -72.05190935)>

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.

.. doctest-skip::

>>> wcsobj.invert(sky)
(<Quantity 1. pix>, <Quantity 2. pix>)

Expand All @@ -240,31 +244,6 @@ Save a WCS object as a pure ASDF file
:ref:`pure_asdf`


Save a WCS object as an ASDF extension in a FITS file
+++++++++++++++++++++++++++++++++++++++++++++++++++++


.. doctest-skip::

>>> from astropy.io import fits
>>> from asdf import fits_embed
>>> hdul = fits.open("example_imaging.fits")
>>> hdul.info()
Filename: example_imaging.fits
No. Name Ver Type Cards Dimensions Format
0 PRIMARY 1 PrimaryHDU 775 ()
1 SCI 1 ImageHDU 71 (600, 550) float32
>>> tree = {"sci": hdul[1].data,
... "wcs": wcsobj}
>>> fa = fits_embed.AsdfInFits(hdul, tree)
>>> fa.write_to("imaging_with_wcs_in_asdf.fits")
>>> fits.info("imaging_with_wcs_in_asdf.fits")
Filename: example_with_wcs.asdf
No. Name Ver Type Cards Dimensions Format
0 PRIMARY 1 PrimaryHDU 775 ()
1 SCI 1 ImageHDU 71 (600, 550) float32
2 ASDF 1 BinTableHDU 11 1R x 1C [5086B]

Reading a WCS object from a file
++++++++++++++++++++++++++++++++

Expand All @@ -280,12 +259,6 @@ from a pure ASDF file or from an ASDF extension in a FITS file.
>>> wcsobj = asdf_file.tree['wcs']


.. doctest-skip::

>>> import asdf
>>> fa = asdf.open("imaging_with_wcs_in_asdf.fits")
>>> wcsobj = fa.tree["wcs"]

Other Examples
--------------

Expand All @@ -309,6 +282,7 @@ Using ``gwcs``
gwcs/pure_asdf.rst
gwcs/wcs_validation.rst
gwcs/points_to_wcs.rst
gwcs/fits_analog.rst


See also
Expand Down