From 223ea23d09f2f29a7aceeae4b56c2a1c5c03109f Mon Sep 17 00:00:00 2001 From: Mike Jarvis Date: Fri, 5 Apr 2024 15:08:20 -0400 Subject: [PATCH] Fix bug in CelestialWCS._local near ra=0 --- galsim/wcs.py | 5 ++++- tests/test_wcs.py | 13 +++++++++++++ 2 files changed, 17 insertions(+), 1 deletion(-) diff --git a/galsim/wcs.py b/galsim/wcs.py index 7d39ef6e8c..014dfa1aa8 100644 --- a/galsim/wcs.py +++ b/galsim/wcs.py @@ -1344,7 +1344,7 @@ def _shiftOrigin(self, origin, world_origin, color): # If the class doesn't define something else, then we can approximate the local Jacobian # from finite differences for the derivatives of ra and dec. Very similar to the - # version for EuclideanWCS, but convert from dra, ddec to du, dv locallat at the given + # version for EuclideanWCS, but convert from dra, ddec to du, dv locally at at the given # position. def _local(self, image_pos, color): @@ -1360,6 +1360,9 @@ def _local(self, image_pos, color): xlist = np.array([ x0, x0+dx, x0-dx, x0, x0 ], dtype=float) ylist = np.array([ y0, y0, y0, y0+dy, y0-dy ], dtype=float) ra, dec = self._radec(xlist,ylist,color) + # Wrap ra to be near ra[0] + ra[ra < ra[0]-np.pi] += 2*np.pi + ra[ra > ra[0]+np.pi] -= 2*np.pi # Note: our convention is that ra increases to the left! # i.e. The u,v plane is the tangent plane as seen from Earth with +v pointing diff --git a/tests/test_wcs.py b/tests/test_wcs.py index bd5e4d3690..920b3ae185 100644 --- a/tests/test_wcs.py +++ b/tests/test_wcs.py @@ -3228,6 +3228,19 @@ def test_razero(): do_celestial_wcs(wcs, 'Astropy file '+file_name) do_wcs_image(wcs, 'Astropy near ra=0') + # Test that the local wcs comes out right where ra crosses 0. + image_pos = wcs.toImage(galsim.CelestialCoord(ra=0.*galsim.degrees, dec=47*galsim.degrees)) + print('0,47 is at ',image_pos) + local_wcs = wcs.local(image_pos) + print('local_wcs = ',local_wcs) + + # Should be very close to the same jacobian a spot a few arcsec over. + image_pos2 = wcs.toImage(galsim.CelestialCoord(ra=3*galsim.arcsec, dec=47*galsim.degrees)) + print('3 arcsec east = ',image_pos2) + local_wcs2 = wcs.local(image_pos2) + print('local_wcs2 = ',local_wcs2) + np.testing.assert_allclose(local_wcs.getMatrix(), local_wcs2.getMatrix(), rtol=1.e-3) + # This file is similar, but adjusted to be located near the south pole. file_name = 'pole.fits' wcs = galsim.AstropyWCS(file_name, dir=dir)