Fix an issue with units in wcs_from_points (#507)

CHANGES.rst

@@ -1,7 +1,9 @@
 0.22.0 (unreleased)
 -------------------
 
-- replace usages of ``copy_arrays`` with ``memmap`` [#503]
+- Replace usages of ``copy_arrays`` with ``memmap`` [#503]
+
+- Fix an issue with units in ``wcs_from_points``. [#507]
 
 0.21.0 (2024-03-10)
 -------------------

docs/gwcs/points_to_wcs.rst

@@ -76,7 +76,10 @@ can be used to convert coordinates from pixel to world.
   >>> gwcs_obj(36.235,642.215)    # doctest: +FLOAT_CMP
   (246.72158004206716, 43.46075091731673)
 
-Or equivalently 
-  >>> gwcs_obj.forward_transform(36.235,642.215)  # doctest: +FLOAT_CMP
+Or using the common WCS API 
+  >>> gwcs_obj.pixel_to_world_values(36.235,642.215)  # doctest: +FLOAT_CMP
   (246.72158004206716, 43.46075091731673)
 
+  >>> gwcs_obj.pixel_to_world(36.235,642.215)  # doctest: +FLOAT_CMP
+  <SkyCoord (ICRS): (ra, dec) in deg
+    (246.7215802, 43.46075103)>

gwcs/tests/test_wcs.py

@@ -423,6 +423,14 @@ def test_wcs_from_points():
     assert_allclose(newra, ra, atol=10**-6)
     assert_allclose(newdec, dec, atol=10**-6)
 
+    newra, newdec = w.pixel_to_world_values(x, y)
+    assert_allclose(newra, ra, atol=10**-6)
+    assert_allclose(newdec, dec, atol=10**-6)
+
+    newsky = w.pixel_to_world(x, y)
+    assert_allclose(newsky.data.lon.deg, ra, atol=10**-6)
+    assert_allclose(newsky.data.lat.deg, dec, atol=10**-6)
+
 
 def test_grid_from_bounding_box_2():
     bb = ((-0.5, 5.5), (-0.5, 4.5))

gwcs/wcstools.py

@@ -247,8 +247,9 @@ def wcs_from_points(xy, world_coords, proj_point='center',
         be passed in.
     projection : `~astropy.modeling.projections.Projection`
         A projection type. One of the projections in
-        `~astropy.modeling.projections.projcodes`. Defaults to TAN projection
-        (`astropy.modeling.projections.Sky2Pix_TAN`).
+        `~astropy.modeling.projections.projcodes`.
+        The direction is from sky to detector.
+        Defaults to TAN projection (`astropy.modeling.projections.Sky2Pix_TAN`).
     poly_degree : int
         Degree of polynomial model to be fit to data. Defaults to 4.
     polynomial_type : str
@@ -302,7 +303,7 @@ def wcs_from_points(xy, world_coords, proj_point='center',
                          "Only one of {} is supported.".format(polynomial_type,
                                                                supported_poly_types.keys()))
 
-    skyrot = models.RotateCelestial2Native(crval[0], crval[1], 180*u.deg)
+    skyrot = models.RotateCelestial2Native(crval[0].deg, crval[1].deg, 180)
     trans = (skyrot | projection)
     projection_x, projection_y = trans(lon, lat)
     poly = supported_poly_types[polynomial_type](poly_degree)