Fixes for solar frames and non-degree units in find_optimal_celestial_wcs

reproject/mosaicking/tests/test_wcs_helpers.py

@@ -287,3 +287,58 @@ def test_input_types(valid_celestial_input_shapes, iterable):
 
  assert_wcs_allclose(wcs_test, wcs_ref)
  assert shape_test == shape_ref
+
+
+SOLAR_HEADER = """
+CRPIX1 = -1374.571094981584 / [pix]
+CRPIX2 = 2081.629159922445 / [pix]
+CRDATE1 = '2017-01-01T00:00:00.000'
+CRDATE2 = '2017-01-01T00:00:00.000'
+CRVAL1 = -619.0078311637853
+CRVAL2 = -407.000970936774
+CDELT1 = 0.01099999994039536
+CDELT2 = 0.01099999994039536
+CUNIT1 = 'arcsec '
+CUNIT2 = 'arcsec '
+CTYPE1 = 'HPLN-TAN'
+CTYPE2 = 'HPLT-TAN'
+PC1_1 = 0.966887196065055
+PC1_2 = -0.01087372434907635
+PC2_1 = 0.01173971407248916
+PC2_2 = 0.9871195868097251
+LONPOLE = 180.0 / [deg]
+DATEREF = '2022-06-02T17:22:53.220'
+OBSGEO-X= -5466045.256954942 / [m]
+OBSGEO-Y= -2404388.737412784 / [m]
+OBSGEO-Z= 2242133.887690042 / [m]
+SPECSYS = 'TOPOCENT'
+VELOSYS = 0.0
+"""
+
+
+@pytest.mark.filterwarnings("ignore::astropy.wcs.wcs.FITSFixedWarning")
+def test_solar_wcs():
+ # Regression test for issues that occurred when trying to find
+ # the optimal WCS for a set of solar WCSes
+
+ pytest.importorskip("sunpy", minversion="2.1.0")
+
+ # Make sure the WCS <-> frame functions are registered
+ import sunpy.coordinates
+
+ wcs_ref = WCS(fits.Header.fromstring(SOLAR_HEADER, sep="\n"))
+
+ wcs1 = wcs_ref.deepcopy()
+ wcs2 = wcs_ref.deepcopy()
+ wcs2.wcs.crpix[0] -= 4096
+
+ wcs, shape = find_optimal_celestial_wcs([((4096, 4096), wcs1), ((4096, 4096), wcs2)])
+
+ wcs.wcs.set()
+
+ assert wcs.wcs.ctype[0] == wcs_ref.wcs.ctype[0]
+ assert wcs.wcs.ctype[1] == wcs_ref.wcs.ctype[1]
+ assert wcs.wcs.cunit[0] == wcs_ref.wcs.cunit[0]
+ assert wcs.wcs.cunit[1] == wcs_ref.wcs.cunit[1]
+
+ assert shape == (4281, 8237)

reproject/mosaicking/wcs_helpers.py

@@ -123,7 +123,7 @@ def find_optimal_celestial_wcs(
 
  # We start off by looping over images, checking that they are indeed
  # celestial images, and building up a list of all corners and all reference
- # coordinates in celestial (ICRS) coordinates.
+ # coordinates in the frame of reference of the first image.
 
  corners = []
  references = []
@@ -162,18 +162,19 @@ def find_optimal_celestial_wcs(
  xc = np.array([-0.5, nx - 0.5, nx - 0.5, -0.5])
  yc = np.array([-0.5, -0.5, ny - 0.5, ny - 0.5])
 
- # We have to do .frame here to make sure that we get an ICRS object
+ # We have to do .frame here to make sure that we get a frame object
  # without any 'hidden' attributes, otherwise the stacking below won't
  # work.
- corners.append(wcs.pixel_to_world(xc, yc).icrs.frame)
+ corners.append(wcs.pixel_to_world(xc, yc).transform_to(frame).frame)
 
  if isinstance(wcs, WCS):
- # We now figure out the reference coordinate for the image in ICRS. The
- # easiest way to do this is actually to use pixel_to_skycoord with the
- # reference position in pixel coordinates. We have to set origin=1
- # because crpix values are 1-based.
+ # We now figure out the reference coordinate for the image in the
+ # frame of the first image. The easiest way to do this is actually
+ # to use pixel_to_skycoord with the reference position in pixel
+ # coordinates. We have to set origin=1 because crpix values are
+ # 1-based.
  xp, yp = wcs.wcs.crpix
- references.append(pixel_to_skycoord(xp, yp, wcs, origin=1).icrs.frame)
+ references.append(pixel_to_skycoord(xp, yp, wcs, origin=1).transform_to(frame).frame)
 
  # Find the pixel scale at the reference position - we take the minimum
  # since we are going to set up a header with 'square' pixels with the
@@ -183,7 +184,7 @@ def find_optimal_celestial_wcs(
 
  else:
  xp, yp = (nx - 1) / 2, (ny - 1) / 2
- references.append(wcs.pixel_to_world(xp, yp).icrs.frame)
+ references.append(wcs.pixel_to_world(xp, yp).transform_to(frame).frame)
 
  xs = np.array([xp, xp, xp + 1])
  ys = np.array([yp, yp + 1, yp])
@@ -192,7 +193,7 @@ def find_optimal_celestial_wcs(
  dy = abs(cs[0].separation(cs[1]).deg)
  resolutions.append(min(dx, dy))
 
- # We now stack the coordinates - however the ICRS class can't do this
+ # We now stack the coordinates - however the frame classes can't do this
  # so we have to use the high-level SkyCoord class.
  corners = SkyCoord(corners)
  references = SkyCoord(references)
@@ -212,15 +213,24 @@ def find_optimal_celestial_wcs(
  if resolution is None:
  resolution = np.min(resolutions) * u.deg
 
- # Determine the resolution in degrees
- cdelt = resolution.to(u.deg).value
-
  # Construct WCS object centered on position
  wcs_final = celestial_frame_to_wcs(frame, projection=projection)
 
+ if wcs_final.wcs.cunit[0] == "":
+ wcs_final.wcs.cunit[0] = "deg"
+
+ if wcs_final.wcs.cunit[1] == "":
+ wcs_final.wcs.cunit[1] = "deg"
+
  rep = reference.represent_as("unitspherical")
- wcs_final.wcs.crval = rep.lon.degree, rep.lat.degree
- wcs_final.wcs.cdelt = -cdelt, cdelt
+ wcs_final.wcs.crval = (
+ rep.lon.to_value(wcs_final.wcs.cunit[0]),
+ rep.lat.to_value(wcs_final.wcs.cunit[1]),
+ )
+ wcs_final.wcs.cdelt = (
+ -resolution.to_value(wcs_final.wcs.cunit[0]),
+ resolution.to_value(wcs_final.wcs.cunit[1]),
+ )
 
  # For now, set crpix to (1, 1) and we'll then figure out where all the
  # images fall in this projection, then we'll adjust crpix.