Merge pull request #330 from mcara/add-num-inv-examples

Add examples of calling numerical_inverse

.travis.yml

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

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/wcs.py

@@ -635,6 +635,73 @@ def numerical_inverse(self, *args, **kwargs):
         ValueError
             Invalid argument values.
 
+        Examples
+        --------
+        >>> from astropy.utils.data import get_pkg_data_filename
+        >>> from gwcs import NoConvergence
+        >>> import asdf
+        >>> import numpy as np
+
+        >>> filename = get_pkg_data_filename('data/nircamwcs.asdf', package='gwcs.tests')
+        >>> w = asdf.open(filename).tree['wcs']
+
+        >>> ra, dec = w([1,2,3], [1,1,1])
+        >>> print(ra)  # doctest: +FLOAT_CMP
+        [5.927628   5.92757069 5.92751337]
+        >>> print(dec)  # doctest: +FLOAT_CMP
+        [-72.01341247 -72.01341273 -72.013413  ]
+
+        >>> x, y = w.numerical_inverse(ra, dec)
+        >>> print(x)  # doctest: +FLOAT_CMP
+        [1.00000005 2.00000005 3.00000006]
+        >>> print(y)  # doctest: +FLOAT_CMP
+        [1.00000004 0.99999979 1.00000015]
+
+        >>> x, y = w.numerical_inverse(ra, dec, maxiter=3, tolerance=1.0e-10, quiet=False)
+        Traceback (most recent call last):
+        ...
+        gwcs.wcs.NoConvergence: 'WCS.numerical_inverse' failed to converge to the
+        requested accuracy after 3 iterations.
+
+        >>> w.numerical_inverse(
+        ...     *w([1, 300000, 3], [2, 1000000, 5], with_bounding_box=False),
+        ...     adaptive=False,
+        ...     detect_divergence=True,
+        ...     quiet=False,
+        ...     with_bounding_box=False
+        ... )
+        Traceback (most recent call last):
+        ...
+        gwcs.wcs.NoConvergence: 'WCS.numerical_inverse' failed to converge to the
+        requested accuracy. After 4 iterations, the solution is diverging at
+        least for one input point.
+
+        >>> # Now try to use some diverging data:
+        >>> divradec = w([1, 300000, 3], [2, 1000000, 5], with_bounding_box=False)
+        >>> print(divradec)  # doctest: +FLOAT_CMP
+        (array([  5.92762673, 148.21600848,   5.92750827]),
+         array([-72.01339464,  -7.80968079, -72.01334172]))
+        >>> try:  # doctest: +SKIP
+        ...     x, y = w.numerical_inverse(*divradec, maxiter=20,
+        ...                                tolerance=1.0e-4, adaptive=True,
+        ...                                detect_divergence=True,
+        ...                                quiet=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: None
+        Indices of poorly converging points: [1]
+        Best solution:
+        [[1.00000040e+00 1.99999841e+00]
+         [6.33507833e+17 3.40118820e+17]
+         [3.00000038e+00 4.99999841e+00]]
+        Achieved accuracy:
+        [[2.75925982e-05 1.18471543e-05]
+         [3.65405005e+04 1.31364188e+04]
+         [2.76552923e-05 1.14789013e-05]]
+
         """
         tolerance = kwargs.get('tolerance', 1e-5)
         maxiter = kwargs.get('maxiter', 50)
@@ -963,14 +1030,14 @@ def correction(pix):
         if (ind is not None or inddiv is not None) and not quiet:
             if inddiv is None:
                 raise NoConvergence(
-                    "'WCS.invert' failed to "
+                    "'WCS.numerical_inverse' failed to "
                     "converge to the requested accuracy after {:d} "
                     "iterations.".format(k), best_solution=pix,
                     accuracy=np.abs(dpix), niter=k,
                     slow_conv=ind, divergent=None)
             else:
                 raise NoConvergence(
-                    "'WCS.invert' failed to "
+                    "'WCS.numerical_inverse' failed to "
                     "converge to the requested accuracy.\n"
                     "After {:d} iterations, the solution is diverging "
                     "at least for one input point."