Merge pull request #340 from spacetelescope/0.16.0x

0.16.0x

.travis.yml

@@ -15,7 +15,7 @@ env:
         # The following versions are the 'default' for tests, unless
         # overidden underneath. They are defined here in order to save having
         # to repeat them for all configurations.
-        - NUMPY_VERSION="1.16"
+        - NUMPY_VERSION="1.17"
         - ASDF_GIT="git+https://github.com/spacetelescope/asdf.git#egg=asdf"
         - ASTROPY_GIT="git+https://github.com/astropy/astropy.git#egg=astropy"
         - PIP_DEPENDENCIES=.[test]
@@ -41,7 +41,7 @@ jobs:
         # Check for sphinx doc build warnings - we do this first because it
         # may run for a long time
         - python: "3.8"
-          env: 
+          env:
             PIP_DEPENDENCIES=".[docs]"
             TEST_COMMAND="make --directory=docs html"
           addons:

CHANGES.rst

@@ -1,3 +1,27 @@
+0.16.0 (2020-12-18)
+-------------------
+New Features
+^^^^^^^^^^^^
+
+- Added an option to `to_fits_sip()` to be able to specify the reference
+  point (``crpix``) of the FITS WCS. [#337]
+
+- Added support for providing custom range of degrees in ``to_fits_sip``. [#339]
+
+Bug Fixes
+^^^^^^^^^
+
+- ``bounding_box`` now works with tuple of ``Quantities``. [#331]
+
+- Fix a formula for estimating ``crpix`` in ``to_fits_sip()`` so that ``crpix``
+  is near the center of the bounding box. [#337]
+
+- Allow sub-pixel sampling of the WCS model when computing SIP approximation in
+  ``to_fits_sip()``. [#338]
+
+- Fixed a bug in ``to_fits_sip`` due to which ``inv_degree`` was ignored. [#339]
+
+
 0.15.0 (2020-11-13)
 -------------------
 New Features

docs/gwcs/imaging_with_distortion.rst

@@ -19,14 +19,15 @@ The imaging example without units:
   >>> from gwcs import wcs
   >>> from gwcs import coordinate_frames as cf
 
-  >>> shift_by_crpix = models.Shift(-2048) & models.Shift(-1024)
+  >>> crpix = (2048, 1024)
+  >>> shift_by_crpix = models.Shift(-crpix[0]) & models.Shift(-crpix[1])
   >>> matrix = np.array([[1.290551569736E-05, 5.9525007864732E-06],
   ...                    [5.0226382102765E-06 , -1.2644844123757E-05]])
   >>> rotation = models.AffineTransformation2D(matrix)
   >>> rotation.inverse = models.AffineTransformation2D(np.linalg.inv(matrix))
   >>> tan = models.Pix2Sky_TAN()
   >>> celestial_rotation =  models.RotateNative2Celestial(5.63056810618, -72.05457184279, 180)
-  >>> det2sky = shift_by_crpix | rotation | tan | celestial_rotation 
+  >>> det2sky = shift_by_crpix | rotation | tan | celestial_rotation
   >>> det2sky.name = "linear_transform"
   >>> detector_frame = cf.Frame2D(name="detector", axes_names=("x", "y"),
   ...                             unit=(u.pix, u.pix))
@@ -46,21 +47,24 @@ First create distortion corrections represented by a polynomial
 model of fourth degree. The example uses the astropy `~astropy.modeling.polynomial.Polynomial2D`
 and `~astropy.modeling.mappings.Mapping` models.
 
-
   >>> poly_x = models.Polynomial2D(4)
-  >>> poly_x.parameters = np.arange(15) * .1
+  >>> poly_x.parameters = [0, 1, 8.55e-06, -4.73e-10, 2.37e-14, 0, -5.20e-06,
+  ...                      -3.98e-11, 1.97e-15, 2.17e-06, -5.23e-10, 3.47e-14,
+  ...                      1.08e-11, -2.46e-14, 1.49e-14]
   >>> poly_y = models.Polynomial2D(4)
-  >>> poly_y.parameters = np.arange(15) * .2
-  >>> distortion = models.Mapping((0, 1, 0, 1)) | poly_x & poly_y
+  >>> poly_y.parameters = [0, 0, -1.75e-06, 8.57e-11, -1.77e-14, 1, 6.18e-06,
+  ...                      -5.09e-10, -3.78e-15, -7.22e-06, -6.17e-11,
+  ...                      -3.66e-14, -4.18e-10, 1.22e-14, -9.96e-15]
+  >>> distortion = ((models.Shift(-crpix[0]) & models.Shift(-crpix[1])) |
+  ...               models.Mapping((0, 1, 0, 1)) | (poly_x & poly_y) |
+  ...               (models.Shift(crpix[0]) & models.Shift(crpix[1])))
   >>> distortion.name = "distortion"
-
-Create an intermediate frame for distortion free coordinates.
 
+Create an intermediate frame for distortion free coordinates.
 
   >>> undistorted_frame = cf.Frame2D(name="undistorted_frame", unit=(u.pix, u.pix),
   ...                                axes_names=("undist_x", "undist_y"))
 
-
 Using the example in :ref:`getting-started`, add the distortion correction to
 the WCS pipeline and initialize the WCS.
 
@@ -70,9 +74,17 @@ the WCS pipeline and initialize the WCS.
   ...             ]
   >>> wcsobj = wcs.WCS(pipeline)
   >>> print(wcsobj)
-         From          Transform    
+         From          Transform
   ----------------- ----------------
            detector       distortion
   undistorted_frame linear_transform
                icrs             None
-
+
+Finally, save this WCS to an ``ASDF`` file:
+
+.. doctest-skip::
+
+  >>> from asdf import AsdfFile
+  >>> tree = {"wcs": wcsobj}
+  >>> wcs_file = AsdfFile(tree)
+  >>> wcs_file.write_to("imaging_wcs_wdist.asdf")

docs/gwcs/using_wcs.rst

@@ -1,95 +1,177 @@
-.. _user_api:
+.. _using_wcs_examples:
 
 Using the WCS object
 ====================
-   
-This section uses the ``imaging_wcs.asdf`` created in :ref:`imaging_example`
+
+This section uses the ``imaging_wcs_wdist.asdf`` created in :ref:`imaging_example`
 to read in a WCS object and demo its methods.
 
 .. doctest-skip::
-   
+
   >>> import asdf
-  >>> asdf_file = asdf.open("imaging_wcs.asdf")
+  >>> asdf_file = asdf.open("imaging_wcs_wdist.asdf")
   >>> wcsobj = asdf_file.tree["wcs"]
   >>> print(wcsobj)    # doctest: +SKIP
-           From          Transform    
+           From          Transform
   ----------------- ----------------
            detector       distortion
   undistorted_frame linear_transform
                icrs             None
 
+
+Inspecting Available Coordinate Frames
+--------------------------------------
+
 To see what frames are defined:
 
 .. doctest-skip::
 
-   >>> print(wcsobj.available_frames)
-       ['detector', 'undistorted_frame', 'icrs']
-   >>> wcsobj.input_frame
-        <Frame2D(name="detector", unit=(Unit("pix"), Unit("pix")), axes_names=('x', 'y'), axes_order=(0, 1))>
-    >>> wcsobj.output_frame
-        <CelestialFrame(name="icrs", unit=(Unit("deg"), Unit("deg")), axes_names=('lon', 'lat'), axes_order=(0, 1), reference_frame=<ICRS Frame>)>
+  >>> print(wcsobj.available_frames)
+  ['detector', 'undistorted_frame', 'icrs']
+  >>> wcsobj.input_frame
+  <Frame2D(name="detector", unit=(Unit("pix"), Unit("pix")), axes_names=('x', 'y'), axes_order=(0, 1))>
+  >>> wcsobj.output_frame
+  <CelestialFrame(name="icrs", unit=(Unit("deg"), Unit("deg")), axes_names=('lon', 'lat'), axes_order=(0, 1), reference_frame=<ICRS 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`.
 
 .. doctest-skip::
-   
+
   >>> skycoord = wcsobj(1, 2, with_units=True)
   >>> print(skycoord)
   <SkyCoord (ICRS): (ra, dec) in deg
-    (5.52886119, -72.05285219)>
+      (5.50090023, -72.04553535)>
   >>> print(skycoord.transform_to("galactic"))
   <SkyCoord (Galactic): (l, b) in deg
-    (306.11346489, -44.89382103)>
+      (306.12713109, -44.8996588)>
+
+Using Bounding Box
+------------------
 
 The WCS object has an attribute :attr:`~gwcs.WCS.bounding_box`
 (default value of ``None``) which describes the range of
 acceptable values for each input axis.
 
 .. doctest-skip::
-   
+
   >>> wcsobj.bounding_box = ((0, 2048), (0, 1000))
   >>> wcsobj((2,3), (1020, 980))
-      array([nan, 133.48248429]), array([nan, -11.24021056])
-      
+  [array([       nan, 5.54527989]), array([         nan, -72.06454341])]
+
 The WCS object accepts a boolean flag called ``with_bounding_box`` with default value of
 ``True``. Output values which are outside the ``bounding_box`` are set to ``NaN``.
 There are cases when this is not desirable and ``with_bounding_box=False`` should be passes.
 
 Calling the :meth:`~gwcs.WCS.footprint` returns the footprint on the sky.
 
 .. doctest-skip::
-   
+
    >>> wcsobj.footprint()
-
+
+
+Manipulating Transforms
+-----------------------
+
 Some methods allow managing the transforms in a more detailed manner.
 
 Transforms between frames can be retrieved and evaluated separately.
 
 .. doctest-skip::
 
-
   >>> dist = wcsobj.get_transform('detector', 'undistorted_frame')
   >>> dist(1, 2)    # doctest: +FLOAT_CMP
-      (47.8, 95.60)
+  (-292.4150238489997, -616.8680129899999)
 
 Transforms in the pipeline can be replaced by new transforms.
 
 .. doctest-skip::
-   
+
   >>> new_transform = models.Shift(1) & models.Shift(1.5) | distortion
-  >>> wcsobj.set_transform('detector', 'focal_frame', new_transform)
+  >>> wcsobj.set_transform('detector', 'undistorted_frame', new_transform)
   >>> wcsobj(1, 2)         # doctest: +FLOAT_CMP
-      (5.5583005430002785, -72.06028278184611)
-
+  (5.501064280097802, -72.04557376712566)
 
 A transform can be inserted before or after a frame in the pipeline.
 
 .. doctest-skip::
-   
+
   >>> scale = models.Scale(2) & models.Scale(1)
   >>> wcsobj.insert_transform('icrs', scale, after=False)
   >>> wcsobj(1, 2)          # doctest: +FLOAT_CMP
-      (11.116601086000557, -72.06028278184611)
+  (11.002128560195604, -72.04557376712566)
+
+
+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.
+
+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
+inverse transformation converts world coordinates back to image coordinates.
+
+.. doctest-skip::
+
+  >>> wcsobj = asdf.open(get_pkg_data_filename('imaging_wcs_wdist.asdf')).tree['wcs']
+
+The most general method available for computing inverse coordinate
+transformation is the `WCS.invert() <gwcs.wcs.WCS.invert>`
+method. This method uses automatic or user-supplied analytical inverses whenever
+available to convert coordinates from the output frame to the input frame.
+When analytical inverse is not available as is the case for the ``wcsobj`` above,
+a numerical solution will be attempted using
+`WCS.numerical_inverse() <gwcs.wcs.WCS.numerical_inverse>`.
+
+Default parameters used by `WCS.numerical_inverse() <gwcs.wcs.WCS.numerical_inverse>`
+or `WCS.invert() <gwcs.wcs.WCS.invert>` methods should be acceptable in
+most situations:
+
+.. doctest-skip::
+
+  >>> world = wcsobj(350, 200)
+  >>> print(wcsobj.invert(*world))  # convert a single point
+  (349.9999994163172, 200.00000017679295)
+  >>> world = wcsobj([2, 350, -5000], [2, 200, 6000])
+  >>> print(wcsobj.invert(*world))  # convert multiple points at once
+  (array([ 2.00000000e+00,  3.49999999e+02, -5.00000000e+03]), array([1.99999972e+00, 2.00000002e+02, 6.00000000e+03])
+
+By default, parameter ``quiet`` is set to `True` in `WCS.numerical_inverse() <gwcs.wcs.WCS.numerical_inverse>`
+and so it will return results "as is" without warning us about possible loss
+of accuracy or about divergence of the iterative process.
+
+In order to catch these kind of errors that can occur during numerical
+inversion, we need to turn off ``quiet`` mode and be prepared to catch
+`gwcs.wcs.WCS.NoConvergence` exceptions. In the next example, let's also add a
+point far away from the image for which numerical inverse fails.
+
+.. doctest-skip::
 
+  >>> from gwcs import NoConvergence
+  >>> world = wcsobj([-85000, 2, 350, 3333, -5000], [-55000, 2, 200, 1111, 6000],
+  ...                with_bounding_box=False)
+  >>> try:
+  ...     x, y = wcsobj.invert(*world, quiet=False, maxiter=40,
+  ...                          detect_divergence=True, with_bounding_box=False)
+  ... except NoConvergence as e:
+  ...     print(f"Indices of diverging points: {e.divergent}")
+  ...     print(f"Indices of poorly converging points: {e.slow_conv}")
+  ...     print(f"Best solution:\n{e.best_solution}")
+  ...     print(f"Achieved accuracy:\n{e.accuracy}")
+  Indices of diverging points: [0]
+  Indices of poorly converging points: [4]
+  Best solution:
+  [[ 1.38600585e+11  6.77595594e+11]
+   [ 2.00000000e+00  1.99999972e+00]
+   [ 3.49999999e+02  2.00000002e+02]
+   [ 3.33300000e+03  1.11100000e+03]
+   [-4.99999985e+03  5.99999985e+03]]
+  Achieved accuracy:
+  [[8.56497375e+02 5.09216089e+03]
+   [6.57962988e-06 3.70445289e-07]
+   [5.31656943e-06 2.72052603e-10]
+   [6.81557583e-06 1.06560533e-06]
+   [3.96365344e-04 6.41822468e-05]]

gwcs/api.py

@@ -296,7 +296,10 @@ def pixel_to_world(self, *pixel_arrays):
         Convert pixel values to world coordinates.
         """
         pixels = self._sanitize_pixel_inputs(*pixel_arrays)
-        return self(*pixels, with_units=True)
+        result = self(*pixels, with_units=True)
+        if self.output_frame.naxes == 1:
+            return result[0]
+        return result
 
     def array_index_to_world(self, *index_arrays):
         """

gwcs/tests/test_wcs.py

@@ -269,6 +269,14 @@ def test_bounding_box():
     with pytest.raises(ValueError):
         w.bounding_box = ((1, 5), (2, 6))
 
+    # Test that bounding_box with quantities can be assigned and evaluates
+    bb = ((1 * u.pix, 5 * u.pix), (2 * u.pix, 6 * u.pix))
+    trans = models.Shift(10 * u .pix) & models.Shift(2 * u.pix)
+    pipeline = [('detector', trans), ('sky', None)]
+    w = wcs.WCS(pipeline)
+    w.bounding_box = bb
+    assert_allclose(w(-1*u.pix, -1*u.pix), (np.nan, np.nan))
+
 
 def test_grid_from_bounding_box():
     bb = ((-1, 9.9), (6.5, 15))
@@ -797,3 +805,26 @@ def test_iter_inv():
     xp, yp = e.value.best_solution.T
     assert np.allclose((x[1:], y[1:]), (xp[1:], yp[1:]))
     assert e.value.divergent[0] == 0
+
+
+def test_tabular_2d_quantity():
+    shape = (3, 3)
+    data = np.arange(np.product(shape)).reshape(shape) * u.m / u.s
+
+    # The integer location is at the centre of the pixel.
+    points_unit = u.pix
+    points = [(np.arange(size) - 0) * points_unit for size in shape]
+
+    kwargs = {
+        'bounds_error': False,
+        'fill_value': np.nan,
+        'method': 'nearest',
+    }
+
+    forward = models.Tabular2D(points, data, **kwargs)
+    input_frame = cf.CoordinateFrame(2, ("PIXEL", "PIXEL"), (0,1), unit=(u.pix, u.pix), name="detector")
+    output_frame = cf.CoordinateFrame(1, "CUSTOM", (0,), unit=(u.m/u.s,))
+    w = wcs.WCS(forward_transform=forward, input_frame=input_frame, output_frame=output_frame)
+
+    bb = w.bounding_box
+    assert all(u.allclose(u.Quantity(b), [0, 2] * u.pix) for b in bb)