Merge pull request #210 from Jammy2211/feature/wcs_coordinates

Feature/wcs coordinates

Author: Jammy2211

autoarray/config/visualize/plots.yaml

@@ -21,7 +21,7 @@ fit:                                       # Settings for plots of all fits (e.g
 fit_imaging: {}                            # Settings for plots of fits to imaging datasets (e.g. FitImagingPlotter).
 inversion:                                 # Settings for plots of inversions (e.g. InversionPlotter).
   subplot_inversion: true                  # Plot subplot of all quantities in each inversion (e.g. reconstrucuted image, reconstruction)?
-  subplot_mappings: true                   # Plot subplot of the image-to-source pixels mappings of each pixelization?
+  subplot_mappings: false                  # Plot subplot of the image-to-source pixels mappings of each pixelization?
   data_subtracted: false                   # Plot individual plots of the data with the other inversion linear objects subtracted?
   reconstruction_noise_map: false          # Plot image of the noise of every mesh-pixel reconstructed value?
   sub_pixels_per_image_pixels: false       # Plot the number of sub pixels per masked data pixels?

autoarray/geometry/geometry_2d.py

@@ -132,10 +132,10 @@ def central_scaled_coordinates(self) -> Tuple[float, float]:
 
     def pixel_coordinates_2d_from(
         self, scaled_coordinates_2d: Tuple[float, float]
-    ) -> Tuple[float, float]:
+    ) -> Tuple[int, int]:
         """
-        Convert a 2D (y,x) scaled coordinate to a 2D (y,x) pixel coordinate, which are returned as floats such that they
-        include the decimal offset from each pixel's top-left corner relative to the input scaled coordinate.
+        Convert a 2D (y,x) scaled coordinate to a 2D (y,x) pixel coordinate, which are returned as integers such that
+        they do not include the decimal offset from each pixel's top-left corner relative to the input scaled coordinate.
 
         The conversion is performed according to the 2D geometry on a uniform grid, where the pixel coordinate origin
         is at the top left corner, such that the pixel [0,0] corresponds to the highest (most positive) y scaled
@@ -190,6 +190,41 @@ def scaled_coordinates_2d_from(
             origins=self.origin,
         )
 
+    def pixel_coordinates_wcs_2d_from(
+        self, scaled_coordinates_2d: Tuple[float, float]
+    ) -> Tuple[float, float]:
+        """
+        Convert a 2D (y,x) scaled coordinate to a 2D (y,x) WCS/FITS-style pixel coordinate.
+
+        The returned pixel coordinates follow the standard WCS convention:
+
+        - Coordinates are 1-based rather than 0-based, so that the centre of the top-left pixel is at (y, x) = (1.0, 1.0).
+        - Coordinates refer to pixel centres, not pixel corners.
+        - Values are continuous floats, so the fractional part encodes the sub-pixel offset from the pixel centre.
+
+        This differs from integer pixel-index conversions (e.g. ``pixel_coordinates_2d_from``), which return 0-based
+        indices associated with pixel corners/top-left positions.
+
+        The mapping from scaled coordinates to WCS pixel coordinates is defined by this geometry's ``origin``: scaled
+        coordinates are first shifted by the specified origin(s) before being converted using the pixel scale and
+        array shape. Changing ``origin`` therefore translates the returned WCS pixel coordinates by a constant offset.
+
+        Parameters
+        ----------
+        scaled_coordinates_2d
+            The 2D (y,x) coordinates in scaled units to be converted to WCS-style pixel coordinates.
+
+        Returns
+        -------
+        A 2D (y,x) WCS pixel coordinate, expressed as 1-based, pixel-centre, floating-point values.
+        """
+        return geometry_util.pixel_coordinates_wcs_2d_from(
+            scaled_coordinates_2d=scaled_coordinates_2d,
+            shape_native=self.shape_native,
+            pixel_scales=self.pixel_scales,
+            origins=self.origin,
+        )
+
     def scaled_coordinate_2d_to_scaled_at_pixel_centre_from(
         self, scaled_coordinate_2d: Tuple[float, float]
     ) -> Tuple[float, float]:

autoarray/geometry/geometry_util.py

@@ -357,6 +357,57 @@ def scaled_coordinates_2d_from(
     return (y_pixel, x_pixel)
 
 
+def pixel_coordinates_wcs_2d_from(
+    scaled_coordinates_2d: Tuple[float, float],
+    shape_native: Tuple[int, int],
+    pixel_scales: ty.PixelScales,
+    origins: Tuple[float, float] = (0.0, 0.0),
+) -> Tuple[float, float]:
+    """
+    Return FITS / WCS pixel coordinates (1-based, pixel-centre convention) as floats.
+
+    This function returns continuous pixel coordinates suitable for Astropy WCS
+    transforms (e.g. ``wcs_pix2world`` with ``origin=1``). Pixel centres lie at
+    integer values; for an image of shape ``(ny, nx)`` the geometric centre is::
+
+        ((ny + 1) / 2, (nx + 1) / 2)
+
+    e.g. ``(100, 100) -> (50.5, 50.5)``.
+
+    Parameters
+    ----------
+    scaled_coordinates_2d
+        The 2D (y, x) coordinates in scaled units which are converted to WCS
+        pixel coordinates.
+    shape_native
+        The (y, x) shape of the 2D array on which the scaled coordinates are
+        defined, used to determine the geometric centre in WCS pixel units.
+    pixel_scales
+        The (y, x) conversion factors from scaled units to pixel units.
+    origins
+        The (y, x) origin in scaled units about which the coordinates are
+        defined. The scaled coordinates are shifted by this origin before being
+        converted to WCS pixel coordinates.
+
+    Returns
+    -------
+    pixel_coordinates_wcs_2d
+        A 2D (y, x) WCS pixel coordinate in the 1-based, pixel-centre
+        convention, returned as floats.
+    """
+    ny, nx = shape_native
+
+    # Geometric centre in WCS pixel coordinates (1-based, pixel centres at integers)
+    ycen_wcs = (ny + 1) / 2.0
+    xcen_wcs = (nx + 1) / 2.0
+
+    # Continuous WCS pixel coordinates (NO int-cast, NO +0.5 binning)
+    y_wcs = (-scaled_coordinates_2d[0] + origins[0]) / pixel_scales[0] + ycen_wcs
+    x_wcs = (scaled_coordinates_2d[1] - origins[1]) / pixel_scales[1] + xcen_wcs
+
+    return (y_wcs, x_wcs)
+
+
 def transform_grid_2d_to_reference_frame(
     grid_2d: np.ndarray, centre: Tuple[float, float], angle: float, xp=np
 ) -> np.ndarray: