From c895e734de90cf869e3df731e80a39c15cd9142a Mon Sep 17 00:00:00 2001 From: Thomas Robitaille Date: Thu, 18 May 2023 12:18:01 +0100 Subject: [PATCH 1/5] Fix find_optimal_celestial_wcs to work with e.g. Solar frames, and make sure that we don't assume the output WCS is in degrees --- reproject/mosaicking/wcs_helpers.py | 34 ++++++++++++++++------------- 1 file changed, 19 insertions(+), 15 deletions(-) diff --git a/reproject/mosaicking/wcs_helpers.py b/reproject/mosaicking/wcs_helpers.py index 327e81477..d0c8551ab 100644 --- a/reproject/mosaicking/wcs_helpers.py +++ b/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,18 @@ 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) 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. From f5069badbccb4874468539a0a5b852e1778bc5f7 Mon Sep 17 00:00:00 2001 From: Thomas Robitaille Date: Thu, 18 May 2023 14:11:25 +0100 Subject: [PATCH 2/5] Added regression test for solar issues --- .../mosaicking/tests/test_wcs_helpers.py | 50 +++++++++++++++++++ 1 file changed, 50 insertions(+) diff --git a/reproject/mosaicking/tests/test_wcs_helpers.py b/reproject/mosaicking/tests/test_wcs_helpers.py index de5ae13e6..f74895231 100644 --- a/reproject/mosaicking/tests/test_wcs_helpers.py +++ b/reproject/mosaicking/tests/test_wcs_helpers.py @@ -287,3 +287,53 @@ 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 +""" + + +def test_solar_wcs(): + + # Regression test for issues that occurred when trying to find + # the optimal WCS for a set of solar WCSes + + import sunpy.coordinates + + wcs_ref = WCS(fits.Header.fromstring(SOLAR_HEADER, sep="\n"), fix=False) + + wcs1 = wcs_ref + wcs2 = wcs_ref.copy() + 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] From 1e093debd1da79bc1a537ada975c09f8e9ce7c22 Mon Sep 17 00:00:00 2001 From: Thomas Robitaille Date: Thu, 18 May 2023 14:22:12 +0100 Subject: [PATCH 3/5] Improve test --- reproject/mosaicking/tests/test_wcs_helpers.py | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/reproject/mosaicking/tests/test_wcs_helpers.py b/reproject/mosaicking/tests/test_wcs_helpers.py index f74895231..bd3ea7c66 100644 --- a/reproject/mosaicking/tests/test_wcs_helpers.py +++ b/reproject/mosaicking/tests/test_wcs_helpers.py @@ -316,17 +316,21 @@ def test_input_types(valid_celestial_input_shapes, iterable): """ +@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"), fix=False) + wcs_ref = WCS(fits.Header.fromstring(SOLAR_HEADER, sep="\n")) - wcs1 = wcs_ref - wcs2 = wcs_ref.copy() + 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)]) @@ -337,3 +341,5 @@ def test_solar_wcs(): 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) From 6da7f6cf09e45118f855b82f6db4ab75c339416e Mon Sep 17 00:00:00 2001 From: Thomas Robitaille Date: Thu, 18 May 2023 14:25:57 +0100 Subject: [PATCH 4/5] Whitespace --- reproject/mosaicking/tests/test_wcs_helpers.py | 1 - 1 file changed, 1 deletion(-) diff --git a/reproject/mosaicking/tests/test_wcs_helpers.py b/reproject/mosaicking/tests/test_wcs_helpers.py index bd3ea7c66..c15b6d79d 100644 --- a/reproject/mosaicking/tests/test_wcs_helpers.py +++ b/reproject/mosaicking/tests/test_wcs_helpers.py @@ -318,7 +318,6 @@ def test_input_types(valid_celestial_input_shapes, iterable): @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 From 5d1b68213939add5f5313cf72c464c6e10368ecd Mon Sep 17 00:00:00 2001 From: Thomas Robitaille Date: Thu, 18 May 2023 14:37:00 +0100 Subject: [PATCH 5/5] Fix issue with WCSes that have no default units --- reproject/mosaicking/wcs_helpers.py | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/reproject/mosaicking/wcs_helpers.py b/reproject/mosaicking/wcs_helpers.py index d0c8551ab..79f314461 100644 --- a/reproject/mosaicking/wcs_helpers.py +++ b/reproject/mosaicking/wcs_helpers.py @@ -216,6 +216,12 @@ def find_optimal_celestial_wcs( # 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.to_value(wcs_final.wcs.cunit[0]),