Fix bug in CelestialWCS._local near ra=0

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

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)