From 2ac22a7ea245847919be02dc5c03670bedf1603d Mon Sep 17 00:00:00 2001 From: Nadia Dencheva Date: Tue, 15 Mar 2022 11:26:09 -0400 Subject: [PATCH 01/12] open for development --- CHANGES.rst | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/CHANGES.rst b/CHANGES.rst index 7a8e1e57..d24c16e6 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -1,3 +1,11 @@ +0.18.2 (Unreleased) +------------------- +Bug Fixes +^^^^^^^^^ + +New Features +^^^^^^^^^^^^ + 0.18.1 (2022-03-15) ------------------- Bug Fixes From 60afaa7c07e4ca3c33e4bc6c36ddd603106f5631 Mon Sep 17 00:00:00 2001 From: William Jamieson Date: Tue, 12 Apr 2022 16:32:27 -0400 Subject: [PATCH 02/12] Pin astropy min version to 5.0.4 --- setup.cfg | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/setup.cfg b/setup.cfg index 52160a49..8c511252 100644 --- a/setup.cfg +++ b/setup.cfg @@ -20,7 +20,7 @@ setup_requires = install_requires = asdf >= 2.8.1 - astropy >= 4.1 + astropy >= 5.0.4 numpy scipy asdf_wcs_schemas From 202be7dfc558c4c45e006ad1a5b718b648e27d30 Mon Sep 17 00:00:00 2001 From: William Jamieson Date: Tue, 12 Apr 2022 16:33:48 -0400 Subject: [PATCH 03/12] Update changes --- CHANGES.rst | 2 ++ 1 file changed, 2 insertions(+) diff --git a/CHANGES.rst b/CHANGES.rst index d24c16e6..13615283 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -3,6 +3,8 @@ Bug Fixes ^^^^^^^^^ +- Pin astropy min version to 5.0.4. [#404] + New Features ^^^^^^^^^^^^ From 790f8cc6e92726cbc5a0db313fedade6b421d516 Mon Sep 17 00:00:00 2001 From: William Jamieson Date: Wed, 8 Jun 2022 14:32:28 -0400 Subject: [PATCH 04/12] Add label to main page for intersphinx link --- docs/index.rst | 2 ++ 1 file changed, 2 insertions(+) diff --git a/docs/index.rst b/docs/index.rst index 031660fc..126a6b82 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -1,3 +1,5 @@ +.. _gwcs: + GWCS Documentation ================== From 99af4a7c7504ea78fdb07d240c66ef6f3c371e9f Mon Sep 17 00:00:00 2001 From: Stuart Mumford Date: Tue, 14 Jun 2022 00:10:38 +0200 Subject: [PATCH 05/12] Fix APE 14 API with 1D CoordinateFrame object (#407) * This is the way * a test * Fix world_to_pixel now inverse dosen't always return tuple * Test inverse as well --- gwcs/api.py | 13 +++++++++---- gwcs/coordinate_frames.py | 2 ++ gwcs/tests/test_api.py | 20 ++++++++++++++++++++ gwcs/wcs.py | 5 +++-- 4 files changed, 34 insertions(+), 6 deletions(-) diff --git a/gwcs/api.py b/gwcs/api.py index c15891d5..9e5012fb 100644 --- a/gwcs/api.py +++ b/gwcs/api.py @@ -312,10 +312,15 @@ def world_to_pixel(self, *world_objects): Convert world coordinates to pixel values. """ result = self.invert(*world_objects, with_units=True) - if not utils.isnumerical(result[0]): - result = [i.value for i in result] - if self.input_frame.naxes == 1: - return result[0] + + if self.input_frame.naxes > 1: + first_res = result[0] + if not utils.isnumerical(first_res): + result = [i.value for i in result] + else: + if not utils.isnumerical(result): + result = result.value + return result def world_to_array_index(self, *world_objects): diff --git a/gwcs/coordinate_frames.py b/gwcs/coordinate_frames.py index dcc60c86..9db35fd6 100644 --- a/gwcs/coordinate_frames.py +++ b/gwcs/coordinate_frames.py @@ -280,6 +280,8 @@ def axes_type(self): def coordinates(self, *args): """ Create world coordinates object""" coo = tuple([arg * un if not hasattr(arg, "to") else arg.to(un) for arg, un in zip(args, self.unit)]) + if self.naxes == 1: + return coo[0] return coo def coordinate_to_quantity(self, *coords): diff --git a/gwcs/tests/test_api.py b/gwcs/tests/test_api.py index 55dedee9..cae9bdd9 100644 --- a/gwcs/tests/test_api.py +++ b/gwcs/tests/test_api.py @@ -10,7 +10,9 @@ from astropy import time from astropy import coordinates as coord from astropy.wcs.wcsapi import HighLevelWCSWrapper +import astropy.modeling.models as m import gwcs.coordinate_frames as cf +import gwcs # Shorthand the name of the 2d gwcs fixture @@ -492,3 +494,21 @@ def test_composite_many_base_frame(): assert len(wao_components) == 2 assert not {c[0] for c in wao_components}.difference({"SPATIAL", "SPATIAL1"}) + + +def test_coordinate_frame_api(): + forward = m.Linear1D(slope=0.1*u.deg/u.pix, intercept=0*u.deg) + + output_frame = cf.CoordinateFrame(1, "SPATIAL", (0,), unit=(u.deg,), name="sepframe") + input_frame = cf.CoordinateFrame(1, "PIXEL", (0,), unit=(u.pix,)) + + wcs = gwcs.WCS(forward_transform=forward, input_frame=input_frame, output_frame=output_frame) + + world = wcs.pixel_to_world(0) + assert isinstance(world, u.Quantity) + + pixel = wcs.world_to_pixel(world) + assert isinstance(pixel, float) + + pixel2 = wcs.invert(world) + assert u.allclose(pixel2, 0*u.pix) diff --git a/gwcs/wcs.py b/gwcs/wcs.py index 730d024f..c1279c4f 100644 --- a/gwcs/wcs.py +++ b/gwcs/wcs.py @@ -480,8 +480,9 @@ def invert(self, *args, **kwargs): Returns ------- - result : tuple - Returns a tuple of scalar or array values for each axis. + result : tuple or value + Returns a tuple of scalar or array values for each axis. Unless + ``input_frame.naxes == 1`` when it shall return the value. """ with_units = kwargs.pop('with_units', False) From b4c472d9916ae73eed3b7caca551d58aa0900bf0 Mon Sep 17 00:00:00 2001 From: Nadia Dencheva Date: Mon, 13 Jun 2022 07:54:07 -0700 Subject: [PATCH 06/12] fix failure with dev version of astropy. --- gwcs/converters/tests/test_wcs.py | 7 ------- 1 file changed, 7 deletions(-) diff --git a/gwcs/converters/tests/test_wcs.py b/gwcs/converters/tests/test_wcs.py index 65488106..d8738864 100644 --- a/gwcs/converters/tests/test_wcs.py +++ b/gwcs/converters/tests/test_wcs.py @@ -2,7 +2,6 @@ # -*- coding: utf-8 -*- import os.path import pytest -import warnings astropy = pytest.importorskip('astropy', minversion='3.0') @@ -103,12 +102,6 @@ def test_composite_frame(tmpdir): def create_test_frames(): """Creates an array of frames to be used for testing.""" - # Suppress warnings from astropy that are caused by having 'dubious' dates - # that are too far in the future. It's not a concern for the purposes of - # unit tests. See issue #5809 on the astropy GitHub for discussion. - from astropy._erfa import ErfaWarning - warnings.simplefilter("ignore", ErfaWarning) - frames = [ cf.CelestialFrame(reference_frame=coord.ICRS()), From f7e62c739de5a74942f2c56bcad4b22904841d63 Mon Sep 17 00:00:00 2001 From: Nadia Dencheva Date: Mon, 13 Jun 2022 09:47:06 -0700 Subject: [PATCH 07/12] fix warnings from tests --- gwcs/tests/test_api_slicing.py | 3 +- gwcs/tests/test_wcs.py | 74 +++++++++++++++++++++++----------- 2 files changed, 53 insertions(+), 24 deletions(-) diff --git a/gwcs/tests/test_api_slicing.py b/gwcs/tests/test_api_slicing.py index caadf048..d4866242 100644 --- a/gwcs/tests/test_api_slicing.py +++ b/gwcs/tests/test_api_slicing.py @@ -1,6 +1,7 @@ + import astropy.units as u from astropy.coordinates import Galactic, SkyCoord, SpectralCoord -from astropy.wcs.wcsapi.sliced_low_level_wcs import SlicedLowLevelWCS +from astropy.wcs.wcsapi.wrappers import SlicedLowLevelWCS from numpy.testing import assert_allclose, assert_equal EXPECTED_ELLIPSIS_REPR = """ diff --git a/gwcs/tests/test_wcs.py b/gwcs/tests/test_wcs.py index 785d1dad..017f0988 100644 --- a/gwcs/tests/test_wcs.py +++ b/gwcs/tests/test_wcs.py @@ -32,9 +32,9 @@ spec = cf.SpectralFrame(name='wave', unit=[u.m, ], axes_order=(2, ), axes_names=('lambda', )) time = cf.TemporalFrame(name='time', unit=[u.s, ], axes_order=(3, ), axes_names=('time', ), reference_frame=Time("2020-01-01")) -pipe = [(detector, m1), - (focal, m2), - (icrs, None) +pipe = [wcs.Step(detector, m1), + wcs.Step(focal, m2), + wcs.Step(icrs, None) ] # Create some data. @@ -75,15 +75,19 @@ def test_init_no_transform(): gw = wcs.WCS(output_frame='icrs') assert len(gw._pipeline) == 2 assert gw.pipeline[0].frame == "detector" - assert gw.pipeline[0][0] == "detector" + with pytest.warns(DeprecationWarning, match="Indexing a WCS.pipeline step is deprecated."): + assert gw.pipeline[0][0] == "detector" assert gw.pipeline[1].frame == "icrs" - assert gw.pipeline[1][0] == "icrs" + with pytest.warns(DeprecationWarning, match="Indexing a WCS.pipeline step is deprecated."): + assert gw.pipeline[1][0] == "icrs" assert np.in1d(gw.available_frames, ['detector', 'icrs']).all() gw = wcs.WCS(output_frame=icrs, input_frame=detector) assert gw._pipeline[0].frame == "detector" - assert gw._pipeline[0][0] == "detector" + with pytest.warns(DeprecationWarning, match="Indexing a WCS.pipeline step is deprecated."): + assert gw._pipeline[0][0] == "detector" assert gw._pipeline[1].frame == "icrs" - assert gw._pipeline[1][0] == "icrs" + with pytest.warns(DeprecationWarning, match="Indexing a WCS.pipeline step is deprecated."): + assert gw._pipeline[1][0] == "icrs" assert np.in1d(gw.available_frames, ['detector', 'icrs']).all() with pytest.raises(NotImplementedError): gw(1, 2) @@ -154,7 +158,7 @@ def test_get_transform(): x, y = 1, 2 fx, fy = tr_forward(1, 2) assert_allclose(w.pipeline[0].transform(x, y), (fx, fy)) - assert_allclose(w.pipeline[0][1](x, y), (fx, fy)) + assert_allclose(w.pipeline[0].transform(x, y), (fx, fy)) assert_allclose((x, y), tr_back(*w(x, y))) assert(w.get_transform('detector', 'detector') is None) @@ -368,9 +372,14 @@ def test_grid_from_bounding_box_step(): def test_wcs_from_points(): np.random.seed(0) hdr = fits.Header.fromtextfile(get_pkg_data_filename("data/acs.hdr"), endcard=False) - with warnings.catch_warnings() as w: - warnings.simplefilter("ignore") + with pytest.warns(astwcs.FITSFixedWarning) as caught_warnings: + # this raises a warning unimportant for this testing the pix2world + # FITSFixedWarning(u'The WCS transformation has more axes (2) than + # the image it is associated with (0)') + # FITSFixedWarning: 'datfix' made the change + # 'Set MJD-OBS to 53436.000000 from DATE-OBS'. [astropy.wcs.wcs] w = astwcs.WCS(hdr) + assert len(caught_warnings) == 2 y, x = np.mgrid[:2046:20j, :4023:10j] ra, dec = w.wcs_pix2world(x, y, 1) fiducial = coord.SkyCoord(ra.mean()*u.deg, dec.mean()*u.deg, frame="icrs") @@ -517,7 +526,14 @@ def test_high_level_api(): class TestImaging(object): def setup_class(self): hdr = fits.Header.fromtextfile(get_pkg_data_filename("data/acs.hdr"), endcard=False) - self.fitsw = astwcs.WCS(hdr) + with pytest.warns(astwcs.FITSFixedWarning) as caught_warnings: + # this raises a warning unimportant for this testing the pix2world + # FITSFixedWarning(u'The WCS transformation has more axes (2) than + # the image it is associated with (0)') + # FITSFixedWarning: 'datfix' made the change + # 'Set MJD-OBS to 53436.000000 from DATE-OBS'. [astropy.wcs.wcs] + self.fitsw = astwcs.WCS(hdr) + assert len(caught_warnings) == 2 a_coeff = hdr['A_*'] a_order = a_coeff.pop('A_ORDER') b_coeff = hdr['B_*'] @@ -544,9 +560,9 @@ def setup_class(self): sky_cs = cf.CelestialFrame(reference_frame=coord.ICRS(), name='sky') det = cf.Frame2D(name='detector') wcs_forward = wcslin | tan | n2c - pipeline = [('detector', distortion), - ('focal', wcs_forward), - (sky_cs, None) + pipeline = [wcs.Step('detector', distortion), + wcs.Step('focal', wcs_forward), + wcs.Step(sky_cs, None) ] self.wcs = wcs.WCS(input_frame=det, @@ -555,7 +571,6 @@ def setup_class(self): self.xv, self.yv = xv, yv - @pytest.mark.filterwarnings('ignore') def test_distortion(self): sipx, sipy = self.fitsw.sip_pix2foc(self.xv, self.yv, 1) sipx = np.array(sipx) + 2048 @@ -759,7 +774,11 @@ def test_to_fits_sip_composite_frame_keep_axis(gwcs_cube_with_separable_spectral assert fw_hdr[f'CTYPE{ra_axis}'] == 'RA---TAN' assert fw_hdr['WCSAXES'] == 2 - fw = astwcs.WCS(fw_hdr) + with pytest.warns(astwcs.FITSFixedWarning, match='The WCS transformation has more axes'): + # this raises a warning unimportant for this testing the pix2world + # FITSFixedWarning(u'The WCS transformation has more axes (3) than + # the image it is associated with (2)') + fw = astwcs.WCS(fw_hdr) gskyval = w(1, 45, 55)[1:] assert np.allclose(gskyval, fw.all_pix2world([[1, 45, 55]], 0)[0][1:]) @@ -805,8 +824,7 @@ def test_to_fits_tab_cube(gwcs_3d_galactic_spectral): assert np.allclose(w(x, y, z), fits_wcs_user_bb.wcs_pix2world(x, y, z, 0), rtol=1e-6, atol=1e-7) - -@pytest.mark.skip(reason="Fails round-trip for -TAB axis 5") +@pytest.mark.filterwarnings('ignore:.*The WCS transformation has more axes.*') def test_to_fits_tab_7d(gwcs_7d_complex_mapping): # gWCS: w = gwcs_7d_complex_mapping @@ -874,7 +892,11 @@ def test_to_fits_no_sip_used(gwcs_spec_cel_time_4d): w = gwcs_spec_cel_time_4d # create FITS headers and -TAB headers - hdr, _ = w.to_fits(degree=3) + with pytest.warns(UserWarning, match='SIP distortion is not supported when the number'): + # UserWarning: SIP distortion is not supported when the number + # of axes in WCS is larger than 2. Setting 'degree' + # to 1 and 'max_inv_pix_error' to None. + hdr, _ = w.to_fits(degree=3) # check that FITS WCS is not using SIP assert not hdr['?_ORDER'] @@ -982,6 +1004,7 @@ def test_to_fits_tab_miri_image(): hdulist = fits.HDUList( [fits.PrimaryHDU(np.ones(w.pixel_n_dim * (2, )), hdr), bt] ) + fits_wcs = astwcs.WCS(hdulist[0].header, hdulist) # test points: @@ -1004,7 +1027,11 @@ def test_to_fits_tab_miri_lrs(): hdulist = fits.HDUList( [fits.PrimaryHDU(np.ones(w.pixel_n_dim * (2, )), hdr), bt[0]] ) - fits_wcs = astwcs.WCS(hdulist[0].header, hdulist) + with pytest.warns(astwcs.FITSFixedWarning, match='The WCS transformation has more axes'): + # this raises a warning unimportant for this testing the pix2world + # FITSFixedWarning(u'The WCS transformation has more axes (3) than + # the image it is associated with (2)') + fits_wcs = astwcs.WCS(hdulist[0].header, hdulist) # test points: (xmin, xmax), (ymin, ymax) = w.bounding_box @@ -1181,12 +1208,11 @@ def test_tabular_2d_quantity(): assert all(u.allclose(u.Quantity(b), [0, 2] * u.pix) for b in bb) -@pytest.mark.filterwarnings("error:.*Indexing a WCS.pipeline step is deprecated.*:DeprecationWarning") def test_initialize_wcs_with_list(): # test that you can initialize a wcs with a pipeline that is a list # containing both Step() and (frame, transform) tuples - # make pipline consisting of tuples and Steps + # make pipeline consisting of tuples and Steps shift1 = models.Shift(10 * u .pix) & models.Shift(2 * u.pix) shift2 = models.Shift(3 * u.pix) pipeline = [('detector', shift1), wcs.Step('extra_step', shift2)] @@ -1195,4 +1221,6 @@ def test_initialize_wcs_with_list(): pipeline.append(extra_step) # make sure no warnings occur when creating wcs with this pipeline - wcs.WCS(pipeline) + with warnings.catch_warnings(): + warnings.simplefilter("error") + wcs.WCS(pipeline) From 719f58a85bd540b7672613dcd1b895b18dcc76df Mon Sep 17 00:00:00 2001 From: Mihai Cara Date: Mon, 6 Dec 2021 22:14:08 -0500 Subject: [PATCH 08/12] Support all projections for lonpole computation --- CHANGES.rst | 4 +++ gwcs/tests/test_utils.py | 4 ++- gwcs/utils.py | 62 ++++++++++++++++++++++++---------------- 3 files changed, 44 insertions(+), 26 deletions(-) diff --git a/CHANGES.rst b/CHANGES.rst index 13615283..f9d531a3 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -23,6 +23,10 @@ Bug Fixes - Updated code in ``region.py`` with latest improvements and bug fixes from ``stsci.skypac.regions.py`` [#382] +- Added support to ``_compute_lon_pole()`` for computation of ``lonpole`` + for all projections from ``astropy.modeling.projections``. This also + extends support for different projections in ``wcs_from_fiducial()``. [#389] + New Features ^^^^^^^^^^^^ diff --git a/gwcs/tests/test_utils.py b/gwcs/tests/test_utils.py index 3bea3934..1d4bbcbd 100644 --- a/gwcs/tests/test_utils.py +++ b/gwcs/tests/test_utils.py @@ -35,13 +35,15 @@ def test_fits_transform(): def test_lon_pole(): tan = models.Pix2Sky_TAN() car = models.Pix2Sky_CAR() + azp = models.Pix2Sky_AZP(mu=-1.35, gamma=25.8458) sky_positive_lat = coord.SkyCoord(3 * u.deg, 1 * u.deg) sky_negative_lat = coord.SkyCoord(3 * u.deg, -1 * u.deg) assert_quantity_allclose(gwutils._compute_lon_pole(sky_positive_lat, tan), 180 * u.deg) assert_quantity_allclose(gwutils._compute_lon_pole(sky_negative_lat, tan), 180 * u.deg) assert_quantity_allclose(gwutils._compute_lon_pole(sky_positive_lat, car), 0 * u.deg) assert_quantity_allclose(gwutils._compute_lon_pole(sky_negative_lat, car), 180 * u.deg) - assert_quantity_allclose(gwutils._compute_lon_pole((0, 34 * u.rad), tan), 180 * u.deg) + assert_quantity_allclose(gwutils._compute_lon_pole((0, 0.34 * u.rad), tan), 180 * u.deg) + assert_quantity_allclose(gwutils._compute_lon_pole((1 * u.rad, 0.34 * u.rad), azp), 180 * u.deg) assert_allclose(gwutils._compute_lon_pole((1, -34), tan), 180) diff --git a/gwcs/utils.py b/gwcs/utils.py index 18d0848c..104558cf 100644 --- a/gwcs/utils.py +++ b/gwcs/utils.py @@ -12,6 +12,7 @@ from astropy import coordinates as coords from astropy import units as u from astropy.time import Time, TimeDelta +from astropy.wcs import Celprm # these ctype values do not include yzLN and yzLT pairs @@ -94,50 +95,61 @@ def get_values(units, *args): def _compute_lon_pole(skycoord, projection): """ - Compute the longitude of the celestial pole of a standard frame in the - native frame. + Compute the longitude of the celestial pole of a standard + frame in the native frame. - This angle then can be used as one of the Euler angles (the other two being skyccord) - to rotate the native frame into the standard frame ``skycoord.frame``. + This angle then can be used as one of the Euler angles + (the other two being skycoord) to rotate the native frame into the + standard frame ``skycoord.frame``. Parameters ---------- skycoord : `astropy.coordinates.SkyCoord`, or sequence of floats or `~astropy.units.Quantity` of length 2 - The fiducial point of the native coordinate system. - If tuple, its length is 2 + The celestial longitude and latitude of the fiducial point - typically + right ascension and declination. These are given by the ``CRVALia`` + keywords in ``FITS``. + projection : `astropy.modeling.projections.Projection` - A Projection instance. + A `~astropy.modeling.projections.Projection` model instance. Returns ------- - lon_pole : float or `~astropy/units.Quantity` - Native longitude of the celestial pole [deg]. + lonpole : float or `~astropy/units.Quantity` + Native longitude of the celestial pole in degrees. - TODO: Implement all projections - Currently this only supports Zenithal and Cylindrical. """ if isinstance(skycoord, coords.SkyCoord): - lat = skycoord.spherical.lat + lon = skycoord.spherical.lon.value + lat = skycoord.spherical.lat.value unit = u.deg + else: lon, lat = skycoord + unit = None + if isinstance(lon, u.Quantity): + lon = lon.to(u.deg).to_value() + unit = u.deg + if isinstance(lat, u.Quantity): + lat = lat.to(u.deg).to_value() unit = u.deg - else: - unit = None - if isinstance(projection, projections.Zenithal): - lon_pole = 180 - elif isinstance(projection, projections.Cylindrical): - if lat >= 0: - lon_pole = 0 - else: - lon_pole = 180 - else: - raise UnsupportedProjectionError("Projection {0} is not supported.".format(projection)) + + cel = Celprm() + cel.ref = [lon, lat] + cel.prj.code = projection.prjprm.code + pvrange = projection.prjprm.pvrange + if pvrange: + i1 = pvrange // 100 + i2 = i1 + (pvrange % 100) + 1 + cel.prj.pv = i1 * [None] + list(projection.prjprm.pv[i1:i2]) + cel.set() + + lonpole = cel.ref[2] if unit is not None: - lon_pole = lon_pole * unit - return lon_pole + lonpole = lonpole * unit + + return lonpole def get_projcode(wcs_info): From 8bcce4771d204125967d2b76a0224cfb4617bdc8 Mon Sep 17 00:00:00 2001 From: Mihai Cara Date: Mon, 25 Jul 2022 00:05:38 -0400 Subject: [PATCH 09/12] Fix a bug in the check for divergence in _fit_2D_poly --- CHANGES.rst | 3 +++ gwcs/wcs.py | 3 +++ 2 files changed, 6 insertions(+) diff --git a/CHANGES.rst b/CHANGES.rst index f9d531a3..ac73f1d1 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -5,6 +5,9 @@ Bug Fixes - Pin astropy min version to 5.0.4. [#404] +- Fixed a bug due to which the check for divercence in ``_fit_2D_poly()`` and + hence in ``to_fits()`` and ``to_fits_sip()`` was ignored. [#414] + New Features ^^^^^^^^^^^^ diff --git a/gwcs/wcs.py b/gwcs/wcs.py index c1279c4f..b57d6d04 100644 --- a/gwcs/wcs.py +++ b/gwcs/wcs.py @@ -2791,6 +2791,8 @@ def _fit_2D_poly(ntransform, npoints, degree, max_error, fit_poly_y(xin, yin)) if max_resid > prev_max_error: raise RuntimeError('Failed to achieve required error tolerance') + prev_max_error = max_resid + if verbose: print(f'Degree = {deg}, max_resid = {max_resid}') if max_resid < max_error: @@ -2804,6 +2806,7 @@ def _fit_2D_poly(ntransform, npoints, degree, max_error, if verbose: print('terminating condition met') break + return fit_poly_x, fit_poly_y, max_resid From 70bca777b815ce4bc6505a8dfe7cba80fc76ccd7 Mon Sep 17 00:00:00 2001 From: Mihai Cara Date: Sun, 24 Jul 2022 23:08:04 -0400 Subject: [PATCH 10/12] Fix how requested SIP fit accuracy and actual residuals are reported --- CHANGES.rst | 3 +++ gwcs/wcs.py | 25 +++++++++++++------------ 2 files changed, 16 insertions(+), 12 deletions(-) diff --git a/CHANGES.rst b/CHANGES.rst index ac73f1d1..021c9842 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -5,6 +5,9 @@ Bug Fixes - Pin astropy min version to 5.0.4. [#404] +- Corrected the reported requested forward SIP accuracy and reported fit + residuals by ``to_fits_sip()`` and ``to_fits()``. [#413] + - Fixed a bug due to which the check for divercence in ``_fit_2D_poly()`` and hence in ``to_fits()`` and ``to_fits_sip()`` was ignored. [#414] diff --git a/gwcs/wcs.py b/gwcs/wcs.py index b57d6d04..8676564e 100644 --- a/gwcs/wcs.py +++ b/gwcs/wcs.py @@ -1827,12 +1827,11 @@ def _to_fits_sip(self, celestial_group, keep_axis_position, yx, yy = ntransform(0, 1) 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, + degree, max_pix_error, plate_scale, u, v, undist_x, undist_y, ud, vd, undist_xd, undist_yd, verbose=verbose @@ -1852,7 +1851,7 @@ def _to_fits_sip(self, celestial_group, keep_axis_position, if max_inv_pix_error: fit_inv_poly_u, fit_inv_poly_v, max_inv_resid = _fit_2D_poly(ntransform, npoints, inv_degree, - max_inv_pix_error, + max_inv_pix_error, 1, U, V, u-U, v-V, Ud, Vd, ud-Ud, vd-Vd, verbose=verbose) @@ -1886,7 +1885,7 @@ def _to_fits_sip(self, celestial_group, keep_axis_position, hdr['B_ORDER'] = fit_poly_x.degree _store_2D_coefficients(hdr, sip_poly_x, 'A') _store_2D_coefficients(hdr, sip_poly_y, 'B') - hdr['sipmxerr'] = (max_resid * plate_scale, 'Max diff from GWCS (equiv pix).') + hdr['sipmxerr'] = (max_resid / plate_scale, 'Max diff from GWCS (equiv pix).') if max_inv_pix_error: hdr['AP_ORDER'] = fit_inv_poly_u.degree @@ -1927,7 +1926,7 @@ def _to_fits_sip(self, celestial_group, keep_axis_position, mat_kind = 'CD' del hdr['CDELT?'] - hdr['sipmxerr'] = (max_resid * plate_scale, 'Max diff from GWCS (equiv pix).') + hdr['sipmxerr'] = (max_resid / plate_scale, 'Max diff from GWCS (equiv pix).') # Construct CD matrix while remapping input axes. # We do not update comments to typical comments for CD matrix elements @@ -2743,7 +2742,7 @@ def _calc_approx_inv(self, max_inv_pix_error=5, inv_degree=None, npoints=16): fit_inv_poly_u, fit_inv_poly_v, max_inv_resid = _fit_2D_poly( ntransform, npoints, None, - max_inv_pix_error, + max_inv_pix_error, 1, undist_x, undist_y, u, v, undist_xd, undist_yd, ud, vd, verbose=True @@ -2755,7 +2754,7 @@ def _calc_approx_inv(self, max_inv_pix_error=5, inv_degree=None, npoints=16): (Shift(crpix[0]) & Shift(crpix[1]))) -def _fit_2D_poly(ntransform, npoints, degree, max_error, +def _fit_2D_poly(ntransform, npoints, degree, max_error, plate_scale, xin, yin, xout, yout, xind, yind, xoutd, youtd, verbose=False): @@ -2780,7 +2779,9 @@ def _fit_2D_poly(ntransform, npoints, degree, max_error, prev_max_error = float(np.inf) if verbose: - print(f'maximum_specified_error: {max_error}') + print(f'Maximum_specified_error: {max_error}') + max_error *= plate_scale + for deg in deglist: poly_x = Polynomial2D(degree=deg) poly_y = Polynomial2D(degree=deg) @@ -2794,20 +2795,20 @@ def _fit_2D_poly(ntransform, npoints, degree, max_error, prev_max_error = max_resid if verbose: - print(f'Degree = {deg}, max_resid = {max_resid}') + print(f'Degree = {deg}, max_resid = {max_resid / plate_scale}') if max_resid < max_error: # Check to see if double sampling meets error requirement. max_resid = _compute_distance_residual(xoutd, youtd, fit_poly_x(xind, yind), fit_poly_y(xind, yind)) if verbose: - print(f'Double sampling check: maximum residual={max_resid}') + print(f'Double sampling check: maximum residual={max_resid / plate_scale}') if max_resid < max_error: if verbose: - print('terminating condition met') + print('Terminating condition met') break - return fit_poly_x, fit_poly_y, max_resid + return fit_poly_x, fit_poly_y, max_resid / plate_scale def _make_sampling_grid(npoints, bounding_box, crpix): From 25341f5dfb84ec10d1a515219d3ceed8ca75d092 Mon Sep 17 00:00:00 2001 From: Mihai Cara Date: Tue, 6 Sep 2022 14:36:01 -0400 Subject: [PATCH 11/12] Fix incorrectly reported sipmxerr value --- CHANGES.rst | 2 +- gwcs/wcs.py | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/CHANGES.rst b/CHANGES.rst index 021c9842..80ddc4fc 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -6,7 +6,7 @@ Bug Fixes - Pin astropy min version to 5.0.4. [#404] - Corrected the reported requested forward SIP accuracy and reported fit - residuals by ``to_fits_sip()`` and ``to_fits()``. [#413] + residuals by ``to_fits_sip()`` and ``to_fits()``. [#413, #419] - Fixed a bug due to which the check for divercence in ``_fit_2D_poly()`` and hence in ``to_fits()`` and ``to_fits_sip()`` was ignored. [#414] diff --git a/gwcs/wcs.py b/gwcs/wcs.py index 8676564e..2c2fdcbb 100644 --- a/gwcs/wcs.py +++ b/gwcs/wcs.py @@ -1885,7 +1885,7 @@ def _to_fits_sip(self, celestial_group, keep_axis_position, hdr['B_ORDER'] = fit_poly_x.degree _store_2D_coefficients(hdr, sip_poly_x, 'A') _store_2D_coefficients(hdr, sip_poly_y, 'B') - hdr['sipmxerr'] = (max_resid / plate_scale, 'Max diff from GWCS (equiv pix).') + hdr['sipmxerr'] = (max_resid, 'Max diff from GWCS (equiv pix).') if max_inv_pix_error: hdr['AP_ORDER'] = fit_inv_poly_u.degree @@ -1926,7 +1926,7 @@ def _to_fits_sip(self, celestial_group, keep_axis_position, mat_kind = 'CD' del hdr['CDELT?'] - hdr['sipmxerr'] = (max_resid / plate_scale, 'Max diff from GWCS (equiv pix).') + hdr['sipmxerr'] = (max_resid, 'Max diff from GWCS (equiv pix).') # Construct CD matrix while remapping input axes. # We do not update comments to typical comments for CD matrix elements From 54ca2aa746fc3e527e08fdc00473b63296abcc96 Mon Sep 17 00:00:00 2001 From: Nadia Dencheva Date: Wed, 7 Sep 2022 12:27:25 -0400 Subject: [PATCH 12/12] prepare 0.18.2 release --- CHANGES.rst | 6 ++---- setup.cfg | 4 ++-- 2 files changed, 4 insertions(+), 6 deletions(-) diff --git a/CHANGES.rst b/CHANGES.rst index 80ddc4fc..c26e530c 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -1,14 +1,12 @@ -0.18.2 (Unreleased) +0.18.2 (2022-09-07) ------------------- Bug Fixes ^^^^^^^^^ -- Pin astropy min version to 5.0.4. [#404] - - Corrected the reported requested forward SIP accuracy and reported fit residuals by ``to_fits_sip()`` and ``to_fits()``. [#413, #419] -- Fixed a bug due to which the check for divercence in ``_fit_2D_poly()`` and +- Fixed a bug due to which the check for divergence in ``_fit_2D_poly()`` and hence in ``to_fits()`` and ``to_fits_sip()`` was ignored. [#414] New Features diff --git a/setup.cfg b/setup.cfg index 8c511252..5293d0ad 100644 --- a/setup.cfg +++ b/setup.cfg @@ -14,13 +14,13 @@ edit_on_github = False [options] zip_safe = False -python_requires = >=3.6 +python_requires = >=3.8 setup_requires = setuptools_scm install_requires = asdf >= 2.8.1 - astropy >= 5.0.4 + astropy >= 5.1 numpy scipy asdf_wcs_schemas