From 7e2dbea7d51f56d99cbfd2610fb3092b55c4c642 Mon Sep 17 00:00:00 2001 From: Nadia Dencheva Date: Fri, 13 Nov 2020 09:37:43 -0500 Subject: [PATCH 01/13] start change log for 0.16.0 --- CHANGES.rst | 3 +++ 1 file changed, 3 insertions(+) diff --git a/CHANGES.rst b/CHANGES.rst index abcbff7d..9e7f04d1 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -1,3 +1,6 @@ +0.16.0 (Unreleased) +------------------- + 0.15.0 (2020-11-13) ------------------- New Features From f198b48354ca49f395985dd418ab6ee5948bc59c Mon Sep 17 00:00:00 2001 From: Stuart Mumford Date: Thu, 19 Nov 2020 12:06:22 +0000 Subject: [PATCH 02/13] APE 14 says return for pixel_to_world should not be a tuple if 1D --- gwcs/api.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/gwcs/api.py b/gwcs/api.py index c15891d5..d07aa566 100644 --- a/gwcs/api.py +++ b/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): """ From ce4bc3ef392722e422996a02fbc58e4ec3e2f025 Mon Sep 17 00:00:00 2001 From: Mihai Cara Date: Wed, 16 Dec 2020 17:40:43 -0500 Subject: [PATCH 03/13] Allow users to specify reference point in to_fits_sip --- CHANGES.rst | 13 ++++++++++++- gwcs/wcs.py | 32 +++++++++++++++++++------------- 2 files changed, 31 insertions(+), 14 deletions(-) diff --git a/CHANGES.rst b/CHANGES.rst index 9e7f04d1..b4732da5 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -1,5 +1,16 @@ -0.16.0 (Unreleased) +0.15.1 (Unreleased) ------------------- +New Features +^^^^^^^^^^^^ + +- Added an option to `to_fits_sip()` to be able to specify the reference + point (``crpix``) of the FITS WCS. [#337] + +Bug Fixes +^^^^^^^^^ + +- Fix a formula for estimating ``crpix`` in ``to_fits_sip()`` so that ``crpix`` + is near the center of the bounding box. [#337] 0.15.0 (2020-11-13) ------------------- diff --git a/gwcs/wcs.py b/gwcs/wcs.py index 9297d839..7d54bcd4 100644 --- a/gwcs/wcs.py +++ b/gwcs/wcs.py @@ -15,7 +15,6 @@ from .api import GWCSAPIMixin from . import coordinate_frames from .utils import CoordinateFrameError -from .utils import _toindex from . import utils from gwcs import coordinate_frames as cf @@ -1310,7 +1309,7 @@ def _order_clockwise(v): vertices = np.array(list(itertools.product(*bb))).T if center: - vertices = _toindex(vertices) + vertices = utils._toindex(vertices) result = np.asarray(self.__call__(*vertices, **{'with_bounding_box': False})) @@ -1366,7 +1365,7 @@ def fix_inputs(self, fixed): def to_fits_sip(self, bounding_box=None, max_pix_error=0.25, degree=None, max_inv_pix_error=0.25, inv_degree=None, - npoints=32, verbose=False): + npoints=32, crpix=None, verbose=False): """ Construct a SIP-based approximation to the WCS in the form of a FITS header @@ -1397,6 +1396,10 @@ def to_fits_sip(self, bounding_box=None, max_pix_error=0.25, degree=None, npoints : int, optional The number of points in each dimension to sample the bounding box for use in the SIP fit. + crpix : list of float, None, optional + Coordinates of the reference point for the new FITS WCS. When not + provided, i.e., when set to `None` (default) the reference pixel + will be chosen near the center of the bounding box. verbose : bool, optional print progress of fits @@ -1432,13 +1435,18 @@ def to_fits_sip(self, bounding_box=None, max_pix_error=0.25, degree=None, bounding_box = self.bounding_box (xmin, xmax), (ymin, ymax) = bounding_box - crpix1 = (xmax - xmin) // 2 - crpix2 = (ymax - ymin) // 2 + if crpix is None: + crpix1 = round((xmax + xmin) / 2, 1) + crpix2 = round((ymax + ymin) / 2, 1) + else: + crpix1 = crpix[0] - 1 + crpix2 = crpix[1] - 1 + crval1, crval2 = transform(crpix1, crpix2) hdr = fits.Header() hdr['naxis'] = 2 - hdr['naxis1'] = xmax - hdr['naxis2'] = ymax + hdr['naxis1'] = int(xmax) + 1 + hdr['naxis2'] = int(ymax) + 1 hdr['ctype1'] = 'RA---TAN-SIP' hdr['ctype2'] = 'DEC--TAN-SIP' hdr['CRPIX1'] = crpix1 + 1 @@ -1678,16 +1686,14 @@ def _calc_approx_inv(self, max_inv_pix_error=5, inv_degree=None, npoints=16): # A bounding_box is needed to proceed. return - (xmin, xmax), (ymin, ymax) = self.bounding_box - crpix1 = (xmax - xmin) // 2 - crpix2 = (ymax - ymin) // 2 - crval1, crval2 = self.forward_transform(crpix1, crpix2) + crpix = np.mean(self.bounding_box, axis=1) + crval1, crval2 = self.forward_transform(*crpix) # Rotate to native system and deproject. Set center of the projection # transformation to the middle of the bounding box ("image") in order # to minimize projection effects across the entire image, # thus the initial shift. - ntransform = ((Shift(crpix1) & Shift(crpix2)) | self.forward_transform + ntransform = ((Shift(crpix[0]) & Shift(crpix[1])) | self.forward_transform | RotateCelestial2Native(crval1, crval2, 180) | Sky2Pix_TAN()) @@ -1710,7 +1716,7 @@ def _calc_approx_inv(self, max_inv_pix_error=5, inv_degree=None, npoints=16): self._approx_inverse = (RotateCelestial2Native(crval1, crval2, 180) | Sky2Pix_TAN() | Mapping((0, 1, 0, 1)) | (fit_inv_poly_u & fit_inv_poly_v) | - (Shift(crpix1) & Shift(crpix2))) + (Shift(crpix[0]) & Shift(crpix[1]))) def _fit_2D_poly(ntransform, npoints, degree, max_error, From 752cfc34c9ff57aa263302a8d02e6d934dda5d2a Mon Sep 17 00:00:00 2001 From: Stuart Mumford Date: Thu, 19 Nov 2020 11:11:25 +0000 Subject: [PATCH 04/13] add test for quantity bounding box --- gwcs/tests/test_wcs.py | 23 +++++++++++++++++++++++ 1 file changed, 23 insertions(+) diff --git a/gwcs/tests/test_wcs.py b/gwcs/tests/test_wcs.py index e72f1137..83a5367d 100644 --- a/gwcs/tests/test_wcs.py +++ b/gwcs/tests/test_wcs.py @@ -797,3 +797,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) From 09f4c2cea4ca6910d312f7b645d4212df81fb600 Mon Sep 17 00:00:00 2001 From: Stuart Mumford Date: Thu, 19 Nov 2020 11:15:58 +0000 Subject: [PATCH 05/13] Do not assume model bounding boxes can be cast to numpy arrays --- gwcs/wcs.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/gwcs/wcs.py b/gwcs/wcs.py index 7d54bcd4..05a092be 100644 --- a/gwcs/wcs.py +++ b/gwcs/wcs.py @@ -313,11 +313,11 @@ def __call__(self, *args, **kwargs): if self.bounding_box is not None: # Currently compound models do not attempt to combine individual model - # bounding boxes. Get the forward transform and assign the ounding_box to it + # bounding boxes. Get the forward transform and assign the bounding_box to it # before evaluating it. The order Model.bounding_box is reversed. axes_ind = self._get_axes_indices() if transform.n_inputs > 1: - transform.bounding_box = np.array(self.bounding_box)[axes_ind][::-1] + transform.bounding_box = [self.bounding_box[ind] for ind in axes_ind][::-1] else: transform.bounding_box = self.bounding_box result = transform(*args, **kwargs) @@ -1204,8 +1204,7 @@ def bounding_box(self): except AttributeError: axes_order = np.arange(transform_0.n_inputs) # Model.bounding_box is in python order, need to reverse it first. - bb = np.array(bb[::-1])[np.array(axes_order)] - return tuple(tuple(item) for item in bb) + return tuple(bb[::-1][i] for i in axes_order) @bounding_box.setter def bounding_box(self, value): From bf30c83bb5370d8c4873c0d04ca641592e5cb94f Mon Sep 17 00:00:00 2001 From: Nadia Dencheva Date: Fri, 18 Dec 2020 09:39:03 -0500 Subject: [PATCH 06/13] fix setting bounding_box with a tuple of quantities --- CHANGES.rst | 2 ++ gwcs/tests/test_wcs.py | 8 ++++++++ gwcs/wcs.py | 4 +++- 3 files changed, 13 insertions(+), 1 deletion(-) diff --git a/CHANGES.rst b/CHANGES.rst index b4732da5..e161da27 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -12,6 +12,8 @@ Bug Fixes - Fix a formula for estimating ``crpix`` in ``to_fits_sip()`` so that ``crpix`` is near the center of the bounding box. [#337] +- ``bounding_box`` now works with tuple of ``Quantities``. [#331] + 0.15.0 (2020-11-13) ------------------- New Features diff --git a/gwcs/tests/test_wcs.py b/gwcs/tests/test_wcs.py index 83a5367d..5e48281e 100644 --- a/gwcs/tests/test_wcs.py +++ b/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)) diff --git a/gwcs/wcs.py b/gwcs/wcs.py index 05a092be..77915939 100644 --- a/gwcs/wcs.py +++ b/gwcs/wcs.py @@ -1219,6 +1219,7 @@ def bounding_box(self, value): value : tuple or None Tuple of tuples with ("low", high") values for the range. """ + print('IN SETTING BB', value) frames = self.available_frames transform_0 = self.get_transform(frames[0], frames[1]) if value is None: @@ -1235,7 +1236,8 @@ def bounding_box(self, value): transform_0.bounding_box = value else: # The axes in bounding_box in modeling follow python order - transform_0.bounding_box = np.array(value)[axes_ind][::-1] + #transform_0.bounding_box = np.array(value)[axes_ind][::-1] + transform_0.bounding_box = [value[ind] for ind in axes_ind][::-1] self.set_transform(frames[0], frames[1], transform_0) def _get_axes_indices(self): From 0e4a5a0c3b462dbecf5bdc4f1fee6e701e353144 Mon Sep 17 00:00:00 2001 From: Nadia Dencheva Date: Fri, 18 Dec 2020 09:58:16 -0500 Subject: [PATCH 07/13] remove print statement --- gwcs/wcs.py | 1 - 1 file changed, 1 deletion(-) diff --git a/gwcs/wcs.py b/gwcs/wcs.py index 77915939..492017a0 100644 --- a/gwcs/wcs.py +++ b/gwcs/wcs.py @@ -1219,7 +1219,6 @@ def bounding_box(self, value): value : tuple or None Tuple of tuples with ("low", high") values for the range. """ - print('IN SETTING BB', value) frames = self.available_frames transform_0 = self.get_transform(frames[0], frames[1]) if value is None: From 6e49b92d851a40c888ab05dc027a405614c20093 Mon Sep 17 00:00:00 2001 From: Nadia Dencheva Date: Fri, 18 Dec 2020 10:07:15 -0500 Subject: [PATCH 08/13] update numpy --- .travis.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.travis.yml b/.travis.yml index e2bae7de..8fdf4b5c 100644 --- a/.travis.yml +++ b/.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] From 0794b41a2dc81a8cb696540092ec5d48c0130613 Mon Sep 17 00:00:00 2001 From: Mihai Cara Date: Wed, 16 Dec 2020 20:40:53 -0500 Subject: [PATCH 09/13] Allow sub-pixel sampling when fitting SIP --- CHANGES.rst | 5 +++- gwcs/wcs.py | 72 ++++++++++++++++++++++++++---------------------- gwcs/wcstools.py | 5 +++- 3 files changed, 47 insertions(+), 35 deletions(-) diff --git a/CHANGES.rst b/CHANGES.rst index e161da27..ae9a2624 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -9,10 +9,13 @@ New Features 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] -- ``bounding_box`` now works with tuple of ``Quantities``. [#331] +- Allow sub-pixel sampling of the WCS model when computing SIP approximation in + ``to_fits_sip()``. [#338] 0.15.0 (2020-11-13) ------------------- diff --git a/gwcs/wcs.py b/gwcs/wcs.py index 492017a0..19e33fbc 100644 --- a/gwcs/wcs.py +++ b/gwcs/wcs.py @@ -13,10 +13,10 @@ import astropy.io.fits as fits from .api import GWCSAPIMixin -from . import coordinate_frames +from . import coordinate_frames as cf from .utils import CoordinateFrameError from . import utils -from gwcs import coordinate_frames as cf +from .wcstools import grid_from_bounding_box try: @@ -251,7 +251,7 @@ def _get_frame_index(self, frame): """ Return the index in the pipeline where this frame is locate. """ - if isinstance(frame, coordinate_frames.CoordinateFrame): + if isinstance(frame, cf.CoordinateFrame): frame = frame.name #frame_names = [getattr(item[0], "name", item[0]) for item in self._pipeline] frame_names = [step.frame if isinstance(step.frame, str) else step.frame.name for step in self._pipeline] @@ -1395,13 +1395,13 @@ def to_fits_sip(self, bounding_box=None, max_pix_error=0.25, degree=None, is ignored. npoints : int, optional The number of points in each dimension to sample the bounding box - for use in the SIP fit. + for use in the SIP fit. Minimum number of points is 3. crpix : list of float, None, optional - Coordinates of the reference point for the new FITS WCS. When not - provided, i.e., when set to `None` (default) the reference pixel - will be chosen near the center of the bounding box. + Coordinates (1-based) of the reference point for the new FITS WCS. + When not provided, i.e., when set to `None` (default) the reference + pixel will be chosen near the center of the bounding box. verbose : bool, optional - print progress of fits + Print progress of fits. Returns ------- @@ -1421,13 +1421,16 @@ def to_fits_sip(self, bounding_box=None, max_pix_error=0.25, degree=None, higher degrees (~7 or higher) will typically fail due floating point problems that arise with high powers. - """ if not isinstance(self.output_frame, cf.CelestialFrame): raise ValueError( "The to_fits_sip method only works with celestial frame transforms") + if npoints < 8: + raise ValueError("Number of sampling points is too small. 'npoints' must be >= 8.") + transform = self.forward_transform + # Determine reference points. if bounding_box is None and self.bounding_box is None: raise ValueError("A bounding_box is needed to proceed.") @@ -1442,6 +1445,10 @@ def to_fits_sip(self, bounding_box=None, max_pix_error=0.25, degree=None, crpix1 = crpix[0] - 1 crpix2 = crpix[1] - 1 + # check that the bounding box has some reasonable size: + if (xmax - xmin) < 1 or (ymax - ymin) < 1: + raise ValueError("Bounding box is too small for fitting a SIP polynomial") + crval1, crval2 = transform(crpix1, crpix2) hdr = fits.Header() hdr['naxis'] = 2 @@ -1463,11 +1470,15 @@ def to_fits_sip(self, bounding_box=None, max_pix_error=0.25, degree=None, ntransform = ((Shift(crpix1) & Shift(crpix2)) | transform | RotateCelestial2Native(crval1, crval2, 180) | Sky2Pix_TAN()) - u, v = _make_sampling_grid(npoints, bounding_box) + + # standard sampling: + u, v = _make_sampling_grid(npoints, bounding_box, crpix=[crpix1, crpix2]) undist_x, undist_y = ntransform(u, v) + # Double sampling to check if sampling is sufficient. - ud, vd = _make_sampling_grid(2 * npoints, bounding_box) + ud, vd = _make_sampling_grid(2 * npoints, bounding_box, crpix=[crpix1, crpix2]) undist_xd, undist_yd = ntransform(ud, vd) + # Determine approximate pixel scale in order to compute error threshold # from the specified pixel error. Computed at the center of the array. x0, y0 = ntransform(0, 0) @@ -1476,12 +1487,16 @@ def to_fits_sip(self, bounding_box=None, max_pix_error=0.25, degree=None, pixarea = np.abs((xx - x0) * (yy - y0) - (xy - y0) * (yx - x0)) plate_scale = np.sqrt(pixarea) max_error = max_pix_error * plate_scale + # The fitting section. - fit_poly_x, fit_poly_y, max_resid = _fit_2D_poly(ntransform, npoints, - degree, max_error, - u, v, undist_x, undist_y, - ud, vd, undist_xd, undist_yd, - verbose=verbose) + fit_poly_x, fit_poly_y, max_resid = _fit_2D_poly( + ntransform, npoints, + degree, max_error, + u, v, undist_x, undist_y, + ud, vd, undist_xd, undist_yd, + verbose=verbose + ) + # The following is necessary to put the fit into the SIP formalism. cdmat, sip_poly_x, sip_poly_y = _reform_poly_coefficients(fit_poly_x, fit_poly_y) # cdmat = np.array([[fit_poly_x.c1_0.value, fit_poly_x.c0_1.value], @@ -1687,6 +1702,7 @@ def _calc_approx_inv(self, max_inv_pix_error=5, inv_degree=None, npoints=16): return crpix = np.mean(self.bounding_box, axis=1) + crval1, crval2 = self.forward_transform(*crpix) # Rotate to native system and deproject. Set center of the projection @@ -1697,11 +1713,12 @@ def _calc_approx_inv(self, max_inv_pix_error=5, inv_degree=None, npoints=16): | RotateCelestial2Native(crval1, crval2, 180) | Sky2Pix_TAN()) - u, v = _make_sampling_grid(npoints, self.bounding_box) + # standard sampling: + u, v = _make_sampling_grid(npoints, self.bounding_box, crpix=crpix) undist_x, undist_y = ntransform(u, v) # Double sampling to check if sampling is sufficient. - ud, vd = _make_sampling_grid(2 * npoints, self.bounding_box) + ud, vd = _make_sampling_grid(2 * npoints, self.bounding_box, crpix=crpix) undist_xd, undist_yd = ntransform(ud, vd) fit_inv_poly_u, fit_inv_poly_v, max_inv_resid = _fit_2D_poly( @@ -1763,21 +1780,10 @@ def _fit_2D_poly(ntransform, npoints, degree, max_error, return fit_poly_x, fit_poly_y, max_resid -def _make_sampling_grid(npoints, bounding_box): - (xmin, xmax), (ymin, ymax) = bounding_box - xsize = xmax - xmin - ysize = ymax - ymin - crpix1 = int(xsize / 2) - crpix2 = int(ysize / 2) - stepsize_x = int(xsize / npoints) - stepsize_y = int(ysize / npoints) - # Ensure last row and column are part of the evaluation grid. - y, x = np.mgrid[: ymax + 1: stepsize_y, : xmax + 1: stepsize_x] - x[:, -1] = xsize - 1 - y[-1, :] = ysize - 1 - u = x - crpix1 - v = y - crpix2 - return u, v +def _make_sampling_grid(npoints, bounding_box, crpix): + step = np.subtract.reduce(bounding_box, axis=1) / (1.0 - npoints) + crpix = np.asanyarray(crpix)[:, None, None] + return grid_from_bounding_box(bounding_box, step=step, center=False) - crpix def _compute_distance_residual(undist_x, undist_y, fit_poly_x, fit_poly_y): diff --git a/gwcs/wcstools.py b/gwcs/wcstools.py index 3041fb9b..6de9d13a 100644 --- a/gwcs/wcstools.py +++ b/gwcs/wcstools.py @@ -8,7 +8,6 @@ from astropy import coordinates as coord from astropy import units as u -from .wcs import WCS from .coordinate_frames import * # noqa from .utils import UnsupportedTransformError, UnsupportedProjectionError from .utils import _compute_lon_pole @@ -53,6 +52,8 @@ def wcs_from_fiducial(fiducial, coordinate_frame=None, projection=None, input_frame : ~gwcs.coordinate_frames.CoordinateFrame` The input coordinate frame. """ + from .wcs import WCS + if transform is not None: if not isinstance(transform, Model): raise UnsupportedTransformError("Expected transform to be an instance" @@ -252,6 +253,8 @@ def wcs_from_points(xy, world_coordinates, fiducial, wcsobj : `~gwcs.wcs.WCS` a WCS object for this observation. """ + from .wcs import WCS + supported_poly_types = {"polynomial": models.Polynomial2D, "chebyshev": models.Chebyshev2D, "legendre": models.Legendre2D From 39a2b76945a9960c2f13d4ec52df0f3511ef0d16 Mon Sep 17 00:00:00 2001 From: Mihai Cara Date: Wed, 18 Nov 2020 00:37:25 -0500 Subject: [PATCH 10/13] Add examples of calling numerical_inverse --- .travis.yml | 2 +- docs/gwcs/imaging_with_distortion.rst | 31 ++++-- docs/gwcs/using_wcs.rst | 138 ++++++++++++++++++++------ gwcs/wcs.py | 71 ++++++++++++- 4 files changed, 201 insertions(+), 41 deletions(-) diff --git a/.travis.yml b/.travis.yml index 8fdf4b5c..e9193461 100644 --- a/.travis.yml +++ b/.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: diff --git a/docs/gwcs/imaging_with_distortion.rst b/docs/gwcs/imaging_with_distortion.rst index 47b950f7..ff66605d 100644 --- a/docs/gwcs/imaging_with_distortion.rst +++ b/docs/gwcs/imaging_with_distortion.rst @@ -26,7 +26,7 @@ The imaging example without units: >>> 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 +46,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 +73,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") diff --git a/docs/gwcs/using_wcs.rst b/docs/gwcs/using_wcs.rst index d6ac1b2d..7762b186 100644 --- a/docs/gwcs/using_wcs.rst +++ b/docs/gwcs/using_wcs.rst @@ -1,58 +1,65 @@ -.. _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 - - >>> wcsobj.output_frame - )> + >>> print(wcsobj.available_frames) + ['detector', 'undistorted_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`. .. doctest-skip:: - + >>> skycoord = wcsobj(1, 2, with_units=True) >>> print(skycoord) + (5.50090023, -72.04553535)> >>> print(skycoord.transform_to("galactic")) + (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. @@ -60,36 +67,111 @@ There are cases when this is not desirable and ``with_bounding_box=False`` shoul 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() ` +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() `. + +Default parameters used by `WCS.numerical_inverse() ` +or `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() ` +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]] diff --git a/gwcs/wcs.py b/gwcs/wcs.py index 492017a0..5049512b 100644 --- a/gwcs/wcs.py +++ b/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." From a5a562101cbcf3ab80568793add1d2eeb77d7d72 Mon Sep 17 00:00:00 2001 From: Mihai Cara Date: Fri, 18 Dec 2020 13:02:50 -0500 Subject: [PATCH 11/13] Define crpix --- docs/gwcs/imaging_with_distortion.rst | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/docs/gwcs/imaging_with_distortion.rst b/docs/gwcs/imaging_with_distortion.rst index ff66605d..a3a870b7 100644 --- a/docs/gwcs/imaging_with_distortion.rst +++ b/docs/gwcs/imaging_with_distortion.rst @@ -19,7 +19,8 @@ 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) From 137280253084a3a970850f371e777cb5b89ef0d3 Mon Sep 17 00:00:00 2001 From: Mihai Cara Date: Fri, 18 Dec 2020 12:39:32 -0500 Subject: [PATCH 12/13] Allow custom degree range for SIP approximation --- CHANGES.rst | 5 +++++ gwcs/wcs.py | 50 +++++++++++++++++++++++++++++++++++++++++--------- 2 files changed, 46 insertions(+), 9 deletions(-) diff --git a/CHANGES.rst b/CHANGES.rst index ae9a2624..4eaf2a5c 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -6,6 +6,8 @@ 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 ^^^^^^^^^ @@ -17,6 +19,9 @@ Bug Fixes - 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 diff --git a/gwcs/wcs.py b/gwcs/wcs.py index 19e33fbc..de5d57b0 100644 --- a/gwcs/wcs.py +++ b/gwcs/wcs.py @@ -1380,26 +1380,50 @@ def to_fits_sip(self, bounding_box=None, max_pix_error=0.25, degree=None, A pair of tuples, each consisting of two numbers Represents the range of pixel values in both dimensions ((xmin, xmax), (ymin, ymax)) + max_pix_error : float, optional Maximum allowed error over the domain of the pixel array. This error is the equivalent pixel error that corresponds to the maximum error in the output coordinate resulting from the fit based on a nominal plate scale. - degree : int, optional - Degree of the SIP polynomial. If supplied, max_pixel_error is ignored. + + degree : int, iterable, None, optional + Degree of the SIP polynomial. Default value `None` indicates that + all allowed degree values (``[1...9]``) will be considered and + the lowest degree that meets accuracy requerements set by + ``max_pix_error`` will be returned. Alternatively, ``degree`` can be + an iterable containing allowed values for the SIP polynomial degree. + This option is similar to default `None` but it allows caller to + restrict the range of allowed SIP degrees used for fitting. + Finally, ``degree`` can be an integer indicating the exact SIP degree + to be fit to the WCS transformation. In this case + ``max_pixel_error`` is ignored. + max_inv_error : float, optional Maximum allowed inverse error over the domain of the pixel array in pixel units. If None, no inverse is generated. - inv_degree : int, optional - Degree of the inverse SIP polynomial. If supplied max_inv_pixel_error - is ignored. + + inv_degree : int, iterable, None, optional + Degree of the SIP polynomial. Default value `None` indicates that + all allowed degree values (``[1...9]``) will be considered and + the lowest degree that meets accuracy requerements set by + ``max_pix_error`` will be returned. Alternatively, ``degree`` can be + an iterable containing allowed values for the SIP polynomial degree. + This option is similar to default `None` but it allows caller to + restrict the range of allowed SIP degrees used for fitting. + Finally, ``degree`` can be an integer indicating the exact SIP degree + to be fit to the WCS transformation. In this case + ``max_inv_pixel_error`` is ignored. + npoints : int, optional The number of points in each dimension to sample the bounding box for use in the SIP fit. Minimum number of points is 3. + crpix : list of float, None, optional Coordinates (1-based) of the reference point for the new FITS WCS. When not provided, i.e., when set to `None` (default) the reference pixel will be chosen near the center of the bounding box. + verbose : bool, optional Print progress of fits. @@ -1510,7 +1534,7 @@ def to_fits_sip(self, bounding_box=None, max_pix_error=0.25, degree=None, if max_inv_pix_error: fit_inv_poly_u, fit_inv_poly_v, max_inv_resid = _fit_2D_poly(ntransform, - npoints, None, + npoints, inv_degree, max_inv_pix_error, U, V, u-U, v-V, Ud, Vd, ud-Ud, vd-Vd, @@ -1747,10 +1771,18 @@ def _fit_2D_poly(ntransform, npoints, degree, max_error, llsqfitter = LinearLSQFitter() # The case of one pass with the specified polynomial degree - if degree: - deglist = [degree] + if degree is None: + deglist = range(1, 10) + elif hasattr(degree, '__iter__'): + deglist = sorted(map(int, degree)) + if set(deglist).difference(range(1, 10)): + raise ValueError("Allowed values for SIP degree are [1...9]") else: - deglist = range(10) + degree = int(degree) + if degree < 1 or degree > 9: + raise ValueError("Allowed values for SIP degree are [1...9]") + deglist = [degree] + prev_max_error = float(np.inf) if verbose: print(f'maximum_specified_error: {max_error}') From 0a54327a1d3f48434b6d58d7f857c27c3088e127 Mon Sep 17 00:00:00 2001 From: Nadia Dencheva Date: Fri, 18 Dec 2020 17:31:20 -0500 Subject: [PATCH 13/13] prepare for release --- CHANGES.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CHANGES.rst b/CHANGES.rst index 4eaf2a5c..0cac0712 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -1,4 +1,4 @@ -0.15.1 (Unreleased) +0.16.0 (2020-12-18) ------------------- New Features ^^^^^^^^^^^^