Fix numerical_inverse API issue (#565)

Author: WilliamJamieson

CHANGES.rst

@@ -1,6 +1,8 @@
 0.25.0 (unreleased)
 -------------------
 
+- Fix API issue with ``wcs.numerical_inverse``. [#565]
+
 0.24.0 (2025-02-04)
 -------------------
 

gwcs/examples.py

@@ -5,8 +5,8 @@
 from astropy.time import Time
 
 from . import coordinate_frames as cf
+from . import geometry, wcs
 from . import spectroscopy as sp
-from . import wcs
 
 # frames
 DETECTOR_1D_FRAME = cf.CoordinateFrame(
@@ -651,3 +651,56 @@ def gwcs_with_pipeline_celestial():
     ]
 
     return wcs.WCS(pipeline)
+
+
+def gwcs_romanisim():
+    targ_pos = coord.SkyCoord(ra=0 * u.deg, dec=0 * u.deg)
+
+    ra_ref = targ_pos.ra.to(u.deg).value
+    dec_ref = targ_pos.dec.to(u.deg).value
+
+    # v2_ref, v3_ref are in arcsec, but RotationSequence3D wants degrees,
+    # so start by scaling by 3600.
+    rot = models.RotationSequence3D([0, 0, 0, dec_ref, -ra_ref], "zyxyz")
+
+    # V2V3 are in arcseconds, while SphericalToCartesian expects degrees,
+    # so again start by scaling by 3600
+    tel2sky = (
+        (models.Scale(1 / 3600) & models.Scale(1 / 3600))
+        | geometry.SphericalToCartesian(wrap_lon_at=180)
+        | rot
+        | geometry.CartesianToSpherical(wrap_lon_at=360)
+    )
+    tel2sky.name = "v23tosky"
+
+    detector = cf.Frame2D(name="detector", axes_order=(0, 1), unit=(u.pix, u.pix))
+    v2v3 = cf.Frame2D(
+        name="v2v3",
+        axes_order=(0, 1),
+        axes_names=("v2", "v3"),
+        unit=(u.arcsec, u.arcsec),
+    )
+    v2v3vacorr = cf.Frame2D(
+        name="v2v3vacorr",
+        axes_order=(0, 1),
+        axes_names=("v2", "v3"),
+        unit=(u.arcsec, u.arcsec),
+    )
+    world = cf.CelestialFrame(reference_frame=coord.ICRS(), name="world")
+
+    va_corr = models.Identity(2)
+    zen2v2v3 = models.EulerAngleRotation(0, -90, 0, "xyz") | (
+        models.Scale(3600) & models.Scale(3600)
+    )
+    tanproj = models.Pix2Sky_Gnomonic()
+    pix2tan = models.Scale(0.11 / 3600) & models.Scale(0.11 / 3600)
+    distortion = pix2tan | tanproj | zen2v2v3
+
+    pipeline = [
+        wcs.Step(detector, distortion),
+        wcs.Step(v2v3, va_corr),
+        wcs.Step(v2v3vacorr, tel2sky),
+        wcs.Step(world, None),
+    ]
+
+    return wcs.WCS(pipeline)

gwcs/tests/conftest.py

@@ -156,3 +156,8 @@ def cart_to_spher():
 @pytest.fixture
 def gwcs_with_pipeline_celestial():
     return examples.gwcs_with_pipeline_celestial()
+
+
+@pytest.fixture
+def gwcs_romanisim():
+    return examples.gwcs_romanisim()

gwcs/tests/test_wcs.py

@@ -1804,3 +1804,12 @@ def test_bounding_box_with_units():
 
     w_gwcs.invert(4 * u.deg, 5 * u.deg)
     w_gwcs.to_fits(bounding_box=([0, 100] * u.pix, [0, 100] * u.pix))
+
+
+def test_direct_numerical_inverse(gwcs_romanisim):
+    xy = (128, 256)
+    coord = gwcs_romanisim(*xy)
+    ra_dec = (np.radians(c) * u.rad for c in coord)
+    out = gwcs_romanisim.numerical_inverse(*ra_dec)
+
+    assert_allclose(xy, out)

gwcs/wcs/_wcs.py

@@ -352,7 +352,7 @@ def _call_backward(
         else:
             # Always strip units for numerical inverse
             args = self._remove_units_input(args, self.output_frame)
-            result = self.numerical_inverse(
+            result = self._numerical_inverse(
                 *args,
                 with_bounding_box=with_bounding_box,
                 fill_value=fill_value,
@@ -678,6 +678,30 @@ def numerical_inverse(
          [2.76552923e-05 1.14789013e-05]]
 
         """  # noqa: E501
+        return self._numerical_inverse(
+            *self._remove_units_input(args, self.output_frame),
+            tolerance=tolerance,
+            maxiter=maxiter,
+            adaptive=adaptive,
+            detect_divergence=detect_divergence,
+            quiet=quiet,
+            with_bounding_box=with_bounding_box,
+            fill_value=fill_value,
+            **kwargs,
+        )
+
+    def _numerical_inverse(
+        self,
+        *args,
+        tolerance=1e-5,
+        maxiter=30,
+        adaptive=True,
+        detect_divergence=True,
+        quiet=True,
+        with_bounding_box=True,
+        fill_value=np.nan,
+        **kwargs,
+    ):
         if kwargs.pop("with_units", False):
             msg = (
                 "Support for with_units in numerical_inverse has been removed, "