gwcs

import numpy as np from astropy.modeling import models from astropy import units as u from astropy import coordinates as coord from gwcs import wcs from gwcs import coordinate_frames as cf # Define the transformations pixel_scale = 0.05 * u.arcsec / u.pixel detector_center_pixel = (500, 500) reference_sky_coord = coord.SkyCoord(ra=20 * u.deg, dec=30 * u.deg, frame='icrs') # 1. Shift pixel coordinates to have (0,0) at detector center pixel_to_intermediate_frame = models.Shift(-detector_center_pixel[0] * u.pixel) & \ models.Shift(-detector_center_pixel[1] * u.pixel) # 2. Apply pixel scale intermediate_to_celestial_angles = models.Scale(pixel_scale) & models.Scale(pixel_scale) # 3. Tangent projection (TAN) from celestial angles to sky coordinates celestial_angles_to_sky = models.Pix2Sky_TAN() # 4. Rotate and shift to final celestial position (simple case for illustration) # In a real scenario, this would be more complex, involving proper rotation models # and matching reference points. Here we directly set the reference_sky_coord. # Define the frames detector_frame = cf.Frame2D(name="detector", axes_names=('x', 'y'), unit=(u.pixel, u.pixel)) celestial_frame = cf.CelestialFrame(reference_frame=reference_sky_coord.frame, name='icrs', unit=(u.deg, u.deg)) # Build the GWCS pipeline pipeline = [ wcs.Step(detector_frame, pixel_to_intermediate_frame | intermediate_to_celestial_angles | celestial_angles_to_sky), wcs.Step(celestial_frame, None) ] # Create the WCS object image_wcs = wcs.WCS(pipeline) # Example transformation pixel_coords = np.array([[100, 200], [500, 500]]) * u.pixel world_coords = image_wcs.pixel_to_world(pixel_coords[:, 0], pixel_coords[:, 1]) print(f"Input pixel coordinates:\n{pixel_coords}") print(f"Output world coordinates:\n{world_coords}")